A MULTI-RESOLUTION PD-SPH COUPLING APPROACH FOR STRUCTURAL FAILURE UNDER FLUID-STRUCTURE INTERACTION
-
摘要: 流固耦合破坏是一类涉及结构变形与破坏以及复杂自由表面现象的强非线性力学问题. 结合近场动力学(peridynamics, PD)与光滑粒子流体动力学(smoothed particle hydrodynamics, SPH)各自的优势并考虑其计算效率问题, 提出一种适用于分析流−固耦合破坏问题的多分辨率PD-SPH混合方法. 分别采用SPH和PD方法以不同的空间和时间分辨率对流体和结构进行离散与求解, 利用具有与流体粒子相同光滑长度的虚粒子处理流−固界面, 以高精度满足界面边界条件. 通过两个经典算例: 液柱静压力下弹性板的变形和溃坝流体冲击弹性闸门的变形问题, 表明提出的多分辨率PD-SPH方法兼具较高的计算精度和计算效率; 对含裂缝的Koyna重力坝水力劈裂问题进行模拟计算, 所得裂缝扩展路径与文献结果吻合, 说明该方法适用于涉及结构破坏的流固耦合问题仿真. 最后尝试采用该方法进行流体冲击作用下含裂纹混凝土板崩塌过程数值仿真, 准确描述混凝土板的断裂破坏和全过程中的流体运动. 多分辨率PD-SPH混合方法或可为流−固耦合作用下的结构损伤破坏仿真提供一种新选择.Abstract: Structural failure under fluid-structure interaction (FSI) is a type of strong nonlinear problem, which involves structural motion, deformation and failure as well as complex free-surface flows. Considering the respective advantages of peridynamics (PD) and smoothed particle hydrodynamics (SPH) as well as their computational efficiency, a multi-resolution PD-SPH coupling approach suitable for solving complicated FSI-concerned structural failure problems was proposed. The fluid and solid are discretized and solved by using SPH and PD approaches with different spatial and temporal resolutions, respectively. To achieve the precise satisfaction of interface boundary conditions, the fluid-structure interface is treated by using virtual particle technology, in which the same smoothing length of virtual particles as fluid particles is adopted. The modeling and analysis for two benchmark tests: large deformation of an elastic plate with hydrostatic pressure, and dam-break flow through an elastic gate, show that the presented multi-resolution PD-SPH coupling strategy and approach is suitable for simulating fluid-structure-interaction problems with satisfactory accuracy and efficiency. Further, the process of hydraulic fracture in Koyna gravity dam with an initial crack is simulated, and the cracking path in the simulation agrees well with available literature results, which indicates that the proposed coupling approach is appropriate for solving FSI-concerned structural failure problems. Finally, the proposed coupling strategy and numerical approach is employed to investigate the collapse process of a concrete slab due to fluid flow impacting, and the whole process of the concrete slab fracture as well as the motion of the fluid are captured with high accuracy. The results show that the proposed multi-resolution PD-SPH coupling approach may provide a potential alternative to simulate the process of structural failure under fluid-structure interaction.
-
引 言
流−固耦合(fluid-structure interaction, FSI)是一种流体与固体相互作用和影响的力学问题, 其广泛存在于诸多工程领域, 如水上降落、液舱晃荡以及海浪砰击舰船与海洋平台等. 由于流体流动、结构变形破坏以及二相间相互作用问题的复杂性, 流−固耦合相关问题, 特别是流−固耦合导致结构破坏问题的数值模拟一直是研究难点. 已有诸多数值方法用于求解FSI问题, 其中, 基于网格的数值方法因计算精度和效率高而被广泛应用, 如有限单元法(finite element method, FEM)[1-2]、有限体积法(finite volume method, FVM)[3]、有限差分法(finite difference method, FDM)[4]及其混合方法[5-6]等. 然而, 基于网格的数值方法在求解FSI问题时往往需要网格重构、界面追踪等特殊数值技术, 增加了算法复杂性, 并且在处理三维裂纹萌生、动态多裂纹扩展等强不连续问题时依然面临挑战. 无网格方法是另一类被广泛使用于FSI问题的数值方法, 其不受网格约束, 能在拉格朗日框架下求解流体域或固体域, 在复杂流动问题模拟以及固体结构的大变形与断裂分析等方面展现出独特的优势, 已成为研究FSI问题的重要手段.
近场动力学(peridynamics, PD)[7-9]与光滑粒子流体动力学(smoothed particle hydrodynamics, SPH)[10]是两种基于非局部假定的无网格方法, 通常将计算域离散为一系列自由移动的粒子, 通过求解控制每个粒子运动的积分方程或偏微分方程获得整个系统的力学行为. PD方法无需预先确定裂纹位置, 可自然地描述裂纹萌生和扩展, 近年来成功应用于固体断裂[11]、冲击破坏[12-13]等不连续力学问题模拟. 针对FSI问题, 部分学者建立了不同的PD模型, 如Zhang等[14]和Zhou等[15]分别建立了用于水力压裂模拟的PD模型, Ren等[16]建立了模拟海冰在海浪作用下力学行为的PD模型. 然而, 当前的PD模型通常不适用于含自由表面流体的FSI问题, 涉及复杂流动的PD相关算法仍有待进一步研究[17-18]. SPH方法适用于流体运动模拟, 且能处理结构动态响应问题. 近年来, Oger等[19]应用弱可压SPH方法对刚性楔形体入水问题进行研究, 获得了较为理想的数值结果, Sun等[20]结合δ+−SPH与完全拉格朗日格式的SPH方法研究了多种流体与弹性结构相互作用问题, Khayyer等[21]建立了不可压缩SPH与传统SPH混合方法, 将其应用于波浪冲击可变形结构问题. 尽管SPH可以单独处理FSI问题, 但目前很少有研究考虑强耦合作用导致的固体破坏.
在PD−SPH混合方法研究方面, Fan等[22-23]针对土壤的爆炸破碎问题, 采用非常规态型PD模型建立了一种固体粒子与流体粒子互为虚粒子的混合方法, Sun等[24]、姚学昊等[25]针对FSI问题分别采用基于接触力的混合方法. 上述方法均使用相同的粒子分辨率离散固体域和流体域, 在求解结构空间尺度远小于流体的问题时需要较高的计算成本. Rahimi等[26]通过一种基于拉格朗日插值的PD−SPH混合方法, 实现了空间多分辨率离散, 但未考虑多时间步长问题.
基于上述研究现状, 本文提出一种用于流−固耦合结构破坏模拟的多时间和空间分辨率PD−SPH混合方法. 该方法采用不同的粒子分辨率和时间步长对固体域和流体域进行计算以降低计算成本. 通过典型的自由表面流体与可变形结构相互作用问题算例验证了方法的精度和效率, 并应用于流体冲击混凝土板破坏过程分析.
1. 数值模拟方法
1.1 PD方法
传统连续介质力学理论中物体内部的相互作用由应力张量的散度
${\nabla} \cdot {{\boldsymbol{\sigma}} }$ 表示[9], 而PD方法中固体域被离散为一系列包含质量m和密度ρ等基本信息的无拓扑关系的粒子Xi, 每个粒子与其近场范围HX (半径为δ的球形区域)内的其他粒子Xj间存在相互作用f$$ {\rho _i}{\ddot {\boldsymbol{u}}_i} = \int_{{H_X}} {{\boldsymbol{f}}({{\boldsymbol{X}}_j},{{\boldsymbol{X}}_i},t)\,} {\text{d}}{V_j} + \tilde {{\boldsymbol{b}}}({{\boldsymbol{X}}_i},t) $$ (1) 式中, i和j为粒子编号;
${\boldsymbol{u}}$ 和$\ddot {\boldsymbol{u}}$ 分别为位移和加速度;$\tilde {{\boldsymbol{b}}}({{\boldsymbol{X}}_i},t)$ 表示t时刻作用在粒子Xi上的外力密度;${\boldsymbol{X}}$ 为PD粒子在初始构形中的位置矢量;$ {V_j} $ 表示粒子Xj的体积.在态型PD理论[7]中, 分别定义参考位置矢量状态
$\underline {\boldsymbol{X }}$ 、形变矢量状态$\underline {\boldsymbol{Y }}$ 、拉伸标量状态$\underline {{e }}$ 及力矢量状态$\underline {\boldsymbol{T}}$ 如下$$ \underline {\boldsymbol{X}} \left\langle {\boldsymbol{\xi}} \right\rangle = {\boldsymbol{\xi}} ,\,\quad \underline {\boldsymbol{Y }} \left\langle {\boldsymbol{\xi}} \right\rangle = {\boldsymbol{\xi}} + {\boldsymbol{\eta}} $$ (2) $$ {\boldsymbol{\xi}} = {{\boldsymbol{X}}_j} - {{\boldsymbol{X}}_i},\quad {\boldsymbol{\eta}} = {{\boldsymbol{u}}_j} - {{\boldsymbol{u}}_i} $$ (3) $$ \underline{e}=\underline{y}-\underline{x}, \quad \underline{y}=|\underline{\boldsymbol{Y}}|, \quad \underline{x}=|\underline{\boldsymbol{X}}| $$ (4) $$ \underline{\boldsymbol{T}}=\underline{\boldsymbol{T}}(\underline{\boldsymbol{Y}}, \varLambda)$$ (5) 式中,
$\varLambda $ 为与本构模型相关的变量. 对于简单弹性材料, 力矢量状态仅与形变矢量状态$\underline{\boldsymbol{Y}}$ 有关, 即$\underline{\boldsymbol{T}} = \underline{\boldsymbol{T}} (\underline{\boldsymbol{Y}} )$ . 同时, 在常规态型PD理论中, 引入力标量状态$\underline{t}$ 使得$$ \underline{\boldsymbol{T}}=\underline{t} \,\underline{\boldsymbol{M}} $$ (6) 式中,
$\underline{\boldsymbol{M}}$ 是变形方向矢量状态, 其定义为$$ \underline{\boldsymbol{M}}=\left\{\begin{array}{cl} \dfrac{\underline{\boldsymbol{Y}}}{|\underline{\boldsymbol{Y}}|}, & |\underline{\boldsymbol{Y}}| \neq 0 \\ 0, & \text { otherwise } \end{array}\right.$$ (7) 定义体积膨胀标量
$ \theta $ $$ \theta=\frac{3}{q_{\mathrm{v}}} \int_{H_{X}}(\underline{\omega} \,\underline{x}) \cdot \underline{e} \mathrm{~d} V_{j} $$ (8) 式中,
$\underline{\omega}$ 为影响函数;$ {q_{\text{v}}} $ 为加权体积因子$$q_{\mathrm{v}}=\int_{H_{X}}(\underline{\omega}\, \underline{x}) \cdot \underline{x} \mathrm{~d} V_{j}$$ (9) 利用体积膨胀标量可将拉伸标量状态分解为球量部分
$\underline{e}^{{\rm{s}}}=\theta \underline{x} / 3$ 和偏量部分$\underline{e}^{\mathrm{d}}=\underline{e}-\underline{e}^{\mathrm{s}}$ . 此时, 借鉴经典理论中弹性应变能密度表达式, 可得到适用于简单弹性材料的常规态型PD变形能密度函数$$ W_{\mathrm{P}}(\underline{\boldsymbol{Y}})=W_{\mathrm{P}}\left(\theta, \underline{e}^{\mathrm{d}}\right)=\frac{k}{2} \theta^{2}+\frac{\alpha}{2}\left(\omega \underline{e}^{\mathrm{d}}\right) \cdot \underline{e}^{\mathrm{d}} $$ (10) 式中, k和α是与体积模量K和剪切模量G相关的PD材料参数, 可根据PD理论中变形能密度与传统理论的应变能密度等效获得[27]. 计算变形能密度函数
$W_{\mathrm{P}}(\underline{\boldsymbol{Y}})$ 关于拉伸标量状态$\underline{e}$ 的弗雷谢(Fréchet)导数则可得到力标量状态$\underline{t}$ $$ \underline{t}=\frac{3 k \theta}{q_{\mathrm{v}}} \underline{\omega} \,\underline{x}+\alpha \underline{\omega} \,\underline{e}^{\mathrm{d}} $$ (11) 此时, f可表示为
$$ \boldsymbol{f}\left(\boldsymbol{X}_{j}, \boldsymbol{X}_{i}, t\right)=\left(\underline{t}\left[\boldsymbol{X}_{j}, t\right]\langle-\boldsymbol{\xi}\rangle+\underline{t}\left[\boldsymbol{X}_{i}, t\right]\langle{\boldsymbol{\xi}}\rangle\right) \underline{\boldsymbol{M}} $$ (12) 针对断裂问题, 假定两个粒子Xi和Xj之间的相对位置满足一定关系时, 粒子对内粒子间的相互作用消失. 可引入间断函数表示
$$ {\mu _{ij}}({{\boldsymbol{X}}_i},{{\boldsymbol{X}}_j},t) = \left\{ {\begin{array}{*{20}{l}} 1,&{{\text{ }}\,{s_{ij}} < {s_0}} \\ 0,&{{\text{ otherwise}}} \end{array}} \right. $$ (13) 式中,
${s_{ij}} = \left( {\left| {{\boldsymbol{\xi}} + {\boldsymbol{\eta}} } \right| - \left| {\boldsymbol{\xi}} \right|} \right)/\left| {\boldsymbol{\xi }} \right|$ 为粒子对的伸长率;${s_0} = \sqrt {5{G_0}/(9 K\delta) }$ 为临界伸长率, G0为传统的临界断裂能密度.根据粒子在某一时刻“键”的断裂数目的统计, 可定义局部损伤$$ {D_i}({{\boldsymbol{X}}_i},t) = 1 - \frac{{\displaystyle\int_{{H_X}} {{\mu _{ij}}({{\boldsymbol{X}}_i},{{\boldsymbol{X}}_j},t){\rm{d}}{V_j}} }}{{\displaystyle\int_{{H_X}} {{\rm{d}}{V_j}} }} $$ (14) 1.2 SPH方法
在SPH方法中, 同样将流体域离散为一系列粒子xi, 通过核函数插值和粒子近似技术可得到离散形式的不可压缩流体控制方程
$$ \left. {\begin{array}{l} {\dfrac{{{\text{D}}{\rho _i}}}{{{\text{D}}t}} = {\rho _i}\displaystyle\sum\limits_j^{} {{{\boldsymbol{v}}_{ij}}{\nabla _i}} {W_{ij}}{V_j} + {\delta _{\text{s}}}h{c_0}\displaystyle\sum\limits_j^{} {{{\boldsymbol{\psi}} _{ij}}{\nabla _i}{W_{ij}}{V_j}} } \\ {{m_i}\dfrac{{{\text{D}}{v_i}}}{{{\text{D}}t}} = - \displaystyle\sum\limits_j^{} {({p_i} + {p_j}){\nabla _i}{W_{ij}}{V_i}{V_j}} + {{\boldsymbol{f}}_{{\Pi }}}{{ + }}{{\boldsymbol{f}}_{{\upsilon }}} + {{\boldsymbol{f}}_{\text{g}}}} \\ {p = c_0^2(\rho - {\rho _0})} \\ {{{\boldsymbol{\psi}} _{ij}} = 2({\rho _j} - {\rho _i})\dfrac{{{{\boldsymbol{x}}_{ji}}}}{{{{| {{{\boldsymbol{x}}_{ji}}} |}^2}}} - (\left\langle {\nabla \rho } \right\rangle _i^{\text{L}} + \left\langle {\nabla \rho } \right\rangle _j^{\text{L}})} \end{array}} \right\} $$ (15) 式中,
${\text{D/D}}t$ 表示物质导数;${{\boldsymbol{x}}_{ji}} = {{\boldsymbol{x}}_j} - {{\boldsymbol{x}}_i}$ 为粒子在当前构形中的相对位置矢量,${{\boldsymbol{v}}_{ij}} = {{\boldsymbol{v}}_i} - {{\boldsymbol{v}}_j}$ 为相对速度矢量; h为光滑长度; p为压力;${{\boldsymbol{f}}_{\text{g}}}$ ,${{\boldsymbol{f}}_{{\Pi }}}$ 和${{\boldsymbol{f}}_{{\upsilon }}}$ 分别为体力、人工黏性力和物理黏性力;$ {c_0} $ 为参考声速, 其值应足够大以满足弱可压缩假定[28]$$ {c_0} \geqslant 10 \max \left( {{{\left|{\boldsymbol{ v}} \right|}_{\max }},\sqrt {{p_{\max }}/{\rho _0}} } \right) $$ (16) 式中,
${\left| {\boldsymbol{v}} \right|_{\max }}$ ,$ {p_{\max }} $ 和$ {\rho _0} $ 分别为流体的最大流速、最大压力以及参考密度.${\nabla _i}$ 表示对粒子xi坐标的空间导数;${W_{ij}} = W({{\boldsymbol{x}}_{ij}},h)$ 为核函数, 本文选用B-样条核函数[10]$$ W({{\boldsymbol{x}}_{ij}},h) = \dfrac{{{\alpha _{\text{w}}}}}{{{h^{dim}}}} \cdot \left\{ {\begin{array}{l} {\dfrac{2}{3} - {q_{\text{s}}} + \dfrac{1}{2}q_{\text{s}}^{\text{3}},\quad 0 \leqslant {q_{\text{s}}} < 1} \\ {\dfrac{1}{6}{{(2 - {q_{\text{s}}})}^3},\quad \;\,\,\,\,\,\,1 \leqslant {q_{\text{s}}} < 2} \\ \,{0,\quad \quad \quad \;\;\;\;\;\;\,\,\quad {q_{\text{s}}} \geqslant 2\quad {\kern 1pt} \,} \end{array}} \right. $$ (17) 式中,
${q_{\text{s}}} = |{{\boldsymbol{x}}_{ij}}|/h$ ; 空间维度dim = 1, 2, 3分别对应${\alpha _{\text{w}}} = {\text{1}}$ ,$15/7{\text{π }}$ ,${\text{3}}/{{2{\text{π }} }}$ . 控制方程组(15)中的连续性方程右侧第二项为人工耗散项[28], 用以抑制因数值噪声产生的压力振荡, 其中$ {\delta _{\text{s}}} $ 为耗散强度控制常数, 常取值为0.1;$ \left\langle {\nabla \rho } \right\rangle _{}^{\text{L}} $ 为修正密度梯度$$ \left. {\begin{array}{l} {\left\langle {\nabla \rho } \right\rangle _i^{\text{L}} = \sum\limits_j {({\rho _j} - {\rho _i}){\boldsymbol{L}_i}{\nabla _i}{W_{ij}}{V_j}} } \\ {{\boldsymbol{L}_i} = {{\left[ {\displaystyle\sum\limits_j {({{\boldsymbol{x}}_j} - {{\boldsymbol{x}}_i}) \otimes {\nabla _i}{W_{ij}}{V_j}} } \right]}^{ - 1}}} \end{array}} \right\} $$ (18) 对于不可压缩流体, 粒子xi受到的物理黏性力可简化为
$$ {{\boldsymbol{f}}_{{\upsilon }}}({{\boldsymbol{x}}_i}) = \sum\limits_j^{} {\frac{{2\mu {{\boldsymbol{x}}_{ij}} \cdot {\nabla _i}{W_{ij}}{V_i}{V_j}}}{{({{| {{{\boldsymbol{x}}_{ij}}} |}^2} + {\varphi ^2})}}{{\boldsymbol{v}}_{ij}}} $$ (19) 式中,
$ \mu $ 为流体动力黏度;$ \varphi = 0.1 h $ 用于防止奇点的出现. 为了提高数值稳定性, Monaghan[29]在动量方程中加入了人工黏性项${{\boldsymbol{f}}_{{\Pi }}}$ 以抑制数值振荡, 其表达式为$$ {{\boldsymbol{f}}_{{\Pi }}}({{\boldsymbol{x}}_i}) = \sum\limits_j {{m_i}{m_j}{\alpha _{\text{π }}}h{{\bar c}_{ij}}\frac{{{{\boldsymbol{v}}_{ij}} \cdot {{\boldsymbol{x}}_{ij}}}}{{{{\bar \rho }_{ij}}\left( {{{| {{{\boldsymbol{x}}_{ij}}} |}^2} + {\varphi ^2}} \right)}}{\nabla _i}{W_{ij}}} $$ (20) 式中,
$ {\alpha _{\text{π }}} $ 为控制参数, 在本文的模拟中取值为0.1;$ {\bar c_{ij}} $ 和$\;{\bar \rho _{ij}}$ 分别为平均声速和平均密度.2. 多分辨率PD−SPH耦合方案
本文采用具有不同初始间距的固体粒子和流体粒子离散计算域, 并通过流体粒子−虚粒子接触耦合方案处理流−固交界面. 为了更显著地提高计算效率, 利用不同的时间步长对固体和流体控制方程进行时间积分.
2.1 多分辨率空间离散
如图1所示, 当SPH流体粒子
$a$ 靠近流−固界面时, PD固体粒子$b$ 和c均位于SPH粒子a的支持域Sa内, 此时粒子$b$ 和c将以虚粒子$b'$ 和$c'$ 的形式加入粒子$a$ 的邻近搜索列表, 从而保证粒子$ {a} $ 的支持域不被流−固界面截断, 实现SPH边界缺陷修正. 考虑到PD粒子与SPH粒子不同的初始间距会导致虚粒子与流体粒子光滑长度不同, 本文中令虚粒子的光滑长度等于SPH粒子的光滑长度h, 这也使得流体粒子$a$ 的支持域Sa的大小与虚粒子$b'$ ,$c'$ 的支持域${S_{b'}}$ ,${S_{c'}}$ 完全相同. 此时, SPH粒子$a$ 的控制方程变为$$ \left. {\begin{array}{l} {\dfrac{{{\text{D}}{\rho _a}}}{{{\text{D}}t}} = {\rho _a}\displaystyle\sum\limits_n^{} {{{\boldsymbol{v}}_{an}}{\nabla _a}} {W_{an}}{V_n} + {\delta _{\text{s}}}h{c_0}\displaystyle\sum\limits_n^{} {{{\boldsymbol{\psi}} _{an}}{\nabla _a}{W_{an}}{V_n}} } \\ {m_a}\dfrac{{{\text{D}}{{\boldsymbol{v}}_a}}}{{{\text{D}}t}} = - \displaystyle\sum\limits_n^{} {({p_a} + {p_n}){\nabla _a}{W_{an}}{V_a}{V_n}} + {{\boldsymbol{f}}_{{\Pi }}}{\text{ + }}{{\boldsymbol{f}}_{{\upsilon }}} + {{\boldsymbol{f}}_{\text{g}}} + {{\boldsymbol{f}}_{\text{r}}} \end{array}} \right\} $$ (21) 式中, n表示与SPH粒子a相邻的其他粒子编号, 包括流体粒子以及由PD粒子
$\kappa $ 伪装的虚粒子$\kappa '$ (含粒子$b'$ 和$c'$ ); 黏性力${{\boldsymbol{f}}_{{\Pi }}}$ 和${{\boldsymbol{f}}_{{\upsilon }}}$ 的计算同样包含虚粒子$\kappa '$ ;${{\boldsymbol{f}}_{\text{r}}}$ 为虚粒子对SPH粒子施加的排斥力, 用以防止粒子穿透流−固界面, 其表达式为$$ {{\boldsymbol{f}}_{\text{r}}}({{\boldsymbol{x}}_a}) = - \sum\limits_{\kappa '} {{K_{\text{r}}}{n_{\text{r}}}\frac{{W{{\left( {{{\boldsymbol{x}}_{a\kappa '}}} \right)}^{{n_{\text{r}}} - 1}}}}{{W{{\left( h \right)}^{{n_{\text{r}}}}}}}{\nabla _a}W\left( {{{\boldsymbol{x}}_{a\kappa '}}} \right){V_a}{V_{\kappa '}}} $$ (22) 式中,
${K_{\text{r}}}$ 和${n_{\text{r}}}$ 为用户自定义参数[10].需要注意的是, 在该耦合方法中, 虚粒子
$\kappa '$ 根据PD粒子在初始构形中的位置分为表面虚粒子(如$b'$ )和内部虚粒子(如$c'$ ), 仅被动地被流体粒子搜索, 自身不进行SPH数值积分, 其位置信息来自对应的PD粒子, 密度等场变量信息则通过插值计算获得[25]. 如图1所示, 对于表面虚粒子$b'$ , 其场变量信息来自其支持域${S_{b'}}$ 内的流体粒子; 内部虚粒子$c'$ 的场变量信息则来自其支持域${S_{c'}}$ 内的流体粒子和表面虚粒子. 虚粒子速度(无滑移边界条件)和密度表达式为$$ {{\boldsymbol{v}}_{\kappa '}} = 2{{\boldsymbol{v}}_\kappa } - {\tilde {\boldsymbol{v}}_{\kappa '}} = 2{{\boldsymbol{v}}_\kappa } - \sum\limits_n {{{\boldsymbol{v}}_n}W_{_{\kappa 'n}}^{{\text{new}}}{V_n}} $$ (23) $$ {\rho _{\kappa '}} = \sum\limits_n {{\rho _n}W_{\kappa 'n}^{{\rm{new}}}{V_n}} $$ (24) 式中,
$W_{\kappa 'n}^{{\rm{new}}} = {W_{\kappa 'n}}/\sum\limits_n {{W_{\kappa 'n}}{V_n}}$ [30];${{\boldsymbol{v}}_\kappa }$ 为PD粒子$\kappa $ 的速度;${\tilde {\boldsymbol{v}}_{\kappa '}}$ 为虚粒子$ \kappa ' $ 支持域内其他粒子的光滑速度[10]. 此时, 耦合方案在流−固界面处满足运动学条件.对于PD部分, 界面处PD粒子受到SPH粒子的作用力, 相当于在耦合界面处对PD粒子施加力边界条件. 为了在流−固界面处满足动力学条件, 保证系统动量守恒, 令SPH粒子
$a$ 对PD粒子$\kappa $ 的作用力${{{\boldsymbol{F}}}_{ato\kappa }}$ $$ {{\boldsymbol{F}}_{ato\kappa }} = - {{\boldsymbol{F}}_{\kappa toa}} = - {{\boldsymbol{F}}_{\kappa 'toa}} $$ (25) 由式(21)可得虚粒子
$\kappa '$ 对SPH粒子$a$ 的作用力$$ {{\boldsymbol{F}}_{\kappa 'toa}} = - \sum\limits_{\kappa '}^{} {({p_a} + {p_{\kappa '}}){\nabla _a}{W_{a\kappa '}}{V_a}{V_{\kappa '}}} + {{\boldsymbol{f}}_{{\Pi }}}{\text{ + }}{{\boldsymbol{f}}_{{\upsilon }}} + {{\boldsymbol{f}}_{\text{r}}} $$ (26) 此时, PD粒子
$\kappa $ 的运动方程可写为$$ {\rho _\kappa }{\ddot {\boldsymbol{u}}_\kappa } = \int_{{H_X}} {{\boldsymbol{f}}({{\boldsymbol{X}}_\kappa },{{\boldsymbol{X}}_j},t)\,} {\text{d}}{V_n} + \tilde {\boldsymbol{b}}({{\boldsymbol{X}}_\kappa },t) + \frac{{{{\boldsymbol{F}}_{ato\kappa }}}}{{{V_\kappa }}} $$ (27) 2.2 多分辨率时间离散
考虑到PD与SPH的完全显式计算, 分别定义结构和流体的计算时间步
$\Delta {t_{\text{S}}}$ 和$\Delta {t_{\text{F}}}$ , 其中, 固体计算的时间步长应满足[31]$$ \Delta {t_{\text{S}}} \propto \frac{\delta }{{c'}} $$ (28) 式中,
$c' = \sqrt {\dfrac{{E(1 - \nu )}}{{(1 + \nu )(1 - 2\nu ){\rho _{\text{S}}}}}}$ 为固体材料声速, E为弹性模量,$\nu $ 为泊松比,${\;\rho _{\text{S}}}$ 为固体材料密度. 流体计算的时间步长则应满足[32]$$ \Delta {t_{\text{F}}} \leqslant 0.2\min [h/({c_0} + {\left| {\boldsymbol{v}} \right|_{\max }}),\sqrt {h/{{\left| {\dot {\boldsymbol{v}}} \right|}_{\max }}} ] $$ (29) 通常, 固体计算的时间步长
$\Delta {t_{\text{S}}}$ 会小于流体计算的时间步长$ \Delta {t_{\text{F}}} $ , 若采用单一时间步长会出现流体计算耗时过多的问题. 本文将采用多时间步长方案, 并通过执行一次流体时间积分同时进行$\ell = \Delta {t_{\text{F}}}/\Delta {t_{\text{S}}}$ 次固体时间积分的方式保证流体与固体计算的时间一致性. 该方案示意图如图2所示.多分辨率PD−SPH耦合计算流程如图3所示.
![]()
图 3 多分辨率PD-SPH耦合计算流程Figure 3. Flow chart for multi-resolution PD-SPH coupling calculation3. 算例验证
为验证多分辨率PD−SPH混合方法的精度和计算效率, 首先对经典的FSI问题进行数值分析. 随后对流固相互作用导致的结构断裂破坏问题进行模拟, 进一步验证耦合方案对于FSI问题的适用性. 本文算例均在Intel i7-7700 CPU (主频
$3.6\,{\text{GHz}}$ )个人计算机上完成, 所取光滑长度和近场范围半径分别为$h = 1.3 d{x_{\text{F}}}$ 和$\delta = 3.0 d{x_{\text{S}}}$ , 其中$d{x_{\text{F}}}$ 和$d{x_{\text{S}}}$ 分别表示SPH粒子与PD粒子的初始间距.3.1 液柱静压力下的铝板变形
首先对液柱静压力作用下的铝板变形问题进行数值模拟. 如图4所示, 长
$L = {\text{1}}{\text{.0}}\,\,{\text{m}}$ 、厚$d{\text{ = }}\,{\text{0}}{\text{.05}}\,\,{\text{m}}$ , 密度${\rho _{\text{S}}} = 2700\,\,{\text{kg/}}{{\text{m}}^3}$ , 弹性模量$E = 67.5\,\,{\text{GPa}}$ , 泊松比$\nu = 0.34$ 的铝板突然受到高$H = 2\,\,{\text{m}}$ , 密度${\;\rho _{\text{F}}} = 1000\;{\text{kg/}}{{\text{m}}^3}$ 的水柱作用. 根据文献[33], 铝板经过初始振荡后达到平衡状态, 其跨中挠度解析解为$ - 6.85 \times {10^{ - 5}}\,\,{\text{m}}$ . 模拟中, 参考声速${c_0} = 60\;{\text{m/s}}$ , 总模拟时长为$5\;{\text{s}}$ .![]()
图 4 液柱静压力下铝板变形的初始条件(单位: m)Figure 4. Initial condition for the aluminum plate under hydrostatic pressure (unit: m)针对当前问题, 考虑到过大的流−固分辨率比会对计算精度和效率会产生一定的影响, 设置了固定结构空间分辨率条件下的3种多分辨率工况, 具体参数设置如表1所示. 图5给出了3种工况下铝板跨中挠度时程, 由图可见, 经过初始振荡, 3种工况下铝板跨中挠度均收敛至约为
$ - {\text{6}}{\text{.70}} \times {\text{1}}{{\text{0}}^{ - {\text{5}}}}{\text{ m}}$ , 与解析解相比误差仅为2.19%. 图6给出了$t = 5\;{\text{s}}$ 时工况III的流体粒子分布与压力场以及铝板挠度场, 可以看出, 此时流体压力场与铝板挠度场均是光滑的. 上述结果表明, 本文多分辨率PD−SPH混合方法计算精度较高.表 1 液柱静压力下铝板变形问题工况布置Table 1. The setup of different cases for aluminum plate deformation with hydrostatic pressureCases $\dfrac{{d{x_{\text{F}}}}}{{d{x_{\text{S}}}}}$ $d{x_{\text{S}}}{\text{/m}}$ $\dfrac{{\Delta {t_{\text{F}}}}}{{\Delta {t_{\text{S}}}}}$ $ \Delta {t_{\text{S}}}{\text{/s}} $ I 1.0 6.25 × 10−3 5.0 1.0 × 10−6 II 2.0 6.25 × 10−3 10.0 1.0 × 10−6 III 4.0 6.25 × 10−3 20.0 1.0 × 10−6 ![]()
图 6 t = 5 s时工况III的流体压力场与板的挠度场Figure 6. Fluid pressure and plate deflection field in case-III at t = 5 s图7为3种工况下FSI系统归一化总能量[21](
$ {E_{\text{T}}}{\text{/}}E_{\text{T}}^{\text{0}} - 1 $ ,$ {E_{\text{T}}} $ 和$E_{\text{T}}^{\text{0}}$ 分别为当前和初始时刻系统总能量)时间历程. 由图可见, 归一化总能量同样存在振动变化过程, 随着FSI系统进入稳定状态, 系统总能量最终保持稳定.$t = 5\;{\text{s}}$ 时, 3种工况的系统总能量损失分别为0.23%, 0.34%以及0.49%, 表明多分辨率混合方法具有较好的能量守恒特性. 综合图5与图7可以发现, 铝板位移初始振荡的持续时间与FSI系统总能量损失有关, 系统能量损伤越少, 其位移振荡时间越长. 该发现与Zhang等[34]所得结论一致, 进一步验证了本文的多分辨率混合方法.![]()
图 7 FSI系统归一化总能量时间历程Figure 7. Time histories for normalized total energy of FSI system图8给出了3种工况下全仿真过程计算时长, 可见, 随着流体粒子间距增大, 流体计算所需时间步长也增大, 总模拟时长显著减少. 工况III的计算速度约为工况I的11.3倍. 由此表明, 多分辨率PD−SPH混合方法能够在保证计算精度的同时有效提高计算效率.
3.2 溃坝流体冲击弹性闸门
Antoci等[35]对溃坝水流冲破弹性闸门问题进行了实验研究, 该问题的初始条件示意图如图9所示. 实验在宽
$L = {\text{0}}{\text{.1}}\;{\text{m}}$ 、内部水深${H_1}{\text{ = }}\;{\text{0}}{\text{.14}}\;{\text{m}}$ 的水箱中进行. 箱体左侧刚性壁面下部放置了上端固定下端自由的弹性橡胶闸门, 其宽度为$s{\text{ = 0}}{\text{.005}}\;{\text{m}}$ 、高度为$H = 0.079\;{\text{m}}$ , 材料参数为: 密度$\;{\rho _{\text{S}}}{\text{ = 11}}00\;{\text{kg/}}{{\text{m}}^3}$ , 弹性模量$E{\text{ = 7}}{\text{.8}}\;{\text{MPa}}$ , 泊松比$\nu = 0.4{\text{7}}$ . 数值模拟中, 参考声速${c_0} = 15\;{\text{m/s}}$ , 水体密度$\;{\rho _{\text{F}}}{\text{ = }}1000\;{\text{kg/}}{{\text{m}}^3}$ , 总模拟时长$0.4\;{\text{s}}$ . 针对该问题的多分辨率工况设置如表2所示.![]()
图 9 溃坝流体冲破弹性闸门的初始条件示意图(单位: mm)Figure 9. Schematic sketch for initial conditions in the dam-break flow through an elastic gate (unit: mm)表 2 溃坝流体冲破弹性闸门问题工况布置Table 2. The setup of different cases for dam-break flow through an elastic gateCases $\dfrac{{d{x_{\text{F}}}}}{{d{x_{\text{S}}}}}$ $d{x_{\text{S}}}{\text{/m}}$ $\dfrac{{\Delta {t_{\text{F}}}}}{{\Delta {t_{\text{S}}}}}$ $ \Delta {t_{\text{S}}}{\text{/s}} $ I 1.0 1.0 × 10−3 1.0 5.0 × 10−6 II 2.0 1.0 × 10−3 1.0 5.0 × 10−6 III 2.0 1.0 × 10−3 4.0 5.0 × 10−6 图10给出了工况III中不同时刻的弹性闸门在时变水压作用下的变形以及自由液面变化过程, 并与实验观察结果[35]进行比较. 由图10可以看出, PD−SPH所得弹性闸门与自由液面形状与实验结果基本吻合. 在
$t = 0.08\,{\text{s}}$ 时, 闸门自由端在水压力作用下向$x$ 轴负方向移动, 形成水流出口. 在$t = 0.16\;{\text{s}}$ 时, 闸门自由端位移与水流出口显著增大. 随后, 由于水位的降低, 水压力减小, 水流出口逐渐减小. 值得注意的是, 由于本文开展的是二维模拟, 实验中的流体飞溅现象不会出现.![]()
图 10 不同时刻弹性闸门变形和自由液面演变Figure 10. Elastic gate deformation and change of fluid free surface图11为弹性闸门自由端水平位移时程曲线. 通过与文献[34]的SPH模拟结果以及文献[35]的实验数据对比可以看出, 3种工况下所得位移结果均吻合较好. 在PD−SPH模拟中, 弹性闸门于
$t = 0.122\;{\text{s}}$ 时刻达到最大变形, 此时其自由端水平位移约为0.043 m. 图12给出了3种工况的计算总时长, 比较可知, 仅采用多分辨空间离散方法(工况II)即可有效提高计算效率; 若采用多时间与空间分辨率, 即工况III, 其计算速度为单一分辨率(工况I)的8.5倍.![]()
图 11 弹性闸门自由端水平位移时程Figure 11. Time histories of horizontal displacements at the free end of elastic gate综上, 本文提出的PD−SPH混合方法在不同分辨率工况下均能有效传递流−固相互作用信息, 可用于涉及复杂流体流动和结构大变形的FSI问题.
3.3 重力坝水力劈裂
为了验证多分辨率PD−SPH混合方法对于FSI致结构破坏问题的适用性, 选取经典的带裂缝Koyna重力坝问题, 分析其在超载作用下的断裂响应. 初始条件如图13所示, 坝高
${\text{103}}\;{\text{m}}$ 、坝顶和坝底宽度分别为$w = 14.8\;{\text{m}}$ 和$L = 70\;{\text{m}}$ . 水平裂缝位于大坝上游面$66.5\;{\text{m}}$ 高程处, 长度为${\text{1}}{\text{.93}}\;{\text{m}}$ . 大坝材料参数为: 密度$\;{\rho _{\text{S}}}{\text{ = 2}}450\;{\text{kg/}}{{\text{m}}^3}$ , 弹性模量$E{\text{ = 25}}\;{\text{GPa}}$ , 泊松比$\nu = 0.2$ . 考虑大坝承受水头超载作用, 根据文献[36], 取上游面超高水头为$10\;{\text{m}}$ , 即大坝上游面水深为113 m. 模拟中, PD粒子间距$d{x_{\text{S}}} = 0.2\;{\text{m}}$ , 时间步长$\Delta {t_{\text{S}}} = 5 \times {10^{ - 6}}{\text{ s}}$ ; 水体密度$\;{\rho _{\text{F}}}{\text{ = }}1000\;{\text{kg/}}{{\text{m}}^3}$ , SPH粒子间距$d{x_{\text{F}}} = 0.4\;{\text{m}}$ , 时间步长$\Delta {t_{\text{F}}} = 5 \times {10^{ - 5}}\;{\text{s}}$ , 参考声速${c_0} = 1500\;{\text{m/s}}$ . 为与文献工况保持一致, 防止水流溢出, 水体高于坝体部分右侧添加了固壁边界; 计算中未考虑上游水渗入裂缝产生的缝面水压.图14给出了Koyna坝在水头超载作用下的裂缝扩展路径. 可以看出, 水力裂缝首先在初始裂缝顶点萌生, 并沿水平方向扩展约3.2 m. 随后, 裂缝开始向右下方扩展, 此时裂缝与水平方向夹角为33°. 本文所得裂缝扩展路径与文献[36]结果对比如图15所示, 可以看出, 两种方法所得裂缝路径基本一致, 表明本文提出的多分辨率PD−SPH混合方法能够较好地再现流体作用下的结构断裂过程.
3.4 水流冲击作用下的混凝土板崩塌问题
通过上述验证, 本节将提出的方法应用于工程常见的水流冲击作用下混凝土板崩塌这一复杂流−固耦合结构破坏问题模拟. 初始条件如图16所示, 水体经由宽为
$w = 0.33\;{\text{m}}$ 的水流入口并以$5\;{\text{m/s}}$ 的速度流入; 混凝土板长$L = 2\;{\text{m}}$ , 高$H = 0.3\;{\text{m}}$ , 其左右两端顶点受到固定约束, 中部含一条长$0.1\;{\text{m}}$ 、宽$0.015\;{\text{m}}$ 的竖向裂缝, 板的具体材料参数为: 密度${\rho _{\text{S}}}{\text{ = 2}}400\;{\text{kg/}}{{\text{m}}^3}$ , 弹性模量$E{\text{ = 35}}\;{\text{GPa}}$ , 泊松比$\nu = 0.2$ . 模拟中, 令PD粒子与SPH粒子的初始间距分别为$d{x_{\text{S}}} = 0.005\;{\text{m}}$ 和$d{x_{\text{F}}} = 0.01\;{\text{m}}$ , 两相计算时间步长为$\Delta {t_{\text{S}}} = 1 \times {10^{ - 6}}{\text{ s}}$ ,$\Delta {t_{\text{F}}} = {\text{5}} \times {10^{ - {\text{6}}}}{\text{ s}}$ ; 水体密度${\rho _{\text{F}}}{\text{ = }}1000 {\text{kg/}}{{\text{m}}^3}$ , 参考声速${c_0} = 60\;{\text{m/s}}$ .![]()
图 16 水流冲击混凝土板模型(单位: m)Figure 16. Model of water impacting on a concrete slab (unit: m)图17为水流冲击作用下混凝土板崩塌过程, 图18则给出了混凝土板上表面中部A点(如图16所示)的竖向位移时程曲线. 在
$t = 0.12\;{\text{s}}$ 时, 部分水体自水流出口流出, 并在重力作用下加速运动;$t = 0.212\;{\text{s}}$ 时, 水流开始冲击混凝土板上部;$t = 0.215\;{\text{s}}$ 时, 初始裂缝开始扩展, 此时A点的竖向位移较小, 仅为$ - 9.18 \times {10^{ - 5}}{\text{ m}}$ . 随着时间推移, 裂缝逐渐向上扩展, 并在$t = 0.222\;{\text{s}}$ 时, 扩展至混凝土板上表面形成贯穿裂缝, A点位移也于此刻产生巨大变化.$t = 0.31\;{\text{s}}$ 时, 左右两块混凝土板在自重及水流冲击作用下分别绕固定点旋转, 使得裂缝左右表面呈倒V形, 此时A点的竖向位移为$ - 0.11\;{\text{m}}$ ; 水流则因混凝土板的阻挡作用形成横向射流. 最终, 水流完全冲开了两块独立运动的混凝土板, 并自由向下运动; A点在$t = 0.45\;{\text{s}}$ 时的竖向位移为$ - 0.57\;{\text{m}}$ . 本文所得混凝土板断裂过程与文献[37]一致, 表明本文提出的多分辨率PD−SPH混合方法适用于复杂FSI问题求解, 能够自然实现流体运动与结构变形破坏全过程连续仿真.![]()
图 17 水流冲击作用下混凝土板崩塌过程Figure 17. Collapse process of concrete slab subjected to impact of water flow4. 小 结
本文针对流−固耦合破坏问题提出一种多时间和空间分辨率的PD−SPH混合方法, 并将其应用于典型问题模拟, 得到以下结论.
(1)该混合方法利用具有不同初始间距的粒子离散固体域和流体域, 采用不同的时间步长对两相的控制方程进行时间积分. 通过将PD粒子设置为具有与流体粒子相同光滑长度的虚粒子处理流−固界面, 可使得界面处满足动力学和运动学条件, 且能够保证系统动量守恒.
(2)通过液柱静压力作用下铝板变形和溃坝流体冲破弹性闸门两个经典考题模拟, 结果表明: 多分辨率PD−SPH混合方法具有较高的计算精度和较好的能量守恒特性, 针对不同问题合理选择流−固分辨率比以及时间步长比可有效提高计算效率.
(3)方法应用于经典Koyna重力坝水力劈裂和水流冲击作用下混凝土板的崩塌问题, 得到了固体结构的运动、变形与断裂过程以及固体结构影响下流体的运动状态. 所得结果验证了多分辨率PD−SPH混合方法在分析流−固耦合破坏问题方面的适用性. 在此基础上, 后续可扩展至三维流−固耦合破坏问题研究. 本文工作希望为分析工程中同时涉及移动自由表面、复杂流动、流−固耦合导致结构大位移和损伤破坏问题的高效数值仿真提供参考.
-
![]()
图 3 多分辨率PD-SPH耦合计算流程
Figure 3. Flow chart for multi-resolution PD-SPH coupling calculation
![]()
图 4 液柱静压力下铝板变形的初始条件(单位: m)
Figure 4. Initial condition for the aluminum plate under hydrostatic pressure (unit: m)
![]()
图 6 t = 5 s时工况III的流体压力场与板的挠度场
Figure 6. Fluid pressure and plate deflection field in case-III at t = 5 s
![]()
图 7 FSI系统归一化总能量时间历程
Figure 7. Time histories for normalized total energy of FSI system
![]()
图 9 溃坝流体冲破弹性闸门的初始条件示意图(单位: mm)
Figure 9. Schematic sketch for initial conditions in the dam-break flow through an elastic gate (unit: mm)
![]()
图 10 不同时刻弹性闸门变形和自由液面演变
Figure 10. Elastic gate deformation and change of fluid free surface
![]()
图 11 弹性闸门自由端水平位移时程
Figure 11. Time histories of horizontal displacements at the free end of elastic gate
![]()
图 16 水流冲击混凝土板模型(单位: m)
Figure 16. Model of water impacting on a concrete slab (unit: m)
![]()
图 17 水流冲击作用下混凝土板崩塌过程
Figure 17. Collapse process of concrete slab subjected to impact of water flow
表 1 液柱静压力下铝板变形问题工况布置
Table 1 The setup of different cases for aluminum plate deformation with hydrostatic pressure
Cases $\dfrac{{d{x_{\text{F}}}}}{{d{x_{\text{S}}}}}$ $d{x_{\text{S}}}{\text{/m}}$ $\dfrac{{\Delta {t_{\text{F}}}}}{{\Delta {t_{\text{S}}}}}$ $ \Delta {t_{\text{S}}}{\text{/s}} $ I 1.0 6.25 × 10−3 5.0 1.0 × 10−6 II 2.0 6.25 × 10−3 10.0 1.0 × 10−6 III 4.0 6.25 × 10−3 20.0 1.0 × 10−6 
下载: 导出CSV
表 2 溃坝流体冲破弹性闸门问题工况布置
Table 2 The setup of different cases for dam-break flow through an elastic gate
Cases $\dfrac{{d{x_{\text{F}}}}}{{d{x_{\text{S}}}}}$ $d{x_{\text{S}}}{\text{/m}}$ $\dfrac{{\Delta {t_{\text{F}}}}}{{\Delta {t_{\text{S}}}}}$ $ \Delta {t_{\text{S}}}{\text{/s}} $ I 1.0 1.0 × 10−3 1.0 5.0 × 10−6 II 2.0 1.0 × 10−3 1.0 5.0 × 10−6 III 2.0 1.0 × 10−3 4.0 5.0 × 10−6
下载: 导出CSV
-
[1] Caraeni D, Casseau V, Habashi WG. Fluid-structure interaction: Extended-FEM approach to solidification. Finite Elements in Analysis and Design, 2020, 177: 103425 doi: 10.1016/j.finel.2020.103425
[2] 汪春辉, 王嘉安, 王超等. 基于S-ALE方法的圆柱体垂直出水破冰研究. 力学学报, 2021, 53(11): 3110-3123 (Wang Chunhui, Wang Jiaan, Wang Chao, et al. Research on vertical movement of cylindrical structure out of water and breaking through ice layer based on S-ALE method. Chinese Journal of Theoretical and Applied Mechanics, 2021, 53(11): 3110-3123 (in Chinese) doi: 10.6052/0459-1879-21-217 [3] Martínez-Ferrer PJ, Qian L, Ma Z, et al. An efficient finite-volume method to study the interaction of two-phase fluid flows with elastic structures. Journal of Fluids and Structures, 2018, 83: 54-71 doi: 10.1016/j.jfluidstructs.2018.08.019
[4] Aus Der Wiesche S. Sloshing dynamics of a viscous liquid in a spinning horizontal cylindrical tank. Aerospace Science and Technology, 2008, 12(6): 448-456 doi: 10.1016/j.ast.2007.10.013
[5] Höllbacher S, Wittum G. Rotational test spaces for a fully-implicit FVM and FEM for the DNS of fluid-particle interaction. Journal of Computational Physics, 2019, 393: 186-213 doi: 10.1016/j.jcp.2019.05.004
[6] Liao KP, Hu CH, Sueyoshi M. Free surface flow impacting on an elastic structure: Experiment versus numerical simulation. Applied Ocean Research, 2015, 50: 192-208 doi: 10.1016/j.apor.2015.02.002
[7] Silling SA, Epton M, Weckner O, et al. Peridynamic states and constitutive modeling. Journal of Elasticity, 2007, 88(2): 151-184 doi: 10.1007/s10659-007-9125-1
[8] 黄丹, 章青, 乔丕忠等. 近场动力学方法及其应用. 力学进展, 2010, 40(4): 448-459 (Huang Dan, Zhang Qing, Qiao Pizhong, et al. A review on peridynamics (PD) method and its applications. Advances in Mechanics, 2010, 40(4): 448-459 (in Chinese) doi: 10.6052/1000-0992-2010-4-J2010-002 [9] 乔丕忠, 张勇, 张恒等. 近场动力学研究进展. 力学季刊, 2017, 38(1): 1-13 (Qiao Pizhong, Zhang Yong, Zhang Heng, et al. A review on advances in peridynamics. Chinese Quarterly of Mechanics, 2017, 38(1): 1-13 (in Chinese) [10] Liu MB, Zhang ZL. Smoothed particle hydrodynamics (SPH) for modeling fluid-structure interactions. Science China: Physics, Mechanics and Astronomy, 2019, 62(8): 984701 doi: 10.1007/s11433-018-9357-0
[11] Niazi S, Chen ZG, Bobaru F. Crack nucleation in brittle and quasi-brittle materials: A peridynamic analysis. Theoretical and Applied Fracture Mechanics, 2021, 112: 102855 doi: 10.1016/j.tafmec.2020.102855
[12] 王涵, 黄丹, 徐业鹏等. 非常规态型近场动力学热黏塑性模型及其应用. 力学学报, 2018, 50(4): 810-819 (Wang Han, Huang Dan, Xu Yepeng, et al. Non-ordinary state-based peridynamic thermal-viscoplastic model and its application. Chinese Journal of Theoretical and Applied Mechanics, 2018, 50(4): 810-819 (in Chinese) doi: 10.6052/0459-1879-18-113 [13] Wu LW, Huang D, Bobaru F. A reformulated rate-dependent visco-elastic model for dynamic deformation and fracture of PMMA with peridynamics. International Journal of Impact Engineering, 2021, 149: 103791 doi: 10.1016/j.ijimpeng.2020.103791
[14] Zhang YB, Huang D, Cai Z, et al. An extended ordinary state-based peridynamic approach for modelling hydraulic fracturing. Engineering Fracture Mechanics, 2020, 234: 107086 doi: 10.1016/j.engfracmech.2020.107086
[15] Zhou XP, Wang YT, Shou YD. Hydromechanical bond-based peridynamic model for pressurized and fluid-driven fracturing processes in fissured porous rocks. International Journal of Rock Mechanics and Mining Sciences, 2020, 132: 104383 doi: 10.1016/j.ijrmms.2020.104383
[16] Ren HF, Zhang CP, Zhao X. Numerical simulations on the fracture of a sea ice floe induced by waves. Applied Ocean Research, 2021, 108: 102527 doi: 10.1016/j.apor.2021.102527
[17] Liu RW, Yan JL, Li SF. Modeling and simulation of ice-water interactions by coupling peridynamics with updated Lagrangian particle hydrodynamics. Computational Particle Mechanics, 2020, 7(2): 241-255 doi: 10.1007/s40571-019-00268-7
[18] Gao Y, Oterkus S. Fluid-elastic structure interaction simulation by using ordinary state-based peridynamics and peridynamic differential operator. Engineering Analysis with Boundary Elements, 2020, 121: 126-142 doi: 10.1016/j.enganabound.2020.09.012
[19] Oger G, Doring M, Alessandrini B, et al. Two-dimensional SPH simulations of wedge water entries. Journal of Computational Physics, 2006, 213(2): 803-822 doi: 10.1016/j.jcp.2005.09.004
[20] Sun PN, Le Touzé D, Oger G, et al. An accurate FSI-SPH modeling of challenging fluid-structure interaction problems in two and three dimensions. Ocean Engineering, 2021, 221: 108552 doi: 10.1016/j.oceaneng.2020.108552
[21] Khayyer A, Gotoh H, Falahaty H, et al. An enhanced ISPH-SPH coupled method for simulation of incompressible fluid-elastic structure interactions. Computer Physics Communications, 2018, 232: 139-164 doi: 10.1016/j.cpc.2018.05.012
[22] Fan HF, Bergel GL, Li SF. A hybrid peridynamics-SPH simulation of soil fragmentation by blast loads of buried explosive. International Journal of Impact Engineering, 2016, 87: 14-27 doi: 10.1016/j.ijimpeng.2015.08.006
[23] Fan HF, Li SF. A peridynamics-SPH modeling and simulation of blast fragmentation of soil under buried explosive loads. Computer Methods in Applied Mechanics and Engineering, 2017, 318: 349-381 doi: 10.1016/j.cma.2017.01.026
[24] Sun WK, Zhang LW, Liew KM. A smoothed particle hydrodynamics-peridynamics coupling strategy for modeling fluid-structure interaction problems. Computer Methods in Applied Mechanics and Engineering, 2020, 371: 113298 doi: 10.1016/j.cma.2020.113298
[25] 姚学昊, 黄丹. 流固耦合问题的PD-SPH建模与分析. 工程力学, doi: 10.6052/j.issn.1000-4750.2021.06.0418 Yao Xuehao, Huang Dan. PD-SPH modelling and analysis for fluid-strucutre interaction problems. Engineering Mechanics, doi: 10.6052/j.issn.1000-4750.2021.06.0418 (in Chinese)
[26] Rahimi MN, Kolukisa DC, Yildiz M, et al. A generalized hybrid smoothed particle hydrodynamics-peridynamics algorithm with a novel Lagrangian mapping for solution and failure analysis of fluid-structure interaction problems. Computer Methods in Applied Mechanics and Engineering, 2022, 389: 114370 doi: 10.1016/j.cma.2021.114370
[27] Le QV, Bobaru F. Surface corrections for peridynamic models in elasticity and fracture. Computational Mechanics, 2018, 61(4): 499-518 doi: 10.1007/s00466-017-1469-1
[28] Zhang AM, Sun PN, Ming FR, et al. Smoothed particle hydrodynamics and its applications in fluid-structure interactions. Journal of Hydrodynamics, 2017, 29(2): 187-216 doi: 10.1016/S1001-6058(16)60730-8
[29] Monaghan JJ. Smoothed particle hydrodynamics. Reports on Progress in Physics, 2005, 68(8): 1703-1759 doi: 10.1088/0034-4885/68/8/R01
[30] 陈飞国, 葛蔚. 多相流动的光滑粒子流体动力学方法研究综述. 力学学报, 2021, 53(9): 2357-2373 (Chen Feiguo, Ge Wei. A review of smoothed particle hydrodynamics family methods for multiphase flow. Chinese Journal of Theoretical and Applied Mechanics, 2021, 53(9): 2357-2373 (in Chinese) doi: 10.6052/0459-1879-21-270 [31] Warren TL, Silling SA, Askari A, et al. A non-ordinary state-based peridynamic method to model solid material deformation and fracture. International Journal of Solids and Structures, 2009, 46(5): 1186-1195 doi: 10.1016/j.ijsolstr.2008.10.029
[32] Peng YX, Zhang AM, Wang SP. Coupling of WCSPH and RKPM for the simulation of incompressible fluid-structure interactions. Journal of Fluids and Structures, 2021, 102: 103254 doi: 10.1016/j.jfluidstructs.2021.103254
[33] Fourey G, Hermange C, Le Touzé D, et al. An efficient FSI coupling strategy between smoothed particle hydrodynamics and finite element methods. Computer Physics Communications, 2017, 217: 66-81 doi: 10.1016/j.cpc.2017.04.005
[34] Zhang C, Rezavand M, Hu XY. A multi-resolution SPH method for fluid-structure interactions. Journal of Computational Physics, 2021, 429: 110028 doi: 10.1016/j.jcp.2020.110028
[35] Antoci C, Gallati M, Sibilla S. Numerical simulation of fluid-structure interaction by SPH. Computers and Structures, 2007, 85(11-14): 879-890 doi: 10.1016/j.compstruc.2007.01.002
[36] 沙莎. 重力坝水力劈裂的数值模拟与坝踵真实应力性态研究. [博士论文]. 北京: 清华大学, 2017 Sha Sha. Numerical simulation of hydraulic fracturing of gravity dams and study on the real stress state of dam heel. [PhD Thesis]. Beijing: Tsinghua University, 2017 (in Chinese)
[37] Cornejo A, Franci A, Zárate F, et al. A fully Lagrangian formulation for fluid-structure interaction problems with free-surface flows and fracturing solids. Computers and Structures, 2021, 250: 106532 doi: 10.1016/j.compstruc.2021.106532
-
期刊类型引用(3)
1. 陆炎洲，李志远，黄丹，姚学昊. 基于近场动力学算子方法的各向异性板自由振动分析. 固体力学学报. 2024(03): 341-351 . 
百度学术
2. 许晓阳，赵雨婷，李家宇，余鹏. 非等温黏弹性复杂流动的改进SPH方法模拟. 力学学报. 2023(05): 1099-1112 . 本站查看
3. 李志远，黄丹，Timon Rabczuk. 近场动力学算子方法. 力学学报. 2023(07): 1593-1603 . 本站查看
其他类型引用(2)
计量
- 文章访问数: 1075
- HTML全文浏览量: 361
- PDF下载量: 203
- 被引次数: 5
