毕业论文范文——基于SPH方法的三维矩形容器液体晃动与冲击数值模拟

西安航空职业学院毕业论文基于SPH方法的三维矩形容器液体晃动与冲击数值模拟姓 名： 专 业： 航空电子 班 级： 完成日期： 指导教师: 摘要：基于SPH方法，采用Pam-crash软件对无阻尼板和装有阻尼板的矩形容器在加速度激励作用下，容器内液体的晃动与冲击进行了三维数值模拟，并将测试点的计算压力及液体晃荡运动模拟与试验进行了对比，吻合较好。与CFD方法相比，SPH方法计算所得结果与试验更为接近，从而验证了SPH方法在求解三维液体大幅晃动问题方面的准确性。分析表明：阻尼板的安装对液体的晃动特性具有显著的影响，迟滞特性所表现的压力双峰特征在大充液比情况下消失。关键词：SPH；自由液面；晃动；冲击；三维数值模拟1 引言带自由液面的液体大幅晃动问题广泛存在于航空航天、石油化工和交通运输等领域，一直受到各国工程研究人员的关注。液体大幅晃动具有高度的非线性，描述晃动的一组方程呈强耦合性，自由液面边界条件未知。对于如此复杂的一个问题，传统的数值计算方法如MAC法、VOF法以及常用的有限元法均存在着一些缺陷：MAC法对硬件要求高；VOF法需要高精度的偏微分方程离散格式；有限元方法在处理液体大变形问题时会出现网格的扭曲缠绕从而导致计算的中止。光滑粒子流体动力学（SPH）方法是由Lucy和Monaghan1提出的。该方法的计算对象是空间运动的粒子，采用拉格朗日方法描述。无网格粒子形式以及拉格朗日性质的SPH方法避免了大变形时网格发生畸变的问题，非常适用于液体大幅晃动问题的研究。自Monaghan2首次将SPH方法应用于自由表面水波的流动模拟以来，一些国内外学者开展了对SPH方法的应用研究。Iglesias等3-4采用SPH方法模拟二维矩形水槽的纵摇，并与试验进行了对比；Yonghwan Kim5利用此方法进行了载液货船液舱液体大幅晃动及与货船运动相耦合的数值模拟。陈正云等6采用SPH方法对二维液舱内液体大幅晃荡问题进行了数值模拟，将计算得到的波面图与试验结果进行了对比，吻合较好；马利7和梁龙河8运用光滑粒子流体动力学耦合有限元的方法进行了碰撞和穿透问题的模拟。目前，相对于二维数值模拟，基于SPH方法的液体大幅晃动三维数值模拟研究比较少。本文采用SPH方法进行三维矩形容器液体晃动与冲击的数值模拟，进一步探讨和验证了采用SPH方法进行带自由液面的液体大幅晃动数值计算的准确性。2 SPH数值方法2.1 SPH方法SPH方法通过构造一个近似场函数来求解一批具有一定质量的粒子在任意时间的速度和能量，表达式如下： (1)式中，和为坐标向量，为光滑核函数，为决定光滑核函数影响域尺寸的光滑长度。一般取Monaghan提出的分段三次B样条函数作为核函数，如式（2）所示： (2)式中，是粒子和之间的距离与光滑长度的比值，为与维数相关的系数，在一维、二维、三维空间中，分别为，。通过对式（1）离散化，可得粒子处的粒子近似式如式（3）所示： (3)式中，和(1, 2, , )为粒子的密度和质量，为粒子影响域里的粒子总数，利用核函数的可微性，可得函数空间导数的粒子近似式： (4)从式（3）和（4）可以看出，粒子近似式将函数及导数的连续积分表达式转换成任意排列的粒子的离散化求和，从而实现了空间的无网格化。2.2 控制方程连续介质流体动力学N-S方程分别由式（5）和式（6）描述： (5) (6)式中，为液体密度，是速度向量，为压强，为液体的动粘系数，为体力。通过对上述方程的转换，可以得到方程的SPH表达式如下9： (7) (8)式中，为人工粘度，用于获得稳定的计算，具体表达式如下：(9)式中，。和为常数，一般取1左右，本文的取值均为1，为声速，为粒子的速度矢量。对于不可压缩流，通过利用人工压缩性10可以将其视为可压缩流，粒子的压强由流体状态方程求得： (10)式中，为液体初始密度，常数，为表征粒子声速的参数，为容器内液面高度。运用SPH方法进行不可压缩流计算时，时间步长与声速成反比。为了使显式时间积分具有稳定性，采用CFL条件进行约束11，即要求时间步长与粒子间距离的最小值成正比，在SPH中即转化为时间步长与最小光滑长度成比例： (11)2.3 固壁边界条件为了阻止粒子穿透固壁边界，在SPH方法中，通过在固壁边界上分布一组虚粒子用于对邻近边界的实粒子作用排斥力，从而阻止邻近边界的这些粒子非物理穿透边界，虚粒子和相邻粒子之间的作用力沿着两粒子的中心连线方向。力的大小如式（12）所示： (12)式中，、和为可调参数。文中的取值为，和一般取12和4。为当前时刻虚粒子和邻近边界实粒子的间距，为截止半径：若太大，会造成初始阶段，实粒子就受到排斥力的作用，扰乱粒子的分布；若太小，粒子可能穿过边界也感受不到排斥力作用，其取值一般与粒子初始间距的大小相近，算例中取其值为粒子的初始间距0.0076m。3 数值算例采用SPH方法，不考虑容器内液体表面所受张力，分别对Craig试验12中无阻尼板和装有阻尼板的矩形容器在平动加速度激励作用下，容器内液体的晃动与冲击进行数值计算，并与试验结果及CFD计算结果12进行了对比。试验中在矩形容器的底面和侧壁布置两个压力测试点P1和P2，容器的长、宽、高分别为500mm、400mm、400mm，容器的主视图及压力测试点布置如图1所示。充液比为0.5，容器固定于一辆轻型卡车上，车由静止加速至40km/h再减速至静止状态，小车的实测加速度曲线及滤波后用于数值计算的曲线如图2所示，表1给出了主要的材料属性参数。表1 材料属性属性名称数值矩形容器及阻尼板密度 / (kg/m3)7830矩形容器及阻尼板杨氏模量 / (Pa)2.07e11矩形容器及阻尼板泊松比0.3液体密度 / (kg/m3)998液体杨氏模量 / (Pa)2.2e9液体粘性系数 / (N.s/m2)0.001 图1 矩形容器示意图 图2 滤波前后的加速度曲线图3比较直观地显示了在加速度激励作用下自由液面晃动的三维模拟及其对应的对称面上的粒子速度矢量图，从图中可以清楚看出液面的翻卷及波浪的破碎，当侧壁处粒子速度方向垂直指向壁面，且运动速度比较大时，会对容器壁产生很大的冲击压力。图4给出了容器运动过程中三种典型的液体晃荡现象驻波、行进波和水跃。图5为有无阻尼板情况下，两种数值计算方法对晃动的模拟与试验的对比，SPH计算较真实地模拟出了试验中自由液面随着容器加速度变化，从平静到大幅晃动的过程，该算例中，相对于CFD方法，SPH方法模拟的结果更接近于试验。图6给出了试验、CFD方法及SPH方法所得的测试点P1处的压力时间历程曲线。从图6（a）中可以看到5s之后测试点的压力所呈现的双峰特征，这是液体晃动过程中迟滞特性的一种典型表现。其中前一个峰值是由于液体对容器壁的冲击而产生的高脉冲压力，时间短，具有比较明显的局部特性。 这可以从图2中的加速度曲线分析出：5-10s这段时间内加速度变化大，容器壁主要承受水跃（图4（a）或者行进波（图4（b）晃动现象引起的瞬时冲击压力。后一个峰值也是由于冲击产生的脉冲压力，但是时间相对较长，而且冲击压力峰值变小，这主要是由于快速并且连续的液体压力作用在没有被完全淹没的容器壁（图4（c），压力特性主要表现为连续性，而不像前一个峰值中所明显表现出的瞬时性。当给容器加装阻尼板后，液体对容器壁的冲击力度降低，由于阻尼板的面积较大，阻尼板的安装一定程度上可以等效为增加了容器的充液比，在大充液比情况下，迟滞现象消失，因此图6（b）中没有明显的双峰特性，且压力峰值减小，晃动现象更多的表现为驻波运动（图4（d）。另外，在无阻尼板工况下，当加速度减至0并基本趋于稳定之后，P1处的压力值开始在初始静态压力值处附近振荡，并呈现出有阻尼的逐渐衰减的趋势；在有阻尼板工况下，P1处压力值则很快收敛于初始静态压力值。通过将图6中的SPH计算结果与试验及CFD计算结果对比可以发现，SPH方法计算所得到的曲线较CFD方法不论在变化趋势、压力作用周期还是压力峰值方面都与试验吻合更好。计算峰值稍低于试验所测峰值，这是由于前面所提及的对小车的加速度激励曲线进行了滤波，滤去了一部分波峰的缘故。对于8-10s时间段，两张图中的计算曲线与试验吻合均不是很理想，从数值计算所得曲线的斜率可以看出压力变化比较突然，这可能是由于试验中低频加速度激励测量误差所致12，该时间段加速度测量值变化程度要比真实情况剧烈，导致数值计算中的液体冲击幅度较大，而在真实试验中，此时间段内容器壁承受的压力主要表现为持续性的压力，因而压力曲线变化比较平缓。 (a) 无阻尼板 (b) 有阻尼板图3 自由液面晃动的三维模拟及对应的对称面粒子速度矢量图 (a) 水跃(无阻尼板) (b) 行进波(无阻尼板) (c) 驻波(无阻尼板) (d) 驻波(有阻尼板)图4 典型晃荡运动 试验 CFD计算 SPH计算 试验 CFD计算 SPH计算 (a) 无阻尼板 (b) 有阻尼板图5 液体晃动的数值模拟与试验对比 (a) 无阻尼板 (b) 有阻尼板图6 SPH数值计算与试验及CFD数值计算压力曲线对比 4 结论采用SPH方法，对矩形容器在平动加速度激励作用下自由液面的晃动进行了三维数值模拟，通过与文献中试验数据的对比，发现计算所得到的压力曲线在变化趋势、压力作用周期和峰值方面都与试验吻合较好，此外，文中算例采用的SPH方法计算所得结果比CFD方法更接近于试验；整个晃动过程模拟也与试验相近，并且比较真实地模拟出了自由液面大幅晃动下采用有限元方法难以模拟的碎波现象。因此，SPH方法进行带自由液面的液体大幅晃动计算是合理的、准确的。本文为求解类似关于液体晃动的问题如飞机油箱燃油晃动、航天器推进剂晃动提供了有效的方法和手段。在本文的计算中，笔者为提高计算精度，所用粒子数达到18万，尽管采用了4个CPU进行并行计算，但是发现并行效率并不高，因此仍花费了大量的计算时间，如何提高并行效率将是下一步所要探讨的问题。另外，文中进行了容器增加阻尼板后的液体晃动计算，阻尼板的数量、形状、弹性和安装位置等因素对容器内液体的晃动及容器壁所受的冲击压力具有比较大的影响，关于容器防晃阻尼板优化这方面的研究，今后将采用SPH方法结合优化算法开展进一步的工作。参考文献：1 R. A. Gingold, J. J. Monaghan. Smoothed Particle Hydrodynamics-Theory and application to non-spherical stars J. Royal Astronomical Society, Monthly Notices, 1977, 181: 375-389.2 J. J. Monaghan. Simulating free surface flow with SPH J. Journal of Computational Physics, 1994, 110: 399-406.3 A. Souto-Iglesias, L. Perez-Rojas, R. Zamora Rodriguez. Simulation of anti-roll tanks and sloshing type problems with smoothed particle hydrodynamics J. Ocean Engineering, 2004, 31: 1169-1192.4 A. Souto-Iglesias, L. Delorme, L. Perez-Rojas et al. Liquid moment amplitude assessment in sloshing type problems with smooth particle hydrodynamics J. Ocean Engineering, 2006, 33: 1462-1484.5 Kim. Yonghwan. Experimental and numerical analyses of sloshing flows J. Journal of Engineering Mathematics, 2007, 58: 191-210.6 陈正云, 朱仁庆, 祁江涛. 基于SPH 法的二维液体大幅晃荡数值模拟 J. 船海工程, 2008, 37(2): 44-47.7 马利, 陶伟明, 郭乙木等. SPH耦合有限元方法的水射流弹塑性碰撞模拟 J. 浙江大学学报(工学版), 2008, 42(2): 259-263.8 梁龙河, 董师舜, 王政. 刚体与SPH耦合的三维数值方法研究 J. 高压物理学报, 2008, 22(3): 275-280.9 L. D. Libersky, A. G. Petschek, T. C. Carney et al. High strain Lagrangian hydrodynamics: A three-dimensional SPH code for dynamic material Response J. Journal Computational Physics, 1993, 109: 67-75.10 J. J. Monaghan. Simulating free surface flows with SPH J. Journal of Computational Physics, 1994, 110(2): 399-406.11 P. W. Randles, L. D. Libersky. Smoothed Particle Hydrodynamics: Some recent improvements and application J. Computer Methods in