Nonlinear finite elements/Newton method for axial bar

From Wikiversity
Jump to navigation Jump to search

Newton's method applied to the axial bar problem[edit]

Let us now return to the axial bar problem and see how the Newton-Raphson scheme can be used to solve the system of equations. The reduced system of equations is

Therefore, the residual is

The terms of the tangent stiffness matrix are given by

For our reduced stiffness matrix, we have

Now, from the stiffness matrix, we see that

The derivatives are,

Therefore, the components of are

In matrix form,

The inverse of (the reduced form shown in equation (8)) is

==== Maple code ====11 A Maple text script for doing these computations symbolically is show below.

> #
> # Setting up the global stiffness matrix and the global
> # tangent stiffness matrix for a two element mesh
> #
> restart;
> with(linalg):
> #
> # Set up mesh
> # 
> nElem := 2;
> xNode := linalg[matrix](3,1,[x_1,x_2,x_3]);
> elem  := linalg[matrix](2,2,[1,2,2,3]);
> xDisp := linalg[matrix](3,1,[u_1,u_2,u_3]);
> #
> # Set up the global stiffness matrix and load vector
> #
> K  := linalg[matrix](3,3,[0,0,0,0,0,0,0,0,0]):
> F  := linalg[matrix](3,1,[0,0,0]):
> T  := linalg[matrix](3,3,[0,0,0,0,0,0,0,0,0]):
> #
> # Compute the global element stiffness, tangent stiffness, and load
> #
> for ii from 1 to nElem do
>   #
>   # Compute the element shape functions and their derivatives
>   #
>   node1 := elem[ii,1]:
>   node2 := elem[ii,2]:
>   x1 := xNode[node1,1]:
>   x2 := xNode[node2,1]:
>   N1 := (x2-x)/(x2-x1):
>   N2 := (x-x1)/(x2-x1):
>   dN1_dx := diff(N1,x):
>   dN2_dx := diff(N2,x):
>   #
>   # The trial function and E
>   #
>   u1 := xDisp[node1,1]:
>   u2 := xDisp[node2,1]:
>   u := u1*N1 + u2*N2:
>   E := E0*(1 - diff(u,x)):
>   #
>   # The stiffness matrix terms
>   #
>   k11 := A*E*dN1_dx*dN1_dx:
>   k12 := A*E*dN1_dx*dN2_dx:
>   k21 := A*E*dN2_dx*dN1_dx:
>   k22 := A*E*dN2_dx*dN2_dx:
>   KK11 := int(k11, x=x1..x2):
>   KK12 := int(k12, x=x1..x2):
>   KK21 := int(k21, x=x1..x2):
>   KK22 := int(k22, x=x1..x2):
>   K11 := simplify(KK11);
>   K12 := simplify(KK12);
>   K21 := simplify(KK21);
>   K22 := simplify(KK22);
>   Ke  := linalg[matrix](2,2,[K11,K12,K21,K22]);
>   #
>   # The load vector terms
>   #
>   f1 := x*N1:
>   f2 := x*N2:
>   F1 := a*int(f1, x=x1..x2):
>   F2 := a*int(f2, x=x1..x2):
>   Fe  := linalg[matrix](2,1,[F1,F2]);
>   #
>   # The derivatives of the stiffness matrix
>   #
>   dK11_u1 := diff(K11,u1):
>   dK11_u2 := diff(K11,u2):
>   dK12_u1 := diff(K12,u1):
>   dK12_u2 := diff(K12,u2):
>   dK21_u1 := diff(K21,u1):
>   dK21_u2 := diff(K21,u2):
>   dK22_u1 := diff(K22,u1):
>   dK22_u2 := diff(K22,u2):
>   #
>   # The tangent stiffness matrix
>   #
>   TT11 := K11 + dK11_u1*u1 + dK12_u1*u2:
>   TT12 := K12 + dK11_u2*u1 + dK12_u2*u2:
>   TT21 := K21 + dK21_u1*u1 + dK22_u1*u2:
>   TT22 := K22 + dK21_u2*u1 + dK22_u2*u2:
>   T11 := simplify(TT11):
>   T12 := simplify(TT12):
>   T21 := simplify(TT21):
>   T22 := simplify(TT22):
>   Te  := linalg[matrix](2,2,[T11,T12,T21,T22]);
>   #
>   # Assemble the global stiffness matrix
>   #
>   jlocal := 0:
>   for jj from node1 to node2 do
>     jlocal := jlocal + 1:
>     klocal := 0:
>     for kk from node1 to node2 do
>       klocal := klocal + 1:
>       K[jj,kk] := K[jj,kk] + Ke[jlocal,klocal]:
>       T[jj,kk] := T[jj,kk] + Te[jlocal,klocal]:
>     end do;
>     F[jj,1] := F[jj,1] + Fe[jlocal,1]:
>   end do:
> end do:
> Ksim  := linalg[matrix](3,3,[0,0,0,0,0,0,0,0,0]):
> Fsim  := linalg[matrix](3,1,[0,0,0]):
> Tsim  := linalg[matrix](3,3,[0,0,0,0,0,0,0,0,0]):
> for ii from 1 to nElem+1 do
>   for jj from 1 to nElem+1 do
>     Ksim[ii,jj] := algsubs(x_2 - x_3 = -h, algsubs(x_1 - x_2 = -h, K[ii,jj])):
>     Tsim[ii,jj] := algsubs(x_2 - x_3 = -h, algsubs(x_1 - x_2 = -h, T[ii,jj])):
>   end do:
>   Fsim[ii,1] := algsubs(x_2 - x_3 = -h, algsubs(x_1 - x_2 = -h, F[ii,1])):
> end do:
> evalm(Ksim);
> evalm(Tsim);
> evalm(Fsim);
> t := linalg[matrix](3,3,[0,0,0,0,0,0,0,0,0]):
> u := xDisp:
> for ii from 1 to nElem+1 do
>   for jj from 1 to nElem+1 do
>     t[ii,jj] := Ksim[ii,jj]:
>     for mm from 1 to nElem+1 do
>       t[ii,jj] := t[ii,jj] + diff(Ksim[ii,mm],u[jj,1])*u[mm,1]:
>     end do:
>   t[ii,jj] := simplify(t[ii,jj]):
>   end do:
> end do:
> evalm(t);
> Tred := linalg[matrix](2,2,[t[2,2],t[2,3],t[3,2],t[3,3]]);
> Tinv := inverse(Tred);

Can the tangent stiffness matrix be assembled too?[edit]

We have derived the tangent stiffness matrix from the assembled global stiffness matrix. If the stiffness matrix is large that does not seem to be the way to go. So the question is, can the tangent stiffness matrix be calculated on an element by element basis and assembled just like the global stiffness matrix?

Let us try to figure that out.

The element stiffness matrix has the form

The tangent stiffness matrix is given by

After going through the algebra, we get

For the two element mesh, the assembled global tangent stiffness matrix is

Applying the boundary condition , we get

This matrix is the same as that shown in equation (8). Therefore, we can actually assemble the global tangent stiffness matrix. You could also show that in more general terms. But we avoid that proof here.

Doing the actual computation[edit]

The Newton iteration scheme is

Let us go through the process. We will have to start with some assumed values of the parameters.

Let m, m, MPa, kN. Since there are two elements, m.

For the first Newton iteration, let the initial guess be

Then, the first Newton iteration gives

The second Newton iteration gives

The third Newton iteration gives

The fourth Newton iteration gives

At this stage we stop the computation because the solution has converged. A plot of the convergence of the Newton iterations in shown in Figure 1.

Figure 1. Convergence of Newton iterations for axially loaded bar.

Once we have computed the nodal displacements, we can compute the element stresses in the usual manner.

We write the strains as

The stresses are computed using

For the axially loaded bar with no body forces, we have a uniform state of strain (and therefore stress). For the numbers that we have used, these stresses are 10 MPa in each element.

A slightly more complicated situation[edit]

At this stage, let us add in the body force terms and work out the solution for another test case.

Let us assume that the input parameters are , , , , and . The stress-strain curve for this material is shown in Figure 2. We use an error tolerance of as a flag to stop the Newton iterations.


Figure 2. Stress-strain relationship for the nonlinear bar.

Mesh with 10 elements[edit]

If we use a mesh containing equally sized elements, we get the solutions shown in Figure 3.

Figure 3(a). FEM displacement solution for nonlinear bar under distributed axial load. We have used 10 linear elements in these calculations.
Figure 3(b). FEM stress solution for nonlinear bar under distributed axial load. We have used 10 linear elements in these calculations.
Figure 3(c). FEM strain solution for nonlinear bar under distributed axial load. We have used 10 linear elements in these calculations.
Figure 3(d). Error norm of FEM solution for nonlinear bar under distributed axial load. We have used 10 linear elements in these calculations.

Mesh with 100 elements[edit]

If we use a mesh containing equally sized elements, we get the solutions shown in Figure 4.

Figure 4(a). FEM displacement solution for nonlinear bar under distributed axial load. We have used 100 linear elements in these calculations.
Figure 4(b). FEM stress solution for nonlinear bar under distributed axial load. We have used 100 linear elements in these calculations.
Figure 4(c). FEM strain solution for nonlinear bar under distributed axial load. We have used 100 linear elements in these calculations.
Figure 4(d). Error norm of FEM solution for nonlinear bar under distributed axial load. We have used 100 linear elements in these calculations.

Even though the displacements and strains are quite different in the linear and the nonlinear cases, the stresses are remarkably exactly the same for the linear and nonlinear cases. This is a result of the form of the modulus that we have chosen to use for this problem. Notice that the error norms show nearly quadratic convergence.

There is however one issue. Due to the shape of the stress-strain curve of the material, Newton iterations do not converge when the applied external load exceeds the amount needed for the strains to exceed 0.5. A method other than Newton-Raphson is needed to solve the problem for such situations.

The Matlab script used for the above calculations is shown below.

Matlab code[edit]

%==================================================================
% Solve nonlinear axial bar with distributed load (FEM)
%==================================================================
function NLAxialBarFEM

  %
  % Input data (with realistic dimensions)
  %
  %E0 = 1.0e8;
  %R = 1.0e5;
  %a = 100.0;
  %W = 0.1;
  %H = 0.1;
  %A = W*H;
  %L = 2;
  %tol = 1.0e-6;

  %
  % Input data (with low compliance)
  %
  E0 = 10;
  R = 2.0;
  a = 1;
  W = 1;
  H = 1;
  A = 1;
  L = 1;
  tol = 1.0e-6;

  %
  % Create mesh
  %
  nElem = 10;
  h = L/nElem;
  nNode = nElem+1;
  for i=1:nNode
    xnode(i) = (i-1)*h;
  end
  for i=1:nElem
    elem(i,:) = [i i+1];
  end
  
  %
  % Initialize u
  %
  u = zeros(nNode,1)

  %
  % Newton iterations
  %
  condition = 1;
  condition_res = 1;
  count = 0;
  while (condition > tol)
    count = count + 1

    %
    % Initialize K, T, u, f
    %
    K = zeros(nNode);
    T = zeros(nNode);
    f = zeros(nNode,1);

    %
    % Create global stiffness matrix and load vector
    %
    for i=1:nElem
      node1 = elem(i,1);
      node2 = elem(i,2);
      Ke = elementStiffness(A, E0, xnode, u, node1, node2);
      Te = elementTangentStiffness(A, E0, xnode, u, node1, node2);
      fe = elementLoad(a, xnode, node1, node2);
      K(node1:node2,node1:node2) = K(node1:node2,node1:node2) + Ke;
      T(node1:node2,node1:node2) = T(node1:node2,node1:node2) + Te;
      f(node1:node2) = f(node1:node2) + fe;
    end

    %
    % Apply Force BCs
    %
    f(nNode) = f(nNode) + R;

    %
    % Apply Displacement BCs
    %
    u0 = 0;
    Kred = K(2:nNode,2:nNode);
    Tred = T(2:nNode,2:nNode);
    fred = f(2:nNode);
    ured = u(2:nNode);

    %
    % Compute residual
    %
    r = Kred*ured - fred;

    %
    % Compute Tinv
    %
    Tinv = inv(Tred);

    %
    % Do newton iteration
    %
    u_new = ured - Tinv*r;
    condition = norm(u_new - ured)/norm(u_new)
    u = [u0; u_new];
    iter(count) = count;
    normres(count) = norm(r);
    normu(count) = condition;
  end
  %
  % Find forces (to check code)
  %
  K = zeros(nNode);
  for i=1:nElem
    node1 = elem(i,1);
    node2 = elem(i,2);
    Ke = elementStiffness(A, E0, xnode, u, node1, node2);
    K(node1:node2,node1:node2) = K(node1:node2,node1:node2) + Ke;
  end
  fsol = K*u;
  sum(fsol);

  %
  % Plot norms
  %
  figure;
  p01 = semilogy(iter, normres, 'ro-'); hold on;
  p02 = semilogy(iter, normu, 'bo-');
  set([p01 p02], 'LineWidth', 2);
  xlabel('Iterations','FontName','palatino','FontSize',18);
  ylabel('Error Norm','FontName','palatino','FontSize',18);
  set(gca, 'LineWidth', 2, 'FontName','palatino','FontSize',18);
  legend([p01 p02],'|r|','|u_{n+1}-u_n|/|u_{n+1}|');
  axis square;
  
  %
  % Plot disp
  %
  figure;
  p0 = plotDisp(E0, A, L, R, a);
  p = plot(xnode, u, 'ro-', 'LineWidth', 2); hold on;
  xlabel('x', 'FontName', 'palatino', 'FontSize', 18);
  ylabel('u(x)', 'FontName', 'palatino', 'FontSize', 18);
  set(gca, 'LineWidth', 3, 'FontName', 'palatino', 'FontSize', 18);
  legend([p0 p], 'Linear (Exact) ','Nonlinear (FEM)');
  axis square;

  %
  % Find stresses
  %
  for i=1:nElem
    node1 = elem(i,1);
    node2 = elem(i,2);
    [eps(i), sig(i)] = elementStrainStress(E0, xnode, u, node1, node2);
  end
  %
  % Plot stress
  %
  figure;
  p0 = plotStress(E0, A, L, R, a);
  for i=1:nElem
    x1 = xnode(elem(i,1));
    x2 = xnode(elem(i,2));
    p1 = plot([x1 x2], [sig(i) sig(i)], 'r-','LineWidth',3); hold on;
  end
  xlabel('x', 'FontName', 'palatino', 'FontSize', 18);
  ylabel('\sigma(x)', 'FontName', 'palatino', 'FontSize', 18);
  set(gca, 'LineWidth', 3, 'FontName', 'palatino', 'FontSize', 18);
  legend([p0 p1], 'Linear (Exact) ','Nonlinear (FEM)');
  axis square;

  %
  % Plot strain
  %
  figure;
  p0 = plotStrain(E0, A, L, R, a);
  for i=1:nElem
    x1 = xnode(elem(i,1));
    x2 = xnode(elem(i,2));
    p1 = plot([x1 x2], [eps(i) eps(i)], 'r-','LineWidth',3); hold on;
  end
  xlabel('x', 'FontName', 'palatino', 'FontSize', 18);
  ylabel('\epsilon(x)', 'FontName', 'palatino', 'FontSize', 18);
  set(gca, 'LineWidth', 3, 'FontName', 'palatino', 'FontSize', 18);
  legend([p0 p1], 'Linear (Exact) ','Nonlinear (FEM)');
  axis square;

%==================================================================
% Calculate element stiffness
%==================================================================
function [Ke] = elementStiffness(A, E0, xnode, u, node1, node2)

  x1 = xnode(node1);
  x2 = xnode(node2);
  h = x2 - x1;
  u1 = u(node1);
  u2 = u(node2);
  C = A*E0/h^2*(h + u1 - u2);
  Ke = C*[[1 -1];[-1 1]];

%==================================================================
% Calculate element tangent stiffness
%==================================================================
function [Te] = elementTangentStiffness(A, E0, xnode, u, node1, node2)

  x1 = xnode(node1);
  x2 = xnode(node2);
  h = x2 - x1;
  u1 = u(node1);
  u2 = u(node2);
  C = A*E0/h^2*(h + 2*u1 - 2*u2);
  Te = C*[[1 -1];[-1 1]];

%==================================================================
% Calculate element load
%==================================================================
function [fe] = elementLoad(a, xnode, node1, node2)

  x1 = xnode(node1);
  x2 = xnode(node2);
  h = x2 - x1;
  fe1 = a*x2/(2*h)*(x2^2-x1^2) - a/(3*h)*(x2^3-x1^3);
  fe2 = -a*x1/(2*h)*(x2^2-x1^2) + a/(3*h)*(x2^3-x1^3);
  fe = [fe1;fe2];

%==================================================================
% Calculate element stress and strain
%==================================================================
function [eps, sig] = elementStrainStress(E0, xnode, u, node1, node2)

  x1 = xnode(node1);
  x2 = xnode(node2);
  h = x2 - x1;
  u1 = u(node1);
  u2 = u(node2);

  B = [-1/h 1/h];
  u = [u1; u2];
  eps = B*u;
  sig = E0*(1-eps)*eps;

%==================================================================
% Plot displacements
%==================================================================
function [p] = plotDisp(E, A, L, R, a)

  dx = 0.01;
  nseg = L/dx;
  for i=1:nseg+1
    x(i) = (i-1)*dx;
    u(i) = 1/(6*A*E)*(-a*x(i)^3 + (6*R + 3*a*L^2)*x(i));
  end
  p = plot(x, u, 'LineWidth', 3); hold on;
  xlabel('x', 'FontName', 'palatino', 'FontSize', 18);
  ylabel('u(x)', 'FontName', 'palatino', 'FontSize', 18);
  set(gca, 'LineWidth', 3, 'FontName', 'palatino', 'FontSize', 18);

%==================================================================
% Plot strains
%==================================================================
function [p] = plotStrain(E, A, L, R, a)

  dx = 0.01;
  nseg = L/dx;
  for i=1:nseg+1
    x(i) = (i-1)*dx;
    eps(i) = 1/(2*A*E)*(-a*x(i)^2 + (2*R + a*L^2));
  end
  p = plot(x, eps, 'LineWidth', 3); hold on;
  xlabel('x', 'FontName', 'palatino', 'FontSize', 18);
  ylabel('\epsilon(x)', 'FontName', 'palatino', 'FontSize', 18);
  set(gca, 'LineWidth', 3, 'FontName', 'palatino', 'FontSize', 18);

%==================================================================
% Plot element stress
%==================================================================
function [p] = plotStress(E, A, L, R, a)

  dx = 0.01;
  nseg = L/dx;
  for i=1:nseg+1
    x(i) = (i-1)*dx;
    sig(i) = 1/(2*A)*(-a*x(i)^2 + (2*R + a*L^2));
  end
  p = plot(x, sig, 'LineWidth', 3); hold on;
  xlabel('x', 'FontName', 'palatino', 'FontSize', 18);
  ylabel('\sigma(x)', 'FontName', 'palatino', 'FontSize', 18);
  set(gca, 'LineWidth', 3, 'FontName', 'palatino', 'FontSize', 18);