计算科学 - 弹性物理学中使用哪些数值方法来模拟变形？ - 吾爱随笔录

弹性物理学中使用哪些数值方法来模拟变形？

计算科学有限元 pde 计算物理学数值建模

2021-12-02 21:10:51

弹性物理学中使用哪些数值方法来模拟变形？例如，下面是 Ansys 中的超弹性变形示例：

也许比超弹性更简单，对于线性弹性，我们有以下方程：

\nabla \cdot σ + F = ρ \ddot{u} ε = \frac{1}{2} [\nabla u + (\nabla u)^{T}] σ = C : ε

$\nabla\cdot\sigma + {F} = \rho\ddot{{u}}\\ {\varepsilon} =\tfrac{1}{2} \left[{\nabla}{u}+({ \nabla}{u})^\mathrm{T}\right]\\ {\sigma} = {C}:{\varepsilon}$

在哪里

$\sigma$ - 应力张量
$\epsilon$ - 应变张量
$u$ - 位移
$C$ - 刚度张量

假设我们在域的网格上应用有限元并结合 Runge-Kutta 方法来处理时间，我们可以求解上述方程并找到表示变形的解 $u$ 。然而，变形似乎意味着某些东西需要移动，直到此时我们有一个静态的、不动的域网格。在这种情况下，什么动作？网格？

更一般地说，用于模拟弹性材料的运动和变形的一般算法是什么，类似于上面的模拟所示？

3个回答

似乎算法的类型有很大差异，具体取决于问题是否是：

准静态弹性或
超弹性

在准静态弹性情况下，一个简单的方法如下：作为每个时间步的第一部分，计算位移场。由于位移现在在网格的每个节点上都是已知的，因此节点可以根据作为每个时间步的第二部分进行位移。 $u$ $u$ $u$

作为deal.ii有限元库的教程案例之一，给出了C++代码示例对此类问题的详细描述：链接

deal.ii 还为超弹性问题提供了一个具体示例。但是，那里的问题描述太复杂了，我不想在这里详细重复。作为一个简短的总结，他们使用三场变分原理和应力张量的一组欧拉-拉格朗日方程。 $u$ $\sigma$

作为深入参考，他们建议：

A. Holzapfel (2001)，非线性固体力学。工程的连续方法，John Wiley & Sons。国际标准书号：0-471-82304-X

如果在那里消除运动，则环在可变载荷下的变形仍然存在。这可以在添加瑞利阻尼的弹性体模型中计算。动画显示了使用 FEM 和 Mathematica 12 计算的橡胶环的变形。

用于计算动画的附录代码（取自教程并针对本案例进行了修改）

Needs["NDSolve`FEM`"]; \[CapitalOmega] = 
 ImplicitRegion[4^2 <= x^2 + y^2 <= 5^2, {x, y}];
mesh = ToElementMesh[\[CapitalOmega], {{-5, 5}, {-5, 5}}, 
   "MaxCellMeasure" -> 0.03];
mesh["Wireframe"]

diffusionCoefficients = 
 "DiffusionCoefficients" -> {{{{-(Y/(1 - \[Nu]^2)), 
       0}, {0, -((Y (1 - \[Nu]))/(2 (1 - \[Nu]^2)))}}, {{0, -((
        Y \[Nu])/(1 - \[Nu]^2))}, {-((Y (1 - \[Nu]))/(
        2 (1 - \[Nu]^2))), 
       0}}}, {{{0, -((Y (1 - \[Nu]))/(2 (1 - \[Nu]^2)))}, {-((
        Y \[Nu])/(1 - \[Nu]^2)), 
       0}}, {{-((Y (1 - \[Nu]))/(2 (1 - \[Nu]^2))), 
       0}, {0, -(Y/(1 - \[Nu]^2))}}}} /. {Y -> 10^2, \[Nu] -> 
    33/100}; massCoefficients = "MassCoefficients" -> {{1, 0}, {0, 1}};
vd = NDSolve`VariableData[{"Time", "DependentVariables", 
     "Space"} -> {t, {u, v}, {x, y}}];
sd = NDSolve`SolutionData[{"Time", "Space"} -> {0., 
     ToNumericalRegion[mesh]}];

methodData = InitializePDEMethodData[vd, sd];

loadCoefficients = "LoadCoefficients" -> {{0}, {0}};

"LoadCoefficients" -> {{0}, {0}};

Subscript[\[CapitalGamma], Nv] = 
 NeumannValue[-3 Sin[Pi t] Sign[y], -1 <= x <= 
   1]; Subscript[\[CapitalGamma], Du] = 
 DirichletCondition[u[x, y] == 0, 
  x == 0]; Subscript[\[CapitalGamma], Nu] = 
 NeumannValue[0, -1 <= x <= 1];
Subscript[\[CapitalGamma], Dv] = 
  DirichletCondition[v[x, y] == 0, y == 0];

initCoeffs = 
 InitializePDECoefficients[vd, 
  sd, {diffusionCoefficients, massCoefficients, loadCoefficients}];

initBCs = 
 InitializeBoundaryConditions[vd, 
  sd, {{Subscript[\[CapitalGamma], Du], Subscript[\[CapitalGamma], 
    Nu]}, {Subscript[\[CapitalGamma], Dv], Subscript[\[CapitalGamma], 
    Nv]}}];

sdpde = DiscretizePDE[initCoeffs, methodData, sd, "Stationary"];
sbcs = DiscretizeBoundaryConditions[initBCs, methodData, sd, 
  "Stationary"];

rhs[t_?NumericQ, uv_, duv_] := 
 Module[{l, s, d, m, tdpde, tbcs, rayleighDamping},
  NDSolve`SetSolutionDataComponent[sd, "Time", t];
  {l, s, d, m} = sdpde["SystemMatrices"];
  tdpde = DiscretizePDE[initCoeffs, methodData, sd, "Transient"];
  tbcs = DiscretizeBoundaryConditions[initBCs, methodData, sd, 
    "Transient"];
  {l, s, d, m} += tdpde["SystemMatrices"];
  rayleighDamping = 0.1*m + 0.04*s;
  DeployBoundaryConditions[{l, s, rayleighDamping, m}, tbcs];
  DeployBoundaryConditions[{l, s, rayleighDamping, m}, sbcs];
  l - s . uv - rayleighDamping . duv
  ]

dof = methodData["DegreesOfFreedom"];
init = dinit = ConstantArray[0, {dof, 1}];

mass = sdpde["MassMatrix"];
stiff = sdpde["StiffnessMatrix"];
rd = 0.1*mass + 0.04*stiff;

sparsity = 
  ArrayFlatten[{{mass["PatternArray"], 
     mass["PatternArray"]}, {rd["PatternArray"], rd["PatternArray"]}}];

Dynamic["time: " <> ToString[CForm[currentTime]]]
AbsoluteTiming[
 tfun = NDSolveValue[{
    mass . uv''[ t] == rhs[t, uv[t], uv'[t]]
    , uv[ 0] == init, uv'[ 0] == dinit}, uv, {t, 0, 2}
   , Method -> {"EquationSimplification" -> "Residual"}
   , Jacobian -> {Automatic, Sparse -> sparsity}
   , EvaluationMonitor :> (currentTime = t;)
   ]]

可视化

split = Span @@@ 
  Transpose[{Most[# + 1], Rest[#]} &[methodData["IncidentOffsets"]]];

graphics = Function[t,
    dmesh = ElementMeshDeformation[mesh, Part[tfun[t], #] & /@ split];
    Show[{
      mesh["Wireframe"["MeshElement" -> "BoundaryElements"]],
      dmesh[
       "Wireframe"[
        "ElementMeshDirective" -> 
         Directive[EdgeForm[Red], FaceForm[]]]]
      }, PlotRange -> {{-6, 6}, {-5.5, 5.5}}]] /@ Range[0, 2, .1];

ListAnimate[graphics, SaveDefinitions -> True]

TL，DR：要么是旧的 Galerkin 有限元方法，要么是无网格/粒子方法。

这里有一些东西要解压。首先，您展示的模拟包括弹性体（环）和硬边界之间的接触，因此约束是非完整的。接触问题比仅具有完整约束的重力下的弹性变形更具挑战性。

在您上面显示的动画中，似乎唯一发生变化的是网格节点的位置。当材料参考配置的变形很小时，这是一种合理的方法，这正是弹性方程为线性时。对于非常大的变形，身体配置中的变换网格可能会变得缠结，正如您可能想象的那样，这是个坏消息。有几种方法可以解决这个问题。

首先，您可以通过定期重新划分网格来更改模拟的内部运作。当三角形变得过于退化时，您可以翻转边缘以恢复网格质量。这种方法可以消除网格缠结的可能性，但很难正确编码。

其次，您可以首先更改您正在模拟的物理问题。如果变形那么大，使用线性方程组不再是描述问题的好方法。非线性弹性，例如参见Neo-Hookean solids，包括当参照体变换的雅可比变得奇异时趋于无穷大的项。这可以防止数学层面的奇怪变形，但也使问题的数值解更具挑战性。

这两个想法都非常符合关于单纯网格的 Galerkin 有限元方法的思想。与这种思维方式更大的不同是使用无网格方法。这包括平滑粒子流体动力学、质点法等。这些通常在计算上更昂贵，因为您需要在每个时间步上进行邻居查询，这很难使缓存友好。但是它们可以更好地适应非常变形或移动的几何形状。

其它你可能感兴趣的问题

上一篇生成随机平滑的二维闭合曲线下一篇Python中具有未知常数的数值积分