CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Wiki > Code: Lid-driven cavity using pressure-free velocity form (2)

Code: Lid-driven cavity using pressure-free velocity form (2)

From CFD-Wiki

(Difference between revisions)
Jump to: navigation, search
(Program script for pressure-free velocity method (LDC4W.m))
(references)
 
(3 intermediate revisions not shown)
Line 43: Line 43:
The output consists of graphic plots of contour levels of the stream function and the pressure levels.  
The output consists of graphic plots of contour levels of the stream function and the pressure levels.  
-
This simplified version for this Wiki resulted from removal of computation of the vorticity, a restart capability, and production of publication-quality plots from one of the research codes used with the paper. It can be compared with a simpler code given in this Wiki to understand coding for higher-order elements. This code can be extended to fifth-order elements (and generalized bicubic elements), which use the same degrees-of-freedom, by simply replacing the functions evaluating the elements.
+
This simplified version for this Wiki resulted from removal of computation of the vorticity, a restart capability, and production of publication-quality plots from one of the research codes used with the paper. It can be compared with a simpler code given in this Wiki to understand coding for higher-order elements. This code can be extended to fifth-order elements (and generalized bicubic elements), which use the same degrees-of-freedom, by simply replacing the class definitions evaluating the elements.
 +
 
 +
This modified version posted 3/15/2013 updates the code to MatLab version R2012b. The element properties are now defined in classes, and the code has been modified to accept class definitions for simple-cubic, quartic and quintic divergence-free elements on quadrilaterals and quadratic and quartic divergence-free elements on triangles.  
===Lid-driven cavity Matlab script===
===Lid-driven cavity Matlab script===
Line 85: Line 87:
</pre>
</pre>
-
===[[PFV program script | Program script for pressure-free velocity method ('''LDC4W.m''')]]===
+
===[[PFV4 program script | Program script for pressure-free velocity method ('''LDC4W.m''')]]===
-
<pre>
+
-
%LDC4W            LID-DRIVEN CAVITY
+
-
% Finite element solution of the 2D Navier-Stokes equation using
+
-
%  4-node, 6-DOF/node, C1,quartic-derived rectangular stream function
+
-
%  elements for the Lid-Driven Cavity.
+
-
%
+
-
% This could be characterized as a VELOCITY-STREAM FUNCTION or
+
-
%  STREAM FUNCTION-VELOCITY method.
+
-
%
+
-
% Reference:  J. T. Holdeman, "A Hermite finite element method
+
-
%  for incompressible fluid flow",
+
-
%  Int. J. Numer. Meth. Fluids, 64, P376-408 (2010).
+
-
%
+
-
% Simplified Wiki version
+
-
%
+
-
% The rectangular problem domain is defined between Cartesian
+
-
%  coordinates Xmin & Xmax and Ymin & Ymax.
+
-
% The computational grid has NumNx nodes in the x-direction
+
-
%  and NumNy nodes in the y-direction.
+
-
% The nodes and elements are numbered column-wise from the 
+
-
%  upper left corner to the lower right corner.
+
-
%
+
-
% This script calls the user-defined functions:
+
-
%    ELS4424r    - class with definitions of velocity shape functions
+
-
%    DMatW        - to evaluate element diffusion matrix
+
-
%    CMatW        - to evaluate element convection matrix
+
-
%    GetPres3W    - to evaluate the pressure
+
-
%    regrade      - to regrade the mesh
+
-
% Uses
+
-
%  ilu          - incomplete LU preconditioner
+
-
%  gmres        - iterative solver
+
-
% Indirectly uses:
+
-
%    Gquad2      - Gauss integraion rules
+
-
%    ELG3412r    - class of pressure basis functions
+
-
%
+
-
% Jonas Holdeman    Updated August 2007, revised March 2013
+
-
clear all;
+
===[[PFV4 diffusion matrix | Diffusion matrix for pressure-free velocity method ('''DMatW.m''')]]===
-
eu = ELS4424r;      % class of functions for velocity
+
-
disp('Lid-driven cavity');
+
===[[PFV4 convection matrix | Convection matrix for pressure-free velocity method ('''CMatW.m''')]]===
-
disp([' Four-node, 24 DOF, ' eu.name ' basis.']);
+
-
% -------------------------------------------------------------
+
===[[PFV4 quartic velocity class| Modified-quartic velocity element class for pressure-free velocity method ('''ELS4424f.m''')]]===
-
  nnd = eu.nnodes;        % nnodes = number of element nodes
+
-
  ndof = eu.nndofs;      % nndofs = number of velocity dofs per node, (6);
+
-
  nd2=ndof*ndof; ND=1:ndof;  % Number of DOF per node - do not change!!
+
-
% -------------------------------------------------------------
+
-
ETstart=clock;
+
-
% Parameters for GMRES solver
+
===[[PFV4 cubic pressure class| Simple-cubic pressure element class for pressure-free velocity method ('''ELG3412f.m''')]]===
-
GM.Tol = 1.e-13;
+
-
GM.MIter = 20;
+
-
GM.MRstrt = 6;
+
-
% parameters for ilu preconditioner
+
-
% Decrease su.droptol if ilu preconditioner fails
+
-
su.type='ilutp';
+
-
su.droptol=1.e-5;
+
-
% Suggested relaxation parameters for given Reynolds number
+
===[[PFV4 get pressure | Consistent pressure for pressure-free velocity method ('''GetPresW.m''')]]===
-
% Re          100    1000    3200  5000  7500  10000  12500
+
-
% ExpCR1                      .210
+
-
% ExpCRO                      .402
+
-
% RelxFacO:    1.0    1.11    .880  .84    .79    .76    .70
+
-
% CritFac:                  1.135
+
-
% Define the problem geometry:
+
===[[PFV4 quadature rules| Gauss quadrature rules on rectangle ('''GQuad2.m''')]]===
-
Xmin = 0.0; Xmax = 1.0;  Ymin = 0.0; Ymax = 1.0;
+
-
% Set mesh grading parameters (set to 1 if no grading).
+
===[[PFV4 mesh regrade | Grade node spacing ('''regrade.m''')]]===
-
% See below for explanation of use of parameters.
+
-
xgrd = .75; ygrd=.75;  % Set xgrd = 1, ygrd=1 for uniform mesh; %
+
-
 
+
-
% Set " RefineBoundary=1 " for additional refinement at boundary,
+
-
%  i.e., split first element along boundary into two.
+
-
RefineBoundary=1;
+
-
 
+
-
%    DEFINE THE MESH
+
-
% Set number of elements in each direction
+
-
NumEx = 18;  NumEy = NumEx;
+
-
 
+
-
% PLEASE CHANGE OR SET NUMBER OF ELEMENTS ABOVE TO CHANGE/SET NUMBER OF NODES!
+
-
NumNx=NumEx+1;    NumNy=NumEy+1;
+
-
 
+
-
%  Define problem parameters:
+
-
% Lid velocity (don't change, vary Re instead unless lid moves to left)
+
-
Vlid=1.;
+
-
 
+
-
% Specify Reynolds number Re, set relaxation factor
+
-
Re=1000.;
+
-
 
+
-
% factor for under/over-relaxation
+
-
RelxFac = 1.;
+
-
 
+
-
% Start with simple non-linear iteration, then switch to Newton method
+
-
%  when non-linear corrections are less than ItThreshold.
+
-
% CAUTION! Large Re may require very small threshold.
+
-
% If ItThreshold = 0 (or sufficiently small) method will never switch.
+
-
ItThreshold = 5.e-2;
+
-
 
+
-
% NUMBER OF NONLINEAR ITERATIONS
+
-
MaxNLit=12; %
+
-
 
+
-
%Selected and labeled contour levels, Ghia, Ghia & Shin for plotting
+
-
clSel=[-.1175,-.115,-.11,-.1,-.09,-.07,-.05,-.03,-.01,-.0001,1e-7,1e-5,...
+
-
        1e-4,5e-4,1e-3,1.5e-3,3e-3];
+
-
LclSel=[1,2,3,4,6,8,10,12,14,16];
+
-
 
+
-
%Contour and labeled levels, Ghia, Ghia & Shin, Re=100, 400, 1000, 3200, ...
+
-
clGGS=[-1.e-7,-1.e-5,-1.e-4,-.01,-.03,-.05,-.07,-.09,-.1,-.11,-.115,-.1175,...
+
-
    1.e-8,1.e-7,1.e-6,1.e-5,5.e-5,1.e-4,2.5e-4,5.e-4,1.e-3,1.5e-3,3.e-3];
+
-
LclGGS=[1,3,4,6,8,10,19,23];
+
-
 
+
-
% Select contour level option
+
-
%CL=clGGS;  CLL=LclGGS;  % GGS
+
-
CL=clSel;  CLL=LclSel;  % GGS (select)
+
-
 
+
-
if (Vlid<0), CL=-CL; end
+
-
 
+
-
% -------------------------------------------------------------
+
-
 
+
-
% Viscosity for specified Reynolds number
+
-
nu=Vlid*(Xmax-Xmin)/Re;
+
-
 
+
-
if (xgrd~=1 || ygrd~=1), meshsp='graded'; else meshsp='uniform'; end
+
-
disp(['Reynolds number = ' num2str(Re) ',  ' num2str(NumEx) 'x' ...
+
-
      num2str(NumEy) ' element ' meshsp ' mesh' ]);
+
-
disp(['Maximum number of nonlinear iterations = ' num2str(MaxNLit)]);
+
-
pause(1);
+
-
 
+
-
% Grade the mesh spacing if desired, call regrade(x,agrd,e).
+
-
% if e=0: refine both sides, 1: refine upper, 2: refine lower
+
-
% if agrd=xgrd|ygrd is the parameter which controls grading, then
+
-
%  if agrd=1 then leave array unaltered.
+
-
%  if agrd<1 then refine (make finer) towards the ends.
+
-
%  if agrd>1 then refine (make finer) towards the center.
+
-
%
+
-
%  Generate equally-spaced nodal coordinates and refine if desired.
+
-
if (RefineBoundary==1)
+
-
  XNc=linspace(Xmin,Xmax,NumNx-2);
+
-
  XNc=[XNc(1),(.62*XNc(1)+.38*XNc(2)),XNc(2:end-1), ...
+
-
      (.38*XNc(end-1)+.62*XNc(end)),XNc(end)];
+
-
  YNc=linspace(Ymax,Ymin,NumNy-2);
+
-
  YNc=[YNc(1),(.62*YNc(1)+.38*YNc(2)),YNc(2:end-1), ...
+
-
      (.38*YNc(end-1)+.62*YNc(end)),YNc(end)];
+
-
else
+
-
  XNc=linspace(Xmin,Xmax,NumNx);
+
-
  YNc=linspace(Ymax,Ymin,NumNy);
+
-
end
+
-
if xgrd ~= 1, XNc=regrade(XNc,xgrd,0); end;
+
-
if ygrd ~= 1, YNc=regrade(YNc,ygrd,0); end;
+
-
 
+
-
[Xgrid,Ygrid]=meshgrid(XNc,YNc);% Generate the x- and y-coordinate meshes.
+
-
 
+
-
% Area-based nodal weights
+
-
wx=zeros(1,NumNx);  wy=zeros(1,NumNy);
+
-
wx(1)=.5*(XNc(2)-XNc(1));
+
-
wx(2:NumNx-1)=.5*(XNc(3:NumNx)-XNc(1:NumNx-2));
+
-
wx(NumNx)=.5*(XNc(NumNx)-XNc(NumNx-1));
+
-
wy(1)=.5*(YNc(1)-YNc(2));
+
-
wy(2:NumNx-1)=.5*(YNc(1:NumNy-2)-YNc(3:NumNy));
+
-
wy(NumNx)=.5*(YNc(NumNy-1)-YNc(NumNy));
+
-
Wa=wy'*wx;
+
-
+
-
% Allocate storage for fields
+
-
psi0=zeros(NumNy,NumNx);
+
-
u0=zeros(NumNy,NumNx);
+
-
v0=zeros(NumNy,NumNx);
+
-
pxx0=zeros(NumNy,NumNx);
+
-
pxy0=zeros(NumNy,NumNx);
+
-
pyy0=zeros(NumNy,NumNx);
+
-
 
+
-
%--------------------Begin grid plot-----------------------
+
-
% ********************** FIGURE 1 *************************
+
-
% Plot the grid
+
-
figure(1);
+
-
clf;
+
-
orient portrait; orient tall;
+
-
subplot(2,2,1);
+
-
hold on;
+
-
plot([Xmax;Xmin],[YNc;YNc],'k');
+
-
plot([XNc;XNc],[Ymax;Ymin],'k');
+
-
hold off;
+
-
axis([Xmin,Xmax,Ymin,Ymax]);
+
-
axis equal;
+
-
axis image;
+
-
title([num2str(NumNx) 'x' num2str(NumNy) ...
+
-
      ' node mesh for Lid-driven cavity']);
+
-
 
+
-
%-------------- End plotting Figure 1 ----------------------
+
-
 
+
-
NumNod=NumNx*NumNy;    % total number of nodes
+
-
MaxDof=ndof*NumNod;        % maximum number of degrees of freedom
+
-
EBC.Mxdof=ndof*NumNod;        % maximum number of degrees of freedom
+
-
 
+
-
NodXY=zeros(NumNod,2);
+
-
nn2nft=zeros(NumNod,2); % node number -> nf & nt
+
-
NodNdx=zeros(NumNod,2);
+
-
% Generate lists of active nodal indices, freedom number & type
+
-
nn=0;  nf=-ndof+1;  nt=1;      %  _______
+
-
for nx=1:NumNx                %  |      |
+
-
  for ny=1:NumNy              %  |      |
+
-
    nn=nn+1;                  %  |_______|
+
-
    NodNdx(nn,:)=[nx,ny];
+
-
    NodXY(nn,:)=[Xgrid(ny,nx),Ygrid(ny,nx)];
+
-
    nf=nf+ndof;              % all nodes have nd (6) dofs
+
-
    nn2nft(nn,:)=[nf,nt];  % dof number & type (all nodes type 1)
+
-
  end;
+
-
end;
+
-
 
+
-
nf2nnt=zeros(MaxDof,2);  % (node, type) associated with dof
+
-
nd=0; dd=1:ndof;
+
-
for nn=1:NumNod
+
-
  for nf=1:ndof
+
-
    nf2nnt(nd+nf,:)=[nn,nf];
+
-
  end
+
-
  nd=nd+ndof;
+
-
end
+
-
 
+
-
NumEl=NumEx*NumEy;
+
-
 
+
-
% Generate element connectivity, from upper left to lower right.
+
-
Elcon=zeros(NumEl,nnd);
+
-
ne=0;  LY=NumNy;
+
-
for nx=1:NumEx
+
-
  for ny=1:NumEy
+
-
    ne=ne+1;
+
-
    Elcon(ne,1)=1+ny+(nx-1)*LY;
+
-
    Elcon(ne,2)=1+ny+nx*LY;
+
-
    Elcon(ne,3)=1+(ny-1)+nx*LY;
+
-
    Elcon(ne,4)=1+(ny-1)+(nx-1)*LY;
+
-
  end  % loop on ny
+
-
end  % loop on nx
+
-
 
+
-
% Begin essential boundary conditions, allocate space
+
-
MaxEBC = ndof*2*(NumNx+NumNy-2);
+
-
EBC.dof=zeros(MaxEBC,1);  % Degree-of-freedom index 
+
-
EBC.typ=zeros(MaxEBC,1);  % Dof type (1,2,3,4,5,6)
+
-
EBC.val=zeros(MaxEBC,1);  % Dof value
+
-
 
+
-
nc=0;
+
-
for nf=1:MaxDof
+
-
  ni=nf2nnt(nf,1);
+
-
  nx=NodNdx(ni,1);  ny=NodNdx(ni,2);
+
-
  x=XNc(nx);        y=YNc(ny);
+
-
  if(x==Xmin || x==Xmax)
+
-
    nt=nf2nnt(nf,2);
+
-
    switch nt;
+
-
    case {1, 2, 3, 5, 6}                                % psi,u,v,pxy,pyy
+
-
      nc=nc+1; EBC.typ(nc)=nt; EBC.dof(nc)=nf; EBC.val(nc)=0;
+
-
    end  % switch (type)
+
-
  elseif (y==Ymin)
+
-
    nt=nf2nnt(nf,2);
+
-
    switch nt;
+
-
    case {1, 2, 3, 4, 5}                                % psi,u,v,pxx,pxy
+
-
      nc=nc+1; EBC.typ(nc)=nt; EBC.dof(nc)=nf; EBC.val(nc)=0;
+
-
    end  % switch (type)
+
-
  elseif (y==Ymax)
+
-
    nt=nf2nnt(nf,2);
+
-
    switch nt;
+
-
    case {1, 3, 4, 5}                                      % psi,v,pxx,pxy
+
-
      nc=nc+1; EBC.typ(nc)=nt; EBC.dof(nc)=nf; EBC.val(nc)=0;
+
-
    case 2
+
-
      nc=nc+1; EBC.typ(nc)=nt; EBC.dof(nc)=nf; EBC.val(nc)=Vlid;  % u
+
-
    end  % switch (type)
+
-
  end  % if (boundary)
+
-
end  % for nf
+
-
EBC.num=nc;
+
-
 
+
-
if (size(EBC.typ,1)>nc)  % Truncate arrays if necessary
+
-
  EBC.typ=EBC.typ(1:nc);
+
-
  EBC.dof=EBC.dof(1:nc);
+
-
  EBC.val=EBC.val(1:nc);
+
-
end    % End ESSENTIAL (Dirichlet) boundary conditions
+
-
 
+
-
% partion out essential (Dirichlet) dofs
+
-
p_vec = (1:EBC.Mxdof)';        % List of all dofs
+
-
EBC.p_vec_undo = zeros(1,EBC.Mxdof);
+
-
% form a list of non-diri dofs
+
-
EBC.ndro = p_vec(~ismember(p_vec, EBC.dof)); % list of non-diri dofs
+
-
% calculate p_vec_undo to restore Q to the original dof ordering
+
-
EBC.p_vec_undo([EBC.ndro;EBC.dof]) = (1:EBC.Mxdof); %p_vec';
+
-
 
+
-
Q=zeros(MaxDof,1); % Allocate space for solution (dof) vector
+
-
 
+
-
% Initialize fields to boundary conditions
+
-
for k=1:EBC.num
+
-
  Q(EBC.dof(k))=EBC.val(k);
+
-
end;
+
-
 
+
-
errpsi=zeros(NumNy,NumNx);  % error correction for iteration
+
-
 
+
-
% Arrays for convergence norm info
+
-
MxNL=max(1,MaxNLit);
+
-
np0=zeros(1,MxNL);    % psi
+
-
nv0=zeros(1,MxNL);    % u & v
+
-
 
+
-
Qs=[];
+
-
 
+
-
% Generate and assemble element matrices
+
-
 
+
-
Dmat = spalloc(MaxDof,MaxDof,56*MaxDof);  % to save the diffusion matrix
+
-
Vdof=zeros(ndof,nnd);
+
-
Xe=zeros(2,nnd);      % coordinates of element corners
+
-
 
+
-
NLitr=0;
+
-
ItType=0;  % Begin with simple iteration
+
-
 
+
-
% <<<<<<<< BEGIN NONLINEAR ITERATION >>>>>>>>>>>>
+
-
while (NLitr<MaxNLit), NLitr=NLitr+1; 
+
-
     
+
-
tclock=clock;  % Start assembly time <<<<<<<<<
+
-
% Generate and assemble element matrices
+
-
Mat=spalloc(MaxDof,MaxDof,56*MaxDof);  % (36)
+
-
RHS=spalloc(MaxDof,1,MaxDof);
+
-
 
+
-
if(ItType==0 && NLitr>1 && nv0(NLitr-1)<ItThreshold && np0(NLitr-1)<ItThreshold)
+
-
  ItType=1; disp(' >>>>> Begin Newton method >>>>>');  % <<<<<<<<<<<<<<<
+
-
end 
+
-
 
+
-
% BEGIN GLOBAL MATRIX ASSEMBLY
+
-
for ne=1:NumEl 
+
-
+
-
  Xe(1:2,1:nnd)=NodXY(Elcon(ne,1:nnd),1:2)';
+
-
     
+
-
  if NLitr == 1   
+
-
%    Fluid element diffusion matrix, save on first iteration   
+
-
    [Emat,Rndx,Cndx] = DMatW(eu,Xe,Elcon(ne,:),nn2nft);
+
-
    Dmat=Dmat+sparse(Rndx,Cndx,Emat,MaxDof,MaxDof);  % Global diff matrix
+
-
  end
+
-
 
+
-
  if (NLitr>1) 
+
-
%    Get stream function and velocities, loop over local nodes of element
+
-
    for n=1:nnd 
+
-
      Vdof(ND,n)=Q((nn2nft(Elcon(ne,n),1)-1)+ND);
+
-
    end
+
-
%    Fluid element convection matrix, first iteration uses Stokes equation.
+
-
    if (ItType==0)
+
-
        % For simple iteration
+
-
      [Emat,Rndx,Cndx] = CMatW(eu,Xe,Elcon(ne,:),nn2nft,Vdof);
+
-
    else
+
-
        % For Newton iteration
+
-
      [Emat,Rndx,Cndx,Rcm,RcNdx] = CMatW(eu,Xe,Elcon(ne,:),nn2nft,Vdof);
+
-
      RHS = RHS + sparse(RcNdx,1,Rcm,MaxDof,1);
+
-
    end
+
-
    Mat=Mat+sparse(Rndx,Cndx,Emat,MaxDof,MaxDof);  % Assemble global
+
-
  end  % end  if(NLitr>1)
+
-
 
+
-
end;  % loop ne over elements
+
-
% END GLOBAL MATRIX ASSEMBLY
+
-
 
+
-
Mat = Mat + nu*Dmat;    % Add in cached/saved global diffusion matrix
+
-
 
+
-
disp(['(' num2str(NLitr) ') Matrix assembly complete, elapsed time = '...
+
-
      num2str(etime(clock,tclock)) ' sec']);  % Assembly time <<<<<<<<<<<
+
-
pause(1);
+
-
 
+
-
Q0 = Q;    % Save dof values
+
-
 
+
-
% Solve system, Start solution time  <<<<<<<<<<<<<<
+
-
tclock=clock;
+
-
 
+
-
RHSr=RHS(EBC.ndro)-Mat(EBC.ndro,EBC.dof)*EBC.val;
+
-
Matr=Mat(EBC.ndro,EBC.ndro);
+
-
Qs=Q(EBC.ndro);
+
-
 
+
-
[Lm,Um] = ilu(Matr,su);      % incomplete LU
+
-
Qr = gmres(Matr,RHSr,GM.MIter,GM.Tol,GM.MRstrt,Lm,Um,Qs); % GMRES
+
-
 
+
-
Q=[Qr;EBC.val];        % Augment active dofs with esential (Dirichlet) dofs
+
-
Q=Q(EBC.p_vec_undo);  % Restore natural order
+
-
 
+
-
stime=etime(clock,tclock); % Solution time <<<<<<<<<<<<<<
+
-
 
+
-
% ********** APPLY RELAXATION FACTOR ***************************
+
-
if(NLitr>1), Q=RelxFac*Q+(1-RelxFac)*Q0; end
+
-
% ****************************************************
+
-
+
-
% Compute change and copy dofs to field arrays
+
-
dsqp=0; dsqv=0;
+
-
for k=1:MaxDof
+
-
  ni=nf2nnt(k,1);
+
-
  nx=NodNdx(ni,1);  ny=NodNdx(ni,2);
+
-
  switch nf2nnt(k,2) % switch on dof type
+
-
    case 1
+
-
      dsqp=dsqp+Wa(ny,nx)*(Q(k)-Q0(k))^2;
+
-
      errpsi(ny,nx)=Q0(k)-Q(k);  psi0(ny,nx)=Q(k);
+
-
    case 2
+
-
      dsqv=dsqv+Wa(ny,nx)*(Q(k)-Q0(k))^2;  u0(ny,nx)=Q(k);
+
-
    case 3
+
-
      dsqv=dsqv+Wa(ny,nx)*(Q(k)-Q0(k))^2;  v0(ny,nx)=Q(k);
+
-
    case 4
+
-
      pxx0(ny,nx)=Q(k);
+
-
    case 5
+
-
      pxy0(ny,nx)=Q(k);
+
-
    case 6
+
-
      pyy0(ny,nx)=Q(k);
+
-
  end  % switch on dof type
+
-
end  % for
+
-
 
+
-
np0(NLitr)=sqrt(dsqp); dP=np0(NLitr);
+
-
nv0(NLitr)=sqrt(dsqv);
+
-
 
+
-
if (np0(NLitr)<=1e-15||nv0(NLitr)<=1e-15)
+
-
  MaxNLit=NLitr; np0=np0(1:MaxNLit); nv0=nv0(1:MaxNLit);  end;
+
-
%  Print solution time
+
-
disp(['(' num2str(NLitr) ') Solution time for linear system = '...
+
-
      num2str(etime(clock,tclock)) ' sec']);
+
-
+
-
%---------- Begin plot of intermediate results ----------
+
-
% ********************** FIGURE 1 *************************
+
-
figure(1);
+
-
 
+
-
% Stream function          (intermediate)
+
-
subplot(2,2,3);
+
-
contour(Xgrid,Ygrid,psi0,8,'k');  % Plot contours (trajectories)
+
-
axis([Xmin,Xmax,Ymin,Ymax]);
+
-
axis equal;  axis image;
+
-
title(['Lid-driven cavity,  Re=' num2str(Re)]);
+
-
 
+
-
% Plot iteration convergence info
+
-
subplot(2,2,2);
+
-
semilogy(1:NLitr,nv0(1:NLitr),'k-+',1:NLitr,np0(1:NLitr),'k-o');
+
-
xlabel('Nonlinear iteration number');
+
-
ylabel('Nonlinear correction');
+
-
axis square;
+
-
title(['Iteration conv.,  Re=' num2str(Re)]);
+
-
legend('U','Psi');
+
-
 
+
-
% Iteration change
+
-
subplot(2,2,4);
+
-
if dP>0
+
-
contour(Xgrid,Ygrid,errpsi,8,'k');  % Plot contours (trajectories)
+
-
axis([Xmin,Xmax,Ymin,Ymax]);
+
-
axis equal;  axis image;
+
-
title('Iteration correction');
+
-
end % dP>0
+
-
% ********************** END FIGURE 1 *************************
+
-
%----------  End plot of intermediate results  ---------
+
-
 
+
-
if (nv0(NLitr)<1e-15), break; end  % Terminate iteration if non-significant
+
-
 
+
-
end;  % <<< END NONLINEAR ITERATION  (while )
+
-
 
+
-
format short g;
+
-
disp('Convergence results by iteration: velocity, stream function');
+
-
disp(['nv0:  ' num2str(nv0)]); disp(['np0:  ' num2str(np0)]);
+
-
 
+
-
% >>>>>>>>>>>>>> BEGIN PRESSURE RECOVERY <<<<<<<<<<<<<<<<<<
+
-
 
+
-
% Essential pressure boundary condition
+
-
% Index of node to apply pressure BC, value at node
+
-
PBCnx=fix((NumNx+1)/2);  % Apply at center of mesh
+
-
PBCny=fix((NumNy+1)/2);
+
-
PBCnod=0;
+
-
for k=1:NumNod
+
-
  if (NodNdx(k,1)==PBCnx && NodNdx(k,2)==PBCny), PBCnod=k; break; end
+
-
end
+
-
if (PBCnod==0), error('Pressure BC node not found');
+
-
else
+
-
  EBCp.nodn = PBCnod;  % Pressure BC node number
+
-
  EBCp.val = 0;  % set P = 0.
+
-
end
+
-
 
+
-
[P,Px,Py] = GetPres3W(eu,NodXY,Elcon,nn2nft,Q,EBCp,0); % nu
+
-
P=reshape(P,NumNy,NumNx);
+
-
Px = reshape(Px,NumNy,NumNx);
+
-
Py = reshape(Py,NumNy,NumNx);
+
-
 
+
-
% >>>>>>>>>>>>> End pressure recovery <<<<<<<<<<<<<
+
-
 
+
-
%--------------Plot of final result-------------
+
-
% ********************** CONTINUE FIGURE 1 *************************
+
-
 
+
-
figure(1);
+
-
textid=['4424X' num2str(NumEx)];
+
-
 
+
-
% Stream function      (final)
+
-
subplot(2,2,3);
+
-
[CT,hn]=contour(Xgrid,Ygrid,psi0,CL,'k');  % Plot contours (trajectories)
+
-
clabel(CT,hn,CL(CLL));
+
-
hold on;
+
-
plot([Xmin,Xmin,Xmax,Xmax,Xmin],[Ymax,Ymin,Ymin,Ymax,Ymax],'k');
+
-
pcolor(Xgrid,Ygrid,psi0);
+
-
shading interp; %flat;  % or interp
+
-
hold off
+
-
axis([Xmin,Xmax,Ymin,Ymax]);
+
-
axis equal;  axis image;
+
-
title(['Stream lines, ' num2str(NumNx) 'x' num2str(NumNy) ...
+
-
    ' mesh,  Re=' num2str(Re)]);
+
-
 
+
-
% Plot pressure contours    (final)
+
-
subplot(2,2,4);
+
-
CPL=[-.002,0,.02,.05,.07,.09,.11,.12,.17,.3];  %
+
-
[CT,hn]=contour(Xgrid,Ygrid,P,CPL,'k');  % Plot pressure contours
+
-
clabel(CT,hn,CPL([3,5,7,10]));
+
-
hold on;
+
-
plot([Xmin,Xmin,Xmax,Xmax,Xmin],[Ymax,Ymin,Ymin,Ymax,Ymax],'k');
+
-
%pcolor(Xgrid,Ygrid,P);
+
-
%shading interp; %flat;
+
-
hold off;
+
-
axis([Xmin,Xmax,Ymin,Ymax]);
+
-
axis equal;  axis image;
+
-
title(['Quartic pressure contours, Re=' num2str(Re)]);
+
-
 
+
-
% ********************* END FIGURE 1 *************************
+
-
% ------------ End final plot --------------
+
-
 
+
-
disp(['Total elapsed time = '...
+
-
    num2str(etime(clock,ETstart)/60) ' min']); % Elapsed time from start <<<
+
-
beep;pause(.25);beep;pause(.25);beep;
+
-
</pre>
+
-
 
+
-
===[[PFV diffusion matrix | Diffusion matrix for pressure-free velocity method ('''DMatW.m''')]]===
+
-
 
+
-
===[[PFV convection matrix | Convection matrix for pressure-free velocity method ('''CMatW.m''')]]===
+
-
 
+
-
===[[PFV quartic velocity class| Modified-quartic velocity element class for pressure-free velocity method ('''ELS4424f.m''')]]===
+
-
 
+
-
===[[PFV cubic pressure class| Simple-cubic pressure element class for pressure-free velocity method ('''ELG3412f.m''')]]===
+
-
 
+
-
===[[PFV get pressure | Consistent pressure for pressure-free velocity method ('''GetPresW.m''')]]===
+
-
 
+
-
===[[PFV quadature rules| Gauss quadrature rules on rectangle ('''GQuad2.m''')]]===
+
-
 
+
-
===[[PFV mesh regrade | Grade node spacing ('''regrade.m''')]]===
+
==references==
==references==
Line 626: Line 116:
[6] {{reference-paper |author = Watkins, D. S. | year = 1976 | title = On the construction of conforming rectangular plate elements | rest = Int. J. Numer. Meth. Fluids, '''10''': 925-933}} }}
[6] {{reference-paper |author = Watkins, D. S. | year = 1976 | title = On the construction of conforming rectangular plate elements | rest = Int. J. Numer. Meth. Fluids, '''10''': 925-933}} }}
 +
 +
[7]{{reference-paper |author = Holdeman, J. T. |year = 2012 | title = A velocity-stream function method for three-dimensional incompressible fluid flow | rest = Comput. Methods Appl. Mech. Engr., '''209-212''': 66-73}}

Latest revision as of 20:07, 15 March 2013

Contents

Lid-driven cavity using pressure-free velocity formulation

This sample code uses four-node quartic Hermite finite elements with simple and Newton iteration.

Theory

The incompressible Navier-Stokes equation is a differential algebraic equation, having the inconvenient feature that there is no explicit mechanism for advancing the pressure in time. Consequently, much effort has been expended to eliminate the pressure from all or part of the computational process. We show a simple, natural way of doing this.

The incompressible Navier-Stokes equation is composite, the sum of two orthogonal equations,

\frac{\partial\mathbf{v}}{\partial t}=\Pi^S(-\mathbf{v}\cdot\nabla\mathbf{v}+\nu\nabla^2\mathbf{v})+\mathbf{f}^S ,
\rho^{-1}\nabla p=\Pi^I(-\mathbf{v}\cdot\nabla\mathbf{v}+\nu\nabla^2\mathbf{v})+\mathbf{f}^I ,

where \Pi^S and \Pi^I are solenoidal and irrotational projection operators satisfying \Pi^S+\Pi^I=1 and \mathbf{f}^S and \mathbf{f}^I are the nonconservative and conservative parts of the body force. This result follows from the Helmholtz Theorem . The first equation is a pressureless governing equation for the velocity, while the second equation for the pressure is a functional of the velocity and is related to the pressure Poisson equation. The explicit functional forms of the projection operator in 2D and 3D are found from the Helmholtz Theorem, showing that these are integro-differential equations, and not particularly convenient for numerical computation.

Equivalent weak or variational forms of the equations, proved to produce the same velocity solution as the Navier-Stokes equation are

(\mathbf{w},\frac{\partial\mathbf{v}}{\partial t})=-(\mathbf{w},\mathbf{v}\cdot\nabla\mathbf{v})-\nu(\nabla\mathbf{w}: \nabla\mathbf{v})+(\mathbf{w},\mathbf{f}^S),
(\mathbf{g}_i,\nabla p)=-(\mathbf{g}_i,\mathbf{v}\cdot\nabla\mathbf{v}_j)-\nu(\nabla\mathbf{g}_i: \nabla\mathbf{v}_j)+(\mathbf{g}_i,\mathbf{f}^I)\,,

for divergence-free test functions \mathbf{w} and irrotational test functions \mathbf{g} satisfying appropriate boundary conditions. Here, the projections are accomplished by the orthogonality of the solenoidal and irrotational function spaces. The discrete form of this is emminently suited to finite element computation of divergence-free flow.

In the discrete case, it is desirable to choose basis functions for the velocity which reflect the essential feature of incompressible flow — the velocity elements must be divergence-free. While the velocity is the variable of interest, the existence of the stream function or vector potential is necessary by the Helmholtz Theorem. Further, to determine fluid flow in the absence of a pressure gradient, one can specify the difference of stream function values across a 2D channel, or the line integral of the tangential component of the vector potential around the channel in 3D, the flow being given by Stokes' Theorem. This leads naturally to the use of Hermite stream function (in 2D) or velocity potential elements (in 3D).

Involving, as it does, both stream function and velocity degrees-of-freedom, the method might be called a velocity-stream function or stream function-velocity method.

We now restrict discussion to 2D continuous Hermite finite elements which have at least first-derivative degrees-of-freedom. With this, one can draw a large number of candidate triangular and rectangular elements from the plate-bending literature. These elements have derivatives as components of the gradient. In 2D, the gradient and curl of a scalar are clearly orthogonal, given by the expressions,

\nabla\phi = \left[\frac{\partial \phi}{\partial x},\,\frac{\partial \phi}{\partial y}\right]^T, \quad
\nabla\times\phi = \left[\frac{\partial \phi}{\partial y},\,-\frac{\partial \phi}{\partial x}\right]^T.

Adopting continuous plate-bending elements, interchanging the derivative degrees-of-freedom and changing the sign of the appropriate one gives many families of stream function elements.

Taking the curl of the scalar stream function elements gives divergence-free velocity elements [3][4]. The requirement that the stream function elements be continuous assures that the normal component of the velocity is continuous across element interfaces, all that is necessary for vanishing divergence on these interfaces.

Boundary conditions are simple to apply. The stream function is constant on no-flow surfaces, with no-slip velocity conditions on surfaces. Stream function differences across open channels determine the flow. No boundary conditions are necessary on open boundaries [3], though consistent values may be used with some problems. These are all Dirichlet conditions.

The algebraic equations to be solved are simple to set up, but of course are non-linear, requiring iteration of the linearized equations.

The finite elements we will use here are a modified form of one due to Gopapacharyulu [1][2] and Watkins [5][6]. These quartic-complete elements have 24 degrees-of-freedom, six degrees-of-freedom at each of the four nodes, and have continuous first derivatives. In the sample code the modified form of the Hermite element was obtained by interchanging first derivatives and the sign of one of them. The degrees-of-freedom of the modified element are the stream function, two components of the solenoidal velocity, and three second derivatives of the stream function. The components of the velocity element, obtained by taking the curl of the stream function element, are continuous at element interfaces. Simple-cubic Hermite elements with three degrees-of-freedom are used for the pressure, though a simple bilinear element could be used as well.

The code implementing the lid-driven cavity problem is written for Matlab. The script below is problem-specific, and calls problem-independent functions to evaluate the element diffusion and convection matricies and evaluate the pressure from the resulting velocity field. These three functions accept general convex quadrilateral elements with straight sides, as well as the rectangular elements used here. Other functions are a GMRES iterative solver using ILU preconditioning and incorporating the essential boundary conditions, and a function to produce non-uniform nodal spacing for the problem mesh.

This "educational code" is a simplified version of code used in [3]. The user interface is the code itself. The user can experiment with changing the mesh, the Reynolds number, and the number of nonlinear iterations performed, as well as the relaxation factor. There are suggestions in the code regarding near-optimum choices for this factor as a function of Reynolds number. Some values are given in the paper as well. For larger Reynolds numbers, a smaller relaxation factor speeds up convergence by smoothing the velocity in the factor \,(\mathbf{v}\cdot\nabla)\, of the convection term, but will impede convergence if made too small. A Newton iteration method can be used as well.

The output consists of graphic plots of contour levels of the stream function and the pressure levels.

This simplified version for this Wiki resulted from removal of computation of the vorticity, a restart capability, and production of publication-quality plots from one of the research codes used with the paper. It can be compared with a simpler code given in this Wiki to understand coding for higher-order elements. This code can be extended to fifth-order elements (and generalized bicubic elements), which use the same degrees-of-freedom, by simply replacing the class definitions evaluating the elements.

This modified version posted 3/15/2013 updates the code to MatLab version R2012b. The element properties are now defined in classes, and the code has been modified to accept class definitions for simple-cubic, quartic and quintic divergence-free elements on quadrilaterals and quadratic and quartic divergence-free elements on triangles.

Lid-driven cavity Matlab script

%LDC4W            LID-DRIVEN CAVITY 
% Finite element solution of the 2D Navier-Stokes equation using 
%   4-node, 6-DOF/node, C1,quartic-derived rectangular stream function 
%   elements for the Lid-Driven Cavity.
% 
% This could be characterized as a VELOCITY-STREAM FUNCTION or 
%   STREAM FUNCTION-VELOCITY method.
%
% Reference:  J. T. Holdeman, "A Hermite finite element method 
%   for incompressible fluid flow", 
%   Int. J. Numer. Meth. Fluids, 64, P376-408 (2010). 
%
% Simplified Wiki version 
%
% The rectangular problem domain is defined between Cartesian 
%   coordinates Xmin & Xmax and Ymin & Ymax.
% The computational grid has NumNx nodes in the x-direction 
%   and NumNy nodes in the y-direction. 
% The nodes and elements are numbered column-wise from the  
%   upper left corner to the lower right corner. 
%
% This script calls the user-defined functions:
%    ELS4424r     - class with definitions of velocity shape functions
%    DMatW        - to evaluate element diffusion matrix 
%    CMatW        - to evaluate element convection matrix
%    GetPres3W    - to evaluate the pressure 
%    regrade      - to regrade the mesh 
% Uses
%   ilu          - incomplete LU preconditioner
%   gmres        - iterative solver
% Indirectly uses:
%    Gquad2      - Gauss integraion rules
%    ELG3412r    - class of pressure basis functions
%
% Jonas Holdeman    Updated August 2007, revised March 2013
 ................... 

Program script for pressure-free velocity method (LDC4W.m)

Diffusion matrix for pressure-free velocity method (DMatW.m)

Convection matrix for pressure-free velocity method (CMatW.m)

Modified-quartic velocity element class for pressure-free velocity method (ELS4424f.m)

Simple-cubic pressure element class for pressure-free velocity method (ELG3412f.m)

Consistent pressure for pressure-free velocity method (GetPresW.m)

Gauss quadrature rules on rectangle (GQuad2.m)

Grade node spacing (regrade.m)

references

[1] Gopalacharyulu, S. (1973), "A higher order conforming rectangular plate element", Int. J. Numer. Meth. Fluids, 6: 305-308.

[2] Gopalacharyulu, S. (1976), "Author's reply to the discussion by Watkins", Int. J. Numer. Meth. Fluids, 10: 472-474.

[3] Holdeman, J. T. (2010), "A Hermite finite element method for incompressible fluid flow", Int. J. Numer. Meth. Fluids, 64: 376-408.

[4] Holdeman, J. T. and Kim, J.W. (2010), "Computation of incompressible thermal flows using Hermite finite elements", Comput. Methods Appl. Mech. Engr., 199: 3297-3304.

[5] Watkins, D. S. (1976), "A comment on Gopalacharyulu's 24 node element", Int. J. Numer. Meth. Fluids, 10: 471-472. }}

[6] Watkins, D. S. (1976), "On the construction of conforming rectangular plate elements", Int. J. Numer. Meth. Fluids, 10: 925-933. }}

[7]Holdeman, J. T. (2012), "A velocity-stream function method for three-dimensional incompressible fluid flow", Comput. Methods Appl. Mech. Engr., 209-212: 66-73.

My wiki