.. # gedit: set fileencoding=utf8 : .. _demo_elastodynamics: ============================================ Time-integration of elastodynamics equation ============================================ This demo is implemented in a single Python file, :download:demo_elastodynamics.py, which contains both the variational forms and the solver. This demo shows how to perform time integration of transient elastodynamics using the generalized-:math:\alpha method [ERL2002]_. In particular it demonstrates how to: * formulate mass and damping forms of the elastodynamics equation * implement the generalized-:math:\alpha method and its influence on the solution * perform efficient computation of stresses using LocalSolver The deformed structure evolution over time along with the axial stress will look as follows: .. image:: elastodynamics-transient_deformation.gif :width: 70% :align: center ----------------------------------------- Introduction and elastodynamics equation ----------------------------------------- The elastodynamics equation combine the balance of linear momentum: .. math:: \nabla \cdot \sigma + \rho b = \rho \ddot{u} where :math:u is the displacement vector field, :math:\ddot{u}=\partial^2 u/\partial t^2 is the acceleration, :math:\rho the material density, :math:b a given body force and :math:\sigma the stress tensor which is related to the displacement through a constitutive equation. In the case of isotropic linearized elasticity, one has: .. math:: \sigma =\lambda \text{tr}(\varepsilon)\mathbb{1} + 2\mu\varepsilon where :math:\varepsilon = (\nabla u + (\nabla u)^T)/2 is the linearized strain tensor, :math:\mathbb{1} is the identity of second-rank tensors and :math:\lambda=\dfrac{E\nu}{(1+\nu)(1-2\nu)},\mu=\dfrac{E}{2(1+\nu)} are the Lame coefficients given as functions of the Young modulus :math:E and the Poisson ratio :math:\nu. The weak form is readily obtained by integrating by part the balance equation using a test function :math:v\in V with :math:V being a suitable function space that satisfies the displacement boundary conditions: .. math:: \int_{\Omega} \rho \ddot{u}\cdot v \, {\rm d} x + \int_{\Omega} \sigma(u):\varepsilon(v) \, {\rm d} x = \int_{\Omega} \rho b \cdot v \, {\rm d} x + \int_{\partial\Omega} (\sigma\cdot n) \cdot v \, {\rm d} s \quad \text{for all } v\in V The previous equation can be written as follows: .. math:: \text{Find }u\in V\text{ such that } m(\ddot{u},v) + k(u,v) = L(v) \quad \text{for all } v\in V where :math:m is the symmetric bilinear form associated with the mass matrix and :math:k the one associated with the stiffness matrix. After introducing the finite element space interpolation, one obtains the corresponding discretized evolution equation: .. math:: \text{Find }\{u\}\in\mathbb{R}^n\text{ such that } \{v\}^T[M]\{\ddot{u}\} + \{v\}^T[K]\{u\} = \{v\}^T\{F\} \quad \text{for all } \{v\}\in\mathbb{R}^n which is a generalized :math:n-dof harmonic oscillator equation. Quite often in structural dynamics, structures do not oscillate perfectly but lose energy through various dissipative mechanisms (friction with air or supports, internal dissipation through plasticity, damage, etc.). Dissipative terms can be introduced at the level of the constitutive equation if these mechanisms are well known but quite often it is not the case. Dissipation can then be modeled by adding an *ad hoc* damping term depending on the structure velocity :math:\dot{u} to the previous evolution equation: .. math:: \text{Find }u\in V\text{ such that } m(\ddot{u},v) + c(\dot{u},v) + k(u,v) = L(v) \quad \text{for all } v\in V The damping form will be considered here as bilinear and symmetric, being therefore associated with a damping matrix :math:[C]. ~~~~~~~~~~~~~~~~~ Rayleigh damping ~~~~~~~~~~~~~~~~~ When little is known about the origin of damping in the structure, a popular choice for the damping matrix, known as *Rayleigh damping*, consists in using a linear combination of the mass and stiffness matrix :math:[C] = \eta_M[M]+\eta_K[K] with two positive parameters :math:\eta_M,\eta_K which can be fitted against experimental measures for instance (usually by measuring the damping ratio of two natural modes of vibration). --------------------------------------------------------------- Time discretization using the generalized-:math:\alpha method --------------------------------------------------------------- We now introduce a time discretization of the interval study :math:[0;T] in :math:N+1 time increments :math:t_0=0,t_1,\ldots,t_N,t_{N+1}=T with :math:\Delta t=T/N denoting the time step (supposed constant). The resolution will make use of the generalized-:math:\alpha method which can be seen as an extension of the widely used Newmark-:math:\beta method in structural dynamics. As an implicit method, it is unconditionally stable for a proper choice of coefficients so that quite large time steps can be used. It also allows for high frequency dissipation and offers a second-order accuracy, i.e. in :math:O(\Delta t^2). The method consists in solving the dynamic evolution equation at intermediate time between :math:t_n and :math:t_{n+1} as follows: .. math:: [M]\{\ddot{u}_{n+1-\alpha_m}\} + [C]\{\dot{u}_{n+1-\alpha_f}\}+[K]\{u_{n+1-\alpha_f}\} = \{F(t_{n+1-\alpha_f})\} with the notation :math:X_{n+1-\alpha} = (1-\alpha)X_{n+1}+\alpha X_{n}. In addition, the following approximation for the displacement and velocity at :math:t_{n+1} are used: .. math:: \begin{align*} \{u_{n+1}\} &= \{u_{n}\}+\Delta t \{\dot{u}_{n}\} + \dfrac{\Delta t^2}{2}\left((1-2\beta)\{\ddot{u}_{n}\}+2\beta\{\ddot{u}_{n+1}\}\right) \\ \{\dot{u}_{n+1}\} &= \{\dot{u}_{n}\} + \Delta t\left((1-\gamma)\{\ddot{u}_{n}\}+\gamma\{\ddot{u}_{n+1}\}\right) \end{align*} It can be seen that these are the relations of the Newmark method. The latter is therefore obtained as a particular case when :math:\alpha_f=\alpha_m=0. The problem can then be formulated in terms of unkown displacement at :math:t_{n+1} with: .. math:: \{\ddot{u}_{n+1}\} = \dfrac{1}{\beta\Delta t^2}\left(\{u_{n+1}\} - \{u_{n}\}-\Delta t \{\dot{u}_{n}\} \right) - \dfrac{1-2\beta}{2\beta}\{\ddot{u}_{n}\} After plugging into the evolution and rearranging the known and unknown terms, one obtains the following system to solve: .. math:: \begin{align*} [\bar{K}]\{u_{n+1}\} &= \{F(t_{n+1-\alpha_f})\} - \alpha_f[K]\{u_n\} \\ &- [C](c_1\{u_n\}+c_2\{\dot{u}_n\}+c_3\{\ddot{u}_n\})-[M](m_1\{u_n\}+m_2\{\dot{u}_n\}+m_3\{\ddot{u}_n\}) \end{align*} where: * :math:[\bar{K}] = [K]+c_1[C]+m_1[M] * :math:c_1 = \dfrac{\gamma(1-\alpha_f)}{\beta\Delta t} * :math:c_2 = 1-\gamma(1-\alpha_f)/\beta * :math:c_3 = \Delta t(1-\alpha_f)(1-\dfrac{\gamma}{2\beta}) * :math:m_1 = \dfrac{(1-\alpha_m)}{\beta\Delta t^2} * :math:m_2 = \dfrac{(1-\alpha_m)}{\beta\Delta t} * :math:m_3 = 1-\dfrac{1-\alpha_m}{2\beta} Once the linear system has been solved for :math:\{u_{n+1}\}, the new velocity and acceleration are computed using the previous formulae. ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Popular choice of parameters ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The most popular choice for the parameters is: :math:\alpha_m,\alpha_f \leq 1/2 and :math:\gamma=\dfrac{1}{2}+\alpha_m-\alpha_f, :math:\beta=\dfrac{1}{4}\left(\gamma+\dfrac{1}{2}\right)^2 which ensures unconditional stability, optimal dissipation and second-order accuracy. --------------- Implementation --------------- We consider a rectangular beam clamped at one end and loaded by a uniform vertical traction at the other end. After importing the relevant modules, the mesh and subdomains for boundary conditions are defined:: from dolfin import * import numpy as np import matplotlib.pyplot as plt # Form compiler options parameters["form_compiler"]["cpp_optimize"] = True parameters["form_compiler"]["optimize"] = True # Define mesh mesh = BoxMesh(Point(0., 0., 0.), Point(1., 0.1, 0.04), 60, 10, 5) # Sub domain for clamp at left end def left(x, on_boundary): return near(x, 0.) and on_boundary # Sub domain for rotation at right end def right(x, on_boundary): return near(x, 1.) and on_boundary Material parameters for the elastic constitutive relation, the material density :math:\rho for the mass matrix and the two parameters defining the Rayleigh damping :math:\eta_M,\eta_K (initially zero damping is considered but this value can be changed) are now defined:: # Elastic parameters E = 1000.0 nu = 0.3 mu = Constant(E / (2.0*(1.0 + nu))) lmbda = Constant(E*nu / ((1.0 + nu)*(1.0 - 2.0*nu))) # Mass density rho = Constant(1.0) # Rayleigh damping coefficients eta_m = Constant(0.) eta_k = Constant(0.) Parameters used for the time discretization scheme are now defined. First, the four parameters used by the generalized-:math:\alpha method are chosen. Here, we used the optimal dissipation and second-order accuracy choice for :math:\beta and :math:\gamma, namely :math:\beta=\dfrac{1}{4}\left(\gamma+\dfrac{1}{2}\right)^2 and :math:\gamma=\dfrac{1}{2}+\alpha_m-\alpha_f with :math:\alpha_m=0.2 and :math:\alpha_f=0.4 ensuring unconditional stability:: # Generalized-alpha method parameters alpha_m = Constant(0.2) alpha_f = Constant(0.4) gamma = Constant(0.5+alpha_f-alpha_m) beta = Constant((gamma+0.5)**2/4.) We also define the final time of the interval, the number of time steps and compute the associated time interval between two steps:: # Time-stepping parameters T = 4.0 Nsteps = 50 dt = Constant(T/Nsteps) We now define the time-dependent loading. Body forces are zero and the imposed loading consists of a uniform vertical traction applied at the right extremity. The loading amplitude will vary linearly from :math:0 to :math:p_0=1 over the time interval :math:[0;T_c=T/5], after :math:T_c the loading is removed. For this purpose, we used the following JIT-compiled Expression. In particular, it uses a conditional syntax using operators ? and : :: p0 = 1. cutoff_Tc = T/5 # Define the loading as an expression depending on t p = Expression(("0", "t <= tc ? p0*t/tc : 0", "0"), t=0, tc=cutoff_Tc, p0=p0, degree=0) A standard vectorial :math:P^1 FunctionSpace is now defined for the displacement, velocity and acceleration fields. We also define a tensorial DG-0 FunctionSpace for saving the stress field evolution:: # Define function space for displacement, velocity and acceleration V = VectorFunctionSpace(mesh, "CG", 1) # Define function space for stresses Vsig = TensorFunctionSpace(mesh, "DG", 0) Test and trial functions are defined and the unkown displacement (corresponding to :math:\{u_{n+1}\} for the current time step) will be represented by the Function u. Displacement, velocity and acceleration fields of the previous increment :math:t_n will respectively be represented by functions u_old, v_old and a_old:: # Test and trial functions du = TrialFunction(V) u_ = TestFunction(V) # Current (unknown) displacement u = Function(V, name="Displacement") # Fields from previous time step (displacement, velocity, acceleration) u_old = Function(V) v_old = Function(V) a_old = Function(V) We now use a MeshFunction for distinguishing the different boundaries and mark the right extremity using an AutoSubDomain. The exterior surface measure ds is then defined using the boundary subdomains. Simple Dirichlet boundary conditions are also defined at the left extremity:: # Create mesh function over the cell facets boundary_subdomains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1) boundary_subdomains.set_all(0) force_boundary = AutoSubDomain(right) force_boundary.mark(boundary_subdomains, 3) # Define measure for boundary condition integral dss = ds(subdomain_data=boundary_subdomains) # Set up boundary condition at left end zero = Constant((0.0, 0.0, 0.0)) bc = DirichletBC(V, zero, left) Python functions are now defined to obtain the elastic stress tensor :math:\sigma (linear isotropic elasticity), the bilinear mass and stiffness forms as well as the damping form obtained as a linear combination of the mass and stiffness forms (Rayleigh damping). The linear form corresponding to the work of external forces is also defined:: # Stress tensor def sigma(r): return 2.0*mu*sym(grad(r)) + lmbda*tr(sym(grad(r)))*Identity(len(r)) # Mass form def m(u, u_): return rho*inner(u, u_)*dx # Elastic stiffness form def k(u, u_): return inner(sigma(u), sym(grad(u_)))*dx # Rayleigh damping form def c(u, u_): return eta_m*m(u, u_) + eta_k*k(u, u_) # Work of external forces def Wext(u_): return dot(u_, p)*dss(3) Functions for implementing the time stepping scheme are also defined. update_a returns :math:\{\ddot{u}_{n+1}\} as a function of the variables at the previous increment and of the new displacement :math:\{u_{n+1}\}. The function accepts a keyword ufl so that the expressions involved can be used with UFL representations if True or with array of values if False (we will make use of both possibilities later). In particular, the time step dt and time-stepping scheme parameters are either Constant or floats depending on the case. Function update_v does the same but for the new velocity :math:\{\dot{u}_{n+1}\} as a function of the previous variables and of the new acceleration. Finally, function update_fields performs the final update at the end of the time step when the new displacement :math:\{u_{n+1}\} has effectively been computed. In this context, the new acceleration and velocities are computed using the vector representation of the different fields. The variables keeping track of the values at the previous increment are now assigned the new values computed for the current increment:: # Update formula for acceleration # a = 1/(2*beta)*((u - u0 - v0*dt)/(0.5*dt*dt) - (1-2*beta)*a0) def update_a(u, u_old, v_old, a_old, ufl=True): if ufl: dt_ = dt beta_ = beta else: dt_ = float(dt) beta_ = float(beta) return (u-u_old-dt_*v_old)/beta_/dt_**2 - (1-2*beta_)/2/beta_*a_old # Update formula for velocity # v = dt * ((1-gamma)*a0 + gamma*a) + v0 def update_v(a, u_old, v_old, a_old, ufl=True): if ufl: dt_ = dt gamma_ = gamma else: dt_ = float(dt) gamma_ = float(gamma) return v_old + dt_*((1-gamma_)*a_old + gamma_*a) def update_fields(u, u_old, v_old, a_old): """Update fields at the end of each time step.""" # Get vectors (references) u_vec, u0_vec = u.vector(), u_old.vector() v0_vec, a0_vec = v_old.vector(), a_old.vector() # use update functions using vector arguments a_vec = update_a(u_vec, u0_vec, v0_vec, a0_vec, ufl=False) v_vec = update_v(a_vec, u0_vec, v0_vec, a0_vec, ufl=False) # Update (u_old <- u) v_old.vector()[:], a_old.vector()[:] = v_vec, a_vec u_old.vector()[:] = u.vector() The system variational form is now built by expressing the new acceleration :math:\{\ddot{u}_{n+1}\} as a function of the TrialFunction du using update_a, which here works as a UFL expression. Using this new acceleration, the same is done for the new velocity using update_v. Intermediate averages using parameters :math:\alpha_m,\alpha_f of the generalized- :math:\alpha method are obtained with a user-defined fuction avg. The weak form evolution equation is then written using all these quantities. Since the problem is linear, we then extract the bilinear and linear parts using rhs and lhs:: def avg(x_old, x_new, alpha): return alpha*x_old + (1-alpha)*x_new # Residual a_new = update_a(du, u_old, v_old, a_old, ufl=True) v_new = update_v(a_new, u_old, v_old, a_old, ufl=True) res = m(avg(a_old, a_new, alpha_m), u_) + c(avg(v_old, v_new, alpha_f), u_) \ + k(avg(u_old, du, alpha_f), u_) - Wext(u_) a_form = lhs(res) L_form = rhs(res) Alternatively, the use of derivative can be made for non-linear problems for instance or one can also directly formulate the system to solve, involving the modified stiffness matrix :math:[\bar{K}] and the various coefficients introduced earlier. Since the system matrix to solve is the same for each time step (constant time step), it is not necessary to factorize the system at each increment. It can be done once and for all and only perform assembly of the varying right-hand side and backsubstitution to obtain the solution much more efficiently. This is done by defining a LUSolver object and asking for reusing the matrix factorization:: # Define solver for reusing factorization