三维非等温黏弹性液滴撞击固壁问题的光滑粒子流体动力学方法模拟

李鸿杰; 许晓阳

doi:10.12404/j.issn.1671-1815.2504853

科学技术与工程 ›› 2026, Vol. 26 ›› Issue (19) :8097-8107 DOI: 10.12404/j.issn.1671-1815.2504853

力学

三维非等温黏弹性液滴撞击固壁问题的光滑粒子流体动力学方法模拟

李鸿杰 ,
许晓阳 ^*

作者信息 +

Numerical Simulation of 3D Non-isothermal Viscoelastic Droplet Impact onto Solid Wall with Smoothed Particle Hydrodynamics Method

LI Hong-jie ,
XU Xiao-yang ^*

Author information +

文章历史 +

摘要

基于光滑粒子流体动力学(smoothed particle hydrodynamics, SPH)方法对三维非等温黏弹性液滴撞击固壁问题进行数值模拟,其中流体的弹性应力通过Phan-Thien-Tanner本构模型进行表征。为解决液滴张力不稳定性问题,引入人工应力项;为提升三维模拟效率,采用GPU并行加速技术。通过模拟二维非等温泊肃叶流动和三维等温液滴撞击固壁问题,并和其他数值方法结果进行比较,验证了所采用SPH方法的有效性。在此基础上,模拟了三维非等温液滴撞击固壁过程,评估了方法的收敛性与并行效率,分析了液滴初始温度、温度依赖系数以及其他物理参数对液滴铺展行为的影响。数值模拟结果表明,该SPH方法能够准确高效地模拟三维非等温黏弹性液滴撞击固壁问题,为相关复杂问题的模拟提供了参考。

Abstract

A three-dimensional non-isothermal viscoelastic droplet impacting onto solid wall was numerically simulated using the smoothed particle hydrodynamics (SPH) method, in which the elastic stress of the fluid was characterized by the Phan-Thien-Tanner (PTT) constitutive model. To suppress tensile instability, an artificial stress term was introduced. GPU parallel computing was employed to improve the computational efficiency of the 3D simulations. The accuracy and validity of the proposed SPH method were verified through simulations of two benchmark cases: a two-dimensional non-isothermal Poiseuille flow and a three-dimensional isothermal droplet impact on a solid wall, with results compared against those from other numerical methods. Based on this validation, the 3D non-isothermal droplet impact process was further investigated. The convergence and parallel performance of the method were assessed, and the effects of initial droplet temperature, temperature-dependent parameters, and other physical properties on the spreading dynamics were analyzed. The numerical results demonstrate that the proposed SPH approach can accurately and efficiently simulate complex 3D non-isothermal viscoelastic droplet impact phenomena, offering valuable insights for related applications.

Graphical abstract

关键词

光滑粒子流体动力学 / PTT模型 / 黏弹性 / 三维流体 / 非等温

Key words

smoothed particle hydrodynamics / PTT model / viscoelastic / three-dimensional fluid / non-isothermal

引用本文

引用格式 ▾

李鸿杰, 许晓阳. 三维非等温黏弹性液滴撞击固壁问题的光滑粒子流体动力学方法模拟[J]. 科学技术与工程, 2026, 26(19): 8097-8107 DOI:10.12404/j.issn.1671-1815.2504853

LI Hong-jie, XU Xiao-yang. Numerical Simulation of 3D Non-isothermal Viscoelastic Droplet Impact onto Solid Wall with Smoothed Particle Hydrodynamics Method[J]. Science Technology and Engineering, 2026, 26(19): 8097-8107 DOI:10.12404/j.issn.1671-1815.2504853

登录浏览全文

4963

注册一个新账户忘记密码

黏弹性液滴撞击固壁广泛存在于自然界与工程应用中,涉及复杂的自由表面流动。在工业领域,其典型应用包括3D打印、喷墨打印、注塑成型等。液滴的动力学行为直接决定了产品的成型精度、表面质量与功能特性^[1]。在3D打印领域,高分子聚合物常以熔融液滴的形式逐层沉积。液滴撞击基板后的形态稳定性与铺展宽度共同决定了打印成型件的几何精度与力学性能^[2]。喷墨打印技术所采用的油墨多为具有黏弹性的非牛顿流体,液滴撞击纸张后的铺展范围直接决定了印刷图案的分辨率与清晰度^[3]。而在注塑成型工艺中,熔体以液滴形式撞击模具壁面,其黏弹性主导的流动变形会影响制品内部应力分布与使用寿命^[4]。同时,上述实际工程过程往往处于非等温环境,壁面的加热或冷却会使流体形成显著温度梯度,进而影响其流变特性与动力学行为^[5]。因此,准确而高效地预测非等温环境下黏弹性流体的流动特性,不仅对理解复杂流体界面动力学具有重要学术意义,也为优化制造工艺、提升产品性能提供理论支撑。

目前,已有研究采用基于网格的数值方法模拟该问题,如有限体积法(Rath等^[6]、庄昕等^[7])、有限差分法(Sharma等^[8])和有限元法(Venkatesan等^[9]、李烨波等^[10])。以上方法发展较为成熟,但在处理涉及自由面剧烈变形问题时,通常需要引入复杂的界面追踪技术[如VOF(volume of fluid)、level set)],其实现流程相对繁琐,且在大变形情况下可能面临网格畸变带来的数值不稳定问题。

光滑粒子流体动力学(smoothed particle hydrodynamics,SPH)方法是一种拉格朗日性质的无网格方法^[11]。与基于网格的数值方法相比,SPH方法在计算过程中不依赖任何网格,能够有效避免因网格变形带来的数值不稳定问题。在模拟自由表面流动以及大变形问题时具有天然优势,无需引入额外的界面捕捉技术。此外,基于粒子的离散化特性,SPH方法更加适合并行化,每个粒子可独立携带多种物理量,便于构建多物理场耦合模型而不依赖复杂的全局网格拓扑。得益于上述特点,SPH方法已在众多科学与工程领域得到广泛应用。在海洋工程领域,Ramezanzadeh等^[12]采用弱可压缩SPH方法模拟了浮体结构在海洋环境下的运动响应,为波浪能装置与离岸结构优化提供了可靠数值工具。在岩土力学领域,张高峰等^[13]提出FEM(finite element method)-SPH耦合算法,揭示了岸坡与桥梁桩基的相互作用致损机理,为库区工程结构的安全防护提供了理论依据。在爆炸物理领域,Kim等^[14]模拟了核反应堆腔室外部蒸汽爆炸的冲击波传播与结构效应,为核反应堆安全评估提供了关键技术支持。在流固耦合领域,Guo等^[15]建立的SPH多相流模型有效缓解了大密度比多相流中的界面非物理间隙问题,并应用于楔形体入水等多相流固耦合场景的精准模拟。

对于黏弹性自由表面流动的SPH模拟,Jiang等^[16]提出一种修正的SPH方法,模拟了Oldroyd-B黏弹性两相流中液滴上升、下落及剪切流中的变形、过冲与破裂等现象。Moinfar等^[17]采用SPH研究了简单剪切流中Giesekus流体包围下液滴的变形行为。Sun等^[18]引入Oldroyd-B本构模型,并结合SPH算法实现了黏弹性血液涌出并扩散的过程。许晓阳等^[19]提出一种改进的SPH方法对液滴撞击固壁问题进行模拟,并与Basilisk开源程序结果进行对比验证了方法的有效性。Li等^[20]将SPH与有限粒子法结合,模拟了Oldroyd-B黏弹性液滴撞击固壁过程,并进一步拓展至双液滴连续撞击的复杂流动场景。上述研究主要基于Oldroyd-B和Giesekus等模型,且大多基于二维等温条件下进行模拟。

目前,SPH模拟三维非等温PTT黏弹性液滴撞击固壁问题的文献较为少见。PTT模型基于聚合物网络理论,具有更强的物理基础与模型灵活性,能够更准确地描述聚合物熔体及浓缩溶液等复杂流体的非线性流变行为^[21]。开展三维数值模拟有助于完整捕捉流动结构的空间演化特征,从而更真实地反映黏弹性流体在实际过程中的复杂动力学行为。此外,引入温度场进一步提升了模型的物理真实性,温度变化不仅直接影响流体的黏度、松弛时间等关键流变参数,同时也参与应力本构关系的演化,从而显著影响液滴撞击后的铺展与回缩特性^[22]。因此,建立三维非等温条件下的PTT黏弹性流动模型,对于准确预测黏弹性流体在实际加工中的非线性响应特性具有重要的理论指导意义。

基于上述背景,在文献[19]的工作基础上,采用SPH方法对三维非等温条件下的PTT液滴撞击固壁过程进行数值模拟研究。针对黏弹性自由表面流动中出现的张力不稳定现象,引入人工应力方法加以抑制;为解决三维问题计算耗时长的问题,采用GPU并行加速技术。通过模拟二维非等温泊肃叶流动并与有限体积法结果对比,以及三维等温液滴撞击固壁过程并与有限差分法结果对比,验证所用SPH方法在处理三维非等温黏弹性流动问题中的有效性。通过改变初始粒子间距进行模拟,评估方法的数值收敛性,同时对比不同粒子总数下CPU与GPU平台的计算性能。此外,还深入分析了温度依赖系数、Weissenberg数等对三维非等温黏弹性液滴流动行为的影响,揭示了温度效应与黏弹性效应对液滴撞击固壁动力学行为的耦合作用机制,为复杂黏弹性流体的撞击、铺展及传热过程提供了理论支撑与数值方法参考,具有重要的理论研究价值与工程应用前景。

1 控制方程

在拉格朗日坐标体系下,非等温黏弹性流动控制方程组可以表示为

(1)$\frac{\mathrm{d} \rho}{\mathrm{~d} t}=-\rho \nabla \cdot \boldsymbol{u}$

(2)$\frac{\mathrm{d} \boldsymbol{u}}{\mathrm{~d} t}=\frac{1}{\rho} \nabla \cdot \boldsymbol{\sigma}+\boldsymbol{F}$

(3)$\rho c_{p} \frac{\mathrm{~d} T}{\mathrm{~d} t}=\nabla \cdot(\kappa \nabla T)$

式中:t为时间;ρ为流体密度;u为流体流速;σ为总应力张量;F为流体所受到的外力;c_p和κ分别为比热容和导热系数;T为温度。

式(2)中的总应力张量σ由各向同性压力p、黏性应力2η_sd和弹性应力张量τ组成,即

(4)σ=-pδ+2η_sd+τ

式(4)中:δ为Kronecker函数;η_s为溶剂黏度;d为形变率张量,计算公式为

(5)$\boldsymbol{d}=\frac{1}{2}\left[\nabla \boldsymbol{u}+(\nabla \boldsymbol{u})^{\mathrm{T}}\right]$

1.1 PTT本构方程

为了求解动量守恒方程[式(2)],引入PTT(Phan-Thien-Tanner)本构方程^[23],即

(6)$f[\operatorname{tr}(\boldsymbol{\tau})] \boldsymbol{\tau}+\lambda_{1} \stackrel{\nabla}{\boldsymbol{\tau}}=2 \eta_{\mathrm{p}} \boldsymbol{d}$

式(6)中:λ₁为流体的松弛时间;η_p为聚合物黏度;$ \stackrel{\nabla}{\boldsymbol{\tau}}$为上随体导数,定义为

(7)$\stackrel{\nabla}{\boldsymbol{\tau}}=\frac{\mathrm{d} \boldsymbol{\tau}}{\mathrm{~d} t}-\boldsymbol{\tau} \nabla \boldsymbol{u}-(\nabla \boldsymbol{u})^{\mathrm{T}} \boldsymbol{\tau}$

f[tr(τ)]为应力系数函数,采用线性形式表示为

(8)f[tr(τ)]=1+$\frac{{\mathit{\lambda }}_{1}\mathit{\epsilon }}{{\mathit{\eta }}_{\mathrm{p}}}$tr(τ)

式(8)中:ε为PTT本构模型中的拉伸参数。溶剂黏度η_s与聚合物黏度η_p满足关系式

(9)η=η_s+η_p

式(9)中:η为流体的总黏度。

采用Reynolds指数模型^[24]表征流体的热流变特性,其中流体的松弛时间和黏度随温度的变化关系可表示为

(10)

$\left\{\begin{array}{l}{\mathit{\lambda }}_{1}\left(\mathit{T}\right)={\mathit{\lambda }}_{1}{\mathrm{e}}^{-\mathit{\varphi }(\mathit{T}-{\mathit{T}}_{0})}\\ {\mathit{\eta }}_{\mathrm{s}}\left(\mathit{T}\right)={\mathit{\eta }}_{\mathrm{s}}{\mathit{e}}^{-\mathit{\varphi }(\mathit{T}-{\mathit{T}}_{0})}\\ {\mathit{\eta }}_{\mathrm{p}}\left(\mathit{T}\right)={\mathit{\eta }}_{\mathrm{p}}{\mathit{e}}^{-\mathit{\varphi }(\mathit{T}-{\mathit{T}}_{0})}\end{array}\right.$

式(10)中:ϕ为温度依赖系数;T₀为基准温度。

为了更好地描述黏弹性流动复杂的流变特性,引入Weissenberg数(Wi)和溶剂黏度比β,即

(11)Wi=$\frac{{\mathit{\lambda }}_{1}\mathit{V}}{\mathit{L}}$

(12)β=$\frac{{\mathit{\eta }}_{\mathrm{s}}}{\mathit{\eta }}$

式中:Wi为黏弹性流体的弹性力与黏性力之比; β为牛顿溶剂黏度与流体总黏度的比值;V为流体的特征速率;L为流体的特征长度。

1.2 状态方程

本文中使用的是弱可压缩SPH方法,根据密度与压力的关系来求解流体的压力p,采取的状态方程为

(13)p=c²(ρ-ρ₀)

式(13)中:ρ₀为流体的初始密度;c为流体中的声速。在SPH模拟中,为了使密度的变化维持在1%以内,人工声速的取值需满足c≥10U_max,其中U_max为流场最大流速^[25]。

2 光滑粒子流体动力学(SPH)方法

2.1 SPH离散

SPH方法是一种无网格粒子法,该方法将流体域离散成一系列均匀分布且数量有限的粒子,随后通过光滑核函数近似求解偏微分方程组,得到SPH数值解。SPH的离散分为积分近似和粒子近似两大步骤。粒子的密度、速度、压力以及温度等物理量通过积分近似定义。变量f(x)的插值近似以积分形式表示为

(14)<f(x)>=∫_Ωf(x')W(x-x',h)dx'

式(14)中:<f(x)>为函数f(x)的近似估计;W为光滑函数;Ω为光滑函数W的作用域;x表示粒子位置;h为光滑长度,决定了光滑函数的作用域。粒子近似是将作用域内的粒子物理量叠加求和。对式(14)进行粒子近似可得

(15)<f(x_i)>=$\stackrel{\mathit{N}}{\sum _{\mathit{j}=1}}\frac{{\mathit{m}}_{\mathit{j}}}{{\mathit{\rho }}_{\mathit{j}}}$f(x_j)W_ij

式(15)中:i和j代表粒子的序号,粒子j在粒子i的支持域内;m_j为粒子质量;W_ij为粒子i与j之间的光滑核函数值。

采用二阶可导的五次样条函数,从而保证SPH数值解的精度与稳定性,其表达式为

(16)W(r,h)=W₀

$\left\{\begin{array}{ll}{(3-\mathit{q})}^{5}-6{(2-\mathit{q})}^{5}+15{(1-\mathit{q})}^{5},& 0\le \mathit{q}<1\\ {(3-\mathit{q})}^{5}-6{(2-\mathit{q})}^{5},& 1\le \mathit{q}<2\\ {(3-\mathit{q})}^{5},& 2\le \mathit{q}<3\\ 0,& \mathit{q}\ge 3\end{array}\right.$

式(16)中:q为相对距离,q=r/h,r=|x-x'|/h为粒子间的相对距离;W₀为五次样条核函数的归一化系数,其取值与空间维数密切相关。

由于二维空间的支持域为圆形区域,而三维空间为球形区域,为保证核函数在各自支持域内的积分为$1\left[\int_{\Omega} W(r, h) \mathrm{d} r=1\right]$,需根据空间维数分别选取对应的归一化系数。在二维情况下,归一化系数为W₀=

$\frac{7}{478\mathrm{\pi }{\mathit{h}}^{2}}$

,在三维情况下为W₀=

$\frac{3}{359\mathrm{\pi }{\mathit{h}}^{3}}$

。

对质量、动量以及温度的守恒方程采取的SPH离散为

(17)$\left\{\begin{array}{l}\left(\frac{\mathrm{d} \rho}{\mathrm{~d} t}\right)_{i}=\sum_{j} m_{j}\left(\boldsymbol{u}_{i}-\boldsymbol{u}_{j}\right) \nabla_{i} W_{i j} \\\left(\frac{\mathrm{~d} \boldsymbol{u}}{\mathrm{~d} t}\right)_{i}=\sum_{j} m_{j}\left(\frac{p_{i}}{\rho_{i}^{2}}+\frac{p_{j}}{\rho_{j}^{2}}\right) \nabla_{i} W_{i j}+\boldsymbol{F} \\\left(\frac{\mathrm{d} T}{\mathrm{~d} t}\right)_{i}=\frac{\kappa}{c_{p}} \sum_{i} \frac{8 m_{j}}{\left(\rho_{i}+\rho_{j}\right)^{2}} \frac{\boldsymbol{r}_{i j} \nabla_{i} W_{i j}}{\left|\boldsymbol{r}_{i j}\right|^{2}}\left(T_{i}-T_{j}\right)\end{array}\right.$

式(17)中:r_ij为粒子i到j的位置矢量。

PTT本构方程的SPH离散主要依赖于速度偏导数的离散,引入速度梯度张量分量

${\mathit{k}}_{\mathit{i}}^{\mathit{\alpha }\mathit{\beta }}$

,即

(18)$k_{i}^{\alpha \beta}=\left(\frac{\nabla u^{\alpha}}{\nabla x^{\beta}}\right)_{i}=\sum_{j} \frac{m_{j}}{\rho_{j}}\left(u_{j}^{\alpha}-u_{i}^{\alpha}\right) \frac{\nabla W_{i j}}{\nabla x_{i}^{\beta}}$

式(18)中:u^α为速度矢量u的第α个分量;上标αβ表示对应张量在(α,β)方向上的分量。

基于此,PTT本构方程的SPH离散可表示为

(19)${\left(\frac{\mathrm{d}{\mathit{\tau }}^{\mathit{\alpha }\mathit{\beta }}}{\mathrm{d}\mathit{t}}\right)}_{\mathit{i}}$=${\mathit{k}}_{\mathit{i}}^{\mathit{\alpha }\mathit{\gamma }}{\mathit{\tau }}_{\mathit{i}}^{\mathit{\gamma }\mathit{\beta }}$+${\mathit{k}}_{\mathit{i}}^{\mathit{\beta }\mathit{\gamma }}{\mathit{\tau }}_{\mathit{i}}^{\mathit{\gamma }\mathit{\alpha }}$-$\frac{\mathit{f}\left[\mathrm{t}\mathrm{r}\right({\mathit{\tau }}_{\mathit{i}}\left)\right]}{{\mathit{\lambda }}_{1}}{\mathit{\tau }}_{\mathit{i}}^{\mathit{\alpha }\mathit{\beta }}$+$\frac{2{\mathit{\eta }}_{\mathrm{p}}}{{\mathit{\lambda }}_{1}}$d^αβ

式(19)中:γ为求和指标,对坐标方向自动求和以完成张量缩并;

${\mathit{k}}_{\mathit{i}}^{\mathit{\alpha }\mathit{\gamma }}{\mathit{\tau }}_{\mathit{i}}^{\mathrm{\gamma }\mathrm{\beta }}$

、

${\mathit{k}}_{\mathit{i}}^{\mathit{\beta }\mathit{\gamma }}{\mathit{\tau }}_{\mathit{i}}^{\mathrm{\gamma }\mathrm{\alpha }}$

分别对应速度梯度张量与弹性应力张量在不同方向分量上的耦合项;d^αβ为形变张量d的(α,β)分量。

2.2 人工应力

张力不稳定是传统SPH方法在固体力学模拟中的常见问题。对于黏弹性流体自由表面问题,由于弹性力的存在也会产生张力不稳定问题,即粒子会出现非物理的聚集,导致应力计算错误,降低计算的稳定性和精度。为防止粒子过于聚集或产生非物理拉伸现象,在动量守恒方程中耦合人工应力项^[26],该方法通过对相互作用的粒子之间施加虚拟的排斥力,其表达式为

(20)${\mathit{S}}_{\mathit{i}\mathit{j}}^{\mathit{\alpha }\mathit{\beta }}$=${\mathit{f}}_{\mathit{i}\mathit{j}}^{\mathit{n}}$(${\mathit{R}}_{\mathit{i}}^{\mathit{\alpha }\mathit{\beta }}$+${\mathit{R}}_{\mathit{j}}^{\mathit{\alpha }\mathit{\beta }}$)

(21)f_ij=$\frac{{\mathit{W}}_{\mathit{i}\mathit{j}}}{\mathit{W}(\mathrm{\Delta }\mathit{x},\mathit{h})}$

(22)n=$\frac{\mathit{W}(0,\mathit{h})}{\mathit{W}(\mathrm{\Delta }\mathit{x},\mathit{h})}$

式中:R^αβ为人工应力张量R的分量;Δx 粒子间的初始间距;f_ij为人工应力系数;n为指数型因子。

在二维空间中,人工应力张量可通过平面直角坐标系的旋转获得,然而在三维空间中,坐标系的旋转计算非常复杂。对此,采取简化的人工应力张量表达式为

(23)R^αβ=$\left\{\begin{array}{ll}\frac{-\mathit{b}{\mathit{\sigma }}^{\mathit{\alpha }\mathit{\beta }}}{{\mathit{\rho }}^{2}},& {\mathit{\sigma }}^{\mathit{\alpha }\mathit{\beta }}>0\\ 0,& {\mathit{\sigma }}^{\mathit{\alpha }\mathit{\beta }}\le 0\end{array}\right.$

式(23)中:b为介于0~1的常量;σ^αβ为应力张量的分量。

数值结果表明,该简化表达式能够在避免三维坐标系复杂旋转计算的同时,保留人工应力抑制张力不稳定性的核心功能,兼顾计算效率与结果准确性^[27]。添加人工应力后,动量守恒的SPH离散格式为

(24)$\left(\frac{\mathrm{d} \boldsymbol{u}}{\mathrm{~d} t}\right)_{i}=\sum_{j} m_{j}\left(\frac{p_{i}}{\rho_{i}^{2}}+\frac{p_{j}}{\rho_{j}^{2}}+S_{i j}^{\alpha \beta}\right) \nabla_{i} W_{i j}+\boldsymbol{F}$

2.3 边界处理

在SPH模拟中,边界处理是极为重要的一环,直接关系到模拟的精度与稳定性。当前应用较为广泛的两种方法是镜像粒子法^[28]与虚拟粒子法^[29]。镜像粒子法通过在边界的对称位置复制流体粒子,从而增强边界附近粒子的支持域,提升边界计算的准确性,但其在处理强剪切或高速冲击问题时容易出现粒子穿透边界的现象。虚拟粒子法则在边界外布置若干层虚拟粒子并赋予其特定物理属性,以避免边界截断对流体粒子支持域的影响。其缺点是获取虚拟粒子的速度、压力等物理量时,需依靠构造伪粒子的手段,利用shepard插值方法实现计算。

采用一种结合实际边界粒子与虚拟粒子的边界处理方法^[19]。如图1所示,首先在边界处按初始粒子间距Δx均匀排列一层边界粒子,并在边界外沿法向扩展三层虚拟粒子,从而防止支持域被截断,提升边界附近的计算准确性。在实现过程中,边界粒子被赋予固定的密度、温度等物理属性,并直接参与控制程求解。边界粒子的弹性应力与应力通过对其邻域内的流体粒子采用shepard插值方法进行计算,表达式为

(25)A_i=$\frac{\sum _{\mathit{j}}\frac{{\mathit{m}}_{\mathit{j}}}{{\mathit{\rho }}_{\mathit{j}}}{\mathit{A}}_{\mathit{j}}{\mathit{W}}_{\mathit{i}\mathit{j}}}{\sum _{\mathit{j}}\frac{{\mathit{m}}_{\mathit{j}}}{{\mathit{\rho }}_{\mathit{j}}}{\mathit{W}}_{\mathit{i}\mathit{j}}}$

式(25)中:A为弹性应力或压力;i为边界粒子;j为i邻域内的流体粒子。

每个虚拟粒子对应有唯一的边界粒子。为满足neumann边界条件,虚拟粒子的弹性应力与压力被设定为其对应边界粒子的值。同时,为实现无滑移边界条件,虚拟粒子的速度被设定为0。

3 基于GPU的并行计算

在SPH方法中,粒子间的相互作用仅依赖邻接粒子物理量,无需复杂耦合计算。在CPU程序中,通过循环依次更新每个粒子,随着粒子数量增加,计算耗时显著增长。由于粒子间无全局依赖性,SPH算法适合于GPU并行计算。对于SPH的三维非等温模拟,其粒子数量大,物理量多,计算量远超二维问题,亟需借助GPU实现高效计算。

如图2所示为SPH方法的GPU并行计算流程。初始化阶段由CPU完成粒子位置、密度、温度等物理量的设定,并分配CPU与GPU所需内存资源。随后,为减少主机与设备间的数据传输开销,将数据一次性传输至GPU。紧接着,在GPU上进行SPH的核心计算包括邻域搜索、物理量变化率求解与更新,结果仅在需要保留时回传至CPU。

在SPH方法中,粒子间的相互作用通常局限于核函数的支持域范围内,表现出强烈的局部性特征。因此,在GPU框架下,若能充分利用这一局部性,将显著提升计算效率。然而,GPU中的多级存储结构(如共享内存和全局内存)在访问速度与容量上存在显著差异^[30]。全局内存容量大但访问延迟高,共享内存访问速度快但容量有限,合理调度这两类内存资源对性能优化至关重要。

引入一种结合粒子支持域特性的GPU内存访问优化策略。如图3所示,左边展示了当前粒子、支持域内的邻接粒子以及支持域外的粒子;右边为GPU计算单元结构,其中每个线程块(block)内的多个线程(thread)可通过共享内存协同工作,并统一访问全局内存。该策略的核心思想是在每个线程块中,利用粒子间交互的局部性特征,将当前粒子支持域内的邻接粒子物理量数据批量加载至共享内存中,从而避免每个线程在计算中频繁从全局内存中读取邻接信息的繁琐。具体而言,首先在每个线程块中通过邻接搜索获取当前粒子所在支持域内的邻接粒子索引信息,设当前粒子i的位置为r_i,核函数的支持半径为h。支持域内的邻接粒子集合定义为

(26)Ω_i={j∣‖r_i-r_j‖≤h}

式(26)中:r_j为粒子j的位置;Ω_i为所有满足距离条件的邻接粒子索引集合。

将这些邻接粒子的物理属性(如位置、速度、应力等)预先批量加载至共享内存中;接着,各线程直接访问共享内存中已加载数据执行相互作用计算。共享内存具有较低访问延迟和线程间数据共享的优势,使得粒子间高频交互数据可在同一线程块内被高效复用。对于部分超出共享内存容量或跨线程块粒子之间的交互,则退回至访问全局内存。该策略实现了高频数据共享、低频数据分离的内存调度机制。一方面显著降低了内存访问延迟,提升了共享数据的利用率;另一方面缓解了全局内存带宽瓶颈,减少了冗余的数据加载。在保证SPH模拟精度的前提下,该策略能够显著提升GPU并行计算的整体性能。

4 数值算例

4.1 有效性验证

4.1.1 二维非等温泊肃叶流动

为了验证SPH方法求解非等温问题的有效性,首先求解二维非等温PTT泊肃叶流动,如图4所示为泊肃叶流动的几何模型,其由两块平行放置的无限长平板组成,板间水平间距为L。初始时刻,流体静止于平板之间,速度为0。流动开始后,流体在平行于x方向的外力F作用下开始流动。针对x轴方向的无限长板采用周期边界条件,计算区域取L_x×L_y=0.5 m×1 m,水平外力F=1.0 N,粒子密度ρ=1 000 kg/m³,无量纲参数取值Wi=1,β=0.1,PPT拉伸参数ε=0.25,其中特征长度L=1,特征速度V=0.625。对于温度条件,流体的初始温度为T_f=443 K,水平长板上下边界温度为T_w=313 K,比热容c_p=1,热导率κ=1,温度依赖系数ϕ=0.002。SPH的粒子初始间距设为Δx=0.012 5,其中流体粒子数量为3 239个,边界粒子为94个,边界外虚拟粒子为282个。五次样条核函数的光滑长度取h=1.5Δx。声速选择为c=10V。为了保持数值稳定性,使用时间步长为Δt=1.0×10^-5。

如图5(a)所示为中心轴上流体速度随时间的变化趋势。在外力与弹性应力的共同作用下,流体速度迅速上升,达到最大值之后出现下降,随后由于弹性产生衰减振荡,最终趋于稳定值。如图5(b)所示为t=0.5时刻流体温度沿y方向的分布情况,由于流体与边界发生热传导,边界附近的流体温度率先开始下降,而中心区域附近的流体温度变化相对滞后。对比有限体积法(finite volume method, FVM)^[31]结果,两者在速度演化及温度分布上吻合良好,验证了本文中所采用SPH方法在非等温黏弹性问题中的准确性。

4.1.2 三维等温液滴撞击固壁

为了进一步验证SPH方法求解三维流动问题的适用性,模拟了三维等温PTT液滴撞击固壁问题。如图6所示,液滴初始形状为半径为R=0.01 m的球体且处于无应力状态下,圆心距下方固壁面的高度为H=0.04 m,初始的下落速度为V=-1.0 m/s。自t>0时刻开始,在重力的驱动下液滴开始加速下落,接触固壁后发生产生明显的铺展形变。重力加速度g=-9.81 m/s²,声速取c=12.5 m/s。液滴的密度为ρ=1 000 kg/m³,松弛时间取λ=0.02 s,无量纲参数Wi=1,β=0.1,PTT拉伸参数取ε=0.1。为避免张力不稳定,人工应力系数取值为b=0.5。相比二维问题,由于三维空间的体积分布特性导致粒子密度需求更高,需要更小的粒子间距以确保足够的邻接粒子和计算精度,因此取粒子初始间距为Δx=0.000 2,对应的总粒子数为885 697个,其中流体粒子数为523 293个,固壁粒子数为90 601,固壁外虚粒子数为271 803,五次样条核函数的光滑长度取h=1.5Δx,使用时间步长Δt=1.0×10^-5。为了便于分析,采取无量纲时间t^*=tV/(2R)来记录液滴的变化过程。

如图7所示为液滴撞击过程中液滴宽度随时间的演化过程。当液滴撞击固壁后,首先会发生水平铺展;待达到最大铺展宽度后,由于弹性力的作用,液滴开始收缩回弹;而当收缩至最小宽度时,因能量耗散的影响,液滴将缓慢进行二次铺展。将本文中SPH方法的计算结果与Tomé等^[32]基于有限差分法(finite difference method, FDM)的计算结果进行对比,两者的变化曲线高度一致。在最大铺展阶段,两种方法所得最大铺展宽度分别约为0.034 7与0.034 6;而在最小收缩阶段,对应的最小宽度分别为0.030 3与0.029 9,二者在关键特征量上的差异极小,验证了SPH方法模拟三维黏弹性液滴撞击固壁问题中的准确性。

通过以上两个算例的模拟分析,本文所采用SPH方法能够准确捕捉黏弹性流体在流动过程中的速度演化特征、热传导行为以及形态变化,具有较高的物理一致性,验证了其在三维非等温黏弹性流动问题中的有效性。

4.2 三维非等温液滴撞击固壁

接下来,将三维等温PTT液滴撞击固壁问题扩展到非等温条件下进行模拟。温度参数设置如下:液滴温度T_f=443 K,冷固壁温度T_w=313 K。液滴比热容c_p=1 J/(kg·K),热导率κ=1 W/(m·K),温度依赖系数ϕ=0.001 K^-1。模拟中,其余物理参数的设定均与等温模拟相同。

如图8所示为液滴撞击固壁过程中六个不同时刻t^*=1.50, 1.75, 2.25, 2.70, 3.00, 5.00下的SPH温度分布粒子结果图。初始阶段(t^*=1.50和t^*=1.75),液滴由于重力加速度和初始速度的共同作用,液滴的高度逐渐减小,宽度持续扩展直至达到最大铺展宽度。由于液滴与固壁之间存在热传导,液滴底部温度开始降低。收缩阶段(t^*=2.25和t^*=2.75),当液滴铺展至最大宽度后由于弹性作用出现了收缩行为。随着时间推移,液滴与壁面之间持续的热交换使得液滴整体温度由下至上逐渐降低,底部区域的温度已接近壁面温度。二次铺展阶段(t^*=3.00和t^*=5.00),液滴表现出缓慢的二次铺展,这是弹性耗散与重力共同作用的结果。最终在t^*=5.00时,液滴温度场已趋于稳定,与冷壁面的温度场基本一致,表明系统热过程已近平衡状态。

为了探究温度对液滴流动状态的影响,对比了等温与非等温条件下液滴铺展宽度随时间的演化过程,结果如图9所示。从图9中可以看出,在初始阶段的铺展行为两者高度相似,表明温度在流动初期对液滴动力学特性的影响比较有限。随着模拟时间的推进,非等温液滴逐渐展现出更明显的收缩趋势,其最小收缩宽度(约0.031 1)明显小于等温情况(约0.031 3)。此外,在后续的二次铺展过程中,非等温液滴的铺展宽度始终略低于等温液滴。

4.2.1 SPH方法收敛性分析

为验证所使用SPH方法的数值收敛性,通过改变流体粒子的初始间距Δx进行对比试验。选取三组粒子间距分别为Δx=0.000 4,0.000 25和0.000 17,对应粒子总数约为15.65万、50.01万和142.54万个,与原始模拟所用间距Δx=0.000 2的结果进行比较。如图10所示为4种不同粒子间距下液滴铺展宽度与中心点温度随时间演化的对比结果。可以看出,尽管粒子间距存在一定差异,但各组模拟在液滴动力学行为与热传导特性方面的演化趋势基本一致,表明本文所采用的SPH方法模拟三维非等温黏弹性液滴撞击固壁问题具备良好的数值收敛性。

4.2.2 GPU性能分析

为评估SPH方法在GPU平台上的加速,选取不同粒子数的三维非等温液滴撞击固壁面算例,在CPU (Intel© Xeon© Gold 6145)与GPU (NVIDIA GeForce RTX 3090 Ti)平台上分别运行,记录其总计算时间,进而计算加速比。如表1所示为不同粒子总数下CPU串行方法与GPU链表法的运行时间比较结果。随着粒子总数从约15万增加至142万, GPU并行始终显著优于CPU串行方法,且在大规模粒子数量的优势愈加明显(最大加速比高达约104倍)。加速比随粒子总数增加而单调上升,表明GPU并行计算在处理三维非等温流动时能够展现出良好适应性和显著的效率优势。

4.2.3 温度的引入对液滴铺展宽度的影响

(1) 液滴初始温度T_f的影响。为了研究温度的引入对液滴铺展宽度的影响,如图11(a)所示为液滴初始温度T_f为443、743、1 043、1 343 K下,液滴铺展宽度随时间t^*的演化过程。随着T_f的升高,液滴的最大铺展宽度逐渐减小。当T_f=1 343 K时,液滴的最大铺展宽度约为0.032 8,小于T_f=1 043 K(约0.033 3)、T_f=743 K(约0.033 7)和T_f=443 K(约0.034 4)条件下的最大铺展宽度。值得注意的是,在液滴达到最小收缩宽度后的二次铺展阶段,初始温度较高的液滴回弹与再铺展行为表现出更弱的趋势。特别地,在t^*>5.0之后,T_f=1 343的铺展速度显著低于其他案例。这一现象可解释如下:初始温度升高,流体黏度降低,使得流动更趋于惯性主导,因而液滴铺展与回缩过程更为缓慢。总体而言,初始温度T_f的升高使得液滴的最大铺展宽度、最小收缩宽度均减小,流动演化显著迟缓。

(2)温度依赖系数ϕ的影响。为了分析温度依赖系数ϕ对液滴流动特性的影响,选取4个不同ϕ值(0.001、0.005、0.01以及0.02)进行计算。如图11(b)所示为不同ϕ下液滴铺展宽度随时间的演化情况。在初始阶段(t^*<2.0),液滴的铺展行为基本一致,但随着时间的推移,不同ϕ下的液滴流动特性逐渐分化。随着ϕ的增大,液滴的最大铺展宽度显著减小。ϕ为0.001、0.005、0.01和0.02下的最大铺展宽度分别约为0.034 6、0.033 9、0.033 3和0.032 3;液滴收缩过程中的最小宽度也呈下降趋势。由于收缩阶段的宽度越小,后续再铺展阶段的动能积累受限,因而最终再铺展宽度也随之减小。该趋势表明,ϕ对温度-黏度耦合效应具有重要调控作用。较大的ϕ会增强液滴黏度对温度变化的敏感性,导致高温区域黏度上升,从而使整体流动行为的弹性恢复作用减弱。在这种情况下,液滴的铺展与收缩行为均受到抑制。

4.2.4 其他物理参数对液滴铺展宽度的影响

(1)Weissenberg数(Wi)的影响。Wi用于衡量流动中的弹性力与黏性力的相对强弱。如图12(a)所示为不同Wi条件下液滴铺展宽度随时间的演化,选取Wi为1、3、5、7、9进行数值模拟。从图12中可以观察到,随着Wi的增大,液滴的最大铺展宽度也随之增加。当Wi=9时,液滴的最大铺展宽度约为0.040 6 m,明显大于Wi=7(约0.040 1 m)、Wi=5(约0.039 3 m)、Wi=3(约0.037 8 m)以及Wi=1(约0.034 4 m)。这一趋势主要源于随着Wi的增大,流体的剪切变稀特性增强,使得液滴在撞击过程中铺展速度加快,从而导致更大的铺展范围。

(2)PTT拉伸参数ε的影响。为进一步分析PTT拉伸参数ε对液滴流动行为的影响,选取ε为0.25、0.50、0.75、1.00共4种情况进行模拟。如图12(b)所示不同ε条件下液滴铺展宽度随时间变化的对比结果。结果表明,ε对液滴的最大铺展宽度影响不显著,但对回弹阶段的收缩行为有明显的影响。具体来讲,随着ε的增大,液滴所达到的收缩极小值增大。当ε=0.25时,液滴的最小收缩宽度约为0.031 1,明显小于ε=0.50、0.75与1.00时的液滴收缩宽度,分别约为0.032 1、0.032 8和0.033 3的最小收缩宽度。这一趋势可归因于ε对流体拉伸黏度的调节作用。当ε增大时,黏弹性流体在相同剪切速率下的拉伸黏度会上升,拉伸作用更加显著。这使得液滴的弹性作用减弱,从而削弱了液滴撞击后因弹性势能释放所产生的收缩行为。

(3)流体总黏度η的影响。如图12(c)所示流体总黏度η对液滴铺展宽度的影响。结果表明,液滴的铺展和收缩行为对黏度变化高度敏感。随着η的增大,液滴的最大铺展宽度明显减小。当η=4、10、15与20时,液滴所达到的最大铺展宽度分别约为0.034 4、0.029 1、0.027 1、0.025 9,呈现出明显的递减趋势。此外,液滴在回弹阶段的收缩行为也随黏度的增加而逐渐减弱。当流体黏度较大时,黏性阻尼效应显著增强,抑制了液滴内弹性能量的释放,从而削弱了液滴的收缩行为。这一现象与拉伸参数ε增大时流体弹性减弱具有类似的物理机制,反映出高黏度下液滴行为逐渐转向以黏性主导的牛顿型流动特征。

(4)溶剂黏度比β的影响。从图12(d)中可看出,随着β的增大,液滴最大铺展宽度呈递减趋势。当β为0.1、0.3、0.5、0.7时,对应的最大铺展宽度约为0.034 4、0.032 0、0.030 3、0.029 1。此外,液滴铺展至稳定形态的过程也受到β的显著影响。特别地,β=0.7的液滴在碰撞后的弹性收缩行为几乎消失。这是由于较大的β值对应更高的牛顿溶剂黏度与更低的聚合物黏度,从而削弱了流体体系中的弹性效应,使液滴的动力学演化趋于黏性主导的牛顿流体特征。

5 结论

应用SPH方法对三维非等温PTT黏弹性液滴撞击固壁过程进行数值模拟。通过引入人工应力项以消除张力不稳定性,结合GPU并行算法以提升计算效率,分析了温度参数和其他物理参数对液滴动力学行为的影响,得到如下结论。

(1)通过模拟二维非等温泊肃叶流动和三维等温液滴撞击固壁问题,并通过与其他数值方法得到的结果进行比较,验证了本文所采用的SPH方法在处理非等温黏弹性流动问题时,具有与传统方法相当的准确性,同时避免了复杂的界面追踪过程。

(2)对于三维非等温液滴撞击固壁,通过不同粒子间距的模拟结果对比,表明本文方法在求解三维非等温黏弹性大变形问题时具有良好的收敛性;当粒子数为142万时,GPU实现了约104倍的加速比,显著提升了计算效率。

(3)利用本文方法易于处理多物理场耦合的优势,系统分析了各参数对液滴动力学行为的影响,揭示了以下物理规律:温度的引入会显著影响黏弹性液滴的流动特性。T_f与ϕ越大,黏性作用越强,铺展范围越小;Wi越大,最大铺展宽度越大;ε、η和β越大,弹性力越弱,收缩行为越不明显。通过系统模拟所总结的物理规律,为深入理解与准确预测实际非等温条件下黏弹性液滴的动力学行为提供理论依据。

参考文献

原文顺序 | 出版日期 | 本文引用

[1]	Zhao D, Zhou H, Wang Y, et al. Drop-on-demand (DOD) inkjet dynamics of printing viscoelastic conductive ink[J]. Additive Manufacturing, 2021, 48: 102451.

[2]	Yu H, Zhang W, Yin B, et al. Modeling extrusion process and layer deformation in 3D concrete printing via smoothed particle hydrodynamics[J]. Computer Methods in Applied Mechanics and Engineering, 2024, 420: 116761.

[3]	Suly P, Sevcik J, Dmonte D J, et al. Inkjet printability assessment of weakly viscoelastic fluid: a semidilute polyvinylpyrrolidone solution ink case study[J]. Langmuir, 2021, 37(28): 8557-8568.

[4]	Chang Y J, Chang R Y, Osswald T A. Modeling and simulation of bulk viscoelasticity for amorphous polymers in injection molding[J]. Physics of Fluids, 2023, 35(5): 053109.

[5]	Akbar S S, Mustafa M. Coupled heat and mass transfer to viscoelastic fluid flow in a rotating frame using series and numerical solutions[J]. International Journal of Heat and Fluid Flow, 2024, 106: 109294.

[6]	Rath S, Mahapatra B. Low Reynolds number pulsatile flow of a viscoelastic fluid through a channel: effects of fluid rheology and pulsation parameters[J]. Journal of Fluids Engineering, 2022, 144(2): 021201.

[7]	庄昕, 刘付军, 孙艳萍, 等. 非等温黏弹性聚合物流体圆柱绕流的高精度数值模拟[J]. 应用数学和力学, 2022, 43(12): 1380-1391.

[8]	Zhuang Xin, Liu Fujun, Sun Yanping, et al. High-accuracy numerical simulation of non-isothermal viscoelastic polymer fluid flow around a cylinder[J]. Applied Mathematics and Mechanics, 2022, 43(12): 1380-1391.

[9]	Sharma A, Koch D L. Finite difference method in prolate spheroidal coordinates for freely suspended spheroidal particles in linear flows of viscous and viscoelastic fluids[J]. Journal of Computational Physics, 2023, 495: 112559.

[10]	Venkatesan J, Ganesan S. Computational modeling of impinging viscoelastic droplets[J]. Journal of Non-Newtonian Fluid Mechanics, 2019, 263: 42-60.

[11]	李烨波, 黄禄渊, 姚瑞, 等. 直立走滑断层凹凸体附近闭锁程度和应力积累的时空演化: 来自三维黏弹性数值模拟的视角[J]. 地质力学学报, 2024, 30(6): 878-892.

[12]	Li Yebo, Huang Luyuan, Yao Rui, et al. Spatiotemporal evolution of locking degree and stress accumulation near asperities on vertical strike-slip faults: insights from 3D viscoelastic numerical simulations[J]. Journal of Geomechanics, 2024, 30(6): 878-892.

[13]	Marrone S, Antuono M, Agresta A, et al. A study on the energy consistency in SPH surface tension modelling[J]. Computer Methods in Applied Mechanics and Engineering, 2025, 433: 117473.

[14]	Ramezanzadeh S, Ozbulut M, Yildiz M. Comprehensive numerical study on the behavior of floating structures under challenging ocean conditions using WCSPH[J]. Journal of Fluids and Structures, 2025, 133: 104273.

[15]	张高峰, 姚国文, 孔国英, 等. 库水位循环下岸坡-桥梁桩基相互作用致损机理[J]. 科学技术与工程, 2023, 23(6): 2317-2325.

[16]	Zhang Gaofeng, Yao Guowen, Kong Guoying, et al. Damage mechanism of slope-bridge pile foundation interaction under reservoir water level cycles[J]. Science Technology and Engineering, 2023, 23(6): 2317-2325.

[17]	Kim J H, Yoo H S, Jo Y B, et al. Simulation of shockwave propagation characteristics in nuclear reactor cavity during external steam explosion using unified SPH[J]. Nuclear Engineering and Technology, 2025, 57(4): 103274.

[18]	Guo C, Zhang H, Qian Z, et al. Smoothed-interface SPH model for multiphase fluid-structure interaction[J]. Journal of Computational Physics, 2024, 518: 113336.

[19]	Jiang T, Li Y, Sun P N, et al. A corrected WCSPH scheme with improved interface treatments for the viscous/viscoelastic two-phase flows[J]. Computational Particle Mechanics, 2022, 9(4): 633-653.

[20]	Moinfar Z, Vahabi S, Vahabi M. Numerical simulation of drop deformation under simple shear flow of Giesekus fluids by SPH[J]. International Journal of Numerical Methods for Heat & Fluid Flow, 2023, 33(1): 263-281.

[21]	Sun P, Liu P X. Bleeding simulation with physical viscoelasticity in smooth particle hydrodynamics[J]. Physics of Fluids, 2024, 36(7): 071909.

[22]	许晓阳, 赵雨婷, 李家宇, 等. 非等温黏弹性复杂流动的改进SPH方法模拟[J]. 力学学报, 2023, 55(5): 1099-1112.

[23]	Xu Xiaoyang, Zhao Yuting, Li Jiayu, et al. Improved SPH simulation of complex non-isothermal viscoelastic flow[J]. Chinese Journal of Theoretical and Applied Mechanics, 2023, 55(5): 1099-1112.

[24]	Li Y, Li Y, Feng Z. A coupled particle model with particle shifting technology for simulating transient viscoelastic fluid flow with free surface[J]. Journal of Computational Physics, 2023, 488: 112213.

[25]	Ribau A M, Ferrás L L, Morgado M L, et al. Analytical study of the annular flow of a generalised Phan-Thien-Tanner fluid[J]. Acta Mechanica, 2024, 235(2): 1307-1317.

[26]	Baranov A V. Influence of temperature and pressure on viscoelastic fluid flow in a plane channel[J]. Journal of Engineering Physics and Thermophysics, 2020, 93(5): 1296-1302.

[27]	Xu J, Yu W, Li C, et al. Influence of viscoelastic properties on induced charge electro-osmosis of Phan-Thien-Tanner fluids around a metal cylinder[J]. Journal of Non-Newtonian Fluid Mechanics, 2025, 337: 105397.

[28]	Wang R, Du Q, Xie Y, et al. Large-eddy simulation of non-isothermal flow and heat transfer in an enclosed rotor-stator cavity[J]. Physics of Fluids, 2023, 35(11): 114112.

[29]	许晓阳, 田凌云. 注塑成型过程的多尺度光滑粒子流体动力学模拟[J]. 清华大学学报(自然科学版), 2024, 64(8): 1401-1413.

[30]	Xu Xiaoyang, Tian Lingyun. Multiscale smoothed particle hydrodynamics simulation of injection molding process[J]. Journal of Tsinghua University (Science and Technology), 2024, 64(8): 1401-1413.

[31]	Xu X, Jiang Y L. SPH simulations of transient non-isothermal viscoelastic flows with free surfaces[J]. Journal of Non-Newtonian Fluid Mechanics, 2022, 309: 104905.

[32]	Feng R, Fourtakas G, Rogers B D, et al. A general smoothed particle hydrodynamics (SPH) formulation for coupled liquid flow and solid deformation in porous media[J]. Computer Methods in Applied Mechanics and Engineering, 2024, 419: 116581.

[33]	Liu X, Xu H, Shao S, et al. An improved incompressible SPH model for simulation of wave-structure interaction[J]. Computers & Fluids, 2013, 71: 113-123.

[34]	Diwakar M K, Mohapatra P K, Tripathi S. Smoothed particle hydrodynamics (SPH) investigation of two-dimensional dam break flows[J]. Environmental Fluid Mechanics, 2025, 25(1): 1-20.

[35]	Seifi N, Al-Mamun A. Optimizing memory access efficiency in CUDA kernel via data layout technique[J]. Journal of Computer and Communications, 2024, 12(5): 124-139.

[36]	Zhuang X, Ouyang J, Li W, et al. Three-dimensional simulations of non-isothermal transient flow and flow-induced stresses during the viscoelastic fluid filling process[J]. International Journal of Heat and Mass Transfer, 2017, 104: 374-391.

[37]	Tomé M F, Merejolli R, Paulo G S, et al. An assessment of the PTT model on the impacting drop problem[J]. Journal of the Brazilian Society of Mechanical Sciences and Engineering, 2018, 40: 1-12.

基金资助

国家自然科学基金(12071367)

陕西省“特支计划”青年拔尖人才项目(289890259)

PDF (6698KB)

访问

被引

详细

导航

Received	Published	Revised
2025-06-29	2026-07-08	2025-12-18
Online Date
2026-08-06

摘要