《毕业论文范文——基于SPH方法的三维矩形容器液体晃动与冲击数值模拟》由会员分享，可在线阅读，更多相关《毕业论文范文——基于SPH方法的三维矩形容器液体晃动与冲击数值模拟（10页珍藏版）》请在金锄头文库上搜索。
1、西安航空职业学院西安航空职业学院毕业论文毕业论文基于基于 SPH 方法的三维矩形容器液体晃动与冲击数值模拟方法的三维矩形容器液体晃动与冲击数值模拟姓姓 名：名： 专专 业：业： 航空电子 班班 级：级： 完成日期：完成日期： 指导教师指导教师: 摘要：摘要：基于 SPH 方法，采用 Pam-crash 软件对无阻尼板和装有阻尼板的矩形容器在加速度激励作用下，容器 内液体的晃动与冲击进行了三维数值模拟，并将测试点的计算压力及液体晃荡运动模拟与试验进行了对比， 吻合较好。与 CFD 方法相比，SPH 方法计算所得结果与试验更为接近，从而验证了 SPH 方法在求解三维液体 大幅晃动问题方面的准确性。
2、分析表明：阻尼板的安装对液体的晃动特性具有显著的影响，迟滞特性所表现 的压力双峰特征在大充液比情况下消失。 关键词：关键词：SPH；自由液面；晃动；冲击；三维数值模拟1引言带自由液面的液体大幅晃动问题广泛存在于航空航天、石油化工和交通运输等领域，一 直受到各国工程研究人员的关注。液体大幅晃动具有高度的非线性，描述晃动的一组方程呈 强耦合性，自由液面边界条件未知。对于如此复杂的一个问题，传统的数值计算方法如MAC 法、VOF法以及常用的有限元法均存在着一些缺陷：MAC法对硬件要求高；VOF法需要高精 度的偏微分方程离散格式；有限元方法在处理液体大变形问题时会出现网格的扭曲缠绕从而 导致计算的中止
3、。 光滑粒子流体动力学（SPH）方法是由 Lucy 和 Monaghan1提出的。该方法的计算对象 是空间运动的粒子，采用拉格朗日方法描述。无网格粒子形式以及拉格朗日性质的 SPH 方法 避免了大变形时网格发生畸变的问题，非常适用于液体大幅晃动问题的研究。自 Monaghan2 首次将 SPH 方法应用于自由表面水波的流动模拟以来，一些国内外学者开展了对 SPH 方法 的应用研究。Iglesias 等3-4采用 SPH 方法模拟二维矩形水槽的纵摇，并与试验进行了对比； Yonghwan Kim5利用此方法进行了载液货船液舱液体大幅晃动及与货船运动相耦合的数值模 拟。陈正云等6采用 SPH 方法
4、对二维液舱内液体大幅晃荡问题进行了数值模拟，将计算得到 的波面图与试验结果进行了对比，吻合较好；马利7和梁龙河8运用光滑粒子流体动力学耦 合有限元的方法进行了碰撞和穿透问题的模拟。目前，相对于二维数值模拟，基于 SPH 方法 的液体大幅晃动三维数值模拟研究比较少。本文采用 SPH 方法进行三维矩形容器液体晃动与 冲击的数值模拟，进一步探讨和验证了采用 SPH 方法进行带自由液面的液体大幅晃动数值计 算的准确性。2SPH 数值方法2.1 SPH 方法 SPH 方法通过构造一个近似场函数来求解一批具有一定质量的粒子在任意时间的( )f x 速度和能量，表达式如下：(1),( )()(, )f xf
5、 x W xx h dx式中，和为坐标向量，为光滑核函数，为决定光滑核函数影响域尺寸的光x,x,(, )W xx hh 滑长度。一般取 Monaghan 提出的分段三次 B 样条函数作为核函数，如式（2）所示：(2)233210132 1( , )(2)12602dRRRW R hRRR 式中，是粒子 和之间的距离与光滑长度的比值，为与维数相关的系数，在一维、二Rijd 维、三维空间中，分别为，。d1/h215/7 h33/2 h通过对式（1）离散化，可得粒子 处的粒子近似式如式（3）所示：i(3)1()(, )( )N j jij jjimf xW xx hf x式中，和(1, 2, , )
6、为粒子的密度和质量，为粒子 影响域里的粒子总数，jjmjNjNi利用核函数的可微性，可得函数空间导数的粒子近似式：(4)1()(, )( )N j jj jjiimf xW xx hf x从式（3）和（4）可以看出，粒子近似式将函数及导数的连续积分表达式转换成任意排 列的粒子的离散化求和，从而实现了空间的无网格化。2.2 控制方程 连续介质流体动力学 N-S 方程分别由式（5）和式（6）描述：(5)dudt (6)21pduuFdt 式中，为液体密度，是速度向量，为压强，为液体的动粘系数，为体力。通过upF 对上述方程的转换，可以得到方程的 SPH 表达式如下9：(7)1N i j jij i
7、j idmdtWux (8)22 1()N ji j jijiji iji ipdupmFdtWx 式中，为人工粘度，用于获得稳定的计算，具体表达式如下：ijijuuuij(9)2 000 , , ijijijijijijijijijcvxvx 式中，22ijijij ijijh vxx 1()2ijijccc 1()2ijij 1()2ijijhhh，。和为常数，一般取 1 左右，本文的取值均为 1，为声速，ijijvvvijijxxxc为粒子的速度矢量。v对于不可压缩流，通过利用人工压缩性10可以将其视为可压缩流，粒子 的压强由流iip 体状态方程求得：(10)0()1)i ipB式中，为
8、液体初始密度，常数，为表征粒子声速的参数，07B0100/cBgH 为容器内液面高度。cH 运用 SPH 方法进行不可压缩流计算时，时间步长与声速成反比。为了使显式时间积分具 有稳定性，采用 CFL 条件进行约束11，即要求时间步长与粒子间距离的最小值成正比，在 SPH 中即转化为时间步长与最小光滑长度成比例：(11)min()ishtC 2.3 固壁边界条件 为了阻止粒子穿透固壁边界，在 SPH 方法中，通过在固壁边界上分布一组虚粒子用于对 邻近边界的实粒子作用排斥力，从而阻止邻近边界的这些粒子非物理穿透边界，虚粒子和相 邻粒子之间的作用力沿着两粒子的中心连线方向。力的大小如式（12）所示：
9、(12)120000 20, 1()(), 1ijijnnijijijij ijPB r rxrrrDrrrr 式中，、和为可调参数。文中的取值为，和一般取 12 和 4。为当D1n2nD10cgH1n2nijr前时刻虚粒子和邻近边界实粒子的间距，为截止半径：若太大，会造成初始阶段，实粒0r子就受到排斥力的作用，扰乱粒子的分布；若太小，粒子可能穿过边界也感受不到排斥力作用，其取值一般与粒子初始间距的大小相近，算例中取其值为粒子的初始间距 0.0076m。3数值算例采用 SPH 方法，不考虑容器内液体表面所受张力，分别对 Craig 试验12中无阻尼板和装 有阻尼板的矩形容器在平动加速度激励作用
10、下，容器内液体的晃动与冲击进行数值计算，并 与试验结果及 CFD 计算结果12进行了对比。试验中在矩形容器的底面和侧壁布置两个压力 测试点 P1 和 P2，容器的长、宽、高分别为 500mm、400mm、400mm，容器的主视图及压力 测试点布置如图 1 所示。充液比为 0.5，容器固定于一辆轻型卡车上，车由静止加速至 40km/h 再减速至静止状态，小车的实测加速度曲线及滤波后用于数值计算的曲线如图 2 所示， 表 1 给出了主要的材料属性参数。表 1 1 材料属性属性名称属性名称数值数值矩形容器及阻尼板密度 / (kg/m3)7830矩形容器及阻尼板杨氏模量 / (Pa)2.07e11矩形
11、容器及阻尼板泊松比0.3液体密度 / (kg/m3)998液体杨氏模量 / (Pa)2.2e9液体粘性系数 / (N.s/m2)0.001图 1 矩形容器示意图 图 2 滤波前后的加速度曲线图 3 比较直观地显示了在加速度激励作用下自由液面晃动的三维模拟及其对应的对称面 上的粒子速度矢量图，从图中可以清楚看出液面的翻卷及波浪的破碎，当侧壁处粒子速度方 向垂直指向壁面，且运动速度比较大时，会对容器壁产生很大的冲击压力。图 4 给出了容器 运动过程中三种典型的液体晃荡现象驻波、行进波和水跃。图 5 为有无阻尼板情况下，两 种数值计算方法对晃动的模拟与试验的对比，SPH 计算较真实地模拟出了试验中自
12、由液面随 着容器加速度变化，从平静到大幅晃动的过程，该算例中，相对于 CFD 方法，SPH 方法模 拟的结果更接近于试验。 图 6 给出了试验、CFD 方法及 SPH 方法所得的测试点 P1 处的压力时间历程曲线。从图 6（a）中可以看到 5s 之后测试点的压力所呈现的双峰特征，这是液体晃动过程中迟滞特性的 一种典型表现。其中前一个峰值是由于液体对容器壁的冲击而产生的高脉冲压力，时间短， 具有比较明显的局部特性。 这可以从图 2 中的加速度曲线分析出：5-10s 这段时间内加速度 变化大，容器壁主要承受水跃（图 4（a） ）或者行进波（图 4（b） ）晃动现象引起的瞬时冲击 压力。后一个峰值也
13、是由于冲击产生的脉冲压力，但是时间相对较长，而且冲击压力峰值变 小，这主要是由于快速并且连续的液体压力作用在没有被完全淹没的容器壁（图 4（c） ） ，压 力特性主要表现为连续性，而不像前一个峰值中所明显表现出的瞬时性。当给容器加装阻尼 板后，液体对容器壁的冲击力度降低，由于阻尼板的面积较大，阻尼板的安装一定程度上可 以等效为增加了容器的充液比，在大充液比情况下，迟滞现象消失，因此图 6（b）中没有明 显的双峰特性，且压力峰值减小，晃动现象更多的表现为驻波运动（图 4（d） ） 。另外，在无 阻尼板工况下，当加速度减至 0 并基本趋于稳定之后，P1 处的压力值开始在初始静态压力值 处附近振荡，
14、并呈现出有阻尼的逐渐衰减的趋势；在有阻尼板工况下，P1 处压力值则很快收 敛于初始静态压力值。 通过将图 6 中的 SPH 计算结果与试验及 CFD 计算结果对比可以发现，SPH 方法计算所 得到的曲线较 CFD 方法不论在变化趋势、压力作用周期还是压力峰值方面都与试验吻合更好。 计算峰值稍低于试验所测峰值，这是由于前面所提及的对小车的加速度激励曲线进行了滤波， 滤去了一部分波峰的缘故。对于 8-10s 时间段，两张图中的计算曲线与试验吻合均不是很理 想，从数值计算所得曲线的斜率可以看出压力变化比较突然，这可能是由于试验中低频加速 度激励测量误差所致12，该时间段加速度测量值变化程度要比真实情
15、况剧烈，导致数值计算 中的液体冲击幅度较大，而在真实试验中，此时间段内容器壁承受的压力主要表现为持续性 的压力，因而压力曲线变化比较平缓。(a) 无阻尼板(b) 有阻尼板图 3 自由液面晃动的三维模拟及对应的对称面粒子速度矢量图 (a) 水跃(无阻尼板) (b) 行进波(无阻尼板) (c) 驻波(无阻尼板) (d) 驻波(有阻尼板)图 4 典型晃荡运动 试验 CFD 计算 SPH 计算 试验 CFD 计算 SPH 计算(a) 无阻尼板 (b) 有阻尼板图 5 液体晃动的数值模拟与试验对比 (a) 无阻尼板 (b) 有阻尼板图 6 SPH 数值计算与试验及 CFD 数值计算压力曲线对比 4结论采用 S
