一种基于无网格流固耦合的降落伞分析方法与流程

本发明涉及一种基于无网格流固耦合的降落伞分析方法，属于航空救生
技术领域：
。
背景技术：
：降落伞作为一种充气织物由于质量轻，易折叠等特点被广泛应用于航空救生、航天减速等领域。降落伞性能是否满足设计要求将直接关系到人员或整个系统的安全。目前，降落伞性能分析主要有试验方法和数值方法。试验方法需要花费大量人力、物力和财力，此外，由于试验数据采集较为困难，试验方法很难揭示其力学机理。而数值方法由于经济性、可重复性等特点已经成为降落伞性能分析的重要手段。目前，工程中使用最多的流固耦合分析方法主要是ALE(ArbitraryEulerian-Lagrangian)方法(可以参考：程涵，余莉，李胜全.基于ALE的降落伞充气过程数值仿真.南京航空航天大学学报,2012,44(3):290-293.或HanCheng,XinghuaZhang,LiYU.StudyofVelocityEffectsonParachuteInflationPerformanceBasedonFluid-StructureInteractionMethod.AppliedMathematicsandMechanics(Englishedition),2014,35(9):1177-1188.)。但ALE分析方法需要在每一个时步，对流场网格进行重构，并进行流场信息的更新，这大大增加了计算分析时间，不利于工程设计的广泛使用。而SPH(SmoothedParticleHydrodynamics)方法作为一种无网格方法具有传统网格方法所没有的优点，由于SPH近似方法的自适应性质，SPH公式构造不受粒子分布的随意性所影响，因此可以很自然地处理一些具有极大变形问题，如流固耦合问题(可以参考：林晓东,卢义玉,汤积仁等.基于SPH-FEM耦合算法的磨料水射流破岩数值模拟.振动与冲击.2014,33(18):170-176.或卢鸣飞,华诚.基于FEM-SPH耦合算法的飞机水上迫降的力学建模与分析.复旦学报(自然科学版).2013,52(5):583-591.)。但SPH方法也存在边界处理困难的问题，很难模拟风洞中恒速气流。技术实现要素：本发明所要解决的技术问题是提供一种基于无网格流固耦合的降落伞分析方法，采用相对压力描述气体及局部计算域的技术，将SPH方法应用于降落伞数值分析，同时解决现有SPH计算无法模拟风洞的难点。本发明为解决上述技术问题采用以下技术方案：本发明提供一种基于无网格流固耦合的降落伞分析方法，包括以下具体步骤：步骤1，根据降落伞实物建立等比折叠状态有限元模型，其中，伞衣采用三角形单元划分，加强带及伞绳采用杆单元划分；步骤2，建立计算域和非计算域，定义计算域与非计算域相连的面为流场域进口，计算域中与流场域进口相对的面为流场域出口；步骤3，在计算域与非计算域内建立流场光滑粒子流体力学SPH粒子模型，其中，SPH粒子模型和降落伞网格相互穿插，且伞衣位于计算域内；步骤4，定义基于相对压力描述的气体状态方程，间接解决流场压力出口边界问题：其中，p为压力，γ为气体比热，ρ为密度，e为单位质量内能，ρ0为初始密度，pe为外界大气压力；步骤5，对非计算域的所有粒子施加速度条件，其值为风洞吹风速度，在流场域入口处获得恒速的粒子流；步骤6，定义计算域中除流场域入口、出口外的4个面为对称边界；步骤7，对计算域中的流场采用基于蛙跳格式的SPH方法计算，伞衣网格采用基于中心时间差分格式的有限元算法进行计算；流体与伞衣结构的耦合计算转化为基于罚函数的接触计算，计算域之外的流体粒子则不进行核函数近似计算；步骤8，根据步骤7计算出的结果对降落伞气动性能进行分析。作为本发明的进一步优化方案，步骤6中对称边界为基于虚粒子的对称边界。作为本发明的进一步优化方案，步骤1中有限元模型在前处理软件LS-PREPOST中建立。作为本发明的进一步优化方案，步骤2中在前处理软件LS-PREPOST中建立计算域和非计算域。作为本发明的进一步优化方案，步骤3中采用前处理软件LS-PREPOST在计算域与非计算域内建立流场光滑粒子流体力学SPH粒子模型。本发明采用以上技术方案与现有技术相比，具有以下技术效果：本发明通过设立局部计算域及改变气体状态方程实现了流场SPH模型与伞衣有限元模型相耦合，大大降低计算时间，且结果可靠。附图说明图1是伞衣有限元模型，其中，(a)是伞衣模型主视图，(b)是伞衣模型右视图。图2是流场域。其中，1-计算域，2-非计算域，3-流场域出口，4-流场域进口，5-降落伞网格，6-SPH粒子。图3是模型装配图。图4是耦合计算流程图。图5是数值结果，其中，(a)是0.01s时伞衣结构变化示意图，(b)是0.01s时流场变化示意图，(c)是0.055s时伞衣结构变化示意图，(d)是0.055s时流场变化示意图，(e)是0.12s时伞衣结构变化示意图，(f)是0.12s时流场变化示意图，(g)是0.25s时伞衣结构变化示意图，(h)是0.25s时流场变化示意图。图6是实验结果和计算结果动载对比。具体实施方式下面结合附图对本发明的技术方案做进一步的详细说明：以某模型伞的充气展开过程为例进行流固耦合分析：首先，在前处理软件LS-PREPOST中，根据该模型伞实物(即如表1所示的伞衣结构和材料参数等)建立有限元模型，其中伞衣采用4646个三角形单元划分，伞绳和加强带采用616个杆单元划分，最终有限元模型可如图1所示。表1伞衣模型参数由于流场采用基于拉格朗日描述的SPH方法计算，因此必须预先建立所有进入计算域的粒子模型。在LS-PREPOST软件建立2.5m×2.5m×15m的区域，并划分为2块区域，即计算域1和非计算域2。其中，计算域1的尺寸为2.5m×2.5m×5m，而非计算域2的尺寸为2.5m×2.5m×10m。如图2所示，计算域1前后两个面分别为流场域出口3和流场域进口4。计算过程中，非计算域2的流体粒子将通过流场域进口4进入计算域，而计算域1中的流体粒子将通过流场域出口3到外部空间，同样外部空间也是非计算域。伞衣位于计算域1，如图2所示。采用前处理软件LS-PREPOST对计算域1和非计算域2建立粒子模型，共1458000个粒子。由于降落伞网格5与SPH粒子6都是基于拉格朗日描述，因此两者相互穿插，最终模型装配图如图3所示。SPH方法是无法定义压力出口边界条件的，因此，本发明通过定义基于相对压力描述的气体状态方程间接解决了这一难题，由于流场域出口4处的SPH粒子和外界不存在压差，不会导致计算错误或失败。因而，无需在流场域出口4进行边界定义。气体状态方程如下所示：其中，p为压力，γ为气体比热，ρ为密度，e为单位质量内能，ρ0为初始密度，pe外界大气压力。同样SPH也是无法定义速度入口边界条件的，因此，对非计算域2的所有粒子施加速度条件，其值为风洞吹风速度。由于非计算域(包括非计算域2和外部空间)的粒子不参与核函数近似计算，因而可以在流场域入口获得恒速的粒子流。而计算域1除去入口和出口外的4个面均采用基于虚粒子的对称边界。对计算域1中的流场和伞衣网格分别采用基于蛙跳格式的SPH算法和中心时间差分格式有限元算法进行计算。在进行SPH核函数近似计算过程中，采用上述公式(1)建立压力和密度之间的关系。流体与伞衣结构的耦合计算则转化为基于罚函数的接触计算。计算域1之外的流体粒子则不进行核函数近似计算，耦合算法可如图4所示，具体算法如下：在SPH算法部分：在每一个时步中，根据更新后的SPH粒子速度和空间位置，判断每个粒子是否在计算域内，若粒子不在，则不进行后续计算，若粒子在计算域内，则进一步按图4中箭头顺序先后计算密度/应变率、压力/应力、粒子力和加速度，而加速度又用于更新下一时步粒子的速度和位置；在有限元算法部分：在每一时步中，根据更新后有限元单元节点位置，按照箭头顺序计算位移、应变/应力，而位移参与第二次单元节点位置更新；该部分的计算采用中心时间差分格式(可参考：BelytschkoT,LiuWK,MoranB.连续体和结构的非线性有限元(庄茁译).北京:清华大学出版社,2002:341-387.)；耦合计算部分：分别根据SPH粒子的速度和空间位置，以及有限单元节点位置，判断伞衣网格流体粒子是否接触，若不接触，则不计算耦合力，若接触则计算耦合力(可参考：程涵.柔性织物折叠建模技术及展开过程数值仿真研究,[博士学位论文].南京:南京航空航天大学,2013.)；计算后得到的耦合力用于参与接触的流体粒子加速度计算，以及参与接触的伞衣网格节点位移的计算。通过上述计算，获得了试验手段较难采集的伞衣结构和流场动态变化，如图5所示。表2为本发明方法与ALE方法耗时比较。同样采用IntelXeonE5-1630处理器及内存32G，在相同粒子(或网格)情况下，本发明计算耗时要少于现有ALE方法，证明了本发明分析的高效性。表2耗时对比本发明方法ALE粒子(或网格)数量14580001458000耗时(h)5362图6为风洞试验和计算动载变化对比，从图中可以发现，风洞试验值与计算值相符，证明了本发明的可行性和准确性。本发明能所采用的全新数值方法能够用于实际工程中降落伞性能分析。以上所述，仅为本发明中的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉该技术的人在本发明所揭露的技术范围内，可理解想到的变换或替换，都应涵盖在本发明的包含范围之内，因此，本发明的保护范围应该以权利要求书的保护范围为准。当前第1页1 2 3