Nonlinear finite elements/Homework 1/Solutions

Problem 1: Mathematical Model Development

Given:

File:NFE Homework1 1.png Figure 1. The motion of a pendulum

Assumptions:

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

Part 1

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

$F_x=-mg\sin \theta .$

The tangential velocity of the bob is

$v_x=\frac{ds}{dt}=\frac{d}{dt}(l\theta )=l\frac{d\theta }{dt}$

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

Therefore, the tangential acceleration of the bob is

$a_x=\frac{d}{dt}\left(l\frac{d\theta }{dt}\right)=l\frac{d^2\theta }{d{t}^2}.$

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

$F_x=m{a}_x⟹-mg\sin \theta =ml\frac{d^2\theta }{d{t}^2}.$

Therefore the equation of motion of the bob is

$\frac{d^2\theta }{d{t}^2}+{\lambda }^2\sin \theta =0,{\lambda }^2=\frac{g}{l}.$

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

Part 2

Why is the equation of motion "nonlinear"?

Let us try the technique we learned in class. Let

$A:=\frac{d^2\theta }{d{t}^2};\lambda =1.$

Then, the equation of motion can be written as

$A+\sin \theta =0.$

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

$y=m\theta .$

Hence the equation of motion is nonlinear.

Part 4

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

From elementary calculus, for small values of $\theta$ we have

$\sin \theta \approx \theta .$

Therefore, for small values of $\theta$, the equation of motion becomes

$\frac{d^2\theta }{d{t}^2}+{\lambda }^2\theta =0.$

Part 5

Why is the small $\theta$ equation of motion "linear"?

In this case we can write the differential equation as

$A+\theta =0.$

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

$y=m\theta .$

Hence the equation of motion is linear.

Part 6

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 $t_1$ to time $t_2$.

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

Let $w(t)$ be a weighting function. After multiplying the ODE by the weighting function and integrating the product over the element, we get

$\underset{t_1}{\overset{t_2}{\int }}(\frac{d^2\theta }{d{t}^2}+{\lambda }^2\theta )w(t)dt=0.$

After integrating the first term by parts, we get

${\left(\frac{d\theta }{dt}w\right)|}_{t_1}^{t_2}-\underset{t_1}{\overset{t_2}{\int }}\frac{d\theta }{dt}\frac{dw}{dt}dt+{\lambda }^2\underset{t_1}{\overset{t_2}{\int }}\theta wdt=0.$

The rearranged weak form is

$\underset{t_1}{\overset{t_2}{\int }}\frac{d\theta }{dt}\frac{dw}{dt}dt-{\lambda }^2\underset{t_1}{\overset{t_2}{\int }}\theta wdt={\left(\frac{d\theta }{dt}w\right)|}_{t_1}^{t_2}.$

Let ${\theta }_h(t)$ be an approximate solution. To get the approximate weak form, we plug ${\theta }_h$ 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 $\omega =d\theta /dt$. Then the approximate weak form is

$\underset{t_1}{\overset{t_2}{\int }}\frac{d{\theta }_h}{dt}\frac{dw}{dt}dt-{\lambda }^2\underset{t_1}{\overset{t_2}{\int }}{\theta }_{h}wdt={\left(\omega w\right)|}_{t_1}^{t_2}.$

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

${\theta }_h(t)={\theta }_1{N}_1(t)+{\theta }_2{N}_2(t)=\sum _{j=1}^2{\theta }_j{N}_j(t)$

where ${\theta }_1$ is the nodal rotation at node $1$ of the element and ${\theta }_2$ is the nodal rotation at node $2$ of the element and $N_1$ and $N_2$ are the corresponding shape functions.

We have seen that instead of thinking of the weighting function as a series of the same form as ${\theta }_h(t)$, we can just consider these functions to have the form $N_i$ - where the index $i$ stands for a row of the finite element system of equations. Therefore, we write

$w=N_i(t).$

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

$\underset{t_1}{\overset{t_2}{\int }}\left(\sum _{j=1}^2{\theta }_j\frac{d{N}_j}{dt}\right)\frac{d{N}_i}{dt}dt-{\lambda }^2\underset{t_1}{\overset{t_2}{\int }}\left(\sum _{j=1}^2{\theta }_j{N}_j\right)N_{i}dt={\left(\omega N_i\right)|}_{t_1}^{t_2}.$

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

$\sum _{j=1}^2\left({\theta }_j\underset{t_1}{\overset{t_2}{\int }}\frac{d{N}_j}{dt}\frac{d{N}_i}{dt}dt\right)-{\lambda }^2\sum _{j=1}^2\left({\theta }_j\underset{t_1}{\overset{t_2}{\int }}{N}_j{N}_{i}dt\right)={\left(\omega N_i\right)|}_{t_1}^{t_2}.$

Collect terms on the left hand side

$\sum _{j=1}^2\left[{\theta }_j\underset{t_1}{\overset{t_2}{\int }}(\frac{d{N}_j}{dt}\frac{d{N}_i}{dt}-{\lambda }^2{N}_j{N}_i)dt\right]={\left(\omega N_i\right)|}_{t_1}^{t_2}.$

Rearrange and get

$\sum _{j=1}^2\left[\underset{t_1}{\overset{t_2}{\int }}(\frac{d{N}_i}{dt}\frac{d{N}_j}{dt}-{\lambda }^2{N}_i{N}_j)dt\right]{\theta }_j={\left(\omega N_i\right)|}_{t_1}^{t_2}.$

The above has the form

$\sum _{j=1}^2{K}_{ij}{\theta }_j=f_i$

where

$K_{ij}=\underset{t_1}{\overset{t_2}{\int }}(\frac{d{N}_i}{dt}\frac{d{N}_j}{dt}-{\lambda }^2{N}_i{N}_j)dt\leftarrow \text{Element stiffness matrix.}$

and

$f_i={\left(\omega N_i\right)|}_{t_1}^{t_2}.\leftarrow \text{Element load vector.}$

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 ${\theta }_i$.

Let us assume that the shape functions are

$N_1(t)=\frac{t_2-t}{\mathrm{\Delta }t}\text{and}{N}_2(t)=\frac{t-t_1}{\mathrm{\Delta }t}$

where $\mathrm{\Delta }t=t_2-t_1$.

Then, the element stiffness matrix terms are

$\begin{array}{cc}{K}_{11}& =\underset{t_1}{\overset{t_2}{\int }}(\frac{d{N}_1}{dt}\frac{d{N}_1}{dt}-{\lambda }^2{N}_1{N}_1)dt=\underset{t_1}{\overset{t_2}{\int }}[\frac{1}{\mathrm{\Delta }{t}^2}-{\lambda }^2{\left(\frac{t_2-t}{\mathrm{\Delta }t}\right)}^2]dt=\frac{1}{3\mathrm{\Delta }t}({\lambda }^2\mathrm{\Delta }{t}^2-3)\\ K_{12}& =\underset{t_1}{\overset{t_2}{\int }}(\frac{d{N}_1}{dt}\frac{d{N}_2}{dt}-{\lambda }^2{N}_1{N}_2)dt=\underset{t_1}{\overset{t_2}{\int }}[-\frac{1}{\mathrm{\Delta }{t}^2}-{\lambda }^2\frac{(t-t_1)(t_2-t)}{\mathrm{\Delta }{t}^2}]dt=\frac{1}{6\mathrm{\Delta }t}({\lambda }^2\mathrm{\Delta }{t}^2+6)\\ K_{22}& =\underset{t_1}{\overset{t_2}{\int }}(\frac{d{N}_2}{dt}\frac{d{N}_2}{dt}-{\lambda }^2{N}_2{N}_2)dt=\underset{t_1}{\overset{t_2}{\int }}[\frac{1}{\mathrm{\Delta }{t}^2}-{\lambda }^2{\left(\frac{t-t_1}{\mathrm{\Delta }t}\right)}^2]dt=\frac{1}{3\mathrm{\Delta }t}({\lambda }^2\mathrm{\Delta }{t}^2-3)\end{array}$

The load vector terms are

$\begin{array}{cc}{f}_1& ={\left(\omega N_1\right)|}_{t_1}^{t_2}=\omega (t_2)N_1(t_2)-\omega (t_1)N_1(t_1)=-\omega (t_1)=:-{\omega }_1\\ f_2& ={\left(\omega N_2\right)|}_{t_1}^{t_2}=\omega (t_2)N_2(t_2)-\omega (t_1)N_2(t_1)=\omega (t_2)=:{\omega }_2\\ \end{array}$

In matrix form:

$\frac{1}{3\mathrm{\Delta }t}\left[\begin{array}{cc}(\lambda \mathrm{\Delta }t{)}^2-3& (\lambda \mathrm{\Delta }t{)}^2/2+3\\ (\lambda \mathrm{\Delta }t{)}^2/2+3& (\lambda \mathrm{\Delta }t{)}^2-3\end{array}\right]\left[\begin{array}{c}{\theta }_1\\ {\theta }_2\end{array}\right]=\left[\begin{array}{c}-{\omega }_1\\ {\omega }_2\end{array}\right]$

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

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 ($t$), mass ($m$), and length ($l$).

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;
         else
           omega0 = 1;
           theta0 = 0.1;
         end

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

         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;
         end

         f(1) = f(1) + omega0;
 
         d1 = theta0;
         for i=1:n
           f(i) = f(i) - K(i,1)*d1;
         end
         Kred = K(2:n,2:n);
         fred = f(2:n);
  
         d = inv(Kred)*fred;
     
         dsol = [d1 d'];
      
         figure;
         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));
         end
         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.

File:PendulumSol1.png Figure 2 (a). Comparisons between FEM and exact solution

File:PendulumSol2.png Figure 2 (b). Comparisons between FEM and exact solution

File:PendulumSol3.png 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 $t=2$ secs. You can see the problem if you increase the time interval to $t=2.2$ secs (see Figure 3).

File:PendulumSol4.png 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

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

In this case, the approximate weak form becomes

$\underset{t_1}{\overset{t_2}{\int }}\frac{d{\theta }_h}{dt}\frac{dw}{dt}dt-{\lambda }^2\underset{t_1}{\overset{t_2}{\int }}\sin ({\theta }_h)wdt={\left(\omega w\right)|}_{t_1}^{t_2}.$

If we plug in the approximate solution, we get

$\sum _{j=1}^2\left(\underset{t_1}{\overset{t_2}{\int }}\frac{d{N}_i}{dt}\frac{d{N}_j}{dt}dt\right){\theta }_j-{\lambda }^2\underset{t_1}{\overset{t_2}{\int }}\sin ({\theta }_h)N_{i}dt={\left(\omega N_i\right)|}_{t_1}^{t_2}.$

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

$\sum _{j=1}^2\left(\underset{t_1}{\overset{t_2}{\int }}\frac{d{N}_i}{dt}\frac{d{N}_j}{dt}dt\right){\theta }_j={\lambda }^2\underset{t_1}{\overset{t_2}{\int }}\sin ({\theta }_h)N_{i}dt+{\left(\omega N_i\right)|}_{t_1}^{t_2}.$

At this stage we have to start with an initial guess for ${\theta }_1$ and ${\theta }_2$, plug them into the $\sin (\theta )$ term, integrate numerically and form a right hand side. We can then solve for the $\theta$ 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

Given:

Mini-centrifuge problem.

File:NFE Homework1 2a.png 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:

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

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

$AE\frac{d^{2}u}{d{x}^2}+q=0.$

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

For a rotating bar, the axial acceleration is

$a(x)={\omega }^{2}x$

where $\omega$ 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

$q(x)=\rho Aa(x)=\rho A{\omega }^{2}x$

where $\rho$ is the density of the material, and $A$ is the area of cross-section.

Therefore, the governing differential equation is

$AE\frac{d^{2}u}{d{x}^2}+\rho A{\omega }^{2}x=0$

or,

$E\frac{d^{2}u}{d{x}^2}+\rho {\omega }^{2}x=0.$

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

$\begin{array}{cc}{K}_{ij}^e& =\underset{x_1}{\overset{x_2}{\int }}\frac{d{N}_j}{dx}E\frac{d{N}_i}{dx}dx\\ f_j^e& =\rho {\omega }^2\underset{x_1}{\overset{x_2}{\int }}x{N}_{j}dx\end{array}$

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^{}$.
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

    /prep7
    !
    ! 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
    *enddo
    !
    ! Create elements
    !
    *do, ii, 1, numElem
        e, ii, ii+1
    *enddo
    !
    ! 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
    *enddo
    
    *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)
    *enddo
    f,numElem+1,fx,fj(numElem)+applF
    !
    ! Apply displacement BCs
    !
    d, 1, all, 0
    finish
    
    !
    ! Solve system of equations
    !
    /solu
    solve
    finish
    
    !
    ! Get the results
    !
    /post1
    ETABLE,saxl,LS,1 
    /output,HW1Prob2,sol
    prnsol,ux 
    pretab,saxl 
    /output,,
    finish

Results from ANSYS

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

$\begin{array}{cc}u(x)& =\frac{\rho {\omega }^2}{6E}(3L^{2}x-x^3)\\ \sigma (x)& =\frac{\rho {\omega }^2}{2}(L^2-x^2)\end{array}$

Node	Location ($x$) (mm)	Nodal Displacement (mm)
1	0	0
2	46.7	0.12
3	93.4	0.22
4	140	0.25

Element	Axial stress (MPa)
1	180
2	139
3	56

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

File:NFE HW1Prob2Sol.png Comparison of FE and exact solutions for axial bar.

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