Nonlinear finite elements/Homework 1/Solutions

From Wikiversity
Jump to navigation Jump to search

Problem 1: Mathematical Model Development[edit | edit source]


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 | edit source]

Part 1[edit | edit source]

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 | edit source]

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 | edit source]

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 | edit source]

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 | edit source]

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 | edit source]

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 | edit source]

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 | edit source]


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 | edit source]

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 | edit source]

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.