【2016】SPH方法进出口数值边界条件研究

第３３卷第５期２０１６年９月计　　算　　物　　理　ＣＨＩＮＥＳＥＪＯＵＲＮＡＬＯＦＣＯＭＰＵＴＡＴＩＯＮＡＬＰＨＹＳＩＣＳ　Ｖｏｌ．３３，Ｎｏ．５Ｓｅｐ．，２０１６文章编号：１００１⁃２４６Ｘ（２０１６）０５⁃０５１６⁃０７ＳＰＨ方法进出口数值边界条件研究周　杰，　徐胜利（清华大学航天航空学院，北京　１０００８４）摘　要：基于虚粒子概念，提出进出口数值边界条件处理方法．在进口边界外设置进口区域，并赋予进口虚粒子相应的物理量；出口边界外设置出口和缓冲区域，出口虚粒子的物理量由计算获得，缓冲区的虚粒子初始物理量是指定的；根据每个时间步的流动情况，改变粒子的区域属性，添加／删除相应的粒子．利用拉格朗日形式的ＳＰＨ方法，通过管内流动问题验证数值进出口边界方法的适用性，研究进出口边界条件在激波管、绕流问题中的应用．提出的进出口边界处理方法，避免了边界附近流体粒子积分截断问题，保证流体粒子能够流出边界，激波能够透射边界．关键词：ＳＰＨ；边界条件；虚粒子；激波管中图分类号：Ｖ２１１􀆰３文献标志码：Ａ收稿日期：２０１５－０８－２６；修回日期：２０１６－０２－１７基金项目：博士后面上基金（２０１５Ｍ５８１０８１）资助项目作者简介：周杰（１９８６－），男，博士，助理研究员，研究方向：超声速弹丸入水现象，Ｅ⁃ｍａｉｌ：[e-mail]０　引言ＳＰＨ方法属于Ｌａｇｒａｎｇｅ类的无网格方法，在计算大变形自由表面和界面等问题时有优势［１－２］．临近计算域边界的粒子，只受到边界内的粒子单边作用，边界条件处理不适当会引起计算结果出错，这涉及到边界粒子积分截断、虚粒子布置及其物理量确定等问题［３］．因此，边界条件处理是ＳＰＨ方法难点之一．计算域边界大致分为固壁和开口两类，进出口属于和添加／删除虚粒子有关的开口边界．目前，ＳＰＨ方法处理进出口边界的基本思想是：专门设置边界区域并填满虚粒子，根据边界物理条件类型，借鉴网格类和解析方法，指定边界区域虚粒子物理量，边界区域虚粒子参与计算，并更新全部或部分物理量．Ｌｅｆｆｅ［４－５］参考网格类处理开口边界的思想，提出了ＳＰＨ中处理开口边界条件的方法．针对不可压缩Ｐｏｉｓｅｕｉｌｌｅ流动和绕流问题，Ｈｏｓｅｉｎｉ［６－７］提出了非齐次压力边界条件，建立数值边界层解决边界附近粒子的求解精度问题．根据明渠上、下游流动的情况，Ｆｅｄｅｒｉｃｏ［８］提出了一种设置进出口边界条件方法计算明渠不同类型的水跃问题．应用特征值方法，Ｌａｓｔｉｗｋａ［９］通过边界内的流体粒子物理量外推方式，指定边界区域虚粒子物理量，结合ＳＰＨ方法，计算并验证了一维喷流和二维绕流流场．根据亚超声速进出口边界条件类型，指定边界区域虚粒子物理量，Ｖａｃｏｎｄｉｏ［１０－１１］计算了激波管和定常流动问题．根据已有的边界区域及其虚粒子概念，本文提出进出口边界的处理方法：①对于进口边界，边界区域被赋值的虚粒子总数保持不变；②对于出口边界，在边界区域下游设置缓冲区域，边界区域虚粒子的物理量由计算获得，缓冲区虚粒子物理量是指定的；③根据每个时间步的流动情况，改变粒子的区域属性，添加／删除相应的粒子．通过管内流动问题验证进出口边界方法的适用性，研究进口边界条件在激波管、绕流问题中的应用．１　模型与方法１􀆰１　进出口边界本文提出的进出口边界条件的设置方法如图１所示．进口区域内的虚粒子流入计算域，虚粒子将变为流体实粒子，需要作以下处理：①补充并布置新的虚粒子，进口区域内的粒子数目保持不变；②将进入计算域的粒子由虚粒子类型转变为流体实粒子，参与后续的流体计算．进口区域的虚粒子在计算过程中遵循以下公式　第５期周　杰等：ＳＰＨ方法进出口数值边界条件研究ｄｖｄｔ＝ｄｅｄｔ＝０，（１）式中，ｅ、ｖ分别表示进口区域虚粒子的内能和速度．值得注意的是，根据计算条件，新添加的进口边界虚粒子初始参量是由计算条件直接给定的．图１　进出口边界附近粒子的处理方法Ｆｉｇ．１　Ｍｅｔｈｏｄｓｏｆｄｅａｌｉｎｇｗｉｔｈｐａｒｔｉｃｌｅｓｎｅａｒｉｎｌｅｔａｎｄｏｕｔｌｅｔｂｏｕｎｄａｒｙ对出口边界，流体实粒子会流进出口区域，将变为出口虚粒子，这时需要作以下处理：①将流体实粒子转变为出口虚粒子；②在虚粒子下游设置缓冲区（或海绵层），计算中需要保证虚粒子和流体粒子计算精度，海绵粒子同样属于虚粒子，并且直接指定其初始参数；③缓冲区虚粒子参与计算（物理量按ＳＰＨ方程更新），但在海绵层下边界由于光滑函数截断产生扰动波并向计算区传播，间隔一定的计算步（扰动波传播到边界区域之前），需要将海绵层粒子重新赋值（一般为初始条件），也就是说，扰动波会影响海绵层粒子计算精度但不会影响到流体粒子的计算精度．计算过程中，出口边界区域虚粒子遵循下列公式：ｄｐｄｔ＝０　或　ｄρｄｔ＝０，（２）式中，ｐ、ρ分别表示出口区域虚粒子的压力和密度．当流体粒子进入出口区域转变为虚粒子时，在计算过程中，出口区域内虚粒子的ｐ或ρ是保持不变的．缓冲区域（理论上宽度大于支持域半径Ｒ即可）中的虚粒子是参与计算的，海绵层下游边界由于积分截断，会产生向内传播的扰动波，需要及时（根据海绵层宽度和扰动波的速度确定）对海绵层内虚粒子的物理量重新赋值．图２是实现进出口边界条件的计算流程图，主要步骤：①ｔｎ时刻，判断流体区域的实粒子是否流入出口区域，并将流入的流体实粒子转化为出口虚粒子；②ｔｎ时刻，判断出口区域的虚粒子是否流出计算区域，并将流出的虚粒子删除，保存其粒子标号；③ｔｎ时刻，判断进口区域的虚粒子是否流入流体区域，并根据流入情况，在进口区域添加新的进口虚粒子，粒子标号应从步骤②中选取（当步骤②中粒子标号不够时，自动添加新的粒子标号，这样处理是为节约内存，提高计算效率）；④根据需要决定是否重新赋值海绵层内虚粒子的参量；⑤完成上述步骤后，进行ＳＰＨ计算，并且开始ｔｎ＋１时刻的循环．固壁边界条件的设置，如图３所示．在每个时间步中当流体粒子到固壁边界的距离小于支持域半径Ｒ，则在边界外对称布置边界虚粒子，赋予边界粒子与对应流体粒子相同的物理量，法向速度相反．为了防止流体粒子非物理穿透固壁边界，边界虚粒子引入与Ｌｅｎｎａｒｄ⁃Ｊｏｎｅｓ势函数类似的排斥力学机制．当流体粒子向边界运动时，与边界粒子的距离ｒｉｊ小于初始间距ｒ０，此时边界粒子会在法向ｎ方向上对流体粒子施加排斥力Ｆｎｉｊ，如下式所示：Ｆｎｉｊ＝Ｄ［（ｒ０／ｒｉｊ）ｌ１－（ｒ０／ｒｉｊ）ｌ２］（ｒｎ／ｒｉｊ），　　ｒ０／ｒｉｊ≥１，（３）式中，Ｄ是由具体问题而定的参数，大小取与速度平方值同一量级，参数ｌ１和ｌ２一般取值为１２和４，可依据具体问题进行修改．１􀆰２　计算模型建立管内流、绕流和激波管问题的计算模型，评估进出口边界条件方法的准确性，如图４所示．计算模型７１５计　　算　　物　　理第３３卷　图２　进出口边界条件计算流程Ｆｉｇ．２　Ｃａｌｃｕｌａｔｉｏｎｐｒｏｃｅｓｓｏｆｉｎｌｅｔａｎｄｏｕｔｌｅｔｂｏｕｎｄａｒｙ图３　固壁边界粒子的布置示意图Ｆｉｇ．３　Ｌａｙｏｕｔｄｉａｇｒａｍｓｏｆｓｏｌｉｄｂｏｕｎｄａｒｙｖｉｒｔｕａｌｐａｒｔｉｃｌｅｓ（ａ）为管内流动，计算区域分为进口区域、流体区域及出口区域，模型的上、下边界为固壁，其中流体区域的长、宽分别Ｌ１＝０􀆰４，Ｌ２＝０􀆰２，进／出口区域的宽度ΔＬ１≥Ｒ（支持域半径），粒子的运动速度Ｖ１＝０􀆰２，监控粒子ｉ和中轴线（Ｙ＝０）上粒子物理量的变化，验证进／出口边界方法的准确性．计算模型（ｂ）为激波管问题，高压区域气体的初始物理量：ｐｈ＝１􀆰０，ｅｈ＝２􀆰５，ρｈ＝１􀆰０；低压区域气体的初始物理量：ρｌ＝０􀆰１７９５，ｅｌ＝１􀆰７９５，ρｌ＝０􀆰２５；其中，ｐ，ｅ，ρ分别为气体的压力、内能和密度，模型中激波管的高、低压区域的长度Ｌ３＝Ｌ４＝０􀆰６，管径Ｌ７＝０􀆰１５，下边界为开口边界，其余为固壁边界，上述物理量均为无量纲．计算模型（ｃ）为绕流问题，在流体区域添加圆形固壁边界（ｒ＝０􀆰２ｍ），因为文中的镜像固壁边界条件处理圆形固壁存在一定困难，故参考Ｖａｃｏｎｄｉｏ［１１］提出的动力学边界条件即在圆柱区域设置气体粒子，粒子的密度、动量及能量等物粒子量按照ＳＰＨ方法进行更新，但是粒子的位置（ｄｘ＝０）保持不变，依靠粒子之间的内力阻止粒子之间的穿透．流体区域的长、宽分别Ｌ６＝２􀆰０ｍ，Ｌ７＝１􀆰８ｍ，气体的运动速度Ｖ２＝５００ｍ·ｓ－１，进／出口区域的宽度ΔＬ２≥Ｒ（支持域半径），模型的上、下边界为固壁，左、右分别为进、出口边界．模拟中使用理想气体状态方程ｐ＝（γ－１）ρｅ，（１）式中，γ＝１．４．采用蛙跳法求解具有拉格朗日性质的ＳＰＨ方法的Ｎ⁃Ｓ方程，引入Ｍｏｎａｇｈａｎ型人工粘度［１２－１３］，采用Ｍｏｎａｇｈａｎ［１４］提出的Ｂ⁃样条函数的光滑函数，用积分表示法来近视某一点的场函数值，时间步长满足ＣＦＬ８１５　第５期周　杰等：ＳＰＨ方法进出口数值边界条件研究图４　计算模型Ｆｉｇ．４　Ｃａｌｃｕｌａｔｉｏｎｍｏｄｅｌｓ条件［１５］．２　结果分析和讨论２􀆰１　算例验证计算模型（ａ）中，粒子的初始参数速度ｖ＝０􀆰２、压力ｐ＝１􀆰０、内能ｅ＝２􀆰５，上、下固壁边界为滑移边界，计算中没有考虑气体粘性．图５为ｔ＝０􀆰２５ｓ时流体粒子的压力云图，由图可知流体区域的粒子压力分布均匀（均为１􀆰０），临近边界（固壁边界、进出口边界）处的流体粒子并没有产生非物理扰动．由此可知文中提出的进出边界条件处理方法是有效的．图５　ｔ＝０􀆰２５ｓ流体粒子的压力云图Ｆｉｇ．５　Ｐｒｅｓｓｕｒｅｃｏｎｔｏｕｒｓｏｆｆｌｕｉｄｐａｒｔｉｃｌｅｓａｔｔ＝０􀆰２５ｓ图６为模型中粒子ｉ的压力ｐｉ、内能ｅｉ随时间的变化，由图可知，粒子ｉ在向下游运动时，其压力ｐｉ、内能ｅｉ均能保持不变．图７为ｔ＝０􀆰２５ｓ时刻，模型中轴线Ｙ＝０上，流体粒子的压力ｐ、内能ｅ随距离的变化也均保持不变．综上所述：本文提出的进出口边界条件，边界粒子删除、添加方法是有效的，边界虚粒子物理量计算准则是合理的，在出口边界处添加缓冲区域，可以有效的避免出口边界由于积分截断而产生非物理扰动．９１５计　　算　　物　　理第３３卷　图６　粒子ｉ的压力、内能随时间的变化Ｆｉｇ．６　Ｐｒｅｓｓｕｒｅａｎｄｉｎｔｅｒｎａｌｅｎｅｒｇｙｏｆｆｌｕｉｄｐａｒｔｉｃｌｅｉｖａｒｙｉｎｇｗｉｔｈｔｉｍｅ图７　ｔ＝０􀆰２５ｓ中轴线Ｙ＝０上粒子的压力、内能Ｆｉｇ．７　ＰｒｅｓｓｕｒｅａｎｄｉｎｔｅｒｎａｌｅｎｅｒｇｙｄｉｓｔｒｉｂｕｔｉｏｎｏｎａｘｉａｌＹ＝０ａｔｔ＝０􀆰２ｓ２􀆰２　方法应用应用文中提出的边界条件处理方法，计算开口激波管问题，如计算模型（ｂ），图８为典型时刻的激波管压力分布云图．其中，ｔ＝０􀆰００ｓ时刻，激波管的初始状态；ｔ＝０􀆰２ｓ、０􀆰４ｓ时刻，激波管膜片破碎后，形成激波向低压区域传播的压力云图；ｔ＝０􀆰５ｓ时刻，激波透射出边界，且出口边界处没有明显的扰动产生．由激波管压力云图可知：（１）设置了海绵层有效地避免了积分截断引起的非物理扰动；（２）对出口区域内虚粒子的处理方法，使得激波能够透射出边界．图８　典型时刻的激波管内压力云图Ｆｉｇ．８　Ｐｒｅｓｓｕｒｅｃｏｎｔｏｕｒｓｉｎｓｈｏｃｋｔｕｂｅａｔｒｅｐｒｅｓｅｎｔａｔｉｖｅｔｉｍｅ图９给出了ｔ＝０􀆰２ｓ时激波管中轴线上气体粒子的压力、密度分布．从图中可以看到，ｘ＝０􀆰３处观察到冲击波，其分辨率约为几个光滑长度；膨胀波位于［－０􀆰２３，－０􀆰０７］之间；不连续接触面位于ｘ＝０􀆰１３处，且出口边界处没有产生非物理扰动．图１０为ｔ＝０􀆰４５ｓ时激波管中轴线上气体粒子的压力、密度分布．由图可知冲击波已经运动出边界，且没有出现反射现象，不连续接触面位于ｘ＝０􀆰１３处．综上可知：文中的出口边界处理方法，可以有效地避免出口边界由于积分截断而产生非物理扰动，同时又能保证激波透射出边界．图１１为圆柱绕流问题的压力云图，其中（ａ）为ＳＰＨ粒子云图，（ｂ）为插值云图．由图可知，气流速度为超声速（５００ｍ·ｓ－１）