报告无网格方法在流体力学中的应用-huangcan

11、概述：传统的CFD数值算法在模拟动边界问题时都需要解决一个问题——网格畸变，无论是采用动网格技术还是采用嵌套网格技术都会影响网格质量，进而影响计算精度，无法真实的捕捉动态气动特性，对于大变形问题，网格可能会出现负体积等现象，使得计算无法进行。由此，人们想到采用一种不需要网格的方法来解决网格畸变问题，SPH(Smoothedparticlehydrodynamics)就是应用于流体问题中的一种无网格方法。1977年，Lucy，Gingold,Monaghan,首次应用SPH方法求解天体物理学问题。目前应用于流体方面的有激波模拟，磁流体动力学，多项流，准不可压流，重力流，热传导，热传递和质量流。本文主要内容包含以下三方面：1理论基础。基于SPH方法，详细的推倒离散N-S方程的过程。2SPH程序。基于离散的N-S方程开发出开放性强的SPH程序，便于为以后发展成为实用强的软件打基础。3算例验证。利用开发出的SPH程序数值模拟简单的流动问题，验证算法的可靠性。2、理论基础2.1、SPH的基本方程SPH方程的构造按照两步进行，第一步是积分表示法(场函数核近似法)；第二步为粒子近似法。2.1.1、函数的积分表示法𝑓(𝑥)的积分表达式为：𝑓(𝑥)≈∫𝑓(𝑥′)𝑊(𝑥−𝑥′,ℎ)𝑑𝑥′Ω其中：𝑊称为光滑核函数(光滑函数，光滑核，核函数)，ℎ是定义光滑函数𝑊的影响区域的光滑长度。光滑函数𝑊需要满足的一些条件有：1)𝑊选用偶函数；2)正则化条件(归一化条件)∫𝑊(𝑥−𝑥′,ℎ)𝑑𝑥′Ω=1；23)当光滑长度趋向于零时具有狄拉克函数性质limℎ→0𝑊(𝑥−𝑥′,ℎ)=𝛿(𝑥−𝑥′)；4)紧支性条件𝑊(𝑥−𝑥′,ℎ)=0，当|𝑥−𝑥′|𝜅ℎ时。κ是与𝑥处光滑函数相关的常数，并确定光滑函数的有效范围；5)在点𝑥上的粒子的支持域内任一点𝑥′处有𝑊(𝑥−𝑥′,ℎ)≥0；6)当粒子间的距离增大时，粒子的光滑函数值应该是单调递减的(衰减性)；7)光滑函数应充分光滑。𝑓(𝑥)导数积分的表达式为：〈∇∙𝑓(𝑥)〉=∫𝑓(𝑥′)𝑊(𝑥−𝑥′,ℎ)𝑠∙𝑛⃗𝑑𝑠−∫𝑓(𝑥′)∙∇𝑊(𝑥−𝑥′,ℎ)𝑑𝑥′Ω如果支持域位于问题域的内部，红色表示的面积项积分为零。若支持域与问题域相互交错，则不为零。2.1.2、粒子近似法在SPH系统中，整个系统是由具有独立的质量、占有独立空间的有限个粒子表示的。粒子叠加求和相对应的离散化过程称为粒子近似法。函数离散化的粒子近似式为：〈𝑓(𝒙)〉=∑𝑚𝑗𝜌𝑗𝑓(𝒙𝒋)𝑊(𝒙−𝒙𝑗,ℎ)𝑁𝑗=1函数空间导数的粒子近似为：〈∇∙𝑓(𝒙)〉=−∑𝑚𝑗𝜌𝑗𝑓(𝒙𝑗)∙∇𝑊(𝒙−𝒙𝑗,ℎ)𝑁𝑗=1粒子𝑖处的函数值和函数导数的表达形式分别为：〈𝑓(𝒙𝑖)〉=∑𝑚𝑗𝜌𝑗𝑓(𝒙𝑗)∙𝑊𝑖𝑗𝑁𝑗=1(1)〈∇∙𝑓(𝒙𝑖)〉=∑𝑚𝑗𝜌𝑗𝑓(𝒙𝑗)∙∇𝑖𝑊𝑖𝑗𝑁𝑗=1(2)其中：𝑊𝑖𝑗=𝑊(𝒙𝑖−𝒙𝑗,ℎ)=𝑊(|𝒙𝑖−𝒙𝑗|,ℎ)，∇𝑖𝑊𝑖𝑗=𝒙𝑖−𝒙𝑗𝑟𝑖𝑗𝜕𝑊𝑖𝑗𝜕𝑟𝑖𝑗=𝒙𝑖𝑗𝑟𝑖𝑗𝜕𝑊𝑖𝑗𝜕𝑟𝑖𝑗2.2、拉格朗日型的Navier-Stokes方程流体动力学的基本控制方程基于的物理守恒定律有：(1)质量守恒3(2)动量守恒(3)能量守恒若使用希腊字母上标∝和𝛽表示坐标方向，则可用指标法来表示方程的叠加。因此这种表示法下的连续性、动量和能量方程的表达式分别为：1、连续性方程𝑑𝜌𝑑𝑡=−𝜌𝜕𝒗𝛽𝜕𝒙𝛽(3)2、动量方程(无外力作用)𝑑𝒗𝛼𝑑𝑡=1𝜌𝜕𝜎𝛼𝛽𝜕𝒙𝛽(4)3、能量方程𝑑𝑒𝑑𝑡=𝜎𝛼𝛽𝜌𝜕𝒗𝛼𝜕𝒙𝛽(5)𝑑𝑒𝑑𝑡=−𝑝𝜌𝜕𝒗𝛽𝜕𝒙𝛽+𝜇2𝜌𝘀𝛼𝛽𝘀𝛼𝛽(6)其中：𝜎𝛼𝛽=−𝑝𝛿𝛼𝛽+𝜏𝛼𝛽，𝜏𝛼𝛽=𝜇𝘀𝛼𝛽，𝘀𝛼𝛽=𝜕𝒗𝛽𝜕𝒙𝛼+𝜕𝒗𝛼𝜕𝒙𝛽−23(𝛻∙𝒗)𝛿𝛼𝛽。σ为总应力张量，𝜏为粘性剪应力，𝘀为剪应变率，𝜇为粘性系数。3、SPH程序结构流程结构图为：4图14基于SPH计算方腔自然对流4.4.1、问题描述模型是一个二维封闭方腔，如图2所示，横轴表示x方向，纵轴表示y方向，方腔的高H=0.1m，方腔的长L=0.1m，左壁面为热壁面，右壁面为冷壁面。初始化模块产生边界虚粒子最近相邻粒子搜索计算光滑函数计算密度变化率密度求和如果使用连续密度法如果使用密度求和计算外力计算人工粘滞力更新位置、速度、能量、密度等计算能量和动量改变率计算内力输出模块更新光滑长度若有需要5图2封闭方腔模型4.4.2、计算结果及对比CFD计算结果SPH计算结果图3t=0.2s，流线图3）时间t=0.4s时，传统CFD的计算结果与SPH计算结果的比较：CFD计算结果SPH计算结果图4t=0.4s，压力P分布云图6CFD计算结果SPH计算结果图5t=0.4s，x方向速度u分布云图CFD计算结果SPH计算结果图6t=0.4s，y方向速度v分布云图CFD计算结果SPH计算结果图7t=0.4s，流线图