|
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.
8.2.1.1.
Calculate the maximum time step for a stable explicit computation.
8.2.1.2.
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.
8.2.1.3.
Repeat 1.2 with
a time step equal to twice the theoretical stable limit.
8.2.1.4.
Repeat 1.2 and
1.3 with Newmark parameters .
8.2.1.5.
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.Â
8.2.2.1.
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)
8.2.2.2.
Use the row-sum method to compute lumped mass matrices for linear and
quadratic quadrilateral elements with sides of unit length and mass density
8.2.2.3.
Modify the finite element code to compute the lumped mass matrix
using the scaled diagonal method instead. Â
Repeat 2.1 and 2.2.
8.2.2.4.
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 .
8.2.3.1.
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.
8.2.3.2.
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)
8.2.3.3.
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.
8.2.3.4.
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.
8.2.5.1.
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)
8.2.5.2.
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.
8.2.5.3.
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
8.2.5.4.
Implementation. Modify the 1D dynamic FEA code to calculate the temperature
variation in a 1D bar.
8.2.5.5.
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
|
