PFVT thermal program script
From CFD-Wiki
Revision as of 20:48, 15 March 2013 by Jonas Holdeman (Talk | contribs)
%TC44SW THERMAL-DRIVEN CAVITY
%
% Finite element solution of the 2D Navier-Stokes equation for the
% bouyancy-driven thermal cavity problem using quartic Hermite elements
% and SEGMENTED SOLUTIONS (T-V split).
%
% This could be characterized as a VELOCITY-STREAM FUNCTION or
% STREAM FUNCTION-VELOCITY method.
%
% Reference: "Computation of incompressible thermal flows using Hermite finite
% elements", Comput. Methods Appl. Mech. Engrg, 199, P3297-3304 (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.
% Segregated version
%
%This script direcly calls the user-defined functions:
% ELS4424r - class with definitions of velocity shape functions
% ELG3412r - class with definitions of temperature & pressure fns
% regrade - to regrade the mesh (optional)
% DMatW - to evaluate element velocity diffusion matrix
% CMatW - to evaluate element velocity convection matrix
% TDMatSW - to evaluate element thermal diffusion matrix
% TCMatSW - to evaluate element thermal convection matrix
% BMatSW - to evaluate element bouyancy matrix
% GetPres3W - to evaluate the pressure
%
% Uses:
% ilu - incomplete LU preconditioner for solver
% gmres - sparse linear solver to solve the sparse system
% Indirectly uses:
% Gquad2 - Gauss integraion rules
% ELG3412r - class of pressure basis functions
%
% Jonas Holdeman, December 2008 revised March 2013
clear all;
eu = ELS4424r; % class of functions for velocity
et = ELG3412r; % class fo functions for emperaure
disp('Bouyancy-driven thermal cavity, T-V split');
disp([' Four-node, 24 dof ' eu.name 'basis, 12 dof cubic thermal basis.']);
% ------------------ Fixed constants --------------------------
nnd = eu.nnodes; % nnodes = number of element nodes
nV = eu.nndofs; % nndofs = number of velocity dofs per node, (6);
nT = et.nndofs; % nndofs = number of temperature dofs per node, (3);
nn = eu.nn; % defines local nodal order, [-1 -1; 1 -1; 1 1; -1 1]
nV2 = nV*nV; % Number velocity DOFs per node squared.
nT2 = nT*nT; % Number temperature DOFs per node squared.
nD = nV+nT; nD2=nD*nD; % Total number of DOFs per node
% -------------------------------------------------------------
% Tolerance parameters for the GMRES iterative sparse solver
GM.VTol = 2.e-12; % convergence tolerance for velociy solutions
GM.TTol = 1.e-13; % convegence tolerance for temperature solutions
GM.MIter = 30; % number of iterations allowed
GM.MRstrt = 6; % number of restarts
% parameters for ilu preconditioner
% Decrease su.droptol if ilu preconditioner fails
su.type='ilutp';
su.droptol=1.e-5;
% -------------------------------------------------------------
% Fixed parameters:
% Temperature on left, right side
TL=1; TR=0;
% Prandtl number
Pr=.71;
% Dimensionless equation form to be used
% Use EquationForm=1 for large Ra, EquationForm=0 for Ra->0
EquationForm=1; % (Choose 0 or 1)
ETstart=clock;
% Define the problem geometry:
Xmin = 0.0; Xmax = 1.0; Ymin = 0.0; Ymax = 1.0;
% Mesh grading parameters
xgrd = 1.0; ygrd=1.0; %
%xgrd = .75; ygrd=.75; %
% Set " RefineBoundary=1 " for additional refinement at boundary,
% i.e., split first element along each boundary into two.
RefineBoundary=0;
% DEFINE THE MESH:
% Number of elements in the x-direction
NumEx = 18;
% Number of elements in the y-direction
NumEy = 18;
NumEl = NumEx*NumEy;
% Number of nodes
NumNx = NumEx+1; NumNy = NumEy+1;
NumNod = NumNx*NumNy;
% --------------------------------------
% Suggested relaxation factors for steady flow
% Ra 10^3 10^4 10^5 10^6
% RelxFac .94 .67 .35 .04
% Rayleigh number Ra
Ra = 10^4;
% Relaxation factor ( <= 1 )
RelxFac=.67;
% Number of nonlinear iterations
MaxNLit=20;
if (xgrd~=1 || ygrd~=1), meshsp='graded'; else meshsp='uniform'; end
disp(['Rayleigh number = ' num2str(Ra) ', ' num2str(NumEx) 'x' ...
num2str(NumEy) ' element ' meshsp ' mesh' ]);
disp(['Maximum number of nonlinear iterations = ' num2str(MaxNLit)]);
disp(' ');
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 towards the ends
% if agrd<1 then refine 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),.5*(YNc(1)+YNc(2)),YNc(2:end-1),(YNc(end-1)+YNc(end))/2.,YNc(end)];
else
XNc=linspace(Xmin,Xmax,NumNx);
YNc=linspace(Ymax,Ymin,NumNy);
end
if xgrd ~= 1, XNc=regrade(XNc,xgrd,0); end; % Refine mesh if desired
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:NumNy-1)=.5*(YNc(1:NumNy-2)-YNc(3:NumNy));
wy(NumNy)=.5*(YNc(NumNy-1)-YNc(NumNy));
Wa=wy'*wx;
% Initial stream function and velocity, temperature & gradient 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);
T0=zeros(NumNy,NumNx);
Tx0=zeros(NumNy,NumNx);
Ty0=zeros(NumNy,NumNx);
%------------------Begin grid plot--------------------
% Plot the grid
figure(1);
clf;
orient portrait; orient tall;
subplot(2,2,2);
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 thermal cavity']);
%-------------------End grid plot---------------------
MaxVdof=nV*NumNod;
MaxTdof=nT*NumNod;
MaxDof=MaxVdof+MaxTdof; % maximum number of degrees of freedom
VBC.Mxdof=MaxVdof; % maximum number V of degrees of freedom
TBC.Mxdof=MaxTdof; % maximum number of T degrees of freedom
NodXY =zeros(NumNod,2);
nn2vft=zeros(NumNod,2); % node number -> nfv & nt
nn2tft=zeros(NumNod,2); % node number -> nft & nt
NodNdx=zeros(NumNod,2);
% Generate lists of active nodal indices, freedom number & type
ni=0; nt=1;
nfv=-nV+1; nft=-nT+1; % ________
for nx=1:NumNx % | |
for ny=1:NumNy % | |
ni=ni+1; % |________|
NodNdx(ni,:)=[nx,ny];
NodXY(ni,:)=[Xgrid(ny,nx),Ygrid(ny,nx)];
nfv=nfv+nV; % all V nodes have 6 dofs
nn2vft(ni,:)=[nfv,nt]; % dof number & type (all nodes type 1)
nft=nft+nT; % all T nodes have 3 dofs
nn2tft(ni,:)=[nft,nt]; % dof number & type (all nodes type 1)
end;
end;
%NumNod=ni; % total number of nodes
nvf2nnt=zeros(MaxVdof,2); % (node, type) associated with dof
ntf2nnt=zeros(MaxTdof,2); % (node, type) associated with dof
nfv=0; nft=0;
for n=1:NumNod
for k=1:nV
nvf2nnt(nfv+k,:)=[n,k];
end
nfv=nfv+nV;
for k=1:nT
ntf2nnt(nft+k,:)=[n,k];
end
nft=nft+nT;
end
% Generate element connectivity, from upper left to lower right.
Elcon=zeros(NumEl,4);
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 (Dirichlet) V boundary conditions
MaxVBC=2*(5*NumNy+5*(NumNx-2)); % Allocate space for VBC
VBC.dof=zeros(MaxVBC,1); % Degree-of-freedom index
VBC.val=zeros(MaxVBC,1); % Dof value
nc=0;
for nf=1:MaxVdof
nt=nvf2nnt(nf,2); ni=nvf2nnt(nf,1);
nx=NodNdx(ni,1); ny=NodNdx(ni,2);
x=XNc(nx); y=YNc(ny);
if((x==Xmin || x==Xmax) && (y==Ymax || y==Ymin)) % corner
nc=nc+1; VBC.dof(nc)=nf; VBC.val(nc)=0; % psi, u, v, pxx, pxy, pyy
elseif(x==Xmin || x==Xmax) % left or right boundary
switch nt;
case {1, 2, 3, 5, 6} % psi, u, v, pxy, pyy (pxx=-v_x free)
nc=nc+1; VBC.dof(nc)=nf; VBC.val(nc)=0;
end % switch (type)
elseif (y==Ymax || y==Ymin) % top, bottom
switch nt
case {1, 2, 3, 4, 5} % psi, u, v, pxx, pxy (pyy=u_y free)
nc=nc+1; VBC.dof(nc)=nf; VBC.val(nc)=0;
end % switch (type)
end % if (boundary)
end % for nf
VBC.num=nc;
if (size(VBC.dof,1)>nc) % Truncate V arrays if necessary
VBC.dof=VBC.dof(1:nc);
VBC.val=VBC.val(1:nc);
end % truncate V
% End ESSENTIAL (Dirichlet) V boundary conditions
% Begin ESSENTIAL (Dirichlet) T boundary conditions
MaxTBC=2*(2*NumNy+(NumNx-2)); % Allocate space for TBC
TBC.dof=zeros(MaxTBC,1); % Degree-of-freedom index
TBC.val=zeros(MaxTBC,1); % Dof value
nc=0;
for nf=1:MaxTdof
nt=ntf2nnt(nf,2); ni=ntf2nnt(nf,1);
nx=NodNdx(ni,1); ny=NodNdx(ni,2);
x=XNc(nx); y=YNc(ny);
if(x==Xmin) % left boundary
switch nt;
case 1
nc=nc+1; TBC.dof(nc)=nf; TBC.val(nc)=TL; % T
case 3
nc=nc+1; TBC.dof(nc)=nf; TBC.val(nc)=0; % Ty
end % switch (type)
elseif (x==Xmax) % right boundary
switch nt;
case 1
nc=nc+1; TBC.dof(nc)=nf; TBC.val(nc)=TR; % T
case 3
nc=nc+1; TBC.dof(nc)=nf; TBC.val(nc)=0; % Ty
end % switch (type)
elseif (y==Ymax || y==Ymin) % top, bottom
switch nt
case 3
nc=nc+1; TBC.dof(nc)=nf; TBC.val(nc)=0; % Ty
end % switch (type)
end % if (boundary)
end % for nf
TBC.num=nc;
if (size(TBC.dof,1)>nc) % Truncate T arrays if necessary
TBC.dof=TBC.dof(1:nc);
TBC.val=TBC.val(1:nc);
end % truncate T
% End ESSENTIAL (Dirichlet) T boundary conditions
% Number active DOFs
ADOFs=MaxDof-VBC.num-TBC.num;
disp(['Number of active DOFs = ' num2str(ADOFs)]);
% partion out essential (Dirichlet) dofs
p_vec = (1:VBC.Mxdof)'; % List of all dofs
VBC.p_vec_undo = zeros(1,VBC.Mxdof);
% form a list of non-diri dofs
VBC.ndro = p_vec(~ismember(p_vec, VBC.dof)); % list of non-diri dofs
% calculate p_vec_undo to restore Q to the original dof ordering
VBC.p_vec_undo([VBC.ndro;VBC.dof]) = (1:VBC.Mxdof); %p_vec';
p_vec = (1:TBC.Mxdof)'; % List of all dofs
TBC.p_vec_undo = zeros(1,TBC.Mxdof);
% form a list of non-diri dofs
TBC.ndro = p_vec(~ismember(p_vec, TBC.dof)); % list of non-diri dofs
% calculate p_vec_undo to restore Q to the original dof ordering
TBC.p_vec_undo([TBC.ndro;TBC.dof]) = (1:TBC.Mxdof); %p_vec';
Qv=zeros(MaxVdof,1); % Allocate space for velocity solution (dof) vector
Qt=zeros(MaxTdof,1); % Allocate space for temperature solution (dof) vector
for k = 1:MaxTdof
nn=ntf2nnt(k,1);
nx=NodNdx(nn,1); ny=NodNdx(nn,2);
switch ntf2nnt(k,2) % switch on DOF type
case 1
Qt(k)=T0(ny,nx);
case 2
Qt(k)=Tx0(ny,nx);
case 3
Qt(k)=Ty0(ny,nx);
end; % switch (type)
end % loop over TDOFs
% Initialize fields to boundary conditions
KB=1:VBC.num;
Qv(VBC.dof(KB))=VBC.val(KB);
KB=1:TBC.num;
Qt(TBC.dof(KB))=TBC.val(KB);
errpsi=zeros(NumNy,NumNx); % error correct for iteration
MxNL=max(1,MaxNLit);
np0=zeros(1,MxNL); % Arrays for convergence info
nv0=zeros(1,MxNL);
nT0=zeros(1,MxNL);
if EquationForm==1
Ctd=1/sqrt(Ra*Pr); Cvd=Ctd*Pr; Cbf=1; % Coefficients for large Ra
else
Ctd=1; Cvd=Pr; Cbf=Pr*Ra; % Coefficients for small Ra
end
DMat = spalloc(MaxVdof,MaxVdof,54*MaxVdof); % for the diffusion matrix
TDMat = spalloc(MaxTdof,MaxTdof,30*MaxTdof); % for thermal diffusion mat.
BMat = spalloc(MaxVdof,MaxTdof,18*MaxVdof); % for the bouyancy matrix
MatV = []; % Velocity submatrix
MatT = []; % Temperature submatrix
Vdof=zeros(6,nnd);
Xe=zeros(2,nnd); % coordinates of element corners
NLitr=0;
NV=1:nV;
ItType=0; % Initialy simple iteration
while NLitr<MaxNLit, NLitr=NLitr+1; % <<< BEGIN NONLINEAR ITERATION
% Generate and assemble element matrices
CMat = spalloc(MaxVdof,MaxVdof,54*MaxVdof); % for fluid convection matrix
TCMat = spalloc(MaxTdof,MaxTdof,30*MaxTdof); % for thermal convection mat.
RHS=spalloc(MaxVdof,1,MaxVdof);
% Copy fields to dof vector
Qv0=Qv; Qt0=Qt;
tclock=clock; % Start assembly time <<<<<<<<<
for ne=1:NumEl % BEGIN GLOBAL MATRIX ASSEMBLY
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,:),nn2vft);
DMat = DMat + sparse(Rndx,Cndx,Emat,MaxVdof,MaxVdof);
% Element thermal diffusion matrix
[Emat,Rndx,Cndx] = TDMatSW(et,Xe,Elcon(ne,:),nn2tft);
TDMat = TDMat + sparse(Rndx,Cndx,Emat,MaxTdof,MaxTdof);
% Global fluid bouyancy matrix
[Emat,Rndx,Cndx] = BMatSW(eu,Xe,Elcon(ne,:),nn2vft,nn2tft);
BMat = BMat + sparse(Rndx,Cndx,Emat,MaxVdof,MaxTdof);
end
if (NLitr>1)
% First iteration uses Stokes equation.
% Get stream function and velocities for linearized Navier-Stokes
for n=1:nnd % Loop over local nodes of element
Vdof(NV,n)=Qv((nn2vft(Elcon(ne,n),1)-1)+NV);
end
% Fluid element convection matrix
[Emat,Rndx,Cndx] = CMatW(eu,Xe,Elcon(ne,:),nn2vft,Vdof);
CMat = CMat + sparse(Rndx,Cndx,Emat,MaxVdof,MaxVdof);
% Thermal convection matrix
[Emat,Rndx,Cndx] = TCMatSW(eu,Xe,Elcon(ne,:),nn2tft,Vdof);
TCMat = TCMat + sparse(Rndx,Cndx,Emat,MaxTdof,MaxTdof);
end; % NLitr>1
end; % END GLOBAL MATRIX ASSEMBLY
MatT = -TCMat-Ctd*TDMat;
MatV = -CMat -Cvd*DMat;
disp(['(' num2str(NLitr) ') Matrix assembly complete, elapsed time = '...
num2str(etime(clock,tclock)) ' sec']); % Assembly time <<<<<<<<<<<
pause(1);
% Solve system
tclock=clock; % Start solution time <<<<<<<<<<<<<<
RHSt=sparse(MaxTdof,1);
RHSt=RHSt(TBC.ndro)-MatT(TBC.ndro,TBC.dof)*TBC.val;
MatT=MatT(TBC.ndro,TBC.ndro);
Qs=Qt(TBC.ndro);
% Solve for temperatures
[Lm,Um] = ilu(MatT,su); % incomplete LU
Qr = gmres(MatT,RHSt,GM.MIter,GM.TTol,GM.MRstrt,Lm,Um,Qs); % GMRES
Qt=[Qr;TBC.val]; % Augment active dofs with esential (Dirichlet) dofs
Qt=Qt(TBC.p_vec_undo); % Restore natural order
RHS=RHS-Cbf*BMat*Qt;
RHS=RHS(VBC.ndro)-MatV(VBC.ndro,VBC.dof)*VBC.val;
MatV=MatV(VBC.ndro,VBC.ndro);
Qs=Qv(VBC.ndro);
% Solve for velocities
[Lm,Um] = ilu(MatV,su); % incomplete LU
Qr = gmres(MatV,RHS,GM.MIter,GM.VTol,GM.MRstrt,Lm,Um,Qs); % GMRES
Qv=[Qr;VBC.val]; % Augment active dofs with esential (Dirichlet) dofs
Qv=Qv(VBC.p_vec_undo); % Restore natural order
% --------------------------------------------
Qv=RelxFac*Qv+(1-RelxFac)*Qv0;
% --------------------------------------------
% Compute change and copy dofs to field arrays
dsqp=0; dsqv=0; dsqt=0;
for k=1:MaxVdof
ni=nvf2nnt(k,1);
nx=NodNdx(ni,1); ny=NodNdx(ni,2);
switch nvf2nnt(k,2) % switch on dof type
case 1
dsqp=dsqp+Wa(ny,nx)*(Qv(k)-Qv0(k))^2; psi0(ny,nx)=Qv(k);
errpsi(ny,nx)=Qv0(k)-Qv(k);
case 2
dsqv=dsqv+Wa(ny,nx)*(Qv(k)-Qv0(k))^2; u0(ny,nx)=Qv(k);
case 3
dsqv=dsqv+Wa(ny,nx)*(Qv(k)-Qv0(k))^2; v0(ny,nx)=Qv(k);
case 4
pxx0(ny,nx)=Qv(k);
case 5
pxy0(ny,nx)=Qv(k);
case 6
pyy0(ny,nx)=Qv(k);
end % switch on dof type
end % for
for k=1:MaxTdof
ni=ntf2nnt(k,1);
nx=NodNdx(ni,1); ny=NodNdx(ni,2);
switch ntf2nnt(k,2) % switch on dof type
case 1
dsqt=dsqt+Wa(ny,nx)*(Qt(k)-Qt0(k))^2; T0(ny,nx)=Qt(k);
case 2
Tx0(ny,nx)=Qt(k);
case 3
Ty0(ny,nx)=Qt(k);
end % switch on dof type
end % for
np0(NLitr)=sqrt(dsqp);
nv0(NLitr)=sqrt(dsqv);
nT0(NLitr)=sqrt(dsqt);
disp(['Solution time for linear system = '...
num2str(etime(clock,tclock)) ' sec, nv0=' num2str(nv0(NLitr)) ...
', nT0=' num2str(nT0(NLitr))]);
if (np0(NLitr)<=1e-15||nv0(NLitr)<=1e-15||nT0(NLitr)<=1e-15)
MaxNLit=NLitr; np0=np0(1:MaxNLit); nv0=nv0(1:MaxNLit); nT0=nT0(1:MaxNLit); end;
%---------- Begin plot of intermediate results ----------
figure(2);
orient portrait;
orient tall;
% Stream function (intermediate)
subplot(2,2,1);
contour(Xgrid,Ygrid,psi0,10,'k'); % Plot contours (trajectories)
hold on;
plot([Xmin,Xmin,Xmax,Xmax,Xmin],[Ymax,Ymin,Ymin,Ymax,Ymax],'k');
hold off;
axis([Xmin,Xmax,Ymin,Ymax]);
axis equal; axis image;
title(['Thermal cavity stream lines; Ra=' num2str(Ra)]);
% Temperature (intermediate)
subplot(2,2,3);
contour(Xgrid,Ygrid,T0,10,'k');
axis([Xmin,Xmax,Ymin,Ymax]);
hold on;
plot([Xmin,Xmin,Xmax,Xmax,Xmin],[Ymax,Ymin,Ymin,Ymax,Ymax],'k');
hold off;
axis equal; axis image;
title(['Thermal cavity isotherms, Ra=' num2str(Ra)]);
% Convergence
subplot(2,2,2);
semilogy(1:MaxNLit,nT0,'k-x',1:MaxNLit,nv0,'k-+',1:MaxNLit,np0,'k-o');
xlabel('Nonlinear iteration number');
ylabel('Nonlinear correction');
axis square;
title(['Iteration conv., Ra=' num2str(Ra)]);
legend('T','U','Psi');
% Error
subplot(2,2,4);
axis([Xmin,Xmax,Ymin,Ymax]);
contour(Xgrid,Ygrid,errpsi,8,'k'); % Plot contours (errors)
hold on;
pcolor(Xgrid,Ygrid,errpsi);
shading interp %flat;
plot([Xmin,Xmin,Xmax,Xmax,Xmin],[Ymax,Ymin,Ymin,Ymax,Ymax],'k');
hold off;
axis([Xmin,Xmax,Ymin,Ymax]);
axis image;
title('Iteration correction, psi');
pause(1);
%---------- End plot of intermediate results ---------
end; % <<< END NONLINEAR ITERATION
format short g;
disp('Convergence results by iteration: temperature, velocity, stream function');
disp(['nT0: ' num2str(nT0)]); 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
% Cubic pressure
[P,Px,Py] = GetPres3W(eu,NodXY,Elcon,nn2vft,Qv,EBCp,Cvd);
P=reshape(P,NumNy,NumNx);
Px = reshape(Px,NumNy,NumNx);
Py = reshape(Py,NumNy,NumNx);
%----------Plot of final result-------------
figure(1);
% Stream function (final)
subplot(2,2,1);
contour(Xgrid,Ygrid,psi0,10,'k'); % Plot contours (trajectories)
hold on;
pcolor(Xgrid,Ygrid,psi0);
shading interp %flat;
plot([Xmin,Xmin,Xmax,Xmax,Xmin],[Ymax,Ymin,Ymin,Ymax,Ymax],'k');
hold off;
axis([Xmin,Xmax,Ymin,Ymax]);
axis equal; axis image;
title(['Stream lines, Ra=' num2str(Ra)]);
% Temperature (final)
subplot(2,2,3);
contour(Xgrid,Ygrid,T0,10,'k');
hold on;
pcolor(Xgrid,Ygrid,T0);
shading interp %flat;
plot([Xmin,Xmin,Xmax,Xmax,Xmin],[Ymax,Ymin,Ymin,Ymax,Ymax],'k');
hold off;
axis([Xmin,Xmax,Ymin,Ymax]);
axis equal; axis image;
title(['Isotherms, Ra=' num2str(Ra)]);...
% Plot pressure contours (final)
subplot(2,2,4);
contour(Xgrid,Ygrid,P,10,'k'); % Plot pressure contours
hold on;
pcolor(Xgrid,Ygrid,P);
shading interp %flat;
plot([Xmin,Xmin,Xmax,Xmax,Xmin],[Ymax,Ymin,Ymin,Ymax,Ymax],'k');
hold off;
axis([Xmin,Xmax,Ymin,Ymax]);
axis equal; axis image;
title(['Pressure contours, Ra=' num2str(Ra)]);
% ************* FIGURE 2 ***********************
figure (2);
clf;
orient tall;
% U-velocity (final)
subplot(2,2,1);
contour(Xgrid,Ygrid,u0,10,'k'); % Plot vector field
hold on;
pcolor(Xgrid,Ygrid,u0);
shading interp %flat;
plot([Xmin,Xmin,Xmax,Xmax,Xmin],[Ymax,Ymin,Ymin,Ymax,Ymax],'k');
hold off;
axis([Xmin,Xmax,Ymin,Ymax]);
axis equal; axis image;
title('U-velocity field');...
% V-velocity (final)
subplot(2,2,3);
contour(Xgrid,Ygrid,v0,10,'k'); % Plot vector field
hold on;
pcolor(Xgrid,Ygrid,v0);
shading interp %flat;
plot([Xmin,Xmin,Xmax,Xmax,Xmin],[Ymax,Ymin,Ymin,Ymax,Ymax],'k');
hold off;
axis([Xmin,Xmax,Ymin,Ymax]);
axis equal; axis image;
title('V-velocity field');...
% Convergence
subplot(2,2,2);
semilogy(1:MaxNLit,nT0,'k-x',1:MaxNLit,nv0,'k-+',1:MaxNLit,np0,'k-o');
xlabel('Nonlinear iteration number');
ylabel('Nonlinear correction');
axis square;
title(['Iteration conv., Ra=' num2str(Ra)]);
legend('T','U','Psi');
subplot(2,2,4);
nc=fix((NumNy+1)/2);
plot(Xgrid(nc,:),T0(nc,:),'k-x');
title('Centerline temperature');...
% ---------- End final plot ------------
disp(['Total elapsed time = '...
num2str(etime(clock,ETstart)/60) ' min']); % Elapsed time from start <<
beep; pause(.5); beep; pause(.5); beep;
