# Time-integration of elastodynamics equation

This demo is implemented in a single Python file, `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-\(\alpha\) method [ERL2002]. In particular it demonstrates how to:

- formulate mass and damping forms of the elastodynamics equation

- implement the generalized-\(\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:

## Introduction and elastodynamics equation

The elastodynamics equation combine the balance of linear momentum:

\[\nabla \cdot \sigma + \rho b = \rho \ddot{u}\]

where \(u\) is the displacement vector field, \(\ddot{u}=\partial^2 u/\partial t^2\) is the acceleration, \(\rho\) the material density, \(b\) a given body force and \(\sigma\) the stress tensor which is related to the displacement through a constitutive equation. In the case of isotropic linearized elasticity, one has:

\[\sigma =\lambda \text{tr}(\varepsilon)\mathbb{1} + 2\mu\varepsilon\]

where \(\varepsilon = (\nabla u + (\nabla u)^T)/2\) is the linearized strain tensor, \(\mathbb{1}\) is the identity of second-rank tensors and \(\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 \(E\) and the Poisson ratio \(\nu\).

The weak form is readily obtained by integrating by part the balance equation using a test function \(v\in V\) with \(V\) being a suitable function space that satisfies the displacement boundary conditions:

\[\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:

\[\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 \(m\) is the symmetric bilinear form associated with the mass matrix and \(k\) the one associated with the stiffness matrix.

After introducing the finite element space interpolation, one obtains the corresponding discretized evolution equation:

\[\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 \(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 \(\dot{u}\) to the previous evolution equation:

\[\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 \([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 \([C] = \eta_M[M]+\eta_K[K]\) with two positive parameters \(\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-\(\alpha\) method

We now introduce a time discretization of the interval study \([0;T]\) in \(N+1\) time increments \(t_0=0,t_1,\ldots,t_N,t_{N+1}=T\) with \(\Delta t=T/N\) denoting the time step (supposed constant). The resolution will make use of the generalized-\(\alpha\) method which can be seen as an extension of the widely used Newmark-\(\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 \(O(\Delta t^2)\).

The method consists in solving the dynamic evolution equation at intermediate time between \(t_n\) and \(t_{n+1}\) as follows:

\[[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 \(X_{n+1-\alpha} = (1-\alpha)X_{n+1}+\alpha X_{n}\). In addition, the following approximation for the displacement and velocity at \(t_{n+1}\) are used:

\[\begin{split}\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*}\end{split}\]

It can be seen that these are the relations of the Newmark method. The latter is therefore obtained as a particular case when \(\alpha_f=\alpha_m=0\).

The problem can then be formulated in terms of unkown displacement at \(t_{n+1}\) with:

\[\{\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:

\[\begin{split}\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*}\end{split}\]

where:

  - \([\bar{K}] = [K]+c_1[C]+m_1[M]\)

  - \(c_1 = \dfrac{\gamma(1-\alpha_f)}{\beta\Delta t}\)

  - \(c_2 = 1-\gamma(1-\alpha_f)/\beta\)

  - \(c_3 = \Delta t(1-\alpha_f)(1-\dfrac{\gamma}{2\beta})\)

  - \(m_1 = \dfrac{(1-\alpha_m)}{\beta\Delta t^2}\)

  - \(m_2 = \dfrac{(1-\alpha_m)}{\beta\Delta t}\)

  - \(m_3 = 1-\dfrac{1-\alpha_m}{2\beta}\)

Once the linear system has been solved for \(\{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: \(\alpha_m,\alpha_f \leq 1/2\) and \(\gamma=\dfrac{1}{2}+\alpha_m-\alpha_f\), \(\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], 0.) and on_boundary

    # Sub domain for rotation at right end
    def right(x, on_boundary):
        return near(x[0], 1.) and on_boundary

Material parameters for the elastic constitutive relation, the material density \(\rho\) for the mass matrix and the two parameters defining the Rayleigh damping \(\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-\(\alpha\) method are chosen. Here, we used the optimal dissipation and second-order accuracy choice for \(\beta\) and \(\gamma\), namely \(\beta=\dfrac{1}{4}\left(\gamma+\dfrac{1}{2}\right)^2\) and \(\gamma=\dfrac{1}{2}+\alpha_m-\alpha_f\) with \(\alpha_m=0.2\) and \(\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 \(0\) to \(p_0=1\) over the time interval \([0;T_c=T/5]\), after \(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 \(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 \(\{u_{n+1}\}\) for the current time step) will be represented by the Function `u`. Displacement, velocity and acceleration fields of the previous increment \(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 \(\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 \(\{\ddot{u}_{n+1}\}\) as a function of the variables at the previous increment and of the new displacement \(\{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 \(\{\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 \(\{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 \(\{\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 \(\alpha_m,\alpha_f\) of the generalized- \(\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 \([\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 while PETSc handles caching factorizations:

    # Define solver for reusing factorization
    K, res = assemble_system(a_form, L_form, bc)
    solver = LUSolver(K, "mumps")
    solver.parameters["symmetric"] = True

We now initiate the time stepping loop. We will keep track of the beam vertical tip displacement over time as well as the different parts of the system total energy. We will also compute the stress field and save it, along with the displacement field, in a `XDMFFile`. The option flush_ouput enables to open the result file before the loop is finished, the `function_share_mesh` option tells that only one mesh is used for all functions of a given time step (displacement and stress) while the `rewrite_function_mesh` enforces that the same mesh is used for all time steps. These two options enables writing the mesh information only once instead of \(2N_{steps}\) times:

    # Time-stepping
    time = np.linspace(0, T, Nsteps+1)
    u_tip = np.zeros((Nsteps+1,))
    energies = np.zeros((Nsteps+1, 4))
    E_damp = 0
    E_ext = 0
    sig = Function(Vsig, name="sigma")
    xdmf_file = XDMFFile("elastodynamics-results.xdmf")
    xdmf_file.parameters["flush_output"] = True
    xdmf_file.parameters["functions_share_mesh"] = True
    xdmf_file.parameters["rewrite_function_mesh"] = False

The time loop is now started, the loading is first evaluated at \(t=t_{n+1-\alpha_f}\). The corresponding system right-hand side is then assembled and the system is solved. The different fields are then updated with the newly computed quantities. Finally, some post-processing is performed: stresses are computed and written to the result file and the tip displacement and the different energies are recorded:

    def local_project(v, V, u=None):
        """Element-wise projection using LocalSolver"""
        dv = TrialFunction(V)
        v_ = TestFunction(V)
        a_proj = inner(dv, v_)*dx
        b_proj = inner(v, v_)*dx
        solver = LocalSolver(a_proj, b_proj)
        solver.factorize()
        if u is None:
            u = Function(V)
            solver.solve_local_rhs(u)
            return u
        else:
            solver.solve_local_rhs(u)
            return

    for (i, dt) in enumerate(np.diff(time)):

        t = time[i+1]
        print("Time: ", t)

        # Forces are evaluated at t_{n+1-alpha_f}=t_{n+1}-alpha_f*dt
        p.t = t-float(alpha_f*dt)

        # Solve for new displacement
        res = assemble(L_form)
        bc.apply(res)
        solver.solve(K, u.vector(), res)

        # Update old fields with new quantities
        update_fields(u, u_old, v_old, a_old)

        # Save solution to XDMF format
        xdmf_file.write(u, t)

        # Compute stresses and save to file
        local_project(sigma(u), Vsig, sig)
        xdmf_file.write(sig, t)

        p.t = t
        # Record tip displacement and compute energies
        # Note: Only works in serial
        if MPI.comm_world.size == 1:
            u_tip[i+1] = u(1., 0.05, 0.)[1]
        E_elas = assemble(0.5*k(u_old, u_old))
        E_kin = assemble(0.5*m(v_old, v_old))
        E_damp += dt*assemble(c(v_old, v_old))
        # E_ext += assemble(Wext(u-u_old))
        E_tot = E_elas+E_kin+E_damp #-E_ext
        energies[i+1, :] = np.array([E_elas, E_kin, E_damp, E_tot])

Note that in the above, the stresses are computed using a `LocalSolver` through the `local_project` function. Since the stress function space is a DG-0 space, the projection on this space can be performed element-wise in a very efficient manner. We therefore take advantage of the `LocalSolver` functionality which is precisely dedicated to such situations. Since this projection is performed at each time step, the savings in terms of computing time can be quite important.

As regards the computation of the various energies, the elastic and kinetic energies are respectively given by:

\[E_{elas} = \int_{\Omega} \dfrac{1}{2}\sigma(u):\varepsilon(u) \, {\rm d} x\]
\[E_{kin} = \int_{\Omega} \dfrac{1}{2}\rho \dot{u}\cdot\dot{u} \, {\rm d} x\]

which are readily computed from the respective stiffness and mass forms `k` and `m` and the current displacement and velocity. The energy related to damping is computed from the corresponding dissipation term \(\mathcal{D}=c(\dot{u},\dot{u})\) and integrated over time:

\[E_{damp} = \int_0^T \mathcal{D} \, {\rm d} t\]

As for the work developed by the external forces, the contribution to the energy is added at each time step. Finally, the total energy of the sytem is given by:

\[E_{tot} = E_{elas}+E_{kin}+E_{damp}-E_{ext}\]

When the time evolution loop is finished, the evolution of the tip displacement as well as the different contributions of the energy are plotted as functions of time:

    if MPI.comm_world.size == 1:
        # Plot tip displacement evolution
        plt.figure()
        plt.plot(time, u_tip)
        plt.xlabel("Time")
        plt.ylabel("Tip displacement")
        plt.ylim(-0.5, 0.5)
        plt.show()

    if (MPI.comm_world.rank == 0):
        # Plot energies evolution
        plt.figure()
        plt.plot(time, energies)
        plt.legend(("elastic", "kinetic", "damping", "total"))
        plt.xlabel("Time")
        plt.ylabel("Energies")
        plt.ylim(0, 0.0011)
        plt.show()

## Analyzing the results

We first consider the case of zero Rayleigh damping \(\eta_M=\eta_K=0\). In this case, it can be observed that the evolution of the total energy depends on the choice of the time-stepping scheme parameters. With \(\alpha_m=\alpha_f=0\), we recover the Newmark-\(\beta\) method with \(\beta=0.25,\gamma=0.5\). This scheme is known for being conservative. This can be observed (figure-left) in the constant total energy for \(t\geq T_c\) when the loading is removed. On the contrary, for non zero alpha parameters, e.g. \(\alpha_m=0.2, \alpha_f=0.4\), it can be observed (figure-right) that the energy is decreasing during this phase, indicating numerical damping. For both cases, the scheme is unconditionally stable. Moreover, these differences vanish when reducing the time step.

Newmark-\(\beta\) method \(\alpha_m=\alpha_f=0\)

Generalized-\(\alpha\) with \(\alpha_m=0.2, \alpha_f=0.4\)

For non-zero Rayleigh damping \(\eta_M=\eta_K=0.01\), the total energy including viscous dissipation tends to oscillate around a constant value, with oscillations vanishing for decreasing time steps.

## References

[ERL2002]

Silvano Erlicher, Luca Bonaventura, Oreste Bursi. The analysis of the Generalized-alpha method for non-linear dynamic problems. Computational Mechanics, Springer Verlag, 2002, 28, pp.83-104, doi:10.1007/s00466-001-0273-z