使用广义α方法(the generalized-α method)求解时变动力学问题

本文讲解如何使用使用广义α方法求解时变动力学问题，附带FEniCS代码

控制方程及其变分形式

首先我们定义我们需要求解的弹性力学方程
$$\nabla \cdot \sigma +\rho b=\rho \ddot{u}\text{\hspace{1em}(1)}$$
其中$\sigma$是应力张量，$b$是体元外力，$\rho$是密度，$\ddot{u}$表示位移$u$对时间的$t$的二阶导数，即加速度。在这个问题中，我们需要求解的是位移$u$。应力张量$\sigma$的表达式如下
$$\begin{gathered} \sigma =\lambda \mathrm{t}\mathrm{r}(ϵ)I+2\mu ϵ \\ ϵ=\frac{\nabla u+(\nabla u{)}^T}{2} \\ \lambda =\frac{E\nu }{(1+\nu )(1-2\nu )},\mu =\frac{E}{2(1+\nu )}\text{\hspace{1em}(2)} \end{gathered}$$
对于式(2)中的参数，杨氏模量$E=1000$，泊松比$\nu =0.3$，密度$\rho =1$，$I$是单位矩阵

上述的过程使用FEniCS实现如下

# 弹性参数
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)))
# 质量密度
rho = Constant(1.0)

# ε
def epsilon(u_):
    return sym(grad(u_))
# 应力张量σ
def sigma(u_):
    return 2.0*mu*epsilon(u_) + lmbda*tr(epsilon(u_))*Identity(len(u_))

对于式(1)我们写出其变分形式
$${\int }_{\Omega }\rho \ddot{u}\cdot vdx-{\int }_{\Omega }(\nabla \cdot \sigma )\cdot vdx={\int }_{\Omega }\rho b\cdot vdx\text{\hspace{1em}(3)}$$
这里我们把待求量$u$放左边，已知量放右边。对于式(3)中的部分我们进一步展开
$$\begin{array}{rl}(\nabla \cdot \sigma )\cdot v& =\frac{\partial {\sigma }_{ij}}{\partial x_i}{v}_j\\ & =\frac{\partial }{\partial x_i}({\sigma }_{ij}{v}_j)-{\sigma }_{ij}\frac{\partial v_j}{\partial x_i}\\ & =\nabla \cdot (\sigma \cdot v)-\sigma \cdot \nabla v\end{array}\text{\hspace{1em}(4)}$$
式(4)中，最后一步的第二项可以写成
σ ⋅ ∇ v = σ ⋅ ( ∇ v ) T =