鞋涉水问题的多体动力学与流-固耦合系统建模及仿真分析

2020-12-30吴晓建祁祺王斌周兵

湖南大学学报（自然科学版） 2020年12期

吴晓建，祁祺†，王斌，周兵

（1.南昌大学机电工程学院，江西，南昌，330031；2.湖南大学汽车车身先进设计制造国家重点实验室，湖南长沙 410082）

CAE（Computer Aided Engineering，CAE）技术在工程分析中发挥着重要的作用，它为工程问题的发现、解决提供了关键性的验证和指导.然而，复杂系统往往同时涉及多体动力学与流-固耦合等不同技术方向，有效而准确的联合仿真已成为工程界共性技术难点.马增辉等[1]采用ALE（Arbitrary Lagrange-Euler，ALE）法进行了水陆两用飞机降落耐波性的研究；江旭东等[2]采用耦合欧拉-拉格朗日CEL（Coupled Euler-Lagrangian，CEL）方法研究了压差驱动式管道机器人的动力特性；雷娟棉等[3]采用CFD（Computational Fluid Dynamics，CFD）的6 自由度动网格法分析了导弹挂载分离过程；杨晓彬等[4]采用CFD 的重叠网格法分析了飞机水上降落的水动力砰击载荷.以上典型复杂工程应用问题侧重于流-固耦合系统建模与仿真，未涉及与多体动力学的联合.本文所研究的行人在雨天行走时的鞋涉水问题则是一个多体多力学与流-固耦合仿真的问题，不同的步态和鞋在涉水过程具有差异化的“涉水-出水”过程，准确地将此过程予以模拟呈现对工程中复杂耦合问题的仿真分析具有良好的参考价值.

鞋的涉水过程指在积水路面行走时，鞋在人体下肢运动的驱动下和水进行接触的过程.整个过程同时涉及结构体、空气和水三相的相互作用，包括流动介质突变引起的载荷突变、自由液面变化、姿态变化等众多的复杂问题.

同时，鞋的运动往往与人体的下肢运动相关，其姿态以及运动形式往往是复杂多变的，进而导致涉水运动成为复杂的多体动力学-流-固耦合问题.

本文将行人行走的鞋涉水问题分析分为两个基本方面，首先建立人体下肢模型，采用多体动力学分析，得到鞋体运动的边界条件；然后建立鞋涉水多体运动的流-固耦合模型.其中，耦合欧拉-拉格朗日CEL 方法，采用基于体积分数的流-固耦合边界追踪算法，自动在结构和流体域间进行载荷、位移、速度等信息的传递，能够解决空化射流冲击[5]、存在大变形的切削问题[6]以及震荡过程的流固耦合问题[7].因而本文将采用CEL 方法对鞋“涉水-出水”过程的流-固耦合问题进行模拟.

最后，基于以上研究方法，对比分析运动鞋和平底鞋的“涉水-出水”过程，通过涉水过程时间和所需能量等，评价两种鞋涉水性能的强弱.本文以生活中常见的行人在湿地行走，其鞋涉水现象为研究对象，开展复杂系统多体动力学-流固耦合仿真技术研究，尽管不是直接的工程问题，但对类似工程问题的建模及分析具有良好的借鉴作用.

1 鞋涉水过程的分析方法

1.1 基于Adams 的人体步态运动

采用Adams 多体动力学仿真软件模拟人体行走步态，进行足部动力学和运动学分析.通过把上身视为一个质量块，下肢部分分成大腿、小腿和足部，以及在各关节之间设立运动副和驱动，建立人体下肢简化模型.这种简化方法可以解决步态仿真模型不易稳定的问题，且较为准确地反映出人体步态运动.

1.2 耦合欧拉-拉格朗日法

应用于传统有限元分析的拉格朗日方法，以物体的坐标为基础、以有限元网格的节点为物体的质点，可以跟踪质点的运动，从而描述物体的运动，但遇到大变形问题时，会由于网格变形扭曲而失去原有的精度，甚至导致计算无法继续.欧拉方法常用于流体动力学计算，以网格为空间，以网格节点为空间点，材料不随单元变形而在网格中流动，因此可以解决包含大变形和流体流动的问题，但在本文鞋涉水过程流固耦合（Fluid-Structure Interaction，FSI）分析中，欧拉分析虽然可有效进行流体流动分析，却在捕捉结构流固交界面上存在一定困难；并且欧拉方法固存的数值耗散问题导致其对计算资源和时间要求较高.由此可知，单纯的拉格朗日算法和单纯的欧拉算法都有各自的缺陷和不足，也有着各自的优势；如果将两者有机地结合起来，可以解决一些只用单一方法所不能解决的问题.耦合欧拉-拉格朗日CEL方法就是基于这一目的最早由Noh[8]提出的，它采用有限差分法求解带有移动边界的二维流体动力学问题.在Noh 的研究中，网格点可以随物质点一起运动，但也可以在空间中固定不动，甚至网格点可以在一个方向上固定，而在另一个方向上随物质一起运动.因此，CEL 算法在解决物体的大位移时，比如碰撞、流体动力学及流体-固体之间相互作用时有强大的优势.

2 基于Adams 的人体步态运动建模

2.1 人体下肢运动的几何模型

2.1.1 动力学建模数据

因本研究只关注人体下肢的运动，故省略人体上肢建模.下肢结构的主要参数如表1 所示，其中，脚直接使用鞋的几何模型.

表1 下肢主要参数Tab.1 Main parameters of human lower limb

2.1.2 模型平衡控制

人体模型是一种非线性且不易稳定的系统，如何调整各关节之角度或速度来保持平衡一直是非常复杂的难点问题.本文提出一种简化方法解决模型运动过程中的平衡问题，即假设人体运动过程中质心位置变化不大，在此条件下胯部的运动和地面成平行状态，因此在胯部设一平衡杆来控制胯部运动实现平衡，如图1 所示.

图1 人体行走几何模型Fig.1 Geometric model of human walking

2.2 人体下肢运动的动力学模型

2.2.1 模型的运动约束

下肢建模是人体能否实现步行及步行质量好坏的关键所在，尤以关节结构为集中体现.人体膝关节由胫骨-股骨和髌骨-股骨关节组成，这样膝关节就成为胫骨和股骨相咬合的关节[9-10]；同时，膝关节和髋关节都只有一个自由度，可以用转动副连接.将地面设置为固定副，胯部和平衡杆之间为移动副.在各转动副设置驱动器，设置地面和脚之间的接触形式为碰撞.创建出符合实际的步态运动模型，如图2 所示.

图2 步态运动仿真模型Fig.2 Gait motion simulation model

2.2.2 模型的运动控制

人的行走是一个周期活动过程，一个步态周期为从脚跟着地到脚跟再次着地的时间.如图3 所示，步态周期分为支撑相和摆动相.其中，支撑相是指脚跟着地到脚离地，每只脚的支撑相又分为首次触地、支撑初期、支撑中期、支撑末期4 个部分，占整个步态周期的62%左右；摆动相是脚尖离地到脚跟着地，分为摆动早期、摆动中期、摆动末期，大约占整个步态周期38%[11].本文使用step 函数对关节角度进行控制，以模拟行人的步态运动.

图3 人体步态周期Fig.3 Human gait cycle

2.2.3 模型的运动仿真

设置仿真总时间为2 s，仿真步数为500，即每个仿真步长为0.004 s，仿真结果如图4 所示.

图4 步态仿真结果Fig.4 Gait simulation results

对比图3 和图4 可以发现，所建立的人体步态仿真模型，其模拟结果与步态周期吻合度高，能够较为准确地反映出行走过程中鞋的运动与姿态.

为进一步开展鞋涉水流-固耦合仿真，提取出鞋参考点的速度和角速度作为有限元仿真的边界条件，参考点选取了髋关节和脚跟连线中点的位置，具体如图5 所示，参考点的速度和角速度值如图6 所示.

图5 参考点的位置Fig.5 Position of reference point

图6 参考点的速度Fig.6 Reference point speed

3 鞋涉水过程流固耦合仿真

3.1 鞋涉水的系统描述

鞋涉水有限元模型主要包括三部分，即欧拉区域、鞋体和刚性地面.其中欧拉区域为欧拉体，鞋体和刚性地面为拉格朗日体.考虑到鞋体-地面-流体之间的复杂的接触关系，通过罚函数法的通用接触算法来描述多体间的运动.

先以运动鞋为例分析其涉水过程.模型如图7所示，运动鞋单只质量为0.29 kg，鞋长250 mm，鞋底最厚处为30 mm、最薄处8 mm.为减少网格数量和计算消耗，暂不考虑鞋底花纹影响，且对鞋体模型进行适当简化.

图7 运动鞋建模Fig.7 Sports shoes model

将鞋体材料简化为橡胶来模拟鞋体的弹性行为.通过Mooney-Rivlin 模型描述橡胶材料鞋体模型的超弹性本构关系，则有：

式中：W 表示应变能密度；C10和C01为材料常数；I1和I2为第一、第二Green 不变量.

Ogden[12]已经证明，在应变能密度和应力、应变之间存在如下关系：

式中：t 为真实应力；λ1为主伸长比.对于单轴拉伸实验，t2=t3=0，三个方向主伸长λ1=λu，λ2=λ3=1/，则应力应变关系为：

对于双轴拉伸实验，t1=t2，t3=0 三个方向主伸长λ3=，λ1=λ2=λB，则应力应变关系为：

对于平面拉伸实验，t3=0，三个方向主伸长λ1=λS，λ2=1，λ3=1/λS则应力应变关系为：

对橡胶材料试件进行单轴拉伸、双轴拉伸和平面拉伸实验，得到橡胶材料的实验数据.由此，根据式（5）、（6）、（7）拟合橡胶材料的拉伸实验数据，结果如图8 所示（各条曲线表示材料拉伸数据和各种本构方案下的拟合曲线）.最终获得Mooney-Rivlin 模型参数C10=0.181 MPa，C01=0.003 84 MPa.

图8 拟合的橡胶应变势能曲线Fig.8 Fitted strain-stress curve

采用CEL 算法时，为防止欧拉区域和刚体地面接触失效，同时保证水花溅起高度的捕捉，需要设置欧拉区域远大于拉格朗日区域.

在Abaqus 中对水的材料属性进行特殊定义，具体通过材料的状态方程实现，在此选用Mie-Gruneisen 方程[13]，如公式（8）所示：

式中：p 为液体压强，Γ0、s 为材料常数；η=1 -ρ0/ρ为名义体积压缩应变；ρ0为材料初始密度；Em为单位质量内能；C0为材料初始声速.

3.2 有限元模型的建立

欧拉区域的主要尺寸为1 000 mm × 500 mm ×150 mm，为了保证欧拉区域和刚性地面的成功接触，设置接触区域尺寸为1 000 mm×500 mm×5 mm.划分鞋体单元为四结点线性四面体单元C3D4，地面单元为四结点三维双线性刚性四边形单元R3D4，欧拉区域为八结点线性欧拉六面体单元EC3D8R.建立有限元模型如图9 所示.

图9 流固耦合过程的有限元模型Fig.9 Finite element model

为防止欧拉材料流出边界，设置欧拉区域边界的所有法向速度为0.地面为解析刚体，约束其所有自由度.将鞋体表面节点耦合到鞋的参考点上，将Adams 得到的运动数据添加到鞋参考点上，以设置鞋体的边界条件，且设置欧拉区域、鞋体和刚性地面的接触为通用接触.

在CEL 算法中，欧拉区域及自由表面通过体积分数描述[14-15].体积分数1 代表欧拉材料充满欧拉单元；0 代表欧拉单元中不存在欧拉材料.初始状态下欧拉区域体积分数都为0，这表示初始状态下欧拉区域没有欧拉材料分布.为了描述水的初始分布，需设定欧拉材料在欧拉区域的初始分布，本文初始分布的体积模型来自于布尔操作，通过体积分数工具指定其初始材料分布.根据鞋体的初始位置，得出欧拉模型和初始水域模型如图10 所示.

图10 欧拉模型和水域初始节点分布Fig.10 Euler model and initial node distribution in waters finite element model

3.3 计算结果分析

经过有限元分析计算，图11 给出了CEL 计算所得鞋涉水过程.选取了几个比较重要的时刻，可以看出，随着鞋的运动，CEL 算法可以清楚地计算水的运动过程.

图11 涉水计算结果Fig.11 Wading results

进一步地，将鞋“涉水-出水”过程分为5 个阶段，如图11（a）～（e），左侧为水的体积分数云图，右侧为鞋体的应力云图.第一阶段鞋体运动使其前端与地面接触，在这一过程中，发生了冲击，使鞋体周围的水有向外排出的趋势，为排水阶段，如图11（a）；第二阶段鞋体开始向前运动，推鞋前端的水向前运动，为推水阶段，如图11（b）；第三阶段，鞋继续向前运动，鞋前的水继续堆积，在鞋前端产生一个水膜，为水膜生成阶段，如图11（c）；第四阶段，鞋继续向前运动，并伴随向上抬起动作，鞋前的水膜开始破裂，为水膜破裂阶段，如图11（d）；第五阶段，鞋前的水膜完全破裂，鞋从水中出来，为完全出水阶段，如图11（e）.

读取水区域动能变化曲线如图12 所示.可以看到出现了能量峰值，峰值和陡升起点分别在0.125 s和0.175 s.

图12 水区域动能变化图Fig.12 Water kinetic energy change chart

这两个时刻的水的体积分数云图如图13（a）所示；水域速度云图如图14（a）所示.可以发现，陡升起点出现在排水阶段，峰值出现在水膜发展阶段，这是由于，排水和推水过程使水体发生大范围运动，使系统动能陡升.在后续阶段中，因鞋体出水，水和鞋的相对冲击减小，导致水域动能减小，但仍出了一个峰值，其起始时间和峰值时间分别在0.335 s 和0.365 s.后两个时刻的体积分数云图如图13（b）所示；水域速度云图如图14（b）所示.可以发现，为了补充因鞋的划水运动而产生的水坑，水开始向水坑运动；并且这段时间内飞溅的水下落到水面上，使系统动能增大.

图13 峰值时刻的体积分数云图Fig.13 Volume fraction at peak time

图14 峰值时刻的水速度云图Fig.14 Water velocity at peak time

4 不同鞋的涉水性能比较

选取鞋长和运动鞋相同的平底鞋作为对比，两者比较明显的特征区别在于平底鞋鞋头前端没有上翘.为确保可比性，设置运动鞋底厚度的最厚和最薄尺寸与运动鞋一致，且平底鞋的前端鞋面高度与运动鞋的大致相当.根据以上原则选择图15 所示平底鞋模型并加以适当简化.

图15 平底鞋建模Fig.15 Flat shoes model

除鞋体模型外，其他模型和行人步态等边界条件均保持不变，通过计算，比较相同时刻水的体积分数云图和速度云图等来分析两种鞋涉水性能.

图16 和图17 分别为两种模型在主要时刻下的体积分数云图和水域速度云图.时刻1 中运动鞋模型处于水膜破裂阶段，平底鞋处于水膜生成阶段，可以发现平底鞋生成的水膜大于运动鞋生成的水膜.时刻2 中运动鞋处于水膜破裂末期，平底鞋处于水膜破裂初期，平底鞋水膜包裹鞋体的面积大于运动鞋.时刻3 中运动鞋处于完全出水阶段，平底鞋处于水膜破裂阶段，平底鞋溅起的水花大于运动鞋溅起的水花.

图18 为运动鞋和平底鞋计算模型的系统动能比较，可以发现平底鞋的水域系统动能大于运动鞋的水域动能，因为步态数据的一致，动能峰值位置基本相同.

图16 体积分数云图结果比较Fig.16 Comparison of volume fraction results

图17 速度云图结果比较Fig.17 Comparison of velocity results

图18 不同鞋水区域系统动能比较Fig.18 Different water kinetic energy comparison

通过对不同时刻下的体积分数云图和水域速度云图分析，可以发现，在相同运动步态下，平底鞋溅起的水花大于运动鞋，因而从涉水到完全出水过程所花费的时间大于运动鞋；通过分析不同模型水系统的动能变化，可以发现平底鞋系统水区域动能大于运动鞋，因此，平底鞋的涉水-出水过程要花费更多的能量.以涉水时间和消耗能量大小为涉水性能评价指标，平底鞋的涉水性能弱于运动鞋.