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

PFV Tconvection matrix 2

From CFD-Wiki

Jump to: navigation, search

Matlab function TCMat4424SW.m for pressure-free simple-cubic thermal convection matrix.

function [TCm,RowNdx,ColNdx]=TCMat4424SW(Xe, Elcon, nn2tft,Vdof)
% TCMat4424SW - Returns the element thermal convection matrix for
%   segregated solution using the simple-cubic Hermite basis Temperature 
%   functions on 4-node rectangular elements with 3 DOF per node using 
%   Gauss quadrature on the reference square. 
% Uses 4-node quartic velocity functions with 6 dof per node. 
% The columns of the array Vdof must contain the three degree-of-freedom 
%   vectors in the nodal order (psi,u,v). 
% The assumed nodal numbering starts with 1 at the lower left corner 
%   of the element and proceeds counter-clockwise around the element. 
%
% Usage:
%   [TCm,RowNdx,ColNdx]=TCMat4424SW(Xe, Elcon, nn2nft,Vdof)
%   Xe(1,:) -  x-coordinates of 4 corner nodes of element.  
%   Xe(2,:) -  y-coordinates of 4 corner nodes of element.  
%   Elcon(4) - connectivity matrix for this element, list of nodes. 
%   nn2tft(1,n) - global freedom number for node n.
%   nn2tft(2,n) - global freedom type for node n.
%   Vdof - (3x4) array of stream function and velocity components 
%          (psi,u,v) to specify the velocity over the element.
%
% Jonas Holdeman, January 2009,  revised July 2011


% Constants and fixed data
nn=[-1,-1; 1,-1; 1,1; -1,1]; % defines local nodal order
nnd = 4;                   % nnd = number of nodes in element
nT = 3;   nfT = nnd*nT;    % nT = number of T dofs per node, nfT = number T dofs. 
nV = 6;   nfV = nnd*nV;    % nV = number of V dofs per node, nfV = number V dofs.
NDT = 1:nT;  NDV = 1:nV; 

% Define 5-point quadrature data once, on first call. 
% Gaussian weights and absissas to integrate 9th degree polynomials exactly. 
global GQ5;
if (isempty(GQ5))   % Has 5-point quadrature data been 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;

global ZS4424tc; global Zg4424tc; global ZG4424tc; %global ZGX3412tc; global ZGY3412tc; 
if (isempty(ZS4424tc)|size(ZS4424tc,2)~=Nq)
% Evaluate and save/cache the set of shape functions at quadrature pts. 
  ZS4424tc=cell(nnd,Nq); Zg4424tc=cell(nnd,Nq); ZG4424tc=cell(nnd,Nq);  
  hh = [2,2];
   for k=1:Nq
     for m=1:nnd
       ZS4424tc{m,k} =Sr(nn(m,:),xa(k),ya(k));
       Zg4424tc{m,k}=gr(nn(m,:),xa(k),ya(k));
       ZG4424tc{m,k}=Gr(nn(m,:),xa(k),ya(k));
     end
   end
end  % if(isempty(*))
% --------------- End fixed data ----------------
 
IsAffine=1;    % Set IsAffine=1 for affine elements, IsAffine=0 for general quad

Ti=cell(nnd);  TBi=cell(nnd);
for m=1:nnd   % Loop over corner nodes 
  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;
  TBi{m} = blkdiag(1,Jt');
end  % Loop m 

% Move Jacobian evaluation inside k-loop for general convex quadrilateral. 
% Jt=[x_q, x_r; y_q, y_r];

Jt=Xe*nn(1:4,:)*.25;             % transpose of Jacobian at (0,0)
J=Jt';
Det=J(1,1)*J(2,2)-J(1,2)*J(2,1);  % Determinant of Jt & J 
Jtd=Jt/Det;
Ji=[J(2,2),-J(1,2); -J(2,1),J(1,1)]/Det;
Jd=Jtd;

% Preallocate arrays
TCm=zeros(nfT,nfT);  S=zeros(2,nfV);  g=zeros(1,nfT);  G=zeros(2,nfT); 

% Begin loop over Gauss-Legendre quadrature points. 
for k=1:Nq  
% Evaluate Ji, Det at (xa,ya) for general quadralateral 
%  Jt=Xe*GBL(n(:,:),xa,ya);               % Jacobian at (xa,ya)
%  Det=Jt(1,1)*Jt(2,2)-Jt(1,2)*Jt(2,1);   % Determinant of J 
%  Jtd=Jt/Det;
%  Ji=[Jt(2,2),-Jt(2,1); -Jt(1,2),Jt(1,1)]/Det;  % Inverse J
   
% Initialize arrays of functions and derivatives at the quadrature point (xa,ya).
mv = 0;  mt = 0;
for m=1:nnd  
  S(:,mv+NDV)= Jtd*ZS4424tc{m,k}*Ti{m};  mv = mv+nV;
  g(1,mt+NDT)= Zg4424tc{m,k}*TBi{m};
  G(:,mt+NDT)= Ji*ZG4424tc{m,k}*TBi{m};  mt = mt+nT;
end  % loop m
   
% Compute the fluid velocity at the quadrature point.
  U = S*Vdof(:);
  
  % Label rows by the test (weight) function index, columns by trial function index? 
% Submatrix ordered by T, Tx, Ty  
   TCm=TCm+(g'*(U(1)*G(1,:)+U(2)*G(2,:)))*(Det*wt(k));
end    % loop k over quadrature points 

gf=zeros(nfT,1);
m=0;
for k=1:nnd      % Loop over element nodes 
  nn=Elcon(k);              % Global node number for this element node
  gf(m+NDT)=(nn2tft(1,nn)-1)+NDT;  % Global freedoms for this node 
  m=m+nT;
end %  loop on k
RowNdx=repmat(gf,1,nfT);
ColNdx=RowNdx';
RowNdx=reshape(RowNdx,nfT*nfT,1);
ColNdx=reshape(ColNdx,nfT*nfT,1); 
TCm=reshape(TCm,nfT*nfT,1);
 
return;

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

function S=Sr(ni,q,r)
%SR  Array of quartic 4-node solenoidal basis functions.
   qi=ni(1); q0=q*ni(1); 
   ri=ni(2); r0=r*ni(2);
  % 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=gr(ni,q,r)
%GR  4-node, 12 DOF, simple cubic Hermite basis functions for temperature.
   qi=ni(1); q0=q*ni(1); 
   ri=ni(2); r0=r*ni(2);
% Scalar cubic Hermite   
g=[1/8*(1+q0)*(1+r0)*(2+q0*(1-q0)+r0*(1-r0)), -1/8/qi*(1+q0)*(1+r0)*(1-q^2), ...
       -1/8/ri*(1+q0)*(1+r0)*(1-r^2)];
return;

function G=Gr(ni,q,r)
%GR  gradient array of scalar simple-cubic basis functions.
   qi=ni(1); q0=q*ni(1);  
   ri=ni(2); r0=r*ni(2); 
   % array of irrotational vectors 
G=[ 1/8*qi*(1+r0)*(r0*(1-r0)+3*(1-q^2)), -1/8*(1+r0)*(1+q0)*(1-3*q0), ...
   -1/8*qi/ri*(1-r^2)*(1+r0); 1/8*ri*(1+q0)*(q0*(1-q0)+3*(1-r^2)), ...
   -1/8/qi*ri*(1-q^2)*(1+q0), -1/8*(1+q0)*(1+r0)*(1-3*r0)];
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;

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