Table of contents
Solve the static linear elasticity PDE using bi-linear elements and an external traction force according to
\[ \begin{cases} -\nabla\cdot\boldsymbol{\sigma}=\boldsymbol{f} & \text{in }\Omega\\ \boldsymbol{\sigma}=2\mu\boldsymbol{\varepsilon}+\lambda\text{tr}\boldsymbol{\varepsilon}\boldsymbol{I} & \text{in }\Omega\\ \boldsymbol{\sigma}\cdot\boldsymbol{n}=\boldsymbol{F} & \text{at }x=3\\ u_{x}=0 & \text{at }x=0\\ u_{y}=0 & \text{at }y=0 \end{cases} \]Let \(E=2200\textrm{MPa}\), \(\nu =0\ldotp 33\), \(t=0\ldotp 2\textrm{mm}\), use plane stress and implement a solver to numerically compute the displacements and stress.
Create a simple geometry that is verifiable by hand. In this case we will create a rectangle on which we can apply the load case shown in the figure above. We want to apply boundary conditions such that the rectangle is in pure tension, meaning we only have stress in the x-direction. All other stress i zero. Thus we can compute the displacement analytically:
\[\delta =\dfrac{F L}{E A}\;\textrm{and}\;\sigma =\dfrac{F}{A}\]where \(F=100N\), \(A=h t\;{\textrm{mm}}^2\) and \(L=3\textrm{mm}\), \(h=1\ldotp 2\textrm{mm}\).
The geometry and mesh is given by:
clear
P = [0, 0
3, 0
3, 1.2
0, 1.2
0.7 0
1.6 0
2.5 0
3 0.5
2.3 1.2
1.5 1.2
0.6 1.2
0 0.5
0.55 0.4
1.45 0.6
2.4 0.45];
nodes = [1 5 13 12
12 13 11 4
5 6 14 13
13 14 10 11
14 15 9 10
6 7 15 14
7 2 8 15
15 8 3 9];
figure; hold on
patch('Vertices', P, 'Faces', nodes,'FaceColor', 'c', ...
'EdgeAlpha',0.3)
axis equal tight
plot(P(:,1),P(:,2),'ok','MarkerFaceColor','k','MarkerSize',14)
text(P(:,1)-0.05,P(:,2),cellstr(num2str([1:size(P,1)]')), ...
'Color','w','FontSize',9)
for iel = 1:size(nodes,1)
inod = nodes(iel,:);
xm = mean(P(inod,1));
ym = mean(P(inod,2));
text(xm,ym,num2str(iel),'FontSize',9)
end
Weak form: Find the displacement field \(\bm u\), with \(\bm u = \bm g\) on the boundary \(\partial \Omega_D\), such that
\[\int_{\Omega } \bm \sigma \left(\bm u\right) : \bm \varepsilon \left(\bm v\right)\;d\Omega =\int_{\Omega } \bm f\cdot \bm v\;d\Omega +\int_{\partial \Omega } \bm F \cdot \bm v\;d\;s\;\;\;\forall v\;\textrm{that}\;\textrm{are}\;0\;\textrm{on}\;\partial \Omega_D \;\;\;\;\;\left(1\right)\]We have (using the Mandel notation), the element stiffness matrix
\[S_K =\int_K \left(2\mu \mathbf B_{\varepsilon }^T \mathbf B_{\varepsilon } +\lambda \mathbf B_D^T \mathbf B_D \right)d K\]the element load vector (zero in our problem)
\[\bm f_K =\int_K \Phi^T \mathbf f \left(\bm x\right)\;d K\]the traction (external force) vector
\[\bm g_E =\int_E \Phi_E^T \bm t \left(\bm x\right)\;d E\]where \(E\) denotes edge (in 2D, where as in 3D it would be a surface). The edge element is one spatial dimension lower, we are formulating the equations on an edge instead of a triangle. For the edge element we have
\[ \def\arraystretch{1.5} \Phi_E \overset{2D\;\textrm{lin}}{=} \left\lbrack \begin{array}{cccc} \varphi_1^E & 0 & \varphi_2^E & 0\\ 0 & \varphi_1^E & 0 & \varphi_2^E \end{array}\right\rbrack , \bm \varphi^E \overset{2D\;\textrm{lin}}{=} \left\lbrack \begin{array}{cc} 1+\xi & \xi \end{array}\right\rbrack \] \[ \def\arraystretch{2.5} \bm B_{\varepsilon } :=\left\lbrack \begin{array}{ccccc} \dfrac{\partial \varphi_1 }{\partial x} & 0 & \dfrac{\partial \varphi_2 }{\partial x} & 0 & \cdots \\ 0 & \dfrac{\partial \varphi_1 }{\partial y} & 0 & \dfrac{\partial \varphi_2 }{\partial y} & \cdots \\ \dfrac{1}{\sqrt{2}}\dfrac{\partial \varphi_1 }{\partial y} & \dfrac{1}{\sqrt{2}}\dfrac{\partial \varphi_1 }{\partial x} & \dfrac{1}{\sqrt{2}}\dfrac{\partial \varphi_2 }{\partial y} & \dfrac{1}{\sqrt{2}}\dfrac{\partial \varphi_2 }{\partial x} & \cdots \end{array}\right\rbrack \] \[ \bm B_{\textrm{div}} :=\left\lbrack \begin{array}{ccccc} \dfrac{\partial \varphi_{\;1} }{\partial x} & \dfrac{\partial \varphi_{\;1} }{\partial y} & \dfrac{\partial \varphi_{\;2} }{\partial x} & \dfrac{\partial \varphi_{\;2} }{\partial y} & \cdots \end{array}\right\rbrack \] \[ \def\arraystretch{1.5} \bm u_K (\bm x)\approx U_K (\bm x)={\underset{\Phi_K (\bm x)}{\underbrace{\left\lbrack \begin{array}{cccccc} \varphi_1 (\bm x) & 0 & \varphi_2 (\bm x) & 0 & \varphi_3 (\bm x) & 0\\ 0 & \varphi_1 (\bm x) & 0 & \varphi_2 (\bm x) & 0 & \varphi_3 (\bm x) \end{array}\right\rbrack } } }_{2\textrm{x6}} {\left\lbrack \begin{array}{c} u_x^1 \\ u_y^1 \\ u_x^2 \\ u_y^2 \\ u_x^3 \\ u_y^3 \end{array}\right\rbrack }_{6\textrm{x1}} \] \[ \bm \varphi =\left\lbrack \begin{array}{c} {\left(\eta -1\right)}\,{\left(\xi -1\right)}\\ -\xi \,{\left(\eta -1\right)}\\ \eta \,\xi \\ -\eta \,{\left(\xi -1\right)} \end{array}\right\rbrack \] \[ \hat{\bm B} =\left\lbrack \begin{array}{cccc} \eta -1 & 1-\eta & \eta & -\eta \\ \xi -1 & -\xi & \xi & 1-\xi \end{array}\right\rbrack \]phi = @(xi,eta)[(eta-1.0).*(xi-1.0),-xi.*(eta-1.0),eta.*xi,-eta.*(xi-1.0)];
Bh = @(xi,eta)[eta-1, 1-eta, eta, -eta;
xi-1, -xi, xi, 1-xi ];
\[
\def\arraystretch{1.5}
\mathit{\mathbf{J}}\left(\xi ,\eta \right)=\left\lbrack \begin{array}{ccc}
\dfrac{\partial \varphi_1 }{\partial \xi } & \cdots & \dfrac{\partial \varphi_4 }{\partial \xi }\\
\dfrac{\partial \varphi_1 }{\partial \eta } & \cdots & \dfrac{\partial \varphi_4 }{\partial \eta }
\end{array}\right\rbrack \left\lbrack \begin{array}{cc}
x_1 & y_1 \\
\vdots & \vdots \\
x_4 & y_4
\end{array}\right\rbrack = \hat{\mathbf B} (\xi ,\eta )\; \mathbf P_e
\]
\[
\def\arraystretch{1.5}
\left\lbrack \begin{array}{c}
\dfrac{\partial \varphi_i }{\partial x}\\
\dfrac{\partial \varphi_i }{\partial y}
\end{array}\right\rbrack = \mathbf J^{-1} \hat{\bf B}
\]
Tasks:
Compute the stiffness matrix
Compute the traction load
Solve the linear system \(\mathit{\mathbf{S}}\;\mathit{\mathbf{u}}=\mathit{\mathbf{f}}\)
Compare displacement error
Compute the stress
E = 2200; %MPa
nu = 0.33;
t = 0.2; %mm
% Plane Stress
lambda = E*nu/(1-nu^2);
mu = E/(2*(1+nu));
[nele, knod] = size(nodes);
[nnod, dof] = size(P);
neq = nnod*dof;
[GP,GW] = GaussQuadrilateral(2);
% Pre-allocate
S = zeros(neq,neq);
sq = 1/sqrt(2);
for iel = 1:nele
% iel
inod = nodes(iel,:);
iP = P(inod,:);
locdof = zeros(1,8);
locdof(1:2:end) = inod*2-1;
locdof(2:2:end) = inod*2;
% figure
% patch(iP(:,1),iP(:,2),'w')
% hold on; axis equal
% text(mean(iP(:,1)), mean(iP(:,2)), sprintf('E_%d',iel), 'BackgroundColor','y')
Se = zeros(8,8);
Beps = zeros(3,8);
Bdiv = zeros(1,8);
% area = 0;
for iG = 1:length(GW)
xi = GP(iG,1); eta = GP(iG,2);
wi = GW(iG);
J = Bh(xi,eta)*iP;
% det(J);
% area = area + det(J)*wi;
B = J\Bh(xi,eta);
Beps(1,1:2:end) = B(1,:);
Beps(2,2:2:end) = B(2,:);
Beps(3,1:2:end) = sq*B(2,:);
Beps(3,2:2:end) = sq*B(1,:);
Bdiv(1:2:end) = B(1,:);
Bdiv(2:2:end) = B(2,:);
Se = Se + t*(2*mu*Beps'*Beps + lambda*Bdiv'*Bdiv)*det(J)*wi;
% xy = phi(xi,eta)*iP;
% text(xy(1),xy(2),sprintf('G_%d',iG))
end
% area
% polyarea(iP(:,1), iP(:,2))
S(locdof, locdof) = S(locdof, locdof) + Se;
end
dy = norm(P(2,:)-P(3,:))
dy = 1.2000
traction = [100/(dy*t),0]'; % traction load N/mm^2
phi=@(xi)[1-xi,xi];
edges = [2 8
8 3];
ieqs = zeros(2*2,1);
PHI = zeros(2,4);
F = zeros(neq,1);
L = 0;
for iedge = 1:size(edges,1)
inod = edges(iedge,:);
ieqs(1:2:end) = 2*inod-1; %x-dofs
ieqs(2:2:end) = 2*inod-0; %y-dofs
xc = P(inod,:);
h = norm(xc(end,:)-xc(1,:));
L = L + h;
xi = 0.5; iw = 1;
PHI(1,1:2:end) = phi(xi);
PHI(2,2:2:end) = phi(xi);
fe = PHI'*traction(:)*t*h*iw;
F(ieqs) = F(ieqs) + fe;
end
% F is returnd as N
sum(F)
ans = 100
LoadNodes = [2,3,8];
F(LoadNodes*2-1)
ans = 3x1
20.8333
29.1667
50.0000
Visualize the load vector
figure; hold on
patch('Vertices', P, 'Faces', nodes,'FaceColor', 'c', ...
'EdgeAlpha',0.1,'DisplayName','mesh')
axis equal tight
quiver(P(LoadNodes,1), P(LoadNodes,2), F(LoadNodes*2-1), F(LoadNodes*2-0),'b', ...
'DisplayName','F=[100,0]N');
u = zeros(neq,1);
Y-Displacements
% Prescribe displacements
YNodes = [1 5 6 7 2];
u(YNodes*2) = 0; %y-displacements
plot(P(YNodes,1),P(YNodes,2),'bd','MarkerFaceColor','b','Displayname','u_y=0')
X-Displacements
XNodes = [1 4 12];
% prescribed displacement
u(XNodes*2-1) = 0; %x-displacements
plot(P(XNodes,1),P(XNodes,2),'rs','MarkerFaceColor','r','Displayname','u_x=0')
hl = legend('show');
set(hl,'Position',[0.60 0.69 0.20 0.20]);
Modify the right-hand side
\[\mathbf{S}\;\mathbf{u}=\mathbf{f}-{\mathbf{S}}_{\mathbf{p}} \;{\mathbf{u}}_{\mathbf{p}}\]presc = unique([XNodes*2-1, YNodes*2]);
free = setdiff(1:neq,presc);
F = F - S(:,presc)*u(presc);
Solve the linear system
u(free) = S(free,free)\F(free);
Visualize the displacement field
U = [u(1:2:end), u(2:2:end)];
UR = sum(U.^2,2).^0.5;
scale = 1;
figure; hold on
patch('Vertices', P, 'Faces', nodes,'FaceColor', 'w','EdgeAlpha',0.5)
patch('Vertices', P+U*scale, 'Faces', nodes,'FaceColor', 'interp', ...
'EdgeAlpha',0.1,'CData',UR)
axis equal tight
title(['Displacement field, scale: ',num2str(scale)])
dx = P(2,1)-P(1,1) %mm
dy = P(3,2)-P(2,2) %mm
A = t*dy %mm^2 Cross sectional area
f = traction(1)*A %N Total external load
deltax = f*dx/(E*A) %mm Total x deformation
epsx = deltax/dx %mm/mm Strain in x-dir
epsy = -nu*(epsx) %mm/mm Strain in y-dir
deltay = epsy*dy %mm Total y deformation
sigma = f/A %N/mm^2 or MPa stress
dx = 3
dy = 1.2000
A = 0.2400
f = 100
deltax = 0.5682
epsx = 0.1894
epsy = -0.0625
deltay = -0.0750
sigma = 416.6667
ux = U(LoadNodes,1)
ux = 3x1
0.5682
0.5682
0.5682
uy = U([3 4 9 10 11],2)
uy = 5x1
-0.0750
-0.0750
-0.0750
-0.0750
-0.0750
figure; hold on
patch('Vertices', P, 'Faces', nodes,'FaceColor', 'w','EdgeAlpha',0.1)
axis equal tight
title('von Mises stress field')
phi = @(xi,eta)[(eta-1.0).*(xi-1.0),-xi.*(eta-1.0),eta.*xi,-eta.*(xi-1.0)];
for iel = 1:nele
inod = nodes(iel,:);
iP = P(inod,:);
locdof = zeros(1,8);
locdof(1:2:end) = inod*2-1;
locdof(2:2:end) = inod*2;
iU = U(inod,:);
dd = iP+iU*scale;
Svm = zeros(4,1);
area = 0;
for iG = 1:length(GW)
xi = GP(iG,1); eta = GP(iG,2); iw = GW(iG);
J = Bh(xi,eta)*iP;
B = J\Bh(xi,eta);
gradU = B*iU;
area = area + det(J)*iw;
epsilon = (gradU'+gradU)/2;
sig = 2*mu*epsilon+lambda*trace(epsilon)*eye(2);
sig_vm = sqrt(sig(1,1)^2-sig(1,1)*sig(2,2)+sig(2,2)^2+3*sig(1,2)^2);
[V,sigmaPrinc] = eig(sig, 'vector');
[~, ind] = sort(abs(sigmaPrinc),'descend');
dir = V(:,ind);
sigmaPrinc = sigmaPrinc(ind);
xy = phi(xi,eta)*dd;
hold on; qscale = 0.1*(sigmaPrinc/norm(sigmaPrinc));
if qscale(1) > 0
quiver(xy(1),xy(2),dir(1,1),dir(2,1),qscale(1),'b','Displayname','1st principal direction')
quiver(xy(1),xy(2),-dir(1,1),-dir(2,1),qscale(1),'b','HandleVisibility','off')
end
if qscale(2) > 0
quiver(xy(1),xy(2),dir(1,2),dir(2,2),qscale(2),'r','Displayname','2nd principal direction')
quiver(xy(1),xy(2),-dir(1,2),-dir(2,2),qscale(2),'r','HandleVisibility','off')
end
Svm(iG) = sig_vm;
end
patch('Faces',1:4,'Vertices',dd,'FaceColor','flat','CData',mean(Svm),'EdgeColor','k','FaceAlpha',0.3)
end
colorbar
function [GP,GW] = GaussQuadrilateral(order)
%GaussQuadrilateral Gauss integration scheme for quadrilateral 2D element
% [x,y] in [0,1]
[gp,gw] = gauss(order,0,1);
[GPx,GPy] = meshgrid(gp,gp);
[Gw1,Gw2] = meshgrid(gw,gw);
GW = [Gw1(:), Gw2(:)];
GW = prod(GW,2);
GP = [GPx(:), GPy(:)];
end
function [x, w, A] = gauss(n, a, b)
%------------------------------------------------------------------------------
% gauss.m
%------------------------------------------------------------------------------
%
% Purpose:
%
% Generates abscissas and weigths on I = [ a, b ] (for Gaussian quadrature).
%
%
% Syntax:
%
% [x, w, A] = gauss(n, a, b);
%
%
% Input:
%
% n integer Number of quadrature points.
% a real Left endpoint of interval.
% b real Right endpoint of interval.
%
%
% Output:
%
% x real Quadrature points.
% w real Weigths.
% A real Area of domain.
%------------------------------------------------------------------------------
% 3-term recurrence coefficients:
n = 1:(n - 1);
beta = 1 ./ sqrt(4 - 1 ./ (n .* n));
% Jacobi matrix:
J = diag(beta, 1) + diag(beta, -1);
% Eigenvalue decomposition:
%
% e-values are used for abscissas, whereas e-vectors determine weights.
%
[V, D] = eig(J);
x = diag(D);
% Size of domain:
A = b - a;
% Change of interval:
%
% Initally we have I0 = [ -1, 1 ]; now one gets I = [ a, b ] instead.
%
% The abscissas are Legendre points.
%
if ~(a == -1 && b == 1)
x = 0.5 * A * x + 0.5 * (b + a);
end
% Weigths:
w = V(1, :) .* V(1, :);
end