文 | 夜寻张

编辑 | 夜寻张

前言

近年来，固体力学在工程领域中的应用日益广泛。为了更好地理解和分析复杂的结构和材料行为，研究人员发展了各种数值模拟方法。

控制方程、固体应力求解器是我们常用的工具，它基于有限元理论和计算力学原理，通过求解力平衡方程来获得结构物内部的应力分布。在实际工程问题中，固体的力学响应往往与周围流体的流动耦合在一起，这就需要考虑流固耦合问题。

为了准确描述物质的动态响应和流体的变形以及相互作用，我们将动态网格演化方法引入到数值模拟中。

为了确保数值模拟结果的可靠性和精度，网格收敛性也是该领域关注的重要问题。

控制方程

控制固体力学的基本定律与控制流体力学的定律相同。对于我们的问题，这些定律是质量守恒定律和动量守恒定律：

质量守恒定律：通过控制体的质量平衡由连续方程表示。

动量守恒定律：控制体内的线动量平衡由动量方程表示。

图1

ρ是密度，v是速度。该方程简化为不可压缩流体的∇v=0，如方程所示。

线性动量守恒：线性动量变化的时间速率等于作用于连续体上的合力。

图2

⊗表示外积，fb是作用在连续体上单位体积的体力（重力、离心力等）的合力，σ是应力张量。

应力张量可以表示为法向应力和剪切应力之和，σ = -pI + T，其中p是压力，T是黏性应力张量，从而得到方程，并加入额外的体积力项。

图3

流体和固体求解器之间的方程的主要区别是T项是如何近似的。

流体：用牛顿流体假设和Stokes假设来确定T为：

图4

固体：除了T项的差异外，固体力学求解器经常在控制方程中使用拉格朗日公式，其中对流项消失，速度的时间导数减少为位移的偏二阶导数。

图5

u为实体的位移，b为单位体积的合力，σ为柯西应力张量。我们需要一个σ和u之间的关系来求解这些方程。线弹性材料假设适用于小应变，假设的应力-应变关系为：

图6

其中D称为材料切线，ε为应力，定义为ε=∂u∂x。考虑到一种各向同性的弹性材料，方程中的张量可以用Voigt符号写成：

图7

与D仅依赖于杨氏模量E和泊松比ν，根据以下关系：

图8

使用STARCCM+ v13.04软件进行流固耦合问题的数值求解。

固体应力求解器

对于翅膀变形的求解，我们使用固体应力求解器，并考虑各向同性的线性弹性定律。如果将一个弹性材料加载然后卸载，它会恢复到原始的非变形的形状。变形不依赖于加载过程或加载速率，而是材料本身不表现出阻尼行为。

在解决柔性拍动翅膀的变形问题时，这种模型似乎是合适的。因为气动和惯性力相对较小，翅膀的预期变形仅在弦长的5%到15%之间。

解决固体问题所需的材料特性包括密度ρs、杨氏模量E和泊松比ν。在STARCCM+中的固体应力模型允许使用有限元方法计算承受载荷下的固体结构的位移。

选择一阶精确的时间离散化方法，使用后向欧拉方法近似计算加速度u n 和速度u n 的时间步n。求解器寻找满足上述控制方程有限元离散化形式的位移场u。

流固耦合

实现分割解法的常用方案是：传统串行交错循环（CSS），在将结构域和流体域都进行同样时间步大小∆t的离散化之后，流体和固体之间的界面Γ处从时间水平n到n+1的过程按照以下步骤进行。

图9：对于一个分区的FSI方法的CSS循环

流固耦合过程通常采用串行交错循环（CSS）算法来实现。下面是CSS算法的步骤：

1、将结构体的速度和位移（预测值）传递到流体边界并更新流体网格。

2、使用新的边界信息解决下一个时间步长的流体子系统问题。

3、从流体场解中提取相关部分，并将其转换为新的加载数据以传递到结构边界。

4、计算结构子系统的响应并演化到下一个时间步长。

在界面上执行的CSS循环过程，G上标表示流体网格节点的速度，Γ下标表示界面边界上的流体速度。

松散耦合算法中，循环只执行一次，仅在时间上达到一阶精度。而紧密耦合算法使用元素求解器的收敛准则。

步骤1、2和3可能会重复执行，直到满足收敛停止准则。这种子迭代方案的实现可以得到二阶时间精度的解。

在STARCCM+中，强耦合物理问题需要使用流固耦合模型同时求解流动和结构，该模型允许在流固界面进行数据交换。

此求解器使用固定点迭代方法对流固界面处的固体位移进行松弛。可以指定位移容差以评估时间步长内的收敛性。

图10：使用CSS算法的接口分辨率的示意图

动态网格演化

在流体中运动的结构需要CFD求解器考虑固体结构的形状和位置的变化。我们使用Morphing来实现。

Morphing是指流体网格的变形。通过移动流体顶点，使其符合固体结构并保持合理质量的流体网格。这个过程通常被称为拓扑恒定操作。所有的单元都保持相同的邻居，单元的形状也可以随时间改变。

流体输运方程被求解，在考虑网格的任意运动时。STARCCM+中的Morpher目前使用基于径向基函数的多重四边形Morphing模型来定义内部顶点的运动。

网格和边界条件

模型中还引入了一个新的固体区域，由两个部分组成：

"Mount"：它由一个固体部分组成，用于支撑机翼并将旋转和旋转运动传递给引擎。除了被定义为无滑移壁面外，其所有表面还都受到固定条件的约束，该部分只执行规定的运动而不发生任何变形。

"Wing"：它由执行振动运动的机翼组成。与"Mount"接触的表面被设置为固体/固体界面，流体域接触的表面被设置为固体/流体界面。所有表面都被定义为无滑移壁面。

这两个部分都使用四面体单元进行网格划分，网格尺寸使得机翼弦长方向有50个单元。这导致固体域由92万个单元组成，流体域后总共约为600万个单元。

固体物理连续体允许我们定义机翼的材料属性，可以修改杨氏模量、泊松比和机翼密度以研究其对振动运动性能的影响。

时间步长

如刚性机翼数值模型中所示，时间离散化更倾向于使用非均匀时间步长，以减少计算时间同时保持描述运动和获得的结果的准确性。由于机翼的变形，机翼与流体之间的相对速度可能高于一般的旋转速度，从而增加了每个时间步长内移动的单元数。

为了保持准确的结果，需要减小全局时间步长的量级。使用相同的非均匀时间离散化，但全局时间步长缩小，当没有加速度发生时的时间步长为dt = T/500，最大迭代次数为每周期1500次。

图11

网格收敛性

遵循与刚性机翼相同的过程，以获得网格收敛性。参考基准案例包含大约6*106个单元，在机翼表面每个弦长有50个单元。

将测试更细和更粗的网格以验证网格收敛性。还测试了不同的时间步长离散化。水代替空气成为环绕机翼的流体，以模拟实验条件。这样可以通过将数值模型提供的力与实验获取的力进行比较来验证数值模型的力。

下表列出了用于网格收敛性测试的不同配置列表。

图12：配置测试的FSI模型的收敛性测试。

图13

配置下模拟的10个周期中，机翼升力系数（CL）的平均演变。“FineFlex”配置由于计算成本较高，仅得到了3个周期的结果。

减小实验时间对结果没有显著影响，但对于具有更高加速度率的关键情况可能存在一些差异，正如刚性模型所示。“Non-Uniform+”时间步长被排除在外，只使用参考或缩小的时间步长。

结果显示随着网格细化，平均CL值减小。刚性模型的研究表明，流体求解器在粗网格配置下提供的结果与真实值之间的误差在5%的范围内。

这表明平均CL值的差异可能来自结构求解器，应该分析数值模型获得的结构变形。

图14：FSI数值模型的网格尺寸收敛性。

四个研究案例展示了下机翼端点位移的情况，粗网格低估了机翼的变形，而趋于与“FineFlex”配置一致。

由于不同案例下机翼的形状不同，这影响了升力产生。随着网格细化，机翼的变形增加，导致升力系数（CL）减小。

在这种参考运动案例过程中，由于机翼的结构特性而引起的变形对升力产生有负面影响。对于其他运动案例或提供更适度灵活性的其他结构特性，情况可能不同。

流体领域的网格尺寸使用了“FineFlex”案例中的相同值进行改变。

我们的目标不是评估机翼的柔性效应，而是验证柔性模型提供的结果。

为了确保流体求解器仍然在误差5%的范围内提供结果，类似于刚性模型，进行了额外的网格配置测试。对于这些新案例，结构领域的网格尺寸保持不变，改变的是流体领域的网格尺寸。

图15：所研究的构型的翼尖位移演变

图16：所研究的构型的瞬时CL演化。

通过使用“CoarseFlex”、“ModerateFlex”和“BaseFlex”案例的值来改变流体领域的网格尺寸，以确定结构求解器是否需要比流体求解器更精细的网格才能收敛到准确的解。下表总结了用于网格收敛性测试的额外配置。

图17：附加配置测试的收敛测试。

根据对额外案例的测试结果，在所有测试的流体领域网格尺寸上，保持结构领域的精细网格明显提高了机翼端点位移的解决方案。

“FineStrucCoarFlow”案例由于领域之间的网格尺寸差异较大，提供的结果不太准确。这种差异在在领域界面处插值解决方案时会产生一些小误差，用来将一个求解器的结果传递给另一个求解器。

其他两个案例提供的结果更接近“FineFlex”案例的结果，显示了流体求解器的网格收敛性。

“FineStrucBaseFlow”或“FineStrucModFlow”案例可能提供比“BaseFlex”配置更准确的结果，也不会显著增加总网格尺寸和计算成本，就像使用“FineFlex”案例那样。

为了确定实验机翼的材料，设定了期望的机翼端点位移为10mm。这种变形被认为对升力系数（CL）产生一些影响：

图18：附加网格配置的翼尖位移演变。

图19：附加网格配置的瞬时CL演变

通过额外测试案例的结果，显示了在结构领域保持精细网格如何显著改善了在各个测试网格尺寸下流体领域机翼端点位移的解决方案。

“FineStrucCoarFlow”案例由于领域之间的网格尺寸差异较大，提供的结果不太准确。会在领域界面处进行解算器之间的结果转换时产生一些小误差。其他两个案例提供的结果与“FineFlex”案例的结果更接近，显示了流体求解器的网格收敛性。

为了验证柔性模型获得的变形，将与相同的追踪器引入到数值机翼中。这使我们能够将实验和数值得到的变形与不同的网格尺寸进行比较。机翼后缘的追踪器的实验结果，可以清楚地看到，随着网格细化，每个步骤的结果都更接近实验测试。

图18

将所有构型的后缘跟踪器与实验变形测量结果进行了比较。

数值跟踪位置和实验跟踪位置之间的最大差异，每个坐标表示为：

图19

每个标准差的差异的平均值也包括在表中，这再次显示了“FineFlex”配置的全局数值结果如何更接近于从测试中获得的实验数据。

为了提供结果的全局视图，我们比较了“FineFlex”跟踪器的位置与实验器的位置，显示了结果很好的一致性。

图20：所有研究构型的数值跟踪和实验跟踪位置的差异。

图21：“FineFlex”跟踪器的位置与实验测量值的比较

通过动量方程法得到的实验测得的柔性振翼的力将用于验证数值模型提供的结果。对于柔性翅膀，只有一个位于翼尖的体积进行了实验测量。其横向尺寸范围从Z=180毫米到Z=260毫米，这意味着实验测得的力只提供了翅膀总力的一部分。

为了比较数值和实验结果，从数值模型解中提取出在测试体积中仅翅尖部分产生的力。展示翅尖处数值力系数演变和实验获得的等效力系数，以及数值所得到的翅膀总力系数。

数值和实验力之间存在小幅差异，但两种情况都提供了相同数量级的结果，并且趋势被很好地捕捉到，从而验证了数值模型。

图22

实验结果（翼尖Exp)、数值结果（尖端Num)和整个机翼的数值结果（总Num)的升降力（底部）和牵引力随翼尖下降行程的变化展示了实验和数值解的涡旋结构对比。

与刚性翅膀类似，两种情况都呈现出典型的振动运动流动结构，具有两个连续的前缘涡（LEV），一个翅尖涡（TP）和一串尾部的后缘涡（TEV）。

通过这一切，我们可以得出结论，柔性翅膀的数值模型经过验证，可以用于更详细地研究柔性对振翼的影响。

图23：比较一个柔性扑翼的实验（左）和数值（右）获得的流动结构的说明。

结语

通过控制方程固体应力求解器、流固耦合、动态网格演化和网格收敛性的研究，我们提出了一系列改进和优化方法，为固体力学数值模拟的发展做出了贡献。

通过探索控制方程固体应力求解器的数值算法和计算效率，我们可以更准确地预测结构物内部的应力分布，并为工程实践提供可靠的设计依据。

对流固耦合问题的了解，使得我们可以更好地理解固体与流体之间的相互作用机制，为解决涉及流固耦合的工程问题提供了思路和方法。动态网格演化方法的应用使得我们能够，进一步提高了数值模拟的准确性。

通过对网格收敛性的研究，我们可以评估数值模拟结果的可靠性和精度，并为模型优化和参数选择提供指导。

相信随着这些研究的不断深入和完善，固体力学数值模拟将在各个工程领域得到更广泛的应用，并推动相关领域的科学研究和工程实践取得新的突破。