Nonlinear finite elements/Homework 1/Solutions
Problem 1: Mathematical Model Development
- The rod and the mass at the end are rigid.
- The mass of the rod is negligible compared to the mass of the bob.
- There is no friction at the pivot.
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.
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.
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
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.
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
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);
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; 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.
- 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).
- 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.
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
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.
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
- The material is 6061-T6 aluminum alloy.
- The density is 2790 kg/m.
- The Young's modulus at 100 C is 68.8 GPa.
- The Poisson's ratio is 0.35.
- The initial yield stress at -80 C is 350 MPa.
- The initial yield stress at 100 C is 250 MPa.
- The fracture strain at -80 C is 0.05.
- The fracture strain at 100 C is 0.50.
- The length is 140 mm.
- The width is 40 mm.
- 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
|Element||Axial stress |
Here's a comparison between the exact solution and the finite element solution
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.