0

I am trying to write some MATLAB code to solve fluid flow problems using the finite element method. I am starting with the relatively-simple steady Stokes problem. I am using the quadrilateral P2P1 element, and following Reddy and Gartling 2010. As far as I understand, the viscous matrix (K(1:18,1:18) in my code below) should be symmetric positive definite for the linear Stokes problem, but mine is symmetric but not positive definite! I checked every single line and cannot still figure out why it's not working!

Could anyone please take a look at the code and please tell me if anything is wrong!

I would definitely appreciate any additional advice/resources, thank you.

clc; clear;

syms r s mu rho x1 x2 x3 x4 x5 x6 x7 x8 x9 y1 y2 y3 y4 y5 y6 y7 y8 y9


% Velocity Shape functions (Reddy and Gartling P65):
% h > ψ, r > ξ, s > η, H > Ψ:
h1 = ( (1-r)(1-s)(-r-s-1)+(1-r^2)(1-s^2) )/4;
h2 = ( (1+r)(1-s)(+r-s-1)+(1-r^2)(1-s^2) )/4;
h3 = ( (1+r)(1+s)(+r+s-1)+(1-r^2)(1-s^2) )/4;
h4 = ( (1-r)(1+s)(-r+s-1)+(1-r^2)(1-s^2) )/4;
h5 = ( 2(1-r^2)(1-s^1)-(1-r^2)(1-s^2) )/4;
h6 = ( 2(1+r^1)(1-s^2)-(1-r^2)(1-s^2) )/4;
h7 = ( 2(1-r^2)(1+s^1)-(1-r^2)(1-s^2) )/4;
h8 = ( 2(1-r^1)(1-s^2)-(1-r^2)(1-s^2) )/4;
h9 = ( 4(1-r^2)(1-s^2) )/4;
H =  [h1; h2; h3; h4; h5; h6; h7; h8; h9];


Hr = diff(H,r); Hs = diff(H,s);


% Coordinate Interpolation (Local to Natural Coor.) (Reddy and Gartling P165)
x_n = [x1; x2; x3; x4; x5; x6; x7; x8; x9];
y_n = [y1; y2; y3; y4; y5; y6; y7; y8; y9];


% The Jacobian operator and its determinant (Reddy and Gartling P69)
J = [Hr.'; Hs.'][x_n y_n];
iJ = J^-1; dJ = det(J);
% Derivatives of H Shape Functions:
Hxy = iJ[Hr.'; Hs.']; Hx = transpose(Hxy(1,:)); Hy = transpose(Hxy(2,:));


% Pressure Basis Functions (Reddy and Gartling P65):
HP = [(1-r)(1-s)/(4); (1+r)(1-s)/(4); (1+r)(1+s)/(4); (1-r)(1+s)/(4)];


% Building the FE Equations (Reddy and Gartling P166):
K11 = mu(HxHx.')dJ;
K22 = mu(HyHy.')dJ;
K12 = muHxHy.'dJ;
K21 = muHyHx.'dJ;
Q1  = HxHP.'dJ;
Q2  = HyHP.'dJ;
Zp  = zeros(max(size(HP)), max(size(HP)));


K = [2K11+K22,  K12,       -Q1; ...
     K21,        K11+2K22, -Q2; ...
    -Q1.',      -Q2.',       Zp];
Mohamed Abdelhamid
  • 105

1 Answers1

0

I realized the shape functions h5, h6, h7, and h8 in the textbook were wrong. They are missing "2" in the second term. They should be as follows:

h5 = (2*(1-r^2)*(1-s^1)-2*(1-r^2)*(1-s^2) )/4;
h6 = (2*(1+r^1)*(1-s^2)-2*(1-r^2)*(1-s^2) )/4;
h7 = (2*(1-r^2)*(1+s^1)-2*(1-r^2)*(1-s^2) )/4;
h8 = (2*(1-r^1)*(1-s^2)-2*(1-r^2)*(1-s^2) )/4;

It makes me sad that such a highly praised book has such a silly error.

Mohamed Abdelhamid
  • 105