Nonlinear finite elements/Homework 1/Solutions

Problem 1: Mathematical Model Development edit


Figure 1. The motion of a pendulum


  1. The rod and the mass at the end are rigid.
  2. The mass of the rod is negligible compared to the mass of the bob.
  3. There is no friction at the pivot.

Solution edit

Part 1 edit

Find the equation of motion of the bob (angular displacement as a function of time). State all other assumptions.

The component gravitational force in the direction tangent to the motion is


The tangential velocity of the bob is


where   is an arc length along the motion of the bob.

Therefore, the tangential acceleration of the bob is


From Newton's second law (balance of momentum) we have


Therefore the equation of motion of the bob is


Additional assumptions are lack of viscosity in the surrounding medium, no wind load, and so on.

Part 2 edit

Why is the equation of motion "nonlinear"?

Let us try the technique we learned in class. Let


Then, the equation of motion can be written as


Without making further assumptions, this equation cannot be expressed as the equation of a line of the form


Hence the equation of motion is nonlinear.

Part 4 edit

Assume that   is "small". Derive the equation of motion for this situation.

From elementary calculus, for small values of   we have


Therefore, for small values of  , the equation of motion becomes


Part 5 edit

Why is the small   equation of motion "linear"?

In this case we can write the differential equation as


This equation can be expressed as the equation of a line of the form


Hence the equation of motion is linear.

Part 6 edit

Derive a finite element formulation for the linear equation of motion.

We know that if we can formulate the finite element system of equations for one element, we can assemble the element equations into a global system.

Let us find the finite element system of equations for one element. In this problem, we discretize the time variable into elements. Let the element extend from time   to time  .

The first step in the process is the derivation of a approximate weak form of the differential equation.

Let   be a weighting function. After multiplying the ODE by the weighting function and integrating the product over the element, we get


After integrating the first term by parts, we get


The rearranged weak form is


Let   be an approximate solution. To get the approximate weak form, we plug   into the exact weak form. Notice that the right hand side of the exact weak form contains a weighted boundary lux term. We do not approximate this term because it is presumed to be specified or known exactly. In this problem, this term represents the angular velocity  . Then the approximate weak form is


At this stage we choose the form of the approximate solution. Let us assume that each element has two nodes and the approximate solution within an element has the form


where   is the nodal rotation at node   of the element and   is the nodal rotation at node   of the element and   and   are the corresponding shape functions.

We have seen that instead of thinking of the weighting function as a series of the same form as  , we can just consider these functions to have the form   - where the index   stands for a row of the finite element system of equations. Therefore, we write


Plug the approximate solution and the weighting function into the weak form to get


Take the summation sign and the constants outside the integral and get


Collect terms on the left hand side


Rearrange and get


The above has the form






Given the element stiffness matrix and the element load vector, we can assemble a global system of equations and solve it for the unknown values of  .

Let us assume that the shape functions are


where  .

Then, the element stiffness matrix terms are


The load vector terms are


In matrix form:


The Maple script for getting the components of the stiffness matrix is shown below.

        N1 := (t2-t)/(t2-t1);
        N2 := (t-t1)/(t2-t1);
        dN1_dt := diff(N1,t);
        dN2_dt := diff(N2,t);
        k11 := dN1_dt*dN1_dt - lambda^2*N1*N1;
        k12 := dN1_dt*dN2_dt - lambda^2*N1*N2;
        k21 := dN2_dt*dN1_dt - lambda^2*N2*N1;
        k22 := dN2_dt*dN2_dt - lambda^2*N2*N2;
        K11 := int(k11, t=t1..t2);
        K12 := int(k12, t=t1..t2);
        K21 := int(k21, t=t1..t2);
        K22 := int(k22, t=t1..t2);
        simplify(K11); simplify(K12); simplify(K21); simplify(K22);

Part 7 edit

Use ANSYS or some other tool of your choice (including your own code) to solve the linear equation of motion via finite elements. Compare with the exact solution. You can use any reasonable values for time ( ), mass ( ), and length ( ).

The Matlab script for solving the problem is given below:

       function PendulumFEM(situation)     
         T = 2;     
         g = 10;    
         l = 1;         
         if (situation == 1)
           omega0 = 0;
           theta0 = 0.1;
         elseif (situation == 2)
           omega0 = 1;
           theta0 = 0;
           omega0 = 1;
           theta0 = 0.1;

         e = 10;
         h = T/e;
         n = e+1;
         for i=1:n
           tnode(i) = (i-1)*h;
         for i=1:e
           elem(i,:) = [i i+1];

         K = zeros(n);
         f = zeros(n,1);
         for i=1:e
           node1 = elem(i,1);
           node2 = elem(i,2);
           Ke = elementStiffness(tnode(node1), tnode(node2), g, l);
           fe = elementLoad(tnode(node1),tnode(node2));
           K(node1:node2,node1:node2) = K(node1:node2,node1:node2) + Ke;
           f(node1:node2) = f(node1:node2) + fe;

         f(1) = f(1) + omega0;
         d1 = theta0;
         for i=1:n
           f(i) = f(i) - K(i,1)*d1;
         Kred = K(2:n,2:n);
         fred = f(2:n);
         d = inv(Kred)*fred;
         dsol = [d1 d'];
         p0 = plotTheta(T, g, l, omega0, theta0);
         p1 = plot(tnode, dsol, 'ro--', 'LineWidth', 3); hold on;
         legend([p0 p1],'Exact','FEM');
         s = sprintf('\\theta(0) = 
         gtext(s, 'FontName', 'palatino', 'FontSize', 18);
       function [p] = plotTheta(T, g, l, omega0, theta0)
         lambda = sqrt(g/l);
         omega0_lambda = omega0/lambda;
         dt = 0.01;
         nseg = T/dt;
         for i=1:nseg+1
           t(i) = (i-1)*dt;
           theta(i) = omega0_lambda*sin(lambda*t(i)) + theta0*cos(lambda*t(i));
         p = plot(t, theta, 'LineWidth', 3); hold on;
         xlabel('t', 'FontName', 'palatino', 'FontSize', 18);
         ylabel('\theta(t)', 'FontName', 'palatino', 'FontSize', 18);
         set(gca, 'LineWidth', 3, 'FontName', 'palatino', 'FontSize', 18);

       function [Ke] = elementStiffness(node1, node2, g, l)
         lambda = sqrt(g/l);
         t1 = node1;
         t2 = node2;
         delT = t2 - t1;
         Ke = zeros(2);
         Ke(1,1) = (lambda*delT)^2 - 3;
         Ke(1,2) = (lambda*delT)^2/2 + 3;
         Ke(2,2) = (lambda*delT)^2 - 3;
         Ke(2,1) = Ke(1,2);

       function [fe] = elementLoad(node1,node2)

         x1 = node1;
         x2 = node2;
         fe1 = 0;
         fe2 = 0;
         fe = [fe1;fe2];

Comparisons are shown below between the FEM solution and the exact solution for three sets of initial conditions.

Figure 2 (a). Comparisons between FEM and exact solution
Figure 2 (b). Comparisons between FEM and exact solution
Figure 2 (c). Comparisons between FEM and exact solution
  • When the initial angular velocity is 0 and the initial angle is nonzero, the FEM solution is close to the exact solution. This is a lucky coincidence. The total time interval has been taken to be 2.0 secs. It turns out that at exactly 2 secs, the angular velocity goes to zero and hence the problem appears as if a angular velocity boundary condition has been applied at time   secs. You can see the problem if you increase the time interval to   secs (see Figure 3).
Figure 3. Comparisons between FEM and exact solution if t = 2.2 secs
  • In problems that involve time, we usually know only the initial state but not the final state. Therefore, integrals over the total time interval cannot be evaluated easily. That is one of the reasons why we don't use finite elements to discretize time. Instead, we use techniques that directly discretize the differential operator, i.e., we work from the strong form.
  • You cannot solve a problem without specifying BCs at all points on the boundary. Applying two BCs at a single point does not work unless they are for different variables.

Part 8 edit

What happens when you try to solve the nonlinear equation of motion with finite elements?

In this case, the approximate weak form becomes


If we plug in the approximate solution, we get


To get a linear system of equations, we can take the   term to the right hand side.


At this stage we have to start with an initial guess for   and  , plug them into the   term, integrate numerically and form a right hand side. We can then solve for the   values and get a new set. We can then plug these back into the RHS and repeat the process. A straightforward linear solution is no longer possible.

Problem 2: Numerical Exercise edit


Mini-centrifuge problem.

Mini centrifuge.

Temperature range: -80 C to 100 C., RPM: 25000 .

Assumptions: Assume that the current material is a light metal (any other material of your choice will also do). The only requirement is that it should be cheap and stiff. Ignore bending and the weight of the sample tubes.


Find the thinnest beam of the material of your choice that can stand the axial load that is generated due to rotation at the new speed. Use ANSYS or your own code to solve the problem. State all simplifying assumptions.

Solution edit

Let us consider one of the spokes and assume that it's a one-dimensional bar. The governing equation equation for a 1-D bar is


where   is a body force which has units of force per unit length.

For a rotating bar, the axial acceleration is


where   is the angular velocity of rotation (assumed constant in this case). The acceleration of the fixed end is assumed to be zero (but may not be zero for the actual centrifuge).

Therefore, the body force per unit length due to the rotation of the bar is


where   is the density of the material, and   is the area of cross-section.

Therefore, the governing differential equation is




Clearly, the solution for the displacement does not depend on the area of cross-section of the material. Hence, the strains do not depend on the area of cross-section and therefore, neither do the stresses. However, the stresses do depend on the length.

If we assume that the width and the length of the spokes have to be maintained constant in order to hold the tubes, the thickness of the spokes cannot be optimized on the basis of inertial forces. We have to consider bending.

If we work out the finite element equations over an element for this problem, we will get


The second equation above shows you how to apply the body forces (assuming ANSYS does not allow that - which is surprising.)

At this stage, an ANSYS solution becomes a formality that may or may not be performed.

For the ANSYS solution, let us assume that

  1. The material is 6061-T6 aluminum alloy.
  2. The density is 2790 kg/m .
  3. The Young's modulus at 100 C is 68.8 GPa.
  4. The Poisson's ratio is 0.35.
  5. The initial yield stress at -80 C is 350 MPa.
  6. The initial yield stress at 100 C is 250 MPa.
  7. The fracture strain at -80 C is 0.05.
  8. The fracture strain at 100 C is 0.50.
  9. The length is 140 mm.
  10. The width is 40 mm.
  11. The thickness is 10 mm.

Here's a sample ANSYS script for the problem

    ! Constants
    pi = 4*atan(1)
    ! Angular velocity
    radians = 2*pi
    second = 60
    omega = 25000*radians/second
    ! Geometry
    length = 0.140
    width = 0.040
    thickness = 0.010
    area = width*thickness
    ! Material
    rho = 2.79e3
    modulus = 68.8e9
    nu = 0.35
    ! External forces
    applF = 0
    bodyF = rho*area*omega**2
    ! Mesh size
    numElem = 3
    elemSize = length/numElem
    ! Element type and real constants
    et, 1, 1
    r, 1, area
    ! Set up material props.
    mp, ex,   1, modulus
    mp, prxy, 1, nu
    ! Create nodes
    *do, ii, 1, numElem+1
        n, ii, elemSize*(ii-1), 0
    ! Create elements
    *do, ii, 1, numElem
        e, ii, ii+1
    ! Create nodal forces that correspond to the body force
    *dim, x,  array, numElem+1, 1
    *dim, fi, array, numElem, 1
    *dim, fj, array, numElem, 1
    x(1) = 0
    *do, ii, 2, numElem+1
         x(ii) = x(ii-1)+elemSize
    *do, e, 1, numElem
        x1 = x(e)
        x2 = x(e+1)
        fi(e) = bodyF/elemSize*(x2*(x2**2-x1**2)/2-(x2**3-x1**3)/3)
        fj(e) = bodyF/elemSize*(-x1*(x2**2-x1**2)/2+(x2**3-x1**3)/3)
    *end do
    ! Apply nodal forces
    f, 1, fx, fi(1)
    *do, n, 2, numElem
        f, n, fx, fi(n)+fj(n-1)
    ! Apply displacement BCs
    d, 1, all, 0
    ! Solve system of equations
    ! Get the results

Results from ANSYS edit

Ideally you should run a convergence study at this stage to see how close your results are to the actual solution and show plots of the results.

Numbers that have five decimal places or are close to machine precision are usually useless when presented in tables.

From the axial bar problem in the learning resources, you know the exact solution for this problem, which is

Node Location ( )
1 0 0
2 46.7 0.12
3 93.4 0.22
4 140 0.25
Element Axial stress
1 180
2 139
3 56

Here's a comparison between the exact solution and the finite element solution

Comparison of FE and exact solutions for axial bar.

The stresses are lower than the yield stress of Aluminum, but they are close. The strains are also relatively small (approximately 0.002). The structure should hold up under the new conditions.