Problems for Chapter 8
Theory and Implementation of the Finite Element Method
8.2. The Finite Element Method For Dynamic Linear Elasticity
8.2.1. Set up the general purpose 2D/3D dynamic finite element code to solve the 1D wave propagation problem shown in the figure. Use material properties (arbitrary units). Mesh the bar with square 4 noded quadrilateral plane strain elements. Assume that the left hand end of the bar is subjected to a constant traction with magnitude p=5 at time t=0.
188.8.131.52. Calculate the maximum time step for a stable explicit computation.
184.108.40.206. Run an explicit dynamics computation (Newmark parameters ) with a time step equal to half the theoretical stable limit. Use the simulation to plot the time variation of displacement, velocity and acceleration at x=0. Use at least 160 time steps. Compare the numerical solution with the exact solution.
220.127.116.11. Repeat 1.2 with a time step equal to twice the theoretical stable limit.
18.104.22.168. Repeat 1.2 and 1.3 with Newmark parameters .
22.214.171.124. Repeat 1.2 and 1.3 with Newmark parameters . Try a time step equal to ten times the theoretical stable limit.
8.2.2. The 2D/3D dynamic finite element code uses the row-sum method to compute lumped mass matrices.
126.96.36.199. Use the row-sum method to compute lumped mass matrices for linear and quadratic triangular elements with sides of unit length and mass density (you can simply print out the matrix from the code)
188.8.131.52. Use the row-sum method to compute lumped mass matrices for linear and quadratic quadrilateral elements with sides of unit length and mass density
184.108.40.206. Modify the finite element code to compute the lumped mass matrix using the scaled diagonal method instead. Repeat 2.1 and 2.2.
220.127.116.11. Repeat problem 1.2 (using Newmark parameters ) using the two versions of the lumped mass matrices, using quadratic quadrilateral elements to mesh the solid. Compare the results with the two versions of the mass matrix.
8.2.3. In this problem you will implement a simple 1D finite element method to calculate the motion of a stretched vibrating string. The string has mass per unit length m, and is stretched by applying a tension T at one end. Assume small deflections. The equation of motion for the string (see Section 10.3.1) is
and the transverse deflection must satisfy at .
18.104.22.168. Weak form of the equation of motion. Let be a kinematically admissible variation of the deflection, satisfying at . Show that if
for all admissible , then w satisfies the equation of motion.
22.214.171.124. Finite element equations Introduce a linear finite element interpolation as illustrated in the figure. Calculate expressions for the element mass and stiffness matrices (you can calculate analytical expressions for the matrices, or evaluate the integrals numerically)
126.96.36.199. Implementation: Write a simple code to compute and integrate the equations of motion using the Newmark time integration procedure described in Section 8.2.5.
188.8.131.52. Testing: Test your code by computing the motion of the string with different initial conditions. Try the following cases:
· , n=1,2…. These are the mode shapes for string, so the vibration should be harmonic. You can calculate the corresponding natural frequency by substituting into the equation of motion and solving for
Use geometric and material parameters L=10, m=1, T=1. Run a series of tests to investigate (i) the effects of mesh size; (ii) the effects of using a lumped or consistent mass matrix; (ii) the effects of time step size; and (iii) the effects of the Newmark parameters on the accuracy of the solution.
8.2.4. Modify the 1D code described in the preceding problem to calculate the natural frequencies and mode shapes for the stretched string. Calculate the first 5 natural frequencies and mode shapes, and compare the numerical results with the exact solution.
8.2.5. The figure below shows a bar with thermal conductivity and heat capacity c. At time t=0 the bar is at uniform temperature, . The sides of the bar are insulated to prevent heat loss; the left hand end is held at fixed temperature , while a flux of heat is applied into to the right hand end of the bar. In addition, heat is generated inside the bar at rate per unit length.
The temperature distribution in the bar is governed by the 1-D heat equation
and the boundary condition at is
In this problem, you will (i) set up a finite element method to solve the heat equation; (ii) show that the finite element stiffness and mass matrix for this problem are identical to the 1-D elasticity problem solved in class; (iii) implement a simple stepping procedure to integrate the temperature distribution with respect to time.
184.108.40.206. Variational statement of the heat equation. Let be an admissible variation of temperature, which must be (a) differentiable and (b) must satisfy on the left hand end of the bar. Begin by showing that if the temperature distribution satisfies the heat equation and boundary condition listed above, then
where denotes the value of at the left hand end of the bar. To show this, use the procedure that was used to derive the principle of virtual work. First, integrate the first integral on the left hand side by parts to turn the into , then use the boundary conditions and governing equation to show that the variational statement is true. (Of course we actually rely on the converse if the variational statement is satisfied for all then the field equations are satisfied)
220.127.116.11. Finite element equations Now, introduce a linear 1D finite element interpolation scheme as described in Section 8.2.5. Show that the diffusion equation can be expressed in the form
where , b=1,2…N denotes the unknown nodal values of temperature, is a heat capacity matrix analogous to the mass matrix defined in 8.2.5, and is a stiffness matrix. Derive expressions for the element heat capacity matrix, and the element stiffness matrix, in terms of relevant geometric and material properties.
18.104.22.168. Time integration scheme Finally, we must devise a way to integrate the temperature distribution in the bar with respect to time. Following the usual FEA procedure, we will use a simple Euler-type time stepping scheme. To start the time stepping, we note that is known at t=0, and note that we can also compute the rate of change of temperature at time t=0 by solving the FEM equations
For a generic time step, our objective is to compute the temperature and its time derivative at time . To compute the temperature at the end of the step, we write
and estimate based on values at the start and end of the time step as
where is an adjustable numerical parameter. The time derivative of temperature at must satisfy the FEA equations
Hence, show that the time derivative of temperature at time can be calculated by solving
whereupon the temperature at follows as
22.214.171.124. Implementation. Modify the 1D dynamic FEA code to calculate the temperature variation in a 1D bar.
126.96.36.199. Testing: Test the code by setting the bar length to 5; x-sect area to 1; heat capacity=100, thermal conductivity=50; set the heat generation in the interior to zero (bodyforce=0) and set the heat input to the right hand end of the bar to 10 (traction=10), then try the following tests:
· To check the code, set the number of elements to L=15; set the time step interval to dt=0.1and the number of steps to nstep=1000. Also, set up the code to use a lumped mass matrix by setting lumpedmass=true and explicit time stepping . You should find that this case runs quickly, and that the temperature in the bar gradually rises until it reaches the expected linear distribution.
· Show that the explicit time stepping scheme is conditionally stable. Try running with dt=0.2 for just 50 or 100 steps.
· Show that the critical stable time step size reduces with element size. Try running with dt=0.1 with 20 elements in the bar.
· Try running with a fully implicit time stepping scheme with , 15 elements and dt=0.2. Use a full consistent mass matrix for this calculation. Show that you can take extremely large time steps with without instability (eg try dt=10 for 10 or 100 steps) in fact the algorithm is unconditionally stable for
(c) A.F. Bower, 2008