ElasticitySolverBiLinearTraction

Table of contents

2D Static Linear Elasticity - bi-linear elements
Solution strategy
Preliminaries
Material model
Compute the stiffness matrix
Traction load
Essential boundary conditions
Verification - Displacements
Strain and Stress computation - Post Processing

2D Static Linear Elasticity - bi-linear elements

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.

Solution strategy

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

Preliminaries

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

Material model

E = 2200; %MPa
nu = 0.33;
t = 0.2; %mm

% Plane Stress
lambda = E*nu/(1-nu^2);
mu = E/(2*(1+nu));

Compute the stiffness matrix

[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

Traction load

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)

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');

Essential boundary conditions

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)])

Verification - Displacements

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)

uy = U([3 4 9 10 11],2)

uy = 5x1    
   -0.0750
   -0.0750
   -0.0750
   -0.0750
   -0.0750

Strain and Stress computation - Post Processing

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