Plane Elasticity: Triangular 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).
Example: Constant rigth traction
Consider the mesh below with 0.5mm of thikness and apply a constant traction Ft=1000N/mm pressure (force along the boundary) on the right boundary.
clear all;
close all;
eval('CircleHolemesh01');
[numNod,ndim]=size(nodes);
numElem=size(elem,1);
numbering=0; %=1 shows nodes and element numbering
plotElements(nodes,elem,numbering);
Select Boundary points
indRight=find(nodes(:,1)> 0.99);
indLeft=find(nodes(:,1)< -0.99);
%indCircle=find(sqrt(nodes(:,1).^2+nodes(:,2).^2)<0.501);
hold on;
plot(nodes(indRight,1),nodes(indRight,2),'mo');
plot(nodes(indLeft,1),nodes(indLeft,2),'mo');
%plot(nodes(indCircle,1),nodes(indCircle,2),'yo');
Material properties
Define material constants
E=10e+6; %Young Modulus (N/mm^2)
mu=0.25; %Poisson's ratio (adimensional)
h=0.5; % thickness (mm)
forceLoad=[1000; 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
h=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=planeElastTriangStiffMatrix(nodes,elem,e,C,h);
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];
colums= rows;
K(rows,colums)=K(rows,colums)+Ke; %assembly
end
Apply Boundary conditions
% Constant traction on Right (natural BC)
nodLoads=indRight'; %nodes where the force is applied (transposed)
Q=applyLoadsTriang(nodes,elem,nodLoads,Q,forceLoad);
% fix displacements (essential BC)
fixedNodes=[indLeft*ndim-1; indLeft*ndim];
u(fixedNodes,1)=0;
Solve the system
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,1)=um; %this way u is a column vector
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
plotPlaneNodElemDespl(nodes, elem, u, 100)
valueToShow=u(1:2:end,:);
titol='Desp. X';
plotContourSolution(nodes,elem,valueToShow,titol,'Jet');
Compute Strain and Stress per element
For triangle elements this only gives you a constant value for each element.
for e=1:numElem
v1=nodes(elem(e,1),:);
v2=nodes(elem(e,2),:);
v3=nodes(elem(e,3),:);
vertices=[v1;v2;v3];
beta=[v2(2)-v3(2),v3(2)-v1(2),v1(2)-v2(2)];
gamma=-[v2(1)-v3(1),v3(1)-v1(1),v1(1)-v2(1)];
%alpha=[v2(1)*v3(2)-v2(2)*v3(1),v3(1)*v1(2)-v3(2)*v1(1),v1(1)*v2(2)-v1(2)*v2(1)];
Area=1/2*det([ones(3,1),vertices]);
B=[beta(1), 0, beta(2), 0, beta(3), 0;
0 , gamma(1),0 , gamma(2), 0, gamma(3);
gamma(1), beta(1),gamma(2), beta(2),gamma(3), beta(3)]/(2*Area);
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];
ue=u(rows,:);
strain(e,:)=B*ue; %it is constant and must be computed for each element
stress(e,:)=C*B*ue; %it must be computed for each element
sTx=stress(e,1);
sTy=stress(e,2);
sTxy=stress(e,3);
VonMisses(e)=sqrt(sTx^2+sTy^2-sTx*sTy+3*sTxy^2);
end
plotStressVM(nodes,elem, VonMisses);
Function planeElastTriangStiffMatrix
function Ke=planeElastTriangStiffMatrix(nodes,elem,e,C,h)
% Computes the stiff matrix for a triangular element using the Plane
% Elasticity problem.
% input:
% nodes and elem as usual
% e --> number of the present element
% C --> Elasticity Matrix
% h --> Thikness
% Output:
% Ke --> Element Stiff Matrix
%
v1=nodes(elem(e,1),:);
v2=nodes(elem(e,2),:);
v3=nodes(elem(e,3),:);
vertices=[v1;v2;v3];
beta=[v2(2)-v3(2),v3(2)-v1(2),v1(2)-v2(2)];
gamma=-[v2(1)-v3(1),v3(1)-v1(1),v1(1)-v2(1)];
%alpha=[v2(1)*v3(2)-v2(2)*v3(1),v3(1)*v1(2)-v3(2)*v1(1),v1(1)*v2(2)-v1(2)*v2(1)];
Area=1/2*det([ones(3,1),vertices]);
B=[beta(1), 0, beta(2), 0, beta(3), 0;
0 , gamma(1),0 , gamma(2), 0, gamma(3);
gamma(1), beta(1),gamma(2), beta(2),gamma(3), beta(3)]/(2*Area);
Ke=(B'*C*B)*h*Area;
end
Function applyLoadsTriang
function Q=applyLoadsTriang(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);
numCov=size(nodLoads,2);
if numCov==1
error('applyLoadTriang: Not unic node allow');
end
for k=1:numElem
aux=[0,0,0]; %initial values (no node is found)
for inod=1:3 %loop for the three nodes
r=find(nodLoads==elem(k,inod)); %find if one node of element k is in the convection node list
if(~isempty(r))
aux(inod)=1; %put 1 on the local position found
end
end
number = aux(1)+2*aux(2)+4*aux(3);
switch (number) %identify the appropriate edge
case 3
ij=[1,2];
case 5
ij=[3,1];
case 6
ij=[2,3];
case 7
error('applyLoadTriang: Corners not allowed !!!!\n');
otherwise, ij=[0,0];
end
if ( ij(1) > 0) %it's an existing edge
n1=elem(k,ij(1));
n2=elem(k,ij(2));
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
(c)Numerical Factory 2019