PFV convection matrix
From CFD-Wiki
(Difference between revisions)
Jonas Holdeman (Talk | contribs)
(Created page with "Matlab function '''CMatW.m''' for pressure-free velocity convection matrix. <pre> function [Cm,RowNdx,ColNdx]=CMatW(Xe,Elcon,nn2nft,Vdof) %CMATW - Returns the affine-mapped ele...")
Newer edit →
(Created page with "Matlab function '''CMatW.m''' for pressure-free velocity convection matrix. <pre> function [Cm,RowNdx,ColNdx]=CMatW(Xe,Elcon,nn2nft,Vdof) %CMATW - Returns the affine-mapped ele...")
Newer edit →
Revision as of 03:52, 29 June 2011
Matlab function CMatW.m for pressure-free velocity convection matrix.
function [Cm,RowNdx,ColNdx]=CMatW(Xe,Elcon,nn2nft,Vdof) %CMATW - Returns the affine-mapped element convection matrix for the simple cubic Hermite % basis functions on 4-node straight-sided quadrilateral elements with 3 DOF % per node using Gauss quadrature on the reference square and row/col indices. % The columns of the array Vdof must contain the 3 nodal degree-of-freedom % vectors in the proper nodal order. % The degrees of freedom in Vdof are the stream function and the two components % u and v of the solenoidal velocity vector. % The assumed nodal numbering starts with 1 at the lower left corner of the % element and proceeds counter-clockwise around the element. % % Usage: % [CM,Rndx,Cndx] = CMatW(Xe,Elcon,nn2nft,Vdof) % Xe(1,:) - x-coordinates of corner nodes of element. % Xe(2,:) - y-coordinates of corner nodes of element. % Elcon - this element connectivity matrix % nn2nft - global number and type of DOF at each node % Vdof - (3x4) array of stream function, velocity components, and second % stream function derivatives to specify the velocity over the element. % % Jonas Holdeman, August 2007, revised June 2011 % Constants and fixed data nd = 3; nd4=4*nd; ND=1:nd; % nd = number of dofs per node, nn=[-1 -1; 1 -1; 1 1; -1 1]; % defines local nodal order % Define 5-point quadrature data once, on first call. % Gaussian weights and absissas to integrate 9th degree polynomials exactly. global GQ5; if (isempty(GQ5)) % 5-point quadrature data defined? If not, define arguments & weights. Aq=[-.906179845938664,-.538469310105683, .0, .538469310105683, .906179845938664]; Hq=[ .236926885056189, .478628670499366, .568888888888889, .478628670499366, .236926885056189]; GQ5.size=25; GQ5.xa=[Aq;Aq;Aq;Aq;Aq]; GQ5.ya=GQ5.xa'; wt=[Hq;Hq;Hq;Hq;Hq]; GQ5.wt=wt.*wt'; end xa=GQ5.xa; ya=GQ5.ya; wt=GQ5.wt; Nq=GQ5.size; % Use GQ5 (5x5) for exact integration % ----------------------------------------------- global Zs3412D2c; global ZS3412c; if (isempty(ZS3412c)|isempty(Zs3412D2c)|size(Zs3412D2c,2)~=Nq) % Evaluate and save/cache the set of shape functions at quadrature pts. Zs3412D2c=cell(4,Nq); ZS3412c=cell(4,Nq); for k=1:Nq for m=1:4 ZS3412c{m,k}= Sr(nn(m,:),xa(k),ya(k)); Zs3412D2c{m,k}=D3s(nn(m,:),xa(k),ya(k)); end end end % if(isempty(*)) % --------------- End fixed data ---------------- Ti=cell(4); % for m=1:4 Jt=Xe*GBL(nn(:,:),nn(m,1),nn(m,2)); JtiD=[Jt(2,2),-Jt(1,2); -Jt(2,1),Jt(1,1)]; % det(J)*inv(J) Ti{m}=blkdiag(1,JtiD); end % Move Jacobian evaluation inside k-loop for general convex quadrilateral. % Jt=[x_q, x_r; y_q, y_r]; Cm=zeros(nd4,nd4); Rcm=zeros(nd4,1); S=zeros(2,nd4); Sx=zeros(2,nd4); Sy=zeros(2,nd4); % Pre-allocate arrays % Begin loop over Gauss-Legendre quadrature points. for k=1:Nq Jt=Xe*GBL(nn(:,:),xa(k),ya(k)); % transpose of Jacobian at (xa,ya) Det=Jt(1,1)*Jt(2,2)-Jt(1,2)*Jt(2,1); % Determinant of Jt & J Jtd=Jt/Det; TL=[Jt(2,2)^2, -2*Jt(2,1)*Jt(2,2), Jt(2,1)^2; ... -Jt(1,2)*Jt(2,2), Jt(1,1)*Jt(2,2)+Jt(2,1)*Jt(1,2), -Jt(1,1)*Jt(2,1); ... Jt(1,2)^2, -2*Jt(1,1)*Jt(1,2), Jt(1,1)^2 ]/Det^2; % Initialize functions and derivatives at the quadrature point (xa,ya). for m=1:4 mm=nd*(m-1); S(:,mm+ND) = Jtd*ZS3412c{m,k}*Ti{m}; Ds = TL*Zs3412D2c{m,k}*Ti{m}; Sx(:,mm+ND) = [Ds(2,:); -Ds(1,:)]; % [Pyx, -Pxx] Sy(:,mm+ND) = [Ds(3,:); -Ds(2,:)]; % [Pyy, -Pxy] end % loop m % Compute the fluid velocity at the quadrature point. U = S*Vdof(:); % Submatrix ordered by psi,u,v Cm = Cm + S'*(U(1)*Sx+U(2)*Sy)*(wt(k)*Det); end % end loop k over quadrature points gf=zeros(nd4,1); m=0; for n=1:4 % Loop over element nodes gf(m+ND)=(nn2nft(1,Elcon(n))-1)+ND; % Get global freedoms m=m+nd; end RowNdx=repmat(gf,1,nd4); % Row indices ColNdx=RowNdx'; % Col indices Cm = reshape(Cm,nd4*nd4,1); RowNdx=reshape(RowNdx,nd4*nd4,1); ColNdx=reshape(ColNdx,nd4*nd4,1); return; % ---------------------------------------------------------------------------- function P2=D3s(ni,q,r) % Second derivatives [Pxx; Pxy; Pyy] of simple cubic stream function. qi=ni(1); q0=q*ni(1); q1=1+q0; ri=ni(2); r0=r*ni(2); r1=1+r0; P2=[-.75*qi^2*(r0+1)*q0, 0, -.25*qi*(r0+1)*(3*q0+1); ... .125*qi*ri*(4-3*(q^2+r^2)), .125*qi*(r0+1)*(3*r0-1), -.125*ri*(q0+1)*(3*q0-1); ... -.75*ri^2*(q0+1)*r0, .25*ri*(q0+1)*(3*r0+1), 0] ; return; function S=Sr(ni,q,r) %S Array of solenoidal basis functions on rectangle. qi=ni(1); q0=q*ni(1); q1=1+q0; ri=ni(2); r0=r*ni(2); r1=1+r0; % array of solenoidal vectors S=[ .125*ri*q1*(q0*(1-q0)+3*(1-r^2)), .125*q1*r1*(3*r0-1), .125*ri/qi*q1^2*(1-q0); ... -.125*qi*r1*(r0*(1-r0)+3*(1-q^2)), .125*qi/ri*r1^2*(1-r0), .125*q1*r1*(3*q0-1)]; return; function G=GBL(ni,q,r) % Transposed gradient (derivatives) of scalar bilinear mapping function. % The parameter ni can be a vector of coordinate pairs. G=[.25*ni(:,1).*(1+ni(:,2).*r), .25*ni(:,2).*(1+ni(:,1).*q)]; return;