CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Wiki > PFV convection matrix 2

PFV convection matrix 2

From CFD-Wiki

Jump to: navigation, search

Matlab function CMat4424W.m for pressure-free velocity convection matrix.

function [Cm,RowNdx,ColNdx,Rcm,RcNdx]=CMat4424W(Xe,Elcon,nn2nft,Vdof,ItType)
%CMAT4424G - Returns the element convection matrix for the quartic Hermite 
%   basis functions on 4-node straight-sided quadrilateral elements with 6 DOF 
%   per node using Gauss quadrature on the reference square. 
% The columns of the array Vdof must contain the six nodal degree-of-freedom 
%   vectors in the proper nodal order. 
% The degrees of freedom in Vdof are the stream function, the two components
%   u and v of the solenoidal velocity vector, and the second derivatives 
%   Pxx, Pxy, Pyy of the stream function.
% 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] = CMat4424W(Xe,Elcon,nn2nft,Vdof)
%   [CM,Rndx,Cndx,Rcm,RcNdx] = CMat4424W(Xe,Elcon,nn2nft,Vdof,ItType)
%   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  - (6x4) array of stream function, velocity components, and second 
%     stream function derivatives to specify the velocity over the element.
%   ItType - Iteration type, 0 for simple, 1 for Newton, default is simple 
%
% Jonas Holdeman, August 2007, revised July 2011 

% Constants and fixed data
nd = 6;  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 8-point quadrature data once, on first call. 
% Gaussian weights and absissas to integrate 15th degree polynomials exactly. 
global GQ8;
if (isempty(GQ8))       % Define 8-point quadrature data once, on first call. 
  Aq=[-.960289856497536,-.796666477413627,-.525532409916329,-.183434642495650, ...
       .183434642495650, .525532409916329, .796666477413627, .960289856497536]; %Abs
  Hq=[ .101228536290376, .222381034453374, .313706645877887, .362683783378362, ...
       .362683783378362, .313706645877887, .222381034453374, .101228536290376]; %Wts
   GQ8.size=64; GQ8.xa=[Aq;Aq;Aq;Aq;Aq;Aq;Aq;Aq]; GQ8.ya=GQ8.xa'; 
   Wt=[Hq;Hq;Hq;Hq;Hq;Hq;Hq;Hq]; GQ8.wt=Wt.*Wt';
end

   xa=GQ8.xa; ya=GQ8.ya; wt=GQ8.wt; Nq=GQ8.size;  % Use GQ8 (8x8) for exact integration

global Zs4424D2c;  global ZS4424c;
if (isempty(ZS4424c)|isempty(Zs4424D2c)|size(ZS4424c,2)~=Nq)
% Evaluate and save/cache the set of shape functions at quadrature pts. 
  Zs4424D2c=cell(4,Nq);  ZS4424c=cell(4,Nq); 
  for k=1:Nq
    for m=1:4
      Zs4424D2c{m,k}=D3s(nn(m,:),xa(k),ya(k));
      ZS4424c{m,k}= Sr(nn(m,:),xa(k),ya(k));
    end
  end
end  % if(isempty(*))

% --------------- End fixed data ----------------
 
if (nargin<5 | isempty(ItType) | nargout<4) ItType=0; end % default iteration type is 'simple'.  
 
Ti=cell(4);
% Jt=[x_q, x_r; y_q, y_r];
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)
  T2=[Jt(1,1)^2, 2*Jt(1,1)*Jt(1,2), Jt(1,2)^2; ...  % alt
     Jt(1,1)*Jt(2,1), Jt(1,1)*Jt(2,2)+Jt(1,2)*Jt(2,1), Jt(1,2)*Jt(2,2); ...
     Jt(2,1)^2, 2*Jt(2,1)*Jt(2,2), Jt(2,2)^2];
  T6=blkdiag(1,JtiD,T2); 
  Bxy=Xe*BLxy(nn(:,:),nn(m,1),nn(m,2));  % Second cross derivatives
  T6(5,2:3)=Bxy([2,1])';
  Ti{m}=T6;
end

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*ZS4424c{m,k}*Ti{m};
    Ds = TL*Zs4424D2c{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); 

  if (ItType~=0)   % iteration type is Newton  SOMETHING WRONG HERE !!
    Ux=Sx*Vdof(:);
    Uy=Sy*Vdof(:);
    Cm = Cm + S'*(Ux*S(1,:)+Uy*S(2,:))*(wt(k)*Det);
    Rcm=Rcm + S'*(U(1)*Ux+U(2)*Uy)*(wt(k)*Det);
  end               % Cm & Rcm for Newton iteration 
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);   
if(ItType==0) Rcm=[]; RcNdx=[]; else RcNdx=gf; end

return;

% ----------------------------------------------------------------------------

function p2=D3s(ni,q,r)
% Second derivatives [Pxx; Pxy; Pyy] of quartic stream function.
 qi=ni(1); q0=q*ni(1); q1=1+q0; 
 ri=ni(2); r0=r*ni(2); r1=1+r0;
 p2=[-3/32*qi^2*q0*(1+r0)^2*(r0*(10*q^2+3*r^2-7)+20-20*q^2-6*r^2), 3/32*qi^2/ri*q0*(1+r0)^3*(1-r0)*(5-3*r0), ...
     -3/16*qi*(1-q^2)*(1+r0)^2*(2-r0)*(1+5*q0), -1/16*(1+q0)*(1+r0)^2*(2-r0)*(1+2*q0-5*q^2), ...
     -1/8*qi/ri*(1+r0)^2*(1-r0)*(1+3*q0), -3/32*qi^2/ri^2*q0*(1+r0)^3*(1-r0)^2; ...
     9/64*qi*ri*(1-q^2)*(1-r^2)*(6-5*q^2-5*r^2), -3/64*qi*(1-q^2)*(1+r0)^2*(1-3*r0)*(7-5*r0), ...
     3/64*ri*(1+q0)^2*(1-r^2)*(1-3*q0)*(7-5*q0), 3/64*ri/qi*(1+q0)^2*(1-r^2)*(1-q0)*(1-5*q0), ...
     1/16*(1+q0)*(1+r0)*(1-3*q0)*(1-3*r0), 3/64*qi/ri*(1-q^2)*(1+r0)^2*(1-r0)*(1-5*r0); ...
     -3/32*ri^2*r0*(1+q0)^2*(q0*(10*r^2+3*q^2-7)+20-20*r^2-6*q^2), 3/16*ri*(1+q0)^2*(1-r^2)*(2-q0)*(1+5*r0), ...
     -3/32*ri^2*r0/qi*(1+q0)^2*(1-q^2)*(5-3*q0), -3/32*ri^2/qi^2*r0*(1+q0)*(1-q^2)^2, ...
     -1/8*ri/qi*(1+q0)*(1-q^2)*(1+3*r0), -1/16*(1+q0)^2*(1+r0)*(2-q0)*(1+2*r0-5*r^2)] ;  
return;

function S=Sr(ni,q,r)
%SR  Array of solenoidal basis functions.
   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=[3/64*ri*(1+q0)^2*(1-r^2)*((1+q0)*(8-9*q0+3*q^2)+(2-q0)*(1-5*r^2)), ... 
   -1/64*(1+q0)^2*(1+r0)^2*(2-q0)*(7-26*r0+15*r^2), ...
    3/64*ri/qi*(1+q0)^3*(1-r^2)*(1-q0)*(5-3*q0), ...
    3/64*ri/qi^2*(1+q0)^3*(1-r^2)*(1-q0)^2, ...
    1/16/qi*(1+q0)^2*(1+r0)*(1-q0)*(1-3*r0), ... 
    1/64/ri*(1+q0)^2*(1+r0)^2*(1-r0)*(1-5*r0)*(2-q0); ... 
   -3/64*qi*(1+r0)^2*(1-q^2)*((1+r0)*(8-9*r0+3*r^2)+(2-r0)*(1-5*q^2)), ... 
    3/64*qi/ri*(1+r0)^3*(1-q^2)*(1-r0)*(5-3*r0), ... 
   -1/64*(1+r0)^2*(1+q0)^2*(2-r0)*(7-26*q0+15*q^2), ... 
   -1/64/qi*(1+r0)^2*(1+q0)^2*(1-q0)*(1-5*q0)*(2-r0), ... 
   -1/16/ri*(1+r0)^2*(1+q0)*(1-r0)*(1-3*q0), ... 
   -3/64*qi/ri^2*(1+r0)^3*(1-q^2)*(1-r0)^2];
return;

function G=GBL(ni,q,r)
% Transposed gradient (derivatives) of scalar bilinear mapping function for straight-sised elements. 
% 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;

function D=BLxy(ni,q,r)
% Transposed second cross-derivative of scalar bilinear mapping function. 
% The parameter ni can be a vector of coordinate pairs. 
  D=[.25*ni(:,1).*ni(:,2)]; 
return;
My wiki