Plane Elasticity: Triangular elements

Plane Elasticity has two different approaches

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 (th=1).
Plane Stress problems: when the thickness is small compared to the section area. It is like study a 2D domain with thickness (th=thickness).

For Linear Triangular Elements we know that the formulas are explicit and they can be stated as K=h*A*B'*C*B

Problem:

Consider a unique triangular element defined in the figure. The two nodes in the base are fixed and edge 2 is pulled by a traction of 20N/cm^2.

The material associated to the triangle has Young Modulus = 30·10^6 N/cm^2, Poisson's Ratio = 0.3 and a thicness th=0.23 cm

Geometry

define nodes and elements

v1=[0,0];

v2=[10,0];

v3=[5,10];

nodes=[v1;v2;v3];

elem=[1,2,3];

numNod=size(nodes,1);

ndim=size(nodes,2);

numElem=size(elem,1);

plotElements(nodes,elem,1);

Material properties

Define material constants

E=30e+6; %Young Modulus (N/cm^2)
mu=0.3; %Poisson's ratio (adimensional)
th=0.23;  % thickness (cm)
t=20; %traction force N/cm^2

Define Plane Elasticity Problem

Plane Stress problem == 1

Plane Strain problem == 2 (thikness no significative)

modelProblem = 1; %Plane Stress
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
%
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)*th*Area;
Fe=zeros(numNod,1); 
K=Ke;

Boundary conditions

edgeVector=(v3-v2)/norm(v3-v2);

normal=[edgeVector(2),-edgeVector(1)];

f=t*normal; %constant traction applied on normal vector direction

%plot traction vector (origin at the edge)

ndiv=10; %number of arrows in the plot

vini=v2;

vfin=v3;

scale=t/100; %scale arrows according to t value

plotEdgeConstantBC(vini,vfin,t,ndiv,scale) %

Q=norm(v3-v2)/2*[0;0;f(1);f(2);f(1);f(2)];

% fix displacements (essential BC)

fixedNodes=[];

nNod=1;

fixedNodes=[fixedNodes, nNod*ndim-1]; %(u1_x=0)

fixedNodes=[fixedNodes, nNod*ndim]; %(u1_y=0)

nNod=2;

fixedNodes=[fixedNodes, nNod*ndim-1]; %(u1_x=0)

fixedNodes=[fixedNodes, nNod*ndim]; %(u1_y=0)

freeNodes = setdiff(1:2*numNod,fixedNodes); %Complementary of fixed nodes

% 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

u(fixedNodes,1)=0;

displacements=[u(1:2:end),u(2:2:end)] %nodes disp x, y

displacements = 3×2

10^-4 ×

         0         0
         0         0
    0.1733    0.0303

% Post Process

strain=B*u %it must be computed for each element

strain = 3×1

10^-5 ×

         0
    0.0303
    0.1733

stress=C*strain %it must be computed for each element

stress = 3×1

0000
0000
0000

sTx=stress(1);

sTy=stress(2);

sTxy=stress(3);

VonMisses=sqrt(sTx^2+sTy^2-sTx*sTy+3*sTxy^2)

VonMisses = 35.7631

plotPlaneNodElemDespl(nodes, elem, u, 1.e+4)

% For triangle elements this only gives you a constant value for each element.

Exercise 1:

Build the function

function Ke=planeElastTriangStiffMatrix(nodes,elem,e,C,th)

it returns the Stiffness Matrix 6x6 for each element

(c)Numerical Factory 2020