mm_viv.m

用无网格法计算VAN DE POL 振子

% THIS IS PROGRAMMED TO RESEARCH THE DYNAMICS OF A MARINE RISER SUBJECTED TO VIV 
% BY HE CHANGKIANG ACCORDING TO THE PAPER 'RESEARCH ON THE DYNAMIC CHARACTERISTICS 
  ...AND RESPONSES OF A MARINE RISER SUBJECTED OT VIV'(Guo Haiyan,Lou Min and Dong 
% Xiaolin). HERE, MESHLESS METHOD "MWLS" INTRODUCE BY ZHANG XIONG AND LIU YAN IN THE BOOK 
% "Meshless Method" IS APPLIED.
% MWLS - MESHLESS WEIGHTED LEAST-SQUARE PROGRAM 
% THE FUNCTIONAL IS CALCULATED BY NODAL SUNMMATION, WHICH AVOIDS NUMERICAL INTEGRATION.
% Date: 3/20/2006
% Revised: 3/20/2006
% Copyright (C) 2006, HE CHANGKIANG, COLLEGE OF CEVIL ENGINEERING,HIT

clear all; clc; warning on
format long;  format compact;

tcputime = cputime;

% SET UP NODE DISTRIBUTION
L = 100;
m = 3;                                 % NUMBER OF BASIS FUNCTION
npoints = 21;                          % NUMBER OF CALCULATION POINTS
nnodes = npoints / 1;                  % NUMBER OF NODES
lx = L / (npoints-1);                  % DISTANCE BETWEEN ADJACENT CALCULATION POINTS
lxi = L / (nnodes-1);                  % DISTANCE BETWEEN ADJACENT NODES
tot = 6000;
tt = zeros(1,tot);                     % TIME CONSUMPTION OF EACH WHILE-LOOP
delta_t = 0.01; 
tseries = delta_t*[0.0:tot]';
x = [0.0 : lx : L];                    % CALCULATION POINTS' COORDINATES
xi = [0.0 : lxi : L];                  % NODAL COORDINATES        

% SET COMPUTATIONAL PARAMETERS, INCLUDING THE DIMENSION OF SUPPORT AND THE VALUE OF PENALTY.
scale = 2.5;                           % SCALE FACTOR
dm = scale*lxi*ones(1,nnodes);         % SUPPORT RADIUS
curpen = 1e7;                          % PENALTY OF TRACTION BOUNDARY
disppen = 1e3*curpen;                  % PENALTY OF DISPLACEMENT BOUNDARY  

% MATRICE INITIALIZATION
K_r = zeros(nnodes,nnodes);            % STIFFNESS MATRIX OF RISER 
K_o = zeros(nnodes,nnodes);            % STIFFNESS MATRIX OF OSSILATOR
P = zeros(nnodes,1);                   % LOAD VECTOR

y = zeros(npoints,tot+1);
q = zeros(npoints,tot+1);
yh = 1e-3.*ones(npoints,3);  % STORE THE LAST THREE RESULTS(INCLUDE THE ONE BEEN CALCULTING)
qh = 1e-3.*ones(npoints,3);

y_try = 1e-3.*ones(npoints,1);
q_try = 1e-3.*ones(npoints,1);
y_ass = zeros(npoints,1);
q_ass = zeros(npoints,1);

% iteration criterion
epsilon_y = 1e-5;
epsilon_q = 1e-5;

% LOOP OVER TATOL TIME SERIES
for i = 1 : tot
    k = i, kk = 0;
    time = i*delta_t;
    
    tic
    while norm(y_try-y_ass) > epsilon_y | norm(q_try-q_ass) > epsilon_q
        
        kk = kk + 1;
        
        y_ass = y_try;  q_ass = q_try;
        
        % EVALUATE TOTAL STIFFNESS AND LOADS MATRICES OF RISER AND OSSILATOR'S
        % DYNAMIC FUNCTION
        [K_r,K_o,F_r,F_o,phi] = KFPRO(m,nnodes,xi,npoints,x,dm,'SPLIN',0,disppen,curpen,delta_t,time,...
                                      y(:,i), yh(:,3), yh(:,2), yh(:,1), q(:,i), qh(:,2), qh(:,1));                     
                                    
        % NODAL FICTIOUS VALUE
        yh(:,3) = K_r\F_r;                            
        qh(:,3) = K_o\F_o;   
        
        % JUST DEBUGGING!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!      
        if kk > 100, error('look out!~iteration cannot converge!!!'); end
        
        % CALCULATE NODAL DISPLACEMENTS AND LIFT COFFICEENT
        y_try = phi * yh(:,3);                     
        q_try = phi * qh(:,3); 
    end
    tt(i) = toc;    
    kk = kk
    
    % CONSERVE NEW FICTITIOUS NODAL VALUES
    yh(:,1) = yh(:,2);  yh(:,2) = yh(:,3);
    qh(:,1) = qh(:,2);  qh(:,2) = qh(:,3);
    
    % CONSERVE RESULTS
    y(:,i+1) = y_ass;  q(:,i+1) = q_ass;
    
    % LET y_try BE DIFFER FROM y_ass
    y_try = 0.9.*y_ass;  q_try = 0.9.*q_ass;
    
end

% DRAW RESULT CURVE
for i = 1 : tot/100
    figure(i);
    plot(y(:,i),x);
end

figure;
subplot(1,2,1);
hu = plot(tseries,y((npoints+1)/2,:));
xlabel('Coordinate','Fontsize',12);
ylabel('Displacement','Fontsize',12);
legend(hu,'Riser');
grid on;

subplot(1,2,2);
hs = plot(tseries,q((npoints+1)/2,:));
xlabel('Coordinate','Fontsize',12);
ylabel('Lift Cofficient','Fontsize',12);
legend(hs,'Ossilator');
grid on;

% Output nodal displacements and stresses
% fid1 = fopen('riserdis.dat','w');
% fid2 = fopen('liftcoef.dat','w');
% 
% fprintf(fid1,'%10s%10s\n', 'x', 'y');
% fprintf(fid2,'%10s%10s\n', 'x', 'q');
% 
% for j = 1 : nnodes
%    fprintf(fid1,'%10.4f%10.4f\n', x(j), y(j));
%    fprintf(fid2,'%10.4f%10.4f\n', x(j), q(j));
% end
%    
% fclose(fid1);
% fclose(fid2);

e = cputime - tcputime
figure
plot(tt);