一种改进的FPM方法.pdf

本发明涉及一种改进的FPM方法，属于计算力学技术领域，该方法包括以下步骤：首先布置流体粒子和边界粒子，其次对粒子的物理属性初始化，接下来利用离散格式公式对任意函数及其导数进行近似，最后按照近似公式进行计算得出数值模拟结果。对比FPM方法，本发明方法能获得对称性良好的系数矩阵，提高对导数的模拟精度，降低对核函数的要求。

权利要求书
1.  一种改进的FPM方法，其特征在于，包括以下步骤：
步骤1、布置实粒子和虚粒子；
步骤2、对粒子的物理属性初始化；
步骤3、利用下述离散格式对任意函数及其导数进行近似：
对于一维问题：
f′fx′=AΣj=1NfWmρΣj=1Nf(x-x′)Wmρ---(1)]]>
A=Σj=1NWmρΣj=1N(x-x′)WmρΣj=1N(x-x′)WmρΣj=1N(x-x′)(x-x′)Wmρ-1---(2)]]>
对于二维问题：
f′fx′fy′=AΣj=1NfWmρΣj=1Nf(x-x′)WmρΣj=1Nf(y-y′)Wmρ---(3)]]>
A=Σj=1NWmρΣj=1N(x-x′)WmρΣj=1N(y-y′)WmρΣj=1N(x-x′)WmρΣj=1N(x-x′)(x-x′)WmρΣj=1N(y-y′)(x-x′)WmρΣj=1N(y-y′)WmρΣj=1N(x-x′)(y-y′)WmρΣj=1N(y-y′)(y-y′)Wmρ-1---(4)]]>
对于三维问题：
f′fx′fy′fz′=AΣj=1NfWmρΣj=1Nf(x-x′)WmρΣj=1Nf(y-y′)WmρΣj=1Nf(z-z′)Wmρ---(5)]]>
A=Σj=1NWmρΣj=1N(x-x′)WmρΣj=1N(y-y′)WmρΣj=1N(z-z′)WmρΣj=1N(x-x′)WmρΣj=1N(x-x′)(x-x′)WmρΣj=1N(y-y′)(x-x′)WmρΣj=1N(z-z′)(x-x′)WmρΣj=1N(y-y′)WmρΣj=1N(x-x′)(y-y′)WmρΣj=1N(y-y′)(y-y′)WmρΣj=1N(z-z′)(y-y′)WmρΣj=1N(z-z′)WmρΣj=1N(x-x′)(z-z′)WmρΣj=1N(y-y′)(z-z′)WmρΣj=1N(z-z′)(z-z′)Wmρ-1---(6)]]>
其中N表示粒子i的相关粒子总数，fˊ表示当前时刻函数在粒子i处的函数值，表示当前时刻函数对x方向的偏导在粒子i处的值，同理，f表示上一时刻函数在粒子j处的函数值，W表示核函数，m表示粒子j的质量，ρ表示粒子j的密度，x-x’、y-y’和z-z’分别表示粒子i与粒子j在x、y和z方向上的位移；
步骤4、按照步骤3中公式计算得出数值模拟结果。

说明书一种改进的FPM方法
技术领域
本发明涉及一种改进的FPM方法，属于计算力学技术领域。
背景技术
流体力学中应用最为广泛的一种无网格方法是SPH(smooth particle hydrodynamics)方法，SPH方法在数值模拟过程中不需要网格，易于解决基于网格的数值方法难以解决的大变形、移动边界等问题，所以SPH方法被应用于自由表面流、多相流动、爆炸与冲击等问题。随着SPH方法的广泛应用，传统SPH方法存在的粒子不一致性、应力不稳定性等缺陷逐渐暴露出来，所以有必要对SPH方法进行修正。
在SPH方法的基础上，2005年的Liu等人提出了FPM方法。FPM方法是对任意函数进行泰勒展开，然后将表达式与核函数及其导数相乘并进行积分，联立求解积分后非线性相关的方程即可得到任意函数及其导数的离散格式。FPM方法能确保粒子近似方案的一阶连续性，是解决传统SPH方法精度不高的一个有效方法。但是FPM方法存在以下缺点：
(1)该方法对函数导数近似的精度较低；
(2)该方法对核函数的要求比较严格。
发明内容
本发明的目的是为了减弱无网格粒子方法对核函数的要求，提高对函数导数的模拟精度，提出了一种改进的FPM方法。
本发明是通过以下技术方案实现的：
一种改进的FPM方法，具体实现步骤如下：
步骤1、布置实粒子即流体粒子和虚粒子即边界粒子；
步骤2、对粒子的物理属性初始化；
步骤3、利用下述离散格式对任意函数及其导数进行近似：
对于一维问题：
f′fx′=AΣj=1NfWmρΣj=1Nf(x-x′)Wmρ---(1)]]>
A=Σj=1NWmρΣj=1N(x-x′)WmρΣj=1N(x-x′)WmρΣj=1N(x-x′)(x-x′)Wmρ-1---(2)]]>
对于二维问题：
f′fx′fy′=AΣj=1NfWmρΣj=1Nf(x-x′)WmρΣj=1Nf(y-y′)Wmρ---(3)]]>
A=Σj=1NWmρΣj=1N(x-x′)WmρΣj=1N(y-y′)WmρΣj=1N(x-x′)WmρΣj=1N(x-x′)(x-x′)WmρΣj=1N(y-y′)(x-x′)WmρΣj=1N(y-y′)WmρΣj=1N(x-x′)(y-y′)WmρΣj=1N(y-y′)(y-y′)Wmρ-1---(4)]]>
对于三维问题：
f′fx′fy′fz′=AΣj=1NfWmρΣj=1Nf(x-x′)WmρΣj=1Nf(y-y′)WmρΣj=1Nf(z-z′)Wmρ---(5)]]>
A=Σj=1NWmρΣj=1N(x-x′)WmρΣj=1N(y-y′)WmρΣj=1N(z-z′)WmρΣj=1N(x-x′)WmρΣj=1N(x-x′)(x-x′)WmρΣj=1N(y-y′)(x-x′)WmρΣj=1N(z-z′)(x-x′)WmρΣj=1N(y-y′)WmρΣj=1N(x-x′)(y-y′)WmρΣj=1N(y-y′)(y-y′)WmρΣj=1N(z-z′)(y-y′)WmρΣj=1N(z-z′)WmρΣj=1N(x-x′)(z-z′)WmρΣj=1N(y-y′)(z-z′)WmρΣj=1N(z-z′)(z-z′)Wmρ-1---(6)]]>
其中N表示粒子i的相关粒子总数，fˊ表示当前时刻函数在粒子i处的函数值，表示当前时刻函数对x方向的偏导在粒子i处的值，同理，f表示上一时刻函数在粒子j处的函数值，W表示核函数，m表示粒子j的质量，ρ表示粒子j的密度；
(4)步骤4、按照步骤3中公式计算得出数值模拟结果。
有益效果
本发明对比已有技术具有如下优点：
(1)MFPM方法的系数矩阵具有良好的对称性；
(2)MFPM方法提高了对导数的近似精度；
(3)MFPM方法降低了对核函数的要求。
附图说明
图1为本发明实施例1初始粒子分布示意图；
图2为本发明实施例1FPM方法和MFPM方法计算函数f＝x5的结果对比；
图3为本发明实施例1FPM方法和MFPM方法计算函数f＝x5一阶导数fx的结果对比；
图4为本发明实施例1FPM方法和MFPM方法计算函数f＝x5二阶导数fxx的结果对比；
图5为本发明实施例1FPM方法和MFPM方法计算函数f＝x5的误差结果对比；
图6为本发明实施例1FPM方法和MFPM方法计算函数f＝x5一阶导数fx的误差结果对比；
图7为本发明实施例1FPM方法和MFPM方法计算函数f＝x5二阶导数fxx的误差结果对比；
图8为本发明实施例1分别采用21个粒子和51个粒子由MFPM方法计算函数f＝x5的结果；
图9为本发明实施例1分别采用21个粒子和51个粒子由MFPM方法计算函数f＝x5一阶导数fx的结果；
图10为本发明实施例1分别采用21个粒子和51个粒子由MFPM方法计算函数f＝x5二阶导数fxx的结果；
图11(a)为本发明实施例2初始粒子分布示意图；
图11(b)为本发明实施例2腔内剪切流动问题边界条件示意图；
图12(a)为本发明实施例2分别采用FPM方法和MFPM方法模拟腔内剪切流得到的方腔水平中心线上的y方向速度分布；
图12(b)为本发明实施例2分别采用FPM方法和MFPM方法模拟腔内剪切流得到的方腔垂直中心线上的x方向速度分布；
图13(a)为本发明实施例2采用MFPM方法模拟腔内剪切流得到的方腔水平中心线上的x方向速度分布；
图13(b)为本发明实施例2采用MFPM方法模拟腔内剪切流得到的方腔垂直中心线上的y方向速度分布；
具体实施方式
为了使本发明的目的、技术方案和优点更加清晰明白，下面结合实施例和附图，对本发明实施例做进一步的详细说明。
实施例1下面以一维线性函数f＝x5来说明本发明的技术方案。
步骤1、布置实粒子即流体粒子和虚粒子即边界粒子：
以一维线性函数f＝x5为例，x的取值范围为-0.5-0.5，粒子分布方式如图1所示，粒子间距为Δ，图中实心圆点为实粒子，空心圆点为虚粒子。
步骤2、对粒子的物理属性初始化：
粒子的密度取为ρ＝1；粒子的质量m＝ρΔ；光滑长度h＝1.2Δ。
步骤3、利用MFPM法的离散格式对函数f＝x5及其导数进行近似；其近似原理如下：
对任意一维函数f(x)在点xˊ处进行泰勒展开，保留至一阶导数可得：
f(x)＝f(x′)+fx(x′)·(x-x′)+o(h2)    (7)
忽略式(7)中的高阶项，在式(7)两边同时乘以核函数W和(x-x′)W并对x进行积分得：
∫f(x)Wdx=f(x′)∫Wdx+▿f(x′)·∫(x-x′)Wdx---(8)]]>
∫f(x)(x-x′)Wdx=f(x′)∫(x-x′)Wdx+fx(x′)·∫(x-x′)(x-x′)Wdx---(9)]]>
将式(8)和式(9)在计算域上利用粒子进行离散，用求和符号代替积分符号，被积函数dx＝m/ρ，则式(8)和式(9)可表达为下述粒子近似式：
Σj=1Nf(x)Wmρ=f(x′)Σj=1NWmρ+fx(x′)·Σj=1N(x-x′)Wmρ---(10)]]>
Σj=1Nf(x)(x-x′)Wmρ=f(x′)Σj=1N(x-x′)Wmρ+fx(x′)·Σj=1N(x-x′)(x-x′)Wmρ---(11)]]>
为了表达简便，下面用f、f′和fx′分别表示f(x)、f(x′)和fx(x′)。
联立求解式(10)和式(11)，即可得到任意函数及其导数的粒子近似式：
f′fx′=AΣj=1NfWmρΣj=1Nf(x-x′)Wmρ---(12)]]>
式中，系数矩阵A的表达式为：
A=a11a12a21a22=Σj=1NWmρΣj=1N(x-x′)WmρΣj=1N(x-x′)WmρΣj=1N(x-x′)(x-x′)Wmρ-1---(13)]]>
利用式(12)给出的粒子近似方案对任意函数及其一阶导数进行近似。
函数的二阶导数通过函数的一阶导数来计算，即将函数的一阶导数视为原函数，通过式(12)对函数的一阶导数再求一次导，则函数二阶导数的近似表达式为：
fxx′=a21Σj=1NfxWmρ+a22Σj=1Nfx(x-x′)Wmρ---(14)]]>
式中，fxxˊ表示函数在粒子j处的二阶导数值。
对于二、三维情况，以此类推，即可得到二、三维函数的离散格式，函数f可替换为速度、温度等任意场函数。
步骤4、利用MFPM方法对函数f＝x5进行数值模拟。
4.1核函数采用应用较为广泛的B-spline函数
W(R,h)=αd23-R2+12R30≤R<116(2-R)31≤R<202≤R---(15)]]>
其中R＝|x-x′|/h，在一维空间下αd为1/h。
图2至图7均为核函数为B-spline函数，粒子间距Δ＝0.05的模拟结果。图2、图3、图4为采用FPM方法与MFPM方法分别对函数f＝x5及其一阶导数、二阶导数的模拟结果。从图2、图3、图4可以看出，FPM方法与MFPM方法所得模拟结果均与解析解基本吻合。图5、图6、图7为采用FPM方法与MFPM方法模拟函数f＝x5及其一阶导数、二阶导数与解析解的误差值。从图5、图6、图7中可以看出，MFPM方法与FPM方法对函数值的模拟精度相同，但MFPM方法对函数f＝x5一阶及二阶导数与解析解的误差比FPM方法的误差小，说明MFPM方法提高了FPM方法对函数导数的模拟精度。
4.2核函数采用常函数
W(R,h)=10≤R<202≤R---(16)]]>
当核函数采用常函数时，在FPM方法中，由于Wx＝0，系数矩阵A不可逆，无法进行数值模拟，因此下面只采用MFPM方法进行数值模拟。图8、图9、图10分别给出了在核函数为常函数的条件下，分别采用21个粒子(Δ＝0.05)和51个粒子(Δ＝0.01)的MFPM方法对函数f＝x5及其一阶和二阶 导数的模拟结果。从图8、图9、图10中可以看出，粒子间距较小时，MFPM方法对函数及其导数的模拟精度较高，当粒子间距为Δ＝0.01时，虽然核函数比较简单，不满足FPM方法对核函数的要求，但MFPM方法所得的模拟结果与解析解基本一致，说明了MFPM方法降低了对核函数的要求。
实施例2下面以腔内剪切流来验证本发明的技术方案。
步骤1、布置实粒子即流体粒子和虚粒子即边界粒子：
以二维腔内剪切流为例，计算域为边长为L的方腔，粒子分布方式如图11(a)所示，图中实心圆点为实粒子，空心圆点为虚粒子。
步骤2、对粒子的物理属性初始化：
方腔的边长L＝0.001m，初始时刻粒子间距dx＝dy＝0.0002m，光滑长度h＝1.2dx，粒子的密度为标况下水的密度ρ＝1000kg/m3，动力粘度系数μ＝10-3kg/(m.s)，粒子的质量m＝ρdxdy。腔内剪切流动的边界条件如图11(b)所示，上壁面粒子x方向的速度vx＝v∞保持不变，y方向的速度vy＝0左右壁面和下壁面的速度为0，其中v∞是参考速度v∞＝0.001m/s。
步骤3、利用MFPM法的二维离散格式即式(3)和式(4)对腔内剪切流的控制方程进行近似：
腔内剪切流的流体力学控制方程为：
dρdt=-ρ(&PartialD;vx&PartialD;x+&PartialD;vy&PartialD;y)dvxdt=-1ρ&PartialD;p&PartialD;x+1ρμ(&PartialD;&epsiv;xx&PartialD;x+&PartialD;&epsiv;xy&PartialD;y)dvydt=-1ρ&PartialD;p&PartialD;y+1ρμ(&PartialD;&epsiv;xy&PartialD;x+&PartialD;&epsiv;yy&PartialD;y)---(17)]]>
其中应变率分量εαβ的表达式为：
&epsiv;xx=43&CenterDot;&PartialD;vx&PartialD;x-23&CenterDot;&PartialD;vy&PartialD;y&epsiv;xy=&epsiv;yx=&PartialD;vy&PartialD;x+&PartialD;vx&PartialD;y&epsiv;yy=43&CenterDot;&PartialD;vy&PartialD;y-23&CenterDot;&PartialD;vx&PartialD;x---(18)]]>
利用式(3)和式(4)对连续方程进行离散可得：
dρidt=-ρia21Σjvx,jWmjρj+a22Σjvx,j(xj-xi)Wmjρj+a23Σjvx,j(yj-yi)Wmjρj+a31Σjvy,jWmjρj+a32Σjvy,j(xj-xi)Wmjρj+a33Σjvy,j(yj-yi)Wmjρj---(19)]]>
其中i和j分别表示i粒子和j粒子，amn表示式(4)中矩阵A的分量，ρi和ρj分别表示上一时刻粒子i和粒子j的密度，vx,j和vy,j表示上一时刻j粒子分别在x方向和y方向的速度分量，mj表示j粒子的质量。
离散后的x方向动量方程为：
dvxidt=-1ρi(a11ΣjpjWmjρj+a12Σjpj(xj-xi)Wmjρj+a13Σjpj(yj-yi)Wmjρj)+μρi(a21Σj&epsiv;xx,jWmjρj+a22Σj&epsiv;xx,j(xj-xi)Wmjρj+a23Σj&epsiv;xx,j(yj-yi)Wmjρj)+μρi(a31Σj&epsiv;xy,jWmjρj+a32Σj&epsiv;xy,j(xj-xi)Wmjρj+a33Σj&epsiv;xy,j(yj-yi)Wmjρj)---(20)]]>
其中pj表示上一时刻粒子j的压力，εαβ,j表示粒子j的应变率分量。同理也可得到y方向的动量方程以及应变率分量的离散形式，这里不再一一给出。
步骤4、利用MFPM方法对腔内剪切流进行数值模拟，得出模拟结果。
图12给出了核函数为式(15)的B-spline函数，采用FPM方法以及MFPM方法与2003年Liu用有限元方法(FEM)所得到的方腔中心线上的速度分布结果。其中图12(a)是方腔水平中心线上的y方向速度分布，图12(b)是方腔垂直中心线上的x方向速度分布，从图12中可以看出，MFPM方法的结果比FPM方法的结果更加接近文献的结果。
图13给出了核函数为式(16)的常函数，采用MFPM方法与2003年Liu所得到的方腔中心线上的速度分布结果。其中图13(a)是方腔水平中心线上的y方向速度分布，图13(b)是方腔垂直中心线上的x方向速度分布，从图13中可以看出，MFPM方法得到的结果与文献基本一致。然而，当核函数为常函数时，FPM方法系数矩阵的行列式为0，因此系数矩阵不可逆，无法对腔内剪切流进行模拟，这证明了MFPM方法比FPM方法对核函数的要求更低。
以上所述仅为本发明的具体实施例，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。