Plane Elasticity: Quadrangular Meshes
Plane Elasticity has two different approaches
- Plane Stress problems: when the thickness is small compared to the section area. It is like study a 2D domain with thickness (h=thickness).
- Plane Strain problems: when the thickness is very large compared to the section area. It is like to study only the 2D section of a 3D domain (h=1).
At the end of this file you will find the new functions needed for quadrilaterals.
Example: Constant rigth traction
Consider the mesh below with 0.05mm of thikness and apply a constant traction Ft=5000N/mm pressure (force along the boundary) on the right boundary. (see the rest material values in the code)
clear;
close all;
eval('meshPlacaForatQuad');
[numNod,ndim]=size(nodes);
numElem=size(elem,1)
numElem = 433
numbering=0;
plotElements(nodes,elem,numbering);
Select Boundary points
indTop=find(nodes(:,2)> 1.99);
indBot=find(nodes(:,2)< 0.01);
indCircle=find(sqrt((nodes(:,1)-1).^2+(nodes(:,2)-1).^2)<0.501);
indRight=find(nodes(:,1)> 4.95);
hold on;
plot(nodes(indCircle,1),nodes(indCircle,2),'yd');
plot(nodes(indRight,1),nodes(indRight,2),'mo');
hold off;
Material properties
Define material constants
E=10e+6; %Young Modulus
mu=0.3; %Poisson's ratio (adimensional)
th=0.05; % thickness
forceLoad=[5000; 0]; %(Fx,Fy) traction force N/mm
Define Plane Elasticity Problem
Plane Stress problem == 1
Plane Strain problem == 2 (thikness no significative)
modelProblem = 1; %Plane Stress (1) or Plane Strain (2)
if (modelProblem == 1)
c11=E/(1-mu^2);
c22=c11;
c12=mu*c11;
c21=c12;
c33=E/(2*(1+mu));
else
th=1;
c11=E*(1-mu)/((1+mu)*(1-2*mu));
c22=c11;
c12=c11*mu/(1-mu);
c21=c12;
c33=E/(2*(1+mu));
end
C=[c11,c12,0; c21,c22,0; 0,0,c33];
%
% compute the shape coefficients
%
K=zeros(numNod*ndim);
F=zeros(numNod*ndim,1);
Q=zeros(numNod*ndim,1);
for e=1:numElem
Ke=planeElastQuadStiffMatrix(nodes,elem,e,C,th);
rows=[elem(e,1)*2-1; elem(e,1)*2; elem(e,2)*2-1; elem(e,2)*2;...
elem(e,3)*2-1; elem(e,3)*2; elem(e,4)*2-1; elem(e,4)*2];
colums= rows;
K(rows,colums)=K(rows,colums)+Ke; %assembly
%Fe=internalConstantForcesQuad(nodes,elem,e,th,internalForce); %if an internal constant force is defined
%F(rows,1)=F(rows,1)+Fe; %assembly
end
Boundary conditions
Constant traction on Right
nodLoads=indRight'; %nodes where the force is applied
Q=applyLoadsQuad(nodes,elem,nodLoads,Q,forceLoad);
% fix displacements (essential BC)
%
u=zeros(numNod*ndim,1); %this way u is a column vector
fixedNodes=[indCircle*ndim-1 indCircle*ndim]; % x and y coordinates
u(fixedNodes)=0;
freeNodes = setdiff(1:2*numNod,fixedNodes); %Complementary of fixed nodes
% Not modify the linear system, this is only valid if BC = 0.
Km=K(freeNodes,freeNodes);
Qm=Q(freeNodes);
% solve the system
um=Km\Qm;
u(freeNodes)=um;
displacements=[u(1:2:end),u(2:2:end)] %nodes disp x, y
% We show only displacements for the last 10 nodes
%displacements(end-9:end,:)
Post Process step
plotPlaneNodElemDespl(nodes, elem, u, 7);
valueToShow=abs(u(1:2:end,:)); % x displacement in each node
titol='Desp. X';
plotContourSolution(nodes,elem,valueToShow,titol,'Jet');
Compute Strain and Stress at each Node
[stress,VonMisses]=computeQuadStrainStressVM(nodes, elem, u, C);
titol='VonMisses Stress';
plotContourSolution(nodes,elem,VonMisses,titol,'Jet')
pause(2) %you can comment this line
applyLoadsQuad function
Apply constant loads on a boundary
function [Q]=applyLoadsQuad(nodes,elem,nodLoads,Q,forceLoad)
%
% (c) Numerical Factory 2019
%
%------------------------------------------------------------------------
% Apply a load BC on a boundary of a quadrilateral meshed domain.
% The behabiour is very similar to the applyConvecQuad funtion used to
% apply convection BC.
%-------------------------------------------------------------------------
numElem=size(elem,1);
[numNod,ndim]=size(nodes);
numCo=size(nodLoads,2);
if (numCo==1)
error('applyLoadQuad: It can''t manage a single node');
end
for k=1:numElem
aux=[0,0,0,0];
for j=1:4
r=find(nodLoads==elem(k,j)); %find if j node of element k is in the boundary nodes list
if(~isempty(r)), aux(j)=1; end
end
%See the Triang version for the explanation of the binary number
number=aux*[1;2;4;8];
% Here we codify both edges and corners on the boundaries
switch (number)
case 3, ij=[1,2,0,0]; % edge 1: aux=[1,1,0,0];
case 6, ij=[2,3,0,0]; % edge 2: aux=[0,1,1,0];
case 7, ij=[1,2,2,3]; % corner 2: aux=[1,1,1,0];
case 9, ij=[4,1,0,0]; % edge 4
case 11, ij=[4,1,1,2]; % corner 1
case 12, ij=[3,4,0,0]; % edge 3
case 13, ij=[3,4,4,1]; % corner 4
case 14, ij=[2,3,3,4]; % corner 3
case 15, error('applyLoadQuad: It can''t manage two corners'); % two corners
otherwise, ij=[0,0,0,0];
end
for j=1:2:3
if ~ij(j)==0
n1=elem(k,ij(j));
n2=elem(k,ij(j+1));
L=norm(nodes(n1,:)-nodes(n2,:));
posit=[n1*ndim-1, n1*ndim];
Q(posit)=Q(posit)+0.5*L*forceLoad; %same (x,y) force on both nodes
posit=[n2*ndim-1, n2*ndim];
Q(posit)=Q(posit)+0.5*L*forceLoad; %same (x,y) force on both nodes
end
end
end
end
planeElastQuadStiffMatrix function
Compute the Quad stiff matrix
function Ke=planeElastQuadStiffMatrix(nodes,elem,e,C,th)
N=2; %number of GaussPoints
[w,ptGaussRef]=gaussValues2DQuad(N);
%
% First compute Ke, Fe for each element
%
v1=nodes(elem(e,1),:);
v2=nodes(elem(e,2),:);
v3=nodes(elem(e,3),:);
v4=nodes(elem(e,4),:);
vertices=[v1;v2;v3;v4];
% Shape functions
Psi1=@(x,y)(1-x).*(1-y)/4;
Psi2=@(x,y)(1+x).*(1-y)/4;
Psi3=@(x,y)(1+x).*(1+y)/4;
Psi4=@(x,y)(1-x).*(1+y)/4;
shapeFunctions = @(x,y)[Psi1(x,y),Psi2(x,y),Psi3(x,y),Psi4(x,y)];
% Shape function derivatives
dPsi11=@(x,y) -(1-y)/4;
dPsi21=@(x,y) (1-y)/4;
dPsi31=@(x,y) (1+y)/4;
dPsi41=@(x,y) -(1+y)/4;
dPsi12=@(x,y) -(1-x)/4;
dPsi22=@(x,y) -(1+x)/4;
dPsi32=@(x,y) (1+x)/4;
dPsi42=@(x,y) (1-x)/4;
% Derivative matrix 2x4
Jacob =@(x,y) [dPsi11(x,y), dPsi21(x,y),dPsi31(x,y),dPsi41(x,y);...
dPsi12(x,y), dPsi22(x,y),dPsi32(x,y),dPsi42(x,y)];
% Compute the corresponding Gaussian points on the domain
% evaluate Shape functions on Gaussian reference points
xx = ptGaussRef(:,1);
yy = ptGaussRef(:,2);
% evaluate Jacobian contribution for each point
% We use a Matlab cell variable in order to load each matrix
numPtG=size(xx,1);
for i=1:numPtG
Jtilde{i}=inv(Jacob(xx(i),yy(i))*vertices);
evalDetJacob(i) = det(Jacob(xx(i),yy(i))*vertices);
Jacobia{i}=Jacob(xx(i),yy(i)); %derivatives of the shape functions
shapeFun{i}=shapeFunctions(xx(i),yy(i));
end
%
% element stiff matrix
%
ik1=1; ik2=1; %index of the derivatives here x-x
Kxx=compKijIntegQuad(ik1,ik2,Jacobia,w,evalDetJacob,Jtilde,numPtG);
ik1=1; ik2=2;
Kxy=compKijIntegQuad(ik1,ik2,Jacobia,w,evalDetJacob,Jtilde,numPtG);
ik1=2; ik2=1;
Kyx=compKijIntegQuad(ik1,ik2,Jacobia,w,evalDetJacob,Jtilde,numPtG);
ik1=2; ik2=2;
Kyy=compKijIntegQuad(ik1,ik2,Jacobia,w,evalDetJacob,Jtilde,numPtG);
KK={};
for i=1:4
for j=1:4
KK{i,j}=[C(1,1)*Kxx(i,j)+C(3,3)*Kyy(i,j), C(1,2)*Kxy(i,j)+C(3,3)*Kyx(i,j);
C(3,3)*Kxy(i,j)+C(2,1)*Kyx(i,j), C(3,3)*Kxx(i,j)+C(2,2)*Kyy(i,j)];
end
end
Ke =th*[KK{1,1}, KK{1,2}, KK{1,3}, KK{1,4};
KK{2,1}, KK{2,2}, KK{2,3}, KK{2,4};
KK{3,1}, KK{3,2}, KK{3,3}, KK{3,4};
KK{4,1}, KK{4,2}, KK{4,3}, KK{4,4}];
end
compKijIntegQuad function
function Kij=compKijIntegQuad(ik1,ik2,Jacobia,w,evalDetJacob,Jtilde,numPtG)
Kij=zeros(4);
for i=1:4
for j=1:4
suma=0;
for ptG=1:numPtG
first=Jacobia{ptG}(1,i)*Jtilde{ptG}(ik1,1)+Jacobia{ptG}(2,i)*Jtilde{ptG}(ik1,2);
second=Jacobia{ptG}(1,j)*Jtilde{ptG}(ik2,1)+Jacobia{ptG}(2,j)*Jtilde{ptG}(ik2,2);
suma=suma+w(ptG)*first*second*evalDetJacob(ptG);
end
Kij(i,j)=suma;
end
end
end
(c)Numerical Factory 2020