删除或更新信息,请邮件至freekaoyan#163.com(#换成@)

最优输运无网格方法及其在液滴表面张力效应模拟中的应用

本站小编 Free考研考试/2021-12-29

摘要:基于网格的数值方法(如有限元、有限体积、有限差分等)在大变形、不连续等问题中遇到挑战, 因此人们提出了多种无网格方法. 最优输运无网格方法是一种新型拉格朗日无网格方法, 但是继承了有限元方法在边界表征、边界处理等方面的优势, 在表面张力模拟中具有较大潜力. 基于拉格朗日方程, 通过将表面张力势能加入拉格朗日函数, 得到的表面张力广义力精确地作用在流体表面, 而且表面张力系数是唯一的输入参数. 给出了最优输运无网格方法轴对称离散格式. 通过对二维/三维泊肃叶流、静止和振动的液滴、液滴变形等典型问题的仿真分析, 验证了最优输运无网格方法在表面张力问题模拟中的精度和收敛性.
关键词: 最优输运无网格方法/
表面张力/
拉格朗日方程/
轴对称

English Abstract


--> --> -->
精确可靠的模拟表面张力对液滴变形与运动的影响在环境工程、微纳米工程以及医药工程等领域有着十分重要的应用. 采用基于网格的数值方法模拟此类问题往往需要复杂的界面追踪算法, 典型的如流体体积函数法[1,2]和水平集方法[3]等. 因此, 越来越多的研究者采用拉格朗日粒子类方法模拟表面张力相关现象, 利用粒子本身实现对边界的自动追踪, 极大地简化了算法, 其中最常用的数值方法是光滑粒子动力学(smoothed particle hydrodynamics, SPH)[4,5]方法. 在SPH方法中, 主要有两种方法来考虑表面张力的影响: 粒子间相互作用力(inter-particle interaction force, IIF)方法[6-10]和连续力方法(continuum surface force, CSF) [11]. IIF方法来源于如下思想, 即表面张力来自于液体表面分子之间和内部分子之间作用力的差异. 受此启发, 在SPH方法中, 可以通过在宏观粒子之间加入额外远程吸引力和近程排斥力来近似模拟表面张力. 这种方法简单灵活, 粒子之间的额外作用力形式非常多. 这种方法的缺点是表面张力系数不是输入参数, 而是需要先模拟某个标准算例, 然后测量与这种宏观粒子之间吸引力等价的表面张力系数. 此外, 这种方法往往存在不收敛的问题. CSF方法通过给不同的流体赋以不同的色值(color)来区分不同的流体, 利用不连续的色函数(color function)求解表面法向以及曲率来计表面张力, 然后将只作用在表面的力转换为连续的体积力. SPH方法中, 界面法向向量及其曲率的计算误差一般较大. 无论哪种表面张力模拟方法, SPH方法都存在如下固有缺点: 1) SPH形函数不能严格满足1阶精度(特别是粒子分布不均匀或者自由界面处), 导致收敛性不好, 常常需要对形函数及其导数进行修正[12,13]; 2)形函数不满足Kronecker delta性质, 因此位移边界条件处理复杂, 一般需要通过各种虚粒子处理边界条件. 因此, 本文尝试采用一种新型无网格方法模拟表面张力问题, 即最优输运无网格(optimized transportation meshfree, OTM)方法[14].
OTM方法采用粒子系统离散连续介质, 利用局部最大熵(local maximum entropy, LME)[15-17]形函数构造连续场, 并结合数学中的最优输运理论计算系统的动能积分. 相比于其他无网格方法, OTM方法在精度、收敛性、边界处理以及数学理论基础上均具有优势. OTM方法具有无网格方法能够处理大变形的优点, 同时又继承了有限元方法的部分优点, 如形函数严格满足一阶精度、形函数在边界处满足弱Kronecker delta性质等. 当前, OTM方法已经成功应用于金属侵彻[18]、超高速碰撞[19-21]、裂纹扩展[22]以及热传导等问题[23], 但是将其应用于表面张力模拟还未见文献报道.
本文从拉格朗日方程出发, 从能量的角度考虑表面张力的影响, 即将表面能加入到拉格朗日函数中, 得到的表面张力广义力精确地作用在流体表面, 而且表面张力系数是唯一的输入参数. 为了减少计算量, 推导了轴对称形式的OTM离散格式, 用于三维算例数值模拟, 并给出了详细的计算流程. 通过模拟泊肃叶流并与解析解对比, 验证了本文基本离散格式和轴对称处理的正确性. 通过模拟静止液滴, 并与Young-Laplace公式对比, 验证了表面张力处理的正确性. 模拟了液滴的周期振荡问题并与解析解对比, 严格验证了液滴中速度分布的空间收敛性. 最后, 模拟了液滴在重力、表面张力以及界面共同作用下的形变问题, 考察了液滴中的压力分布并与理论解对比, 进一步验证了本文OTM离散格式.
2
2.1.连续介质离散
-->OTM方法中的连续介质离散介于拉格朗日有限元与SPH方法之间. OTM方法将有限元网格中每个小网格转换为1个粒子, 位于小网格形心. 对于三角形, 形心为三条中线的交点. 粒子存储所有的物理量, 与SPH方法类似. 同时, 保留所有的网格节点, 其中存储了速度和位置, 如图1所示. 初始时刻, 边界上只有节点, 且节点随着流场运动. 边界节点很好地表征了连续介质的边界位置, 与拉格朗日有限元方法类似. 此外, 由于局部最大熵形函数满足弱Kronecker delta性质, 边界条件处理简单. 抛弃有限元方法中单元和节点之间的固定连接关系, 而是采用近邻粒子搜索方式动态确定, 从而可以处理大变形, 这与SPH方法类似.
图 1 利用有限元网格将连续介质离散为粒子和节点 (a)有限元网格; (b)用粒子与节点离散连续介质
Figure1. Continuum discretized by particles (red symbols) and nodes (white symbols) in the OTM method by means of finite element mesh: (a) Finite element mesh; (b) continuum discretized by particles and nodes.

2
2.2.LME形函数
-->最优输运无网格方法与更新拉格朗日有限元方法最大的1个区别是采用LME形函数, 因此不需要网格, 属于无网格方法. 给定1个节点集$ {\boldsymbol{x}}_{a}(a= $$ 1, 2, \cdots , n) $, 每个节点的形函数记为$ {N_a}({\boldsymbol{x}}) $. 利用形函数可以从离散数值插值出1个连续场:
$ {u^h}({\boldsymbol{x}}) = \sum\limits_a {{N_a}({\boldsymbol{x}}){u_a}} , $
其中, $ {u^h}({\boldsymbol{x}}) $是插值出来的连续场, $ {u_a} $是物理量u在节点a处的值. 在有限元和SPH等方法中, 形函数及其导数都具有简单解析表达式, 但是LME形函数没有显式的解析表达式, 而是需要数值求解.
如果要求(1)式可以精确重构常数场, 那么有$\displaystyle \sum\nolimits_a {{N_a}({\boldsymbol{x}}) = 1}$, 即形函数满足归一化条件. 此外, 对于固定的x, 如果形函数大于零, 那么形函数$ {N}_{a}(a=1, 2, \cdots , n) $可以看成是某个事件的概率分布, 其熵为
$ H = - \sum\limits_a {{N_a}\lg ({N_a})} , $
其中, 由于x是固定值所以将其省略. 为了减少计算量, 形函数一般都是局域的, 即当$ \left| {{\boldsymbol{x}} - {{\boldsymbol{x}}_a}} \right| $大于某一临界值时, 形函数$ {N_a}({\boldsymbol{x}}) $为零. 因此, 可以定义如下形函数影响域平均大小:
$ U = \sum\limits_a {{N_a}{{({\boldsymbol{x}} - {{\boldsymbol{x}}_a})}^2}} . $
形函数影响域越大, 计算量越大. LME形函数的基本思路是让形函数熵最大, 同时影响域又尽量小, 因此, 构造了如下优化问题:
$ \mathop {\min }\limits_{{N_a}} (f) = \beta U - H, \tag{4a}$
$ {N_a} \geqslant 0\quad a = 1,2, \cdots ,n,\tag{4b} $
$ \sum\limits_a {{N_a} = 1} , \tag{4c}$
$ \sum\limits_a {{N_a}} {{\boldsymbol{x}}_a} = {\boldsymbol{x}}, \tag{4d}$
其中, β是1个可调参数, 用于调节两个目标的权重. 最后1个约束条件表示形函数严格满足一阶精度. 上述优化问题解的存在性、唯一性以及求解方法文献[15-17]中有详细论述, 本文只给出最后结果. LME形函数可以写成如下形式:
$\begin{split}{N_a}({\boldsymbol{x}},{\lambda ^*}({\boldsymbol{x}})) =\;&\frac{1}{{Z({\boldsymbol{x}},{\lambda ^*}({\boldsymbol{x}}))}}\exp [ - \beta {({\boldsymbol{x}} - {{\boldsymbol{x}}_a})^2}\\&+ {{\boldsymbol{\lambda }}^*}{({\boldsymbol{x}})^{\text{T}}} \cdot ({\boldsymbol{x}} - {{\boldsymbol{x}}_a})],\end{split} $
其中Z函数的定义为
$ Z({\boldsymbol{x}},{\boldsymbol{\lambda }}) = \sum\limits_a {\exp \big[ { - \beta {{({\boldsymbol{x}} - {{\boldsymbol{x}}_a})}^2} + {{\boldsymbol{\lambda }}^{\text{T}}} \cdot ({\boldsymbol{x}} - {{\boldsymbol{x}}_a})} \big],} $
$ {{\boldsymbol{\lambda }}^*}({\boldsymbol{x}}) $是以下无约束优化问题的解:
$ {{\boldsymbol{\lambda }}^*}({\boldsymbol{x}}) = \arg \,\mathop {\min }\limits_{\boldsymbol{\lambda }} \;\lg Z({\boldsymbol{x}},{\boldsymbol{\lambda }}). $
定义:
$\begin{split} & {\boldsymbol{r}}({\boldsymbol{x}},{\boldsymbol{\lambda }}) \equiv {\partial _{\boldsymbol{\lambda }}}\lg Z({\boldsymbol{x}},{\boldsymbol{\lambda }}) \\=& \sum\limits_a {{N_a}({\boldsymbol{x}},{\boldsymbol{\lambda }})} ({\boldsymbol{x}} - {{\boldsymbol{x}}_a}), \end{split}\tag{8a}$
$\begin{split}\;& {\boldsymbol{J}}({\boldsymbol{x}},{\boldsymbol{\lambda }}) \equiv {\partial _{\boldsymbol{\lambda }}}{\partial _{\boldsymbol{\lambda }}}\lg Z({\boldsymbol{x}},{\boldsymbol{\lambda }})\\=\;& \sum\limits_a {{N_a}({\boldsymbol{x}},{\boldsymbol{\lambda }})} (x - {x_a}) \otimes ({\boldsymbol{x}} - {{\boldsymbol{x}}_a})\\&- {\boldsymbol{r}}({\boldsymbol{x}},{\boldsymbol{\lambda }}) \otimes {\boldsymbol{r}}({\boldsymbol{x}},{\boldsymbol{\lambda }}). \end{split}\tag{8b}$
那么, 可以采用如下迭代方法求解$ {{\boldsymbol{\lambda }}^*}({\boldsymbol{x}}) $:
$ {{\boldsymbol{\lambda }}^{(k + 1)}}({\boldsymbol{x}}) = {{\boldsymbol{\lambda }}^{(k)}}({\boldsymbol{x}}) - {\boldsymbol{J}}{(x,{{\boldsymbol{\lambda }}^{(k)}}({\boldsymbol{x}}))^{ - 1}}*{\boldsymbol{r}}({\boldsymbol{x}},{{\boldsymbol{\lambda }}^{(k)}}({\boldsymbol{x}})). $
初始时刻, $ {{\boldsymbol{\lambda }}^{(0)}} $一般取零. 流场发生变化后, $ {{\boldsymbol{\lambda }}^{(0)}} $一般取上1个时刻的$ {{\boldsymbol{\lambda }}^*} $. 当β等于零时, 得到的最大熵形函数是全局的, 不适合数值计算. 当β趋于无限大时, 即不考虑熵最大, 只要求形函数影响域最小, 那么LME形函数可以连续过渡到有限元线性形函数. LME形函数具有以下特点: 1)满足一阶精度; 2)不依赖网格; 3)边界处满足弱Kronecker delta性质; 4)形函数非负.
得到$ {{\boldsymbol{\lambda }}^*}({\boldsymbol{x}}) $与形函数后, 形函数导数为
$ \nabla {N_a}({\boldsymbol{x}}) = - {N_a}({\boldsymbol{x}}){\boldsymbol{J}}{({\boldsymbol{x}},{{\boldsymbol{\lambda }}^*}({\boldsymbol{x}}))^{ - 1}}({\boldsymbol{x}} - {{\boldsymbol{x}}_a}). $

2
2.3.基本离散方程
-->首先考虑无黏、可压缩以及等温流体. 当初始状态已知, 节点的位置描述了流场的形变的流动, 实际上完全描述了流场, 因此, 节点的位置可以当成系统的广义坐标. 本文下标p表示粒子, 用下标a, b, c表示节点. 由于无黏, 这是1个保守系统, 拉格朗日方程为
$ \frac{{\text{d}}}{{{\text{d}}t}}\left(\frac{{\partial E}}{{\partial {{\boldsymbol{v}}_a}}} \right) = - \frac{{\partial V}}{{\partial {{\boldsymbol{x}}_a}}} = {{\boldsymbol{f}}_a},\;a = 1,2, \cdots N, $
其中, E是系统动能, V是系统势能, $ {{\boldsymbol{f}}_a} $是广义节点力, $ {{\boldsymbol{v}}_a} $是节点速度, $ {{\boldsymbol{x}}_a} $是节点位置, N是节点总数. 系统动能和势能为所有粒子动能和内能之和:
$ E = \sum\limits_{p} {\frac{1}{2}{m_{p}}} {\boldsymbol{v}}_{p}^2, \tag{12a}$
$ V = \sum\limits_{p} {{m_{p}}{e_{p}}} ({\rho _{p}}),\tag{12b} $
其中$ {m_p} $表示粒子质量, $ {{\boldsymbol{v}}_p} $表示粒子速度, $ {e_p} $表示粒子单位质量内能, $ {\rho _p} $表示粒子密度.
粒子内能和压力的关系通过物态方程(equation of state, EOS)描述:
$ p({\rho _{p}}) = \rho _{p}^2\frac{{\partial e({\rho _{p}})}}{{\partial {\rho _{p}}}}. \tag{13}$
与有限元一样, 粒子位置和速度可以通过节点位置和速度插值得到:
$ {{\boldsymbol{x}}_{p}} = \sum\limits_a {{{\boldsymbol{x}}_a}{N_a}} ({{\boldsymbol{x}}_{p}}),\tag{14a} $
$ {\text{δ}} {{\boldsymbol{x}}_{p}} = \sum\limits_a {{\text{δ}} {{\boldsymbol{x}}_a}{N_a}} ({{\boldsymbol{x}}_{p}}), \tag{14b}$
$ {{\boldsymbol{v}}_{p}} = \sum\limits_a {{{\boldsymbol{v}}_a}{N_a}} ({{\boldsymbol{x}}_{p}}), \tag{14c}$
$ {\text{δ}} {{\boldsymbol{v}}_{p}} = \sum\limits_a {{\text{δ}} {{\boldsymbol{v}}_a}{N_a}} ({{\boldsymbol{x}}_{p}}), \tag{14d}$
$ \nabla \cdot {\text{δ}} {{\boldsymbol{x}}_{p}} = \sum\limits_a {{\text{δ}} {{\boldsymbol{x}}_a} \cdot \nabla {N_a}} ({{\boldsymbol{x}}_{p}}), \tag{14e}$
$ \nabla {\text{δ}} {{\boldsymbol{x}}_{p}} = \sum\limits_a {{\text{δ}} {{\boldsymbol{x}}_a} \otimes \nabla {N_a}} ({{\boldsymbol{x}}_{p}}),\tag{14f} $
其中$ {\text{δ}} {{\boldsymbol{x}}_{p}} $$ {\text{δ}} {{\boldsymbol{x}}_a} $分别为粒子和节点的虚位移, $ {N_a}({{\boldsymbol{x}}_{p}}) $表示节点a处形函数在粒子p处的值, $ \nabla {N_a}({{\boldsymbol{x}}_p}) $表示形函数梯度, 求和是对粒子的所有近邻节点. 粒子密度变化和位置变化的关系由连续性条件给出:
$ {\text{δ}} {\rho _p} = - {\rho _p}\nabla \cdot {\text{δ}} {{\boldsymbol{x}}_p}. $
将(14e)式代入(15)式, 得到
$ \frac{{\partial {\rho _p}}}{{\partial {{\boldsymbol{x}}_a}}} = - {\rho _p}\nabla {N_{a,p}}, $
利用求导链式规则, 并且利用(13)式和(16)式, 得到
$ \begin{split} {{\boldsymbol{f}}_a} =\;& - \frac{{\partial V}}{{\partial {{\boldsymbol{x}}_a}}} = - \sum\limits_p {{m_p}\frac{{{\text{d}}{e_p}}}{{{\text{d}}{\rho _p}}}\frac{{\partial {\rho _p}}}{{\partial {{\boldsymbol{x}}_a}}}} \hfill \\ = \;&\sum\limits_p^{} {\frac{{{m_p}}}{{{\rho _p}}}p({\rho _p})\nabla {N_{a,p}}} . \end{split} $
利用(12a)式、(14c)式和(14d)式得到
$ \frac{{\partial E}}{{\partial {{\boldsymbol{v}}_a}}} = \sum\limits_p {{m_p}{{\boldsymbol{v}}_p}\frac{{\partial {{\boldsymbol{v}}_p}}}{{\partial {{\boldsymbol{v}}_a}}}} = \sum\limits_p {{m_p}\left(\sum\limits_b {{{\boldsymbol{v}}_b}{N_{b,p}}} \right){N_{a,p}}} . $
在(18)式中交换求和顺序, 拉格朗日方程变成:
$ \sum\limits_b {{{\boldsymbol{M}}_{a,b}}} {{\boldsymbol{\dot v}}_b} = {{\boldsymbol{f}}_a}, $
其中, M为相容质量矩阵
$ {{\boldsymbol{M}}_{a,b}} = \sum\limits_p {{m_p}{N_{a,p}}{N_{b,p}}} . $
(19)式与更新拉格朗日有限元格式有两点区别: 1)有限元格式中的单元积分对应了OTM中的粒子积分; 2)有限元中使用基于网格的形函数, 而OTM方法中使用无网格的局部最大熵形函数. 注意到局部最大熵形函数可以连续过渡到有限元线性形函数, 因此, OTM方法可以看成是一种特殊的单点积分拉格朗日有限元方法.
在有限元中方法, 为避免矩阵求逆, 常常将相容质量矩阵按行求和, 得到集中质量矩阵, 大大减少了计算量和存储量, OTM方法中同样如此. 利用形函数的单位分解性质, 可以将粒子的质量分配给其相邻的节点, 得到如下节点质量定义:
$ {m_a} = \sum\limits_p {{m_p}} {N_{a,p}}. $
注意到粒子质量和局部最大熵形函数都是严格大于零的, 因此, 节点质量也是严格大于零的. 这样连续介质既可以看成是粒子系统, 也可以看成是节点系统. 很明显, 这两种系统的总质量和总动量都是严格相等的.
$ \sum\limits_a {{m_a}} = \sum\limits_p {{m_p}} \sum\limits_a {{N_{a,p}}} = \sum\limits_p {{m_p}} ,\tag{22a} $
$ \begin{split}\;&\sum\limits_a {{m_a}{{\boldsymbol{v}}_a}} = \sum\limits_a {\left(\sum\limits_p {{m_p}} {N_{a,p}}\right){{\boldsymbol{v}}_a}} \\=\;& \sum\limits_p {{m_p}\left(\sum\limits_a {{{\boldsymbol{v}}_a}{N_{a,p}}} \right)} = \sum\limits_p {{m_p}{{\boldsymbol{v}}_p}} .\\\end{split} \tag{22b}$
系统的总势能通过粒子系统计算. 系统的总动能既可以通过粒子系统计算, 也可以通过节点系统近似计算:
$ E = \sum\limits_p {\frac{1}{2}{m_p}{{\boldsymbol{v}}_p}^2} \approx \sum\limits_a {\frac{1}{2}{m_a}{{\boldsymbol{v}}_a}^2} . $
将(23)式代入(11)式, 得到如下离散方程:
$ {m_a}{{\boldsymbol{\dot v}}_a} = {{\boldsymbol{f}}_a}. $
从(20)式和(21)式可知, 节点质量刚好就是相容质量矩阵按行求和. 利用节点广义力和节点质量可以更新节点的速度和位置:
$ { \boldsymbol{v}}_a^{k + 1} = { \boldsymbol{v}}_a^k + {\boldsymbol{f}}_a^k/ m_a^k {\text{δ}} t, \tag{25a}$
$ {\boldsymbol{x}}_a^{k + 1} = {{\boldsymbol{x}}_a^k} + {{\boldsymbol{v}}_a^{k + 1}}{\text{δ}} t. \tag{25b}$
用(25b)式更新节点位置后, 粒子位置根据(14a)式更新. 注意到节点位置变化代表了流场的形变, 因此两个相邻时刻的连续映射$\varphi :{{\boldsymbol{x}}^k} \to {{\boldsymbol{x}}^{k + 1}}$可以通过更新的节点位置插值得到
$ {\varphi _{k \to k + 1}}({\boldsymbol{x}}) = \sum\limits_a {{{\boldsymbol{x}}_a^{k + 1}}{N_a}} ({\boldsymbol{x}}). $
上述映射的梯度即为变形梯度张量:
$\begin{split}&{{\boldsymbol{F}}_{k \to k + 1}}({\boldsymbol{x}}) = \nabla {\varphi _{k \to k + 1}}({\boldsymbol{x}}) \\=\;& \sum\limits_a {{{\boldsymbol{x}}_{a,k + 1}} \otimes \nabla {N_a}({\boldsymbol{x}})} .\end{split} $
变形梯度张量的行列式在粒子位置处的值表征了粒子的密度变化
$ \rho _p^{k + 1} = \rho _p ^k/\det ({F_{k \to k + 1}}({x_p})) . $
(28)式可以用于更新粒子密度. 粒子密度更新后, 根据物态方程更新粒子压力.
3
2.3.1.流体黏性的处理
-->以上推导没有考虑黏性, 因而是保守系统. 处理非保守力的基本想法将其看成外力, 然后利用虚功原理推导对应的广义力. 假设$ {{\boldsymbol{f}}_p} $是作用在粒子p上的外力, 那么此力的虚功为
$\begin{split}{\text{δ}} W =\;& \sum\limits_p {{{\boldsymbol{f}}_p} \cdot {\text{δ}} {{\boldsymbol{x}}_p} = } \sum\limits_p {{{\boldsymbol{f}}_p} \cdot \left( {\sum\limits_a {{\text{δ}} {{\boldsymbol{x}}_a}{N_{a,p}}} } \right)} \\=\;& \sum\limits_a {\left( {\sum\limits_p {{{\boldsymbol{f}}_p}{N_{a,p}}} } \right)} \cdot {\text{δ}} {{\boldsymbol{x}}_a}. \end{split}$
由虚位移的独立性, 得到相应的广义力为
$ {{\boldsymbol{f}}_a} = \sum\limits_p {{{\boldsymbol{f}}_p}{N_{a,p}}} . $
例如, 作用在1个粒子上的重力$ {{\boldsymbol{f}}_p} = {m_p}{\boldsymbol{g}} $. 根据(30)式得到广义力${{\boldsymbol{f}}_a} = \displaystyle \sum\nolimits_p {{m_p}{\boldsymbol{g}}{N_{a, p}}}$. 考虑到重力实际上是保守力, 因此也可以将重力势能加入拉格朗日函数中, 得到重力对应的广义力. 简单计算可以证明两者是一致的.
材料的很多行为(如弹塑性, 黏弹性等)都可以统一用1个应力张量描述, 而应力张量的影响也可以看成外力, 如下式中的动量守恒方程:
$ \frac{{D{\boldsymbol{v}}}}{{Dt}} = \frac{1}{\rho }\nabla \cdot {\boldsymbol{\sigma }}. $
应力张量散度的物理意义为单位体积受到的力, 因此, 应力张量虚功为
$ \begin{split} {\text{δ}} W =\;& \int {\nabla \cdot {\boldsymbol{\sigma }} \cdot } {\text{δ}} {\boldsymbol{x}}{\text{d}}V = \int {[\nabla \cdot ({\boldsymbol{\sigma }} \cdot } {\text{δ}} {\boldsymbol{x}}) - {\boldsymbol{\sigma }}:\nabla {\text{δ}} {\boldsymbol{x}}]{\text{d}}V \hfill \\ = \;&- \int {{\boldsymbol{\sigma }}:\nabla } {\text{δ}} {\boldsymbol{x}}{\text{d}}V \approx - \sum\limits_p {{{\boldsymbol{\sigma }}_p}:} \nabla {\text{δ}} {{\boldsymbol{x}}_p} \dfrac{m_p}{\rho _p} \hfill \\ = \;&- \sum\limits_p {{{\boldsymbol{\sigma }}_p}:} \left(\sum\limits_a {{\text{δ}} {{\boldsymbol{x}}_a} \otimes \nabla {N_{a,p}}} \right) \dfrac{m_p}{\rho _p} \hfill \\ = \;&- \sum\limits_a {\left(\sum\limits_p {{{\boldsymbol{\sigma }}_p} \cdot \nabla {N_{a,p}}\dfrac{m_p}{\rho _p}} \right) \cdot {\text{δ}} {{\boldsymbol{x}}_a}} , \\[-15pt]\end{split} $
其中用到了高斯积分公式, 并且假设了边界上虚位移或者应力张量为零, 这在固壁边界以及自由表面边界上都成立. 同理, 得到广义力为
$ {{\boldsymbol{f}}_a} = - \sum\limits_p {\frac{{{m_p}}}{{{\rho _p}}}{{\boldsymbol{\sigma }}_p} \cdot \nabla {N_{a,p}}} . $
如果只考虑牛顿流体的黏性, 那么黏性应力张量为
$ {\boldsymbol{\sigma }} = 2\mu (\nabla {\boldsymbol{v}} + \nabla {{\boldsymbol{v}}^{\text{T}}}), $
其中μ为黏性系数.
同理, (31)式也能通过势能和拉格朗日方程得到. 从(33)式可以看到, 1个粒子在无限短时间内受到的力只和应力张量本身有关, 而和应力应变关系无关. 因此, 可以临时假设1个线性应力应变关系, 从而得到简单的弹性势能. 对弹性势能求变分, 就得到了相应的广义力, 详细证明见附录A. 两种方法的结果是一致的.
3
2.3.2.表面张力处理
-->考虑二维情况, 液体边界被边界节点表征, 如图2所示.
图 2 二维条件下的边界节点(白点)
Figure2. Boundary nodes (white symbols) in two-dimensional coordinate.

表面势能正比于表面积:
$ V = \gamma ({r_{ab}} + {r_{ac}} + \cdots ) \text{, } $
其中$ {r_{ab}} = \left| {{{\boldsymbol{x}}_a} - {{\boldsymbol{x}}_b}} \right| $, $ {r_{ac}} = \left| {{{\boldsymbol{x}}_a} - {{\boldsymbol{x}}_c}} \right| $为节点之间的间距. γ是表面张力系数. 对(35)式求变分, 得到广义力
$ {{\boldsymbol{f}}_a} = - \frac{{\partial V}}{{\partial {{\boldsymbol{x}}_a}}} = - \gamma \left(\frac{{{{\boldsymbol{x}}_a} - {{\boldsymbol{x}}_b}}}{{{r_{ab}}}} + \frac{{{{\boldsymbol{x}}_a} - {{\boldsymbol{x}}_c}}}{{{r_{ac}}}}\right). $
可以看到, (36)式非常简洁, 可以解释为节点a同时受到两个相邻节点bc的吸引力, 吸引力大小刚好等于表面张力系数. (36)式中没有可调参数, 表面张力系数是唯一的输入参数. 此外, 表面张力对应的广义力只作用在表面节点上, 即精确的作用在流体表面.
3
2.3.3.轴对称处理
-->记(r, θ, h)为柱坐标系, 并且假设环向没有位移, 而且所有物理量都和环向坐标无关, 即$ {\text{δ}} {u_\theta } = 0 $, $ \partial /\partial \theta = 0 $. 定义:
$ {{\boldsymbol{x}}_a} = ({r_a},{h_a}), \tag{37a}$
$ {{\boldsymbol{v}}_a} = ({\dot r_a},{\dot h_a}), \tag{37b}$
$ \nabla f = (\partial f/\partial r,\partial f/\partial h) \text{, }\tag{37c} $
$ \nabla \cdot {\boldsymbol{f}} = \partial {f_r}/\partial r + \partial {f_h}/\partial h.\tag{37d} $
那么, (11)—(14)式依然成立. 速度矢量散度在柱坐标下的表达式为
$ \nabla \cdot ({\boldsymbol{v}}) = \frac{{\partial {v_r}}}{{\partial r}} + \frac{{{v_r}}}{r} + \frac{{\partial {v_h}}}{{\partial h}} = \nabla \cdot {\boldsymbol{v}} + \frac{{{v_r}}}{r}, $
因此, 粒子密度变化和粒子位置变化关系式(15)需要修改为
$ \begin{split} {\text{δ}} {\rho _p} = \;&- {\rho _p}\left(\sum\limits_a {{\text{δ}} {{\boldsymbol{x}}_a}} \cdot \nabla {N_{a,p}} + \frac{{{\text{δ}} {r_p}}}{{{r_p}}}\right) \hfill \\ =\; &- {\rho _p}\sum\limits_a {\left({\text{δ}} {{\boldsymbol{x}}_a} \cdot \nabla {N_{a,p}} + \frac{{{N_{a,p}}{\text{δ}} {r_a}}}{{{r_p}}}\right)} , \end{split} $
于是
$ \frac{{\partial {\rho _p}}}{{\partial {{\boldsymbol{x}}_a}}} = - {\rho _p}\left(\nabla {N_{a,p}} + \frac{{{N_{a,p}}}}{{{r_p}}}{{\boldsymbol{n}}_r}\right), $
其中$ {{\boldsymbol{n}}_r} $是径向单位向量. 于是, 轴对称条件下的广义力为
${{\boldsymbol{f}}_a} = \sum\limits_p {\frac{{{m_p}}}{{{\rho _p}}}p({\rho _p})\left( {\nabla {N_{a,p}} + \frac{{{N_{a,p}}}}{{{r_p}}}{{\boldsymbol{n}}_r}} \right)} , $
(18)—(26)式仍然成立. 注意(27)式中的变形梯度张量的行列式表征的是粒子在(r, h)平面上的面积的变化:
$ {S^{k + 1}} = {S^k}\det \left({{\boldsymbol{F}}_{k \to k + 1}}({{\boldsymbol{x}}_p})\right), $
因此, 可以先更新粒子在(r, h)平面上的面积$ {S^{k + 1}} $, 然后再更新密度:
$ {\rho ^{k + 1}} = \frac{{{m_p}}}{{2{\text{π }}{r_p}^{k + 1}{S^{k + 1}}}}. $
很明显, 外力对应的广义力(30)式仍然成立. 黏性力对应的广义力也仍然成立, 详细证明见附录B.
三维轴对称条件下, 表面势能为
$ V = {\text{π }}\gamma ({\tilde r_{ab}}{r_{ab}} + {\tilde r_{ac}}{r_{ac}} + \cdots ), $
其中, $ {\tilde r_{ab}} = {r_a} + {r_b} $, $ {\tilde r_{ac}} = {r_a} + {r_c} $. $ {r_k}(k = a, b, c) $为节点的径向坐标. 对(44)式求变分, 得到表面张力对应的广义力:
$\begin{split}\;&{{\boldsymbol{f}}_a} = - \frac{{\partial V}}{{\partial {{\boldsymbol{x}}_a}}}= - {\text{π }}\gamma \Bigg({\tilde r_{ab}}\frac{{{{\boldsymbol{x}}_a} - {{\boldsymbol{x}}_b}}}{{{r_{ab}}}} \\& \qquad + {\tilde r_{ac}}\frac{{{{\boldsymbol{x}}_a} - {{\boldsymbol{x}}_c}}}{{{r_{ac}}}} + {r_{ab}}{{\boldsymbol{n}}_r} + {r_{ac}}{{\boldsymbol{n}}_r}\Bigg), \\[-15pt]\end{split}$
式中, 前两项为节点b和节点c对节点a的吸引力, 后两项为对称轴对边界节点的吸引力. 可见, 笛卡尔坐标系下的OTM计算格式只需要经过微小改动, 就能计算三维轴对称问题.
2
2.4.OTM方法详细步骤
-->笛卡尔坐标系下, 采用OTM方法模拟表面张力的详细步骤总结如下.
1)初始化: 给定节点的位置和速度, 给定粒子的位置和其他物理量.
2)近邻粒子搜索: 找到所有粒子-节点对.
3)对于所有粒子-节点对, 计算形函数$ {N_{a, p}} $和形函数导数$ \nabla {N_{a, p}} $.
4)计算节点质量、粒子速度梯度、黏性应力张量和广义节点力:
$ {m_a} = \sum\limits_b {{{\boldsymbol{M}}_{a,b}}} = \sum\limits_p {{m_p}} {N_{a,p}}, $
$ \nabla {{\boldsymbol{v}}_p} = \sum\limits_a {{{\boldsymbol{v}}_a} \otimes \nabla {N_{a,p}}} , $
$ {{\boldsymbol{\sigma }}_p} = 2\mu (\nabla {{\boldsymbol{v}}_p} + \nabla {\boldsymbol{v}}_p^{\text{T}}), $
$ \begin{split} {{\boldsymbol{f}}_a} = \;&\sum\limits_p {\frac{{{m_p}}}{{{\rho _p}}}p({\rho _p})\nabla {N_{a,p}}} - \sum\limits_p {\frac{{{m_p}}}{{{\rho _p}}}{{\boldsymbol{\sigma }}_p} \cdot \nabla {N_{a,p}}} \\ &+ \sum\limits_p {{m_p}{\boldsymbol{g}}{N_{a,p}}} - \gamma \left(\frac{{{{\boldsymbol{x}}_a} - {{\boldsymbol{x}}_b}}}{{{r_{ab}}}} + \frac{{{{\boldsymbol{x}}_a} - {{\boldsymbol{x}}_c}}}{{{r_{ac}}}}\right). \end{split} $
5)更新节点速度和位置:
$ \begin{split}& {\boldsymbol{v}}_a ^{k + 1} = \boldsymbol{v}_a^k + {\boldsymbol{f}}_a ^k/ m_a^k{\text{δ}} t, \\ & {\boldsymbol{x}}_a^{k + 1} = {\boldsymbol{x}}_a^k + {\boldsymbol{v}}_a^{k + 1}{\text{δ}} t. \end{split}$
6)根据边界条件调整节点位置, 然后更新粒子位置:
$ {\boldsymbol{x}}_p ^{k+1} = \sum\limits_a {{\boldsymbol{x}}_a^{k + 1}{N_{a,p}}} . $
7)计算变形梯度张量并更新粒子密度和压力:
$ {{\boldsymbol{F}}_{k \to k + 1}}({\boldsymbol{x}}) = \sum\limits_a {{{\boldsymbol{x}}_{a,k + 1}} \otimes \nabla {N_a}({\boldsymbol{x}})} , $
$ \rho _p^{k + 1} = \rho_p^k/\det ({{\boldsymbol{F}}_{k \to k + 1}}({{\boldsymbol{x}}_p})), $
$ p_p^{k + 1} = p(\rho_p^{k + 1}) . $
8)完成1个时间步长, 回到步骤2).
所有的初始化借助于有限元网格. 网格节点变成节点, 网格形心放置粒子, 粒子体积即为网格大小. 物态方程如下:
$ p = {\rho _0}{c_0^2}\big[{{{\left( { {\rho }/{{{\rho _0}}}} \right)}^7} - 1} \big], $
其中, $ {\rho _0} $表示材料初始密度, ρ表示当前密度, $ {c_0} $表示声速, p表示压力.
2
3.1.泊肃叶流
-->两块无限长的二维平板中间充满了静止液体, 液体在平行于平板方向的体积力作用下开始运动, 称之为泊肃叶(Poiseuille)流. 如果是无限长的三维圆管, 那么称之为Hagen-Poiseuille流. 对于二维泊肃叶, 参数和文献[24]一样, 解析解见文献[24]. 将平板距离变成圆管半径, 即为Hagen-Poiseuille流, 解析解见文献[25]. OTM模拟结果见图3.
图 3 速度分布的OTM模拟与解析解对比 (a) 二维泊肃叶流; (b) 三维轴对称泊肃叶流
Figure3. Comparison of OTM (symbols) and analytic solutions (solid curves) for velocity profile: (a) Two-dimensional Poiseuilleflow; (b) axisymmetric Hagen-Poiseuille flow.

OTM模拟结果与理论解吻合很好. 为了模拟无限长管道, 用到了周期边界条件.
2
3.2.静止液滴
-->二维液滴的半径R = 0.2 m, 密度ρ = 1 kg/m3, 表面张力系数γ = 1 Pa·m, 黏性系数η = 0.5 Pa·s, 声速c0 = 50 m/s. 静止状态下液滴的理论压力值可以根据Young-Laplace 公式得到ρ = γ/R = 5 Pa. 如果是三维液滴, 那么其理论压力值为ρ = 2γ/R = 10 Pa. OTM模拟的平均压力分别为 5.004 Pa和 9.967 Pa. 压力分布如图4所示.
图 4 静止液滴的压力场 (a) 二维液滴; (b) 轴对称三维液滴
Figure4. Pressure field: (a) Two-dimensional rod; (b) axisymmetric three dimensional drop.

由于弱可压假设, 压力场中有误差. 但是液滴中的平均压力非常接近理论值(误差0.5%以内), 证明了本方法的精度.
2
3.3.液滴的周期振荡
-->在上1个算例的基础上, 将黏性减小到$ \eta = $$ 5 \times {10^{ - 3}} $ Pa·s, 并且附加1个散度为0的初始速度场$ {v_x} = x $, $ {v_y} = - y $给所有节点. 二维液滴的理论振荡周期为$ T = 2{\text{π }}\sqrt {\rho {R^3}/6\gamma } $ [26]. 为了验证OTM方法的收敛性, 采用不同的粒子总数模拟这个问题, 并且检查两条相互垂直线$ x = 0 $$ y = 0 $上的物理量分布, 时刻$ t = T/2 $, 如图5所示.
图 5 t = T/2时刻的压力和速度分布 (a) x = 0上的压力分布; (b) x = 0上的速度分布; (c) y = 0上的压力分布; (d) y = 0上的速度分布
Figure5. Pressure and velocity profile at t = T/2: (a) Pressure profile at x = 0; (b) velocity profile at x = 0; (c) pressure profile at y = 0; (d) velocity profile at y = 0.

可以看到, 压力分布的收敛不明显, 但是速度分布的收敛很明显. 改变表面张力系数的大小, 得到的周期与理论值的对比如图6所示. 周期是根据液滴右上部分质心轨迹测量得到, 与文献[26]一样.
图 6 二维液滴振荡 (a) 振荡周期理论解和数值解得对比; (b) 表面张力系数为$ \gamma = 1 $时液滴右上部分质心的轨迹
Figure6. Two-dimensional rod oscillation: (a) Comparison of period between the numerical (symbols) and analytical (solid curve) results; (b) center of mass position history when $ \gamma = 1 $.

对于轴对称的三维液滴, 散度为零的初始速度场取为$ {v_{\text{r}}} = r $, $ {v_h} = - 2 h $, 其他物理量不变. 振荡周期的理论解为$ T = 2\pi \sqrt {\rho {R^3}/8\gamma } $[27]. 模拟结果与理论的对比如图7所示. 模拟的周期根据对称轴端点的位置轨迹测量. 可以看到, 两种情况下OTM计算结果和理论非常吻合.
图 7 轴对称三维液滴振荡周期与理论的对比
Figure7. Three-dimensional droplet oscillation, comparing of period between the numerical (symbols) and analytical (solid curve) results.

2
3.4.液滴形变
-->将液滴半径扩大到 1 m, 表面张力系数扩大到10 Pa·m, 加上大小为10 m/s2的重力. y = –1 m处放置1个光滑的固定平板. 液滴在重力、表面张力以及光滑平板的共同作用下产生变形, 如图8所示.
图 8 压力场 (a) t = 0 s; (b) t = 0.02 s; (c) t = 0.4 s; (d) t = 4 s
Figure8. Pressure field at several typical times: (a) t = 0 s; (b) t = 0.02 s; (c) t = 0.4 s; (d) t = 4 s.

图8(d)中, 压力等值线与重力垂直, 和理论相符. 液滴内部的压力场需要满足两个条件, 一是压力沿Y轴的线性分布规律$ {\text{δ}} p = \rho g{\text{δ}} Y $, 另一个是自由表面附近的压力满足Young-Laplace方程$ p = \gamma /R $. 这两个条件实际上决定了液滴的形状, 即液滴自由表面部分的微分方程为
$ R = \dfrac{{{{(1 + y{'^2})}^{3/2}}}}{{\left| {y'} \right|}} = \dfrac{\gamma }{{{p_0} + \rho g(H - y)}},$
其中Hp0分别为液滴最高点的Y坐标和压力. 上述理论曲率半径在实际计算时误差太大, 因此, 本文采用如下方法计算曲率半径: 对于每个边界节点, 找到与其相邻的4个表面节点, 然后用1个圆来拟合这5个节点, 即$\mathop {\min }\limits_{{x_0}, r} \displaystyle \sum\nolimits_{i = 1}^5 {{{(\left| {{x_i} - {x_0}} \right| - r)}^2}}$, 自由变量为圆的中心坐标和半径. 圆的半径即为此节点的曲率半径.
液体中自由边界的右边部分如图9(a)所示. 采用曲率半径计算的边界压力和根据EOS计算的压力比较如图9(b)所示.
图 9 边界位置和压力 (a) 边界位置; (b) 边界压力与Y轴坐标的关系, 压力根据表面曲率和EOS计算
Figure9. Boundary position and pressure profile: (a) Boundary position; (b) boundary pressure versus Y coordinate, computed by Young-Laplace equation (circles) and EOS (points).

两种方法得到压力值大致相等, 说明模拟的自由表面形状和理论相符. 靠近固壁处误差较大, 经检查是因为固壁附近节点的分布更加无序, 导致搜索算法对初始参数敏感. 压力在Y轴方向的线性分布和理论相符, 拟合直线的斜率也和理论$ {\text{d}}p/{\text{d}}Y = - \rho g = - 10 $ Pa/m相符.
本文基于拉格朗日方程, 通过将表面张力势能加入拉格朗日函数, 得到了OTM方法框架下的表面张力效应表达式, 并且给出了轴对称的处理方式. 模拟了泊肃叶流, 静止和振荡的液滴, 液滴在重力、表面张力以及光滑平板共同作用下的形变等典型算例, 通过与解析解对比, 验证了OTM方法在模拟表面张力现象中具有如下优势: 1)表面边界节点精确的表征了边界的位置, 保证了表面张力势能的准确计算, 而且表面张力对应的广义力精确的作用在流体表面; 2)表面张力对应的广义力形式简洁, 具有直观的物理意义, 而且表面张力系数是唯一的输入参数, 不需要任何可调参数; 3)可以保证速度分布的空间收敛性.
本文在连续介质离散时, 节点和节点、节点和粒子之间没有任何固定联系, 而是通过近邻粒子搜索动态确定. 但是, 在计算表面张力势能时, 为了方便计算表面积, 假设了表面节点之间的连接关系不变. 因此, 所有表面节点可以看成是1个表面有限元网格. 表面有限元网格与内部节点在拉格朗日方程框架下相互作用和协调运动. 下一步, 将考虑表面拓扑结构发生变化的情况, 如多个液滴的融合过程.
下面给出(33)式的势能推导方法. 将应力和应变张量写成向量形式${:} $
$ {{\boldsymbol{\varepsilon}} ^{\text{T}}} = \{ {\varepsilon _{xx}},{\varepsilon _{yy}},{\varepsilon _{zz}},{\varepsilon _{yz}},{\varepsilon _{xz}},{\varepsilon _{xy}}\} , \tag{A1}$
$ {{\boldsymbol{\sigma}} ^{\text{T}}} = \{ {\sigma _{xx}},{\sigma _{yy}},{\sigma _{zz}},{\sigma _{yz}},{\sigma _{xz}},{\sigma _{xy}}\} .\tag{A2} $
应变与位移的关系为
$ {\boldsymbol{\varepsilon}} = {\boldsymbol{Lu}}, \tag{A3}$
其中L是微分矩阵
$ {\boldsymbol{L}} = \left[ {\begin{array}{*{20}{c}} {\partial /\partial x}&0&0 \\ 0&{\partial /\partial y}&0 \\ 0&0&{\partial /\partial z} \\ 0&{\partial /\partial z}&{\partial /\partial y} \\ {\partial /\partial z}&0&{\partial /\partial x} \\ {\partial /\partial y}&{\partial /\partial x}&0 \end{array}} \right], \tag{A4}$
u是位移. 广义胡克定律为
$ {\boldsymbol{\sigma}} = {\boldsymbol{c\varepsilon}}, \tag{A5}$
其中c是1个对称矩阵:
$ {\boldsymbol{c}} = {{\boldsymbol{c}}^{\text{T}}}. \tag{A6}$
以下4个等式成立:
$ {{\boldsymbol{\varepsilon }}_p} = {\boldsymbol{L}}{u_p} = \sum\limits_a {{\boldsymbol{L}}{N_{a,p}}} {u_a},\quad \frac{{\partial {{\boldsymbol{\varepsilon }}_p}}}{{\partial {u_a}}} = {\boldsymbol{L}}{N_{a,p}},\tag{A7} $
$ {{\boldsymbol{\sigma }}_p} = {\boldsymbol{cL}}{u_p} = \sum\limits_a {{\boldsymbol{cL}}{N_{a,p}}} {u_a},\quad \frac{{\partial {{\boldsymbol{\sigma }}_p}}}{{\partial {u_a}}} = {\boldsymbol{cL}}{N_{a,p}}. \tag{A8}$
将(16)式写成矩阵形式:
$ {\left(\frac{{\partial {\rho _p}}}{{\partial {x_a}}}\right)^{\text{T}}} = - {\rho _p}\nabla {N_{a,p}}, \tag{A9}$
其中约定了形函数导数$ \nabla {N_{a, p}} $为列向量. 对于固定的应力张量$ {{\boldsymbol{\sigma}}_p} $, 由于c可以取得足够大, 使得应变为小变形, 于是弹性势能为
$ V = \sum\limits_p {\frac{{{m_p}}}{{2{\rho _p}}}{\boldsymbol{\sigma }}_p^{\text{T}}{{\boldsymbol{\varepsilon }}_p}} = \sum\limits_p {\frac{{{m_p}}}{{2{\rho _p}}}{\boldsymbol{\varepsilon }}_p^{\text{T}}{{\boldsymbol{\sigma }}_p}} . \tag{A10}$
对弹性势能求变分, 得到
$ \begin{split}{f_a} =\;& - {\left( {\frac{{\partial V}}{{\partial {u_a}}}} \right)^{\rm{T}}}\\= \;& - {\sum\limits_p {\frac{{{m_p}}}{{2{\rho _p}}}\left( {{\boldsymbol{\sigma }}_p^{\rm{T}}\frac{{\partial {{\boldsymbol{\varepsilon }}_p}}}{{\partial {u_a}}} + {\boldsymbol{\varepsilon }}_p^{\rm{T}}\frac{{\partial {{\boldsymbol{\sigma }}_p}}}{{\partial {u_a}}} - {\boldsymbol{\sigma }}_p^{\rm{T}}{{\boldsymbol{\varepsilon }}_p}\frac{1}{{{\rho _p}}}\frac{{\partial {\rho _p}}}{{\partial {u_a}}}} \right)} ^{\rm{T}}}\\ =\;& - \sum\limits_p {\frac{{{m_p}}}{{2{\rho _p}}}({{\boldsymbol{L}}^{\rm{T}}}{N_{a,p}}{{\boldsymbol{\sigma }}_p} + } {({\boldsymbol{cL}}{N_{a,p}})^{\rm{T}}}{{\boldsymbol{\varepsilon }}_p} + {\boldsymbol{\sigma }}_p^{\rm{T}}{{\boldsymbol{\varepsilon }}_p}\nabla {N_{a,p}})\\ = \;&- \sum\limits_p {\frac{{{m_p}}}{{2{\rho _p}}}({{\boldsymbol{L}}^{\rm{T}}}{N_{a,p}}{{\boldsymbol{\sigma }}_p} + {{\boldsymbol{L}}^{^{\rm{T}}}}{N_{a,p}}{{\boldsymbol{c}}^{\rm{T}}}{{\boldsymbol{\varepsilon }}_p}} + {\boldsymbol{\sigma }}_p^{\rm{T}}{{\boldsymbol{\varepsilon }}_p}\nabla {N_{a,p}})\\ = \;& - \sum\limits_p {\frac{{{m_p}}}{{{\rho _p}}}{{\boldsymbol{L}}^{\rm{T}}}{N_{a,p}}} {{\boldsymbol{\sigma }}_p} - \sum\limits_p {\frac{{{m_p}}}{{2{\rho _p}}}} {\boldsymbol{\sigma }}_p^{\rm{T}}{{\boldsymbol{\varepsilon }}_p}\nabla {N_{a,p}}.\end{split} \tag{A11}$
将向量形式改写为常用的张量形式, 得到
$ {f_a} = - \sum\limits_p {\frac{{{m_p}}}{{{\rho _p}}}{{\boldsymbol{\sigma }}_p} \cdot \nabla {N_{a,p}}} - \sum\limits_p {\frac{{{m_p}}}{{2{\rho _p}}}} ({{\boldsymbol{\sigma }}_p}:{{\boldsymbol{\varepsilon }}_p})\nabla {N_{a,p}} . \tag{A12}$
(A12)式第二项中有应变, 因此是小量, 可以忽略
$ {f_a} = - \sum\limits_p {\frac{{{m_p}}}{{{\rho _p}}}{{\boldsymbol{\sigma }}_p} \cdot \nabla {N_{a,p}}} .\tag{A13} $
(A13)式和 (33)式一致.
两个张量的双点乘是1个标量, 因此, 可以在柱坐标下求解
$\begin{split}\;&{\boldsymbol{\sigma }}:\nabla {\text{δ}} x = \\\;&\left( {\begin{array}{*{20}{c}} {{\sigma _{rr}}}&{{\sigma _{r\theta }}}&{{\sigma _{rh}}} \\ {{\sigma _{\theta r}}}&{{\sigma _{\theta \theta }}}&{{\sigma _{\theta h}}} \\ {{\sigma _{hr}}}&{{\sigma _{h\theta }}}&{{\sigma _{hh}}} \end{array}} \right) : \left( {\begin{array}{*{20}{c}} {{\text{δ}} {x_{r,r}}}&{{\text{δ}} {x_{\theta ,r}}}&{{\text{δ}} {x_{h,r}}} \\ {{\text{δ}} {x_{r,\theta }}}&{{\text{δ}} {x_{\theta ,\theta }}}&{{\text{δ}} {x_{h,\theta }}} \\ {{\text{δ}} {x_{r,h}}}&{{\text{δ}} {x_{\theta ,h}}}&{{\text{δ}} {x_{h,h}}} \end{array}} \right).\\ \end{split} \tag{B1}$
考虑到轴对称, 得到
$ \begin{split} \sigma :\nabla {\text{δ}} x =\;& \left( {\begin{array}{*{20}{c}} {{{\boldsymbol{\sigma}}_{rr}}}&{{\sigma _{r\theta }}}&{{\sigma _{rh}}} \\ {{\sigma _{\theta r}}}&{{\sigma _{\theta \theta }}}&{{\sigma _{\theta h}}} \\ {{\sigma _{hr}}}&{{\sigma _{h\theta }}}&{{\sigma _{hh}}} \end{array}} \right):\left( {\begin{array}{*{20}{c}} {{\text{δ}} {x_{r,r}}}&0&{{\text{δ}} {x_{h,r}}} \\ 0&0&0 \\ {{\text{δ}} {x_{r,h}}}&0&{{\text{δ}} {x_{h,h}}} \end{array}} \right) \hfill \\ =\;& \left( {\begin{array}{*{20}{c}} {{\sigma _{rr}}}&{{\sigma _{rh}}} \\ {{\sigma _{hr}}}&{{\sigma _{hh}}} \end{array}} \right):\left( {\begin{array}{*{20}{c}} {{\text{δ}} {x_{r,r}}}&{{\text{δ}} {x_{h,r}}} \\ {{\text{δ}} {x_{r,h}}}&{{\text{δ}} {x_{h,h}}} \end{array}} \right) \end{split} . \tag{B2}$
利用定义(37a)式和(37c)式, 虚功为
$ \begin{split}{\text{δ}} W =\;& - \int {\boldsymbol{\sigma }(r,h):\nabla } {\text{δ}} x(r,h)r{\rm{d}}r{\rm{d}}\theta {\rm{d}}h\\ \approx \;& - \sum\limits_p {{\boldsymbol{\sigma }_p}:} \nabla {\text{δ}} {x_p}{m_p}/{\rho _p}.\end{split}\tag{B3}$
(B3)式与(32)式相同.
相关话题/计算 系统 质量 文献 节点

  • 领限时大额优惠券,享本站正版考研考试资料!
    大额优惠券
    优惠券领取后72小时内有效,10万种最新考研考试考证类电子打印资料任你选。涵盖全国500余所院校考研专业课、200多种职业资格考试、1100多种经典教材,产品类型包含电子书、题库、全套资料以及视频,无论您是考研复习、考证刷题,还是考前冲刺等,不同类型的产品可满足您学习上的不同需求。 ...
    本站小编 Free壹佰分学习网 2022-09-19
  • 一种基于摄动理论的不连续系统Lyapunov指数算法
    摘要:Lyapunov指数是识别系统非线性动力学特征的重要标志,但是目前的算法通用性不足且计算流程复杂.本文在经典的Lyapunov指数算法的基础上,基于摄动理论提出了一种适用于不连续系统的Lyapunov指数计算方法.首先,以系统状态参数初始值和沿相空间每个基本矢量的扰动量为初始条件,确定相轨迹. ...
    本站小编 Free考研考试 2021-12-29
  • 时间尺度上非迁移Birkhoff系统的Mei对称性定理
    摘要:研究并证明时间尺度上非迁移Birkhoff系统的Mei对称性定理.首先,建立任意时间尺度上Pfaff-Birkhoff原理和广义Pfaff-Birkhoff原理,由此导出时间尺度上非迁移Birkhoff系统(包括自由Birkhoff系统、广义Birkhoff系统和约束Birkhoff系统)的动 ...
    本站小编 Free考研考试 2021-12-29
  • 周期驱动系统的非平衡热输运与热力学几何
    摘要:随着对微纳尺度系统的深入理解和实验技术的进步,发生在这些小系统中的热输运和能量转换近期吸引了大量研究.不同于依赖静态热力学力(如温差、电势差等)的非平衡稳态调控手段,受时间驱动的非平衡非稳态小系统具有特有的高可调性和普遍性,其研究同时具有基础价值和应用潜力.本文从几何这一基本概念出发,分析了热 ...
    本站小编 Free考研考试 2021-12-29
  • 新的具有宽参数范围的五维保守超混沌系统的动力学研究
    摘要:保守系统因为没有吸引子,与常见的耗散系统相比,它的遍历性更好,伪随机性更强,安全性更高,更适合应用于混沌保密通信等领域.基于此,设计了一个新的具有宽参数范围的五维保守超混沌系统.首先,进行Hamilton能量和Casimir能量分析,证明了新系统满足Hamilton能量保守且能够产生混沌.然后 ...
    本站小编 Free考研考试 2021-12-29
  • 基本费米子质量和代问题
    摘要:研究了基本费米子的质量分布,并找到一组描述基本费米子质量在特定分布模式下的经验关系式.这启发我们对基本费米子质量等级和基本费米子具有三代的根源进行深入的思考,提出了一种理论模型,解释了基本费米子为什么具有三代,并讨论了基本费米子质量等级和自旋的起源.关键词:质量等级/基本费米子/代/自旋Eng ...
    本站小编 Free考研考试 2021-12-29
  • 电子束对ZnO和TiO<sub>2</sub>辐照损伤的模拟计算
    摘要:电子辐照在材料中产生的缺陷主要是相互独立的空位-间隙原子对,由于不同靶原子的离位阈能不同,通过改变电子束的能量可以调控在材料中产生的缺陷类型,同时,电子的注量又可以决定电子辐照产生的缺陷的浓度.ZnO和TiO2的磁光电特性受Zn空位、Ti空位、O空位、Zn间隙原子、Ti间隙原子等缺陷的影响,因 ...
    本站小编 Free考研考试 2021-12-29
  • 球差对高功率激光上行大气传输光束质量的影响
    摘要:地基激光空间碎片清除等激光烧蚀推进在太空中的应用中,激光功率已远超过大气非线性自聚焦临界功率,因此自聚焦效应是影响光束质量的重要因素.此外,由于高功率激光产生过程中的非线性效应,光束常伴有球差.本文采用数值模拟方法,研究了球差对高功率激光上行大气传输光束质量的影响.研究表明:对于大尺寸(光束发 ...
    本站小编 Free考研考试 2021-12-29
  • 高阶耦合相振子系统的同步动力学
    摘要:由大量耦合相振子组成的Kuramoto模型是研究各种自持续振荡系统同步相变和集体动力学的重要模型.近些年,高阶耦合Kuramoto模型引起了广泛的研究兴趣,尤其高阶耦合结构在模拟编码和信息存储的动力学方面起到重要作用.为了研究高阶耦合的影响,本文通过考虑频率与耦合之间的关联对高阶耦合的Kura ...
    本站小编 Free考研考试 2021-12-29
  • 高阶效应下对称三量子点系统中光孤子稳定性研究
    摘要:利用多重尺度法解析地研究了窄脉冲探测光激发下半导体三量子点分子系统中高阶效应对光孤子稳定性的影响.结果表明,由标准非线性薛定谔方程所描述的光孤子在传播的过程中会出现较大衰减,而由高阶非线性薛定谔方程所描述的光孤子却有着较为良好的稳定性.此外,数值模拟光孤子间的相互作用发现,由标准非线性薛定谔方 ...
    本站小编 Free考研考试 2021-12-29
  • 宇宙线高能粒子对测试质量充电机制
    摘要:测试质量是空间引力波测量的核心传感器,宇宙线高能粒子能够穿透航天器屏蔽对其造成电荷注入,进而产生库仑力和洛伦兹力噪声对引力波科学探测造成严重影响.本文采用蒙特卡洛仿真方法,探究了不同宇宙线高能粒子对测试质量的充电过程和机制.研究结果表明,在同一能谱下随着截止能量的降低充电速率逐步增大,充电速率 ...
    本站小编 Free考研考试 2021-12-29