www.pudn.com > 87361026FEA.rar > incmp_flow2d.m, change:2006-02-28,size:7609b


% incmp_flow2d.m
%
% by Jack Chessa
% Northwestern University
%
% A 2D program that solves incompressible flow using a projection method
% and characteristic stabilization. Assumes only triangular elements.
clear
colordef black
state = 0;

%******************************************************************************
%***                             I N P U T                                  ***
%******************************************************************************

% HERE WE DEFINE SOME SYSTEM PARPMETERS

% GMSH FILE PARPMETERS
gmshFileName = 'step.msh';  % name of gmsh input file
domainID     = 13;          % zone id for domain mesh
inletBndyID   = 9;          % zone id for displacement boundary mesh
noSlipBndyID   = 11;        % zone id for displacement boundary mesh
pressBndyID   = 10;         % zone id for traction boundary mesh

% MATERIAL PROPERTIES
rho   = 10;    % density
mu    = 0.3;   % poisson ratio
gravity = [0 0]; % gravity vector

% TIME STEPPING PARAMETERS
tstart = 0;       % initial time
tfin = 0.12;      % final time
deltInit = .001;  % initial time step
Cr = 0.75;        % C.F.L. number

% INLET CONDITION
tss=0.5;
uinss=1;

%******************************************************************************
%***                    P R E - P R O C E S S I N G                         ***
%******************************************************************************

% READ GMSH FILE
[node,elements,elemType]=msh2mlab(gmshFileName);

% GET NODE AND CONNECTIVITY MATRICIES
node=node(:,1:2);                    % a nodal coordinate matix (reference config)
element  = elements{domainID};       % connectivity matrix for interior
inletBndy = elements{inletBndyID};  % connectivity matrix for displacement boundary
noSlipBndy = elements{noSlipBndyID}; % connectivity matrix for displacement boundary
pressBndy = elements{pressBndyID};   % connectivity matrix for traction boundary

% CHECK FOR BAD JACOBIANS
tricheck(node,element);

% GET NODES ON NO-SLIP BOUNDARY
noSlipNodes=unique(noSlipBndy);
inletNodes=unique(inletBndy);
pressNodes=unique(pressBndy);

numnode=size(node,1);    % number of nodes
numelem=size(element,1); % number of elements

% DEFINE SYSTEM DATA STRUCTURES
xs=1:numnode;                % x portion of u and v vectors
ys=(numnode+1):2*numnode;    % y portion of u and v vectors
u=zeros(2*numnode,1);        % nodal velocity vector
us=zeros(2*numnode,1);       % nodal fractional velocity vector
p=zeros(numnode,1);          % nodal pressure vector
x=node;                      % nodal coordinates in current configureation
fgrav=zeros(2*numnode,1);    % gravity force vector
M=sparse(2*numnode,2*numnode);        % mass matrix
Kp=sparse(numnode,numnode);           % pressure stiffness matrix
K=sparse(2*numnode,2*numnode);        % viscous stiffness matrix
G=sparse(2*numnode,numnode);          % gradient matirx

% DEFINE INITIAL CONDITIONS
u(:)=0;
p(:)=0;

% DEFINE BOUNDARY CONDITIONS


%******************************************************************************
%***                          P R O C E S S I N G                           ***
%******************************************************************************

%%%%%%%%%%%%%%%%%%% CALCULATE SYTEM MATRICIES
[W,Q]=quadrature( 2, 'TRIANGULAR', 2 );  % define quadrature rule
for e=1:numelem                          % start of element loop
  
  sctr=element(e,:);
  sctrx=sctr;
  sctry=sctr+numnode;
  
  for q=1:size(W,1)                      % quadrature loop
    pt=Q(q,:);                           % quadrature point
    wt=W(q);                             % quadrature weight
    [N,dNdxi]=lagrange_basis('T3',pt);   % element shape functions
    jac=x(sctr,:)'*dNdxi;  % element Jacobian matrix w.r.t ref configuration    
    dNdx=dNdxi*inv(jac);
    detj=det(jac);
    
    M(sctrx,sctrx)=M(sctrx,sctrx)+N*N'*detj*wt;
    M(sctry,sctry)=M(sctry,sctry)+N*N'*detj*wt;
    
    Kp(sctr,sctr)=Kp(sctr,sctr)+dNdx*dNdx'*detj*wt;
    
    K(sctrx,sctrx)=K(sctrx,sctrx)+(mu/rho)*dNdx*dNdx'*detj*wt;
    
    G(sctrx,sctr)=G(sctrx,sctr)+N*dNdx(:,1)'*detj*wt;
    G(sctry,sctr)=G(sctry,sctr)+N*dNdx(:,2)'*detj*wt;
    
    fgrav(sctrx)=fgrav(sctrx)+N*gravity(1)*detj*wt;
    fgrav(sctry)=fgrav(sctry)+N*gravity(2)*detj*wt;
    
  end                                    % of quadrature loop
end   % of element loop
M=full(sum(M)');  % lump the mass matrix

% ***************************** TEMPORAL LOOP *********************************
tn=tstart;
delt=deltInit;
step=0;

while ( tn <= tfin ) 
  tn=tn+delt;
  step=step+1;
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % STEP 1- get fractional time step
  f=delt*(-K*u+fgrav);
  [W,Q]=quadrature( 2, 'TRIANGULAR', 2 );  % define quadrature rule
  for e=1:numelem                          % start of element loop
    
    sctr=element(e,:);
    sctrx=sctr;
    sctry=sctr+numnode;
    
    for q=1:size(W,1)                      % quadrature loop
      pt=Q(q,:);                           % quadrature point
      wt=W(q);                             % quadrature weight
      [N,dNdxi]=lagrange_basis('T3',pt);   % element shape functions
      jac=x(sctr,:)'*dNdxi;  % element Jacobian matrix w.r.t ref configuration    
      dNdx=dNdxi*inv(jac);
      detj=det(jac);
      
      upt=[N'*u(sctrx); N'*u(sctry)];
      Lpt=[dNdx'*u(sctrx), dNdx'*u(sctry)]';
      advect=Lpt*upt;
      
      f(sctrx)=f(sctrx)-delt*N*advect(1)*detj*wt;  % advection terms
      f(sctry)=f(sctry)-delt*N*advect(2)*detj*wt;
      
      % CG stabilization terms
      f(sctrx)=f(sctrx)-delt^2/2*(dNdx*upt)*(advect(1)+gravity(1))*detj*wt;  
      f(sctry)=f(sctry)-delt^2/2*(dNdx*upt)*(advect(2)+gravity(2))*detj*wt;
      
    end                                    % of quadrature loop
  end   % of element loop

  % solve for fractional velocity
  us=u+f./M;
  
  % apply boundary conditions
  us([noSlipNodes noSlipNodes+numnode])=0;
  if ( tn<tss )
    uxin=tn/tss*uinss;
  else
    uxin=uinss;
  end
   
  uyin=0;
  unin=uxin;
  us(inletNodes)=uxin;
  us(inletNodes+numnode)=uyin;

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % STEP 2 - solve Poisson pressure equation
  
  % get influx force
  f=G'*us;
  %f=zeros(numnode,1);
  [W,Q]=quadrature( 1, 'GAUSS', 1 );  % define quadrature rule
  for e=1:size(inletBndy,1)
    
    sctr=inletBndy(e,:);
    
    for q=1:size(W,1)                      % quadrature loop
      pt=Q(q,:);                           % quadrature point
      wt=W(q);                             % quadrature weight
      N=lagrange_basis('L2',pt);           % element shape functions
      jac=norm( x(sctr(2),:)-x(sctr(1),:) )/2;
      detj=det(jac);
      
      f(sctr)=f(sctr)-N*(-1)*unin*detj*wt;
      
    end % of quadrature loop
  end  % of element loop
  
  % apply bcs
  Kbc=Kp;
  bcwt=sum(diag(Kbc));
  for i=1:size(pressNodes,1)
    dof=pressNodes(i);
    Kbc(dof,:)=0;  
    Kbc(:,dof)=0;   
    Kbc(dof,dof)=bcwt;  
    f(dof)=0;
  end
  
  % solve
  p=(rho/delt)*( Kbc\f );

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % STEP 3 - solve for velocity correction
  du=-(delt/rho)*(G*p)./M;
  du(inletNodes)=0;
  du(inletNodes+numnode)=0;
  du(noSlipNodes)=0;
  du(noSlipNodes+numnode)=0;
  
  u=us+du;
  
  % PLOT RESULTS
  clf
  %trisurf(element,x(:,1),x(:,2),0*x(:,1),u(ys))
  trisurf(element,x(:,1),x(:,2),0*x(:,1),us(xs))
  %trisurf(element,x(:,1),x(:,2),0*x(:,1),p)
  shading interp
  axis equal
  view(2)
  colorbar
  
  pause(0.1)
  
end % of time loop
% *****************************************************************************
state=1;