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

原创

已于 2022-03-04 15:35:44 修改 · 3.4k 阅读

12 ·

CC 4.0 BY-SA版权

文章标签：

#数学 #物理学

于 2021-04-17 17:27:59 首次发布

本文介绍了如何使用广义α方法解决时变动力学问题，详细阐述了控制方程的变分形式，时间离散过程，并提供了数值算例和FEniCS代码实现。通过这种方法，可以更稳定地处理动态问题，避免直接震荡，确保数值计算的准确性。

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

控制方程及其变分形式

首先我们定义我们需要求解的弹性力学方程
$\nabla\cdot\textcolor{red}{\sigma}+\rho b = \rho \ddot{u} \tag{1}$
其中 $σ\sigma$ 是应力张量， $b$ 是体元外力， $ρ\rho$ 是密度， $u¨\ddot{u}$ 表示位移 $u$ 对时间的 $t$ 的二阶导数，即加速度。在这个问题中，我们需要求解的是位移 $u$ 。应力张量 $σ\sigma$ 的表达式如下
$\sigma = \lambda\mathrm{tr}(\epsilon)I+2\mu\epsilon \\ \epsilon=\frac{\nabla u+(\nabla u)^T}{2} \\ \lambda=\frac{E\nu}{(1+\nu)(1-2\nu)},\mu=\frac{E}{2(1+\nu)} \tag{2}$
对于式(2)中的参数，杨氏模量 $E = 1000$ ，泊松比 $ν=0.3\nu=0.3$ ，密度 $ρ=1\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 v dx-\int_\Omega\textcolor{red}{(\nabla\cdot\sigma)\cdot v} dx=\int_\Omega\rho b\cdot v dx \tag{3}$
这里我们把待求量 $u$ 放左边，已知量放右边。对于式(3)中的部分我们进一步展开
$\begin{aligned} (\nabla\cdot\sigma)\cdot v&=\frac{\partial \sigma_{ij}}{\partial x_i}v_j\\ &=\frac{\partial}{\partial x_i}(\sigma_{ij}v_j)-\textcolor{red}{\sigma_{ij}\frac{\partial v_j}{\partial x_i}}\\ &=\nabla \cdot(\sigma\cdot v)-\sigma \cdot \nabla v \end{aligned} \tag{4}$
式(4)中，最后一步的第二项可以写成
$\sigma \cdot \nabla v=\sigma \cdot (\nabla v)^T=\sigma(u) \cdot \epsilon(v) \tag{5}$

最低0.47元/天解锁文章