CN110532727B - 可用于常见非牛顿流体的数值模拟方法 - Google Patents

可用于常见非牛顿流体的数值模拟方法 Download PDF

Info

Publication number
CN110532727B
CN110532727B CN201910847006.5A CN201910847006A CN110532727B CN 110532727 B CN110532727 B CN 110532727B CN 201910847006 A CN201910847006 A CN 201910847006A CN 110532727 B CN110532727 B CN 110532727B
Authority
CN
China
Prior art keywords
fluid
external force
newtonian
common non
calculating
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201910847006.5A
Other languages
English (en)
Other versions
CN110532727A (zh
Inventor
吴伟伟
王洲舟
俞亮
丁爽
宋爱平
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Yangzhou University
Original Assignee
Yangzhou University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Yangzhou University filed Critical Yangzhou University
Priority to CN201910847006.5A priority Critical patent/CN110532727B/zh
Publication of CN110532727A publication Critical patent/CN110532727A/zh
Application granted granted Critical
Publication of CN110532727B publication Critical patent/CN110532727B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Abstract

本发明涉及流体动力学计算机模拟运算领域内的一种可用于常见非牛顿流体的数值模拟方法，首先基于格子玻尔兹曼方法建立含外力项的多松弛参数格子玻尔兹曼方法；之后根据三种常见非牛顿流体的流变方程，将它们的非牛顿特征转换成离散外力项，代入计算后，完成数值模拟过程。该方法可以有效改善应用格子玻尔兹曼方法模拟非牛顿流体时的稳定性和准确性，结合介观格子玻尔兹曼方法原理简单、计算方便、易于实现的特点，可指导非牛顿流体的流动仿真。

Description

可用于常见非牛顿流体的数值模拟方法
技术领域
本发明涉及流体动力学计算机模拟运算领域，特别涉及一种非牛顿流体的数值模拟方法。
背景技术
随着流体动力学的不断发展，针对复杂的流体流动，涌现出一些不同于传统有限元计算的方法，如光滑粒子流方法、分子动力学方法以及格子玻尔兹曼方法等，从微观或介观角度更细致地诠释流体的流动过程，解决了利用传统有限元软件仿真，求解复杂微分方程组容易陷入局部解而无法求得最终解的难题，其中格子玻尔兹曼方法具有原理简单、计算方便、易于实现的特点，因此越来越广泛地用于解决复杂的流体流动。
中国专利数据库中，公开了一种移动粒子半隐式算法中自由表面流动模型的构建方法，其申请号：CN201210349290.1，申请日：20120919，公开号：CN102867094A，公开日：20130109，该方法引入的湍流模型包括静态Smagorinsky模型和拉格朗日形式的动态Smagorinsky模型；采用变黏度牛顿流体模型处理本构方程形如μ＝f(|γ|)的非牛顿流体；引入三次样条核函数，采用光滑粒子流体动力学方法的散度离散方案离散非牛顿流体的剪切应力，处理非牛顿流体自由表面流动，其可用于非牛顿流体的计算。
现有技术的不足之处在于：对于非牛顿流体仿真，在利用普通格子玻尔兹曼方法时，会出现稳定性较弱，误差偏大的情况。
发明内容
本发明的目的是提供一种可用于常见非牛顿流体的数值模拟方法，可以有效改善格子玻尔兹曼方法模拟非牛顿流体流动时稳定性弱和误差偏大的不足。
本发明的目的是这样实现的：一种可用于常见非牛顿流体的数值模拟方法，包括步骤：
(1)基于格子玻尔兹曼方法，建立含有外力项的多松弛参数格子玻尔兹曼方法；
(2)根据三种常见非牛顿流体的流变方程，将它们的非牛顿特征转换成离散外力项，代入(1)中计算，完成数值模拟过程。
所述步骤(1)包括如下分步骤：
S2.1定义各物理量初值，及划分网格区域；
S2.2计算平衡态分布函数
Figure BDA0002195598030000021
式中feq(r,t)-流体在时刻t，位置r时平衡态分布函数；
ωi-权重系数，具体定义为
Figure BDA0002195598030000022
ρ-流体的密度；
ei-离散速度，具体定义为
Figure BDA0002195598030000023
u-流体的流动速度；
cs-格子声速，具体定义为
Figure BDA0002195598030000024
c-格子速度，具体定义为c＝δx/δt；
δx-格子步长；
δt-时间步长；
S2.3计算碰撞步，具体公式为
Figure BDA0002195598030000031
式中f+(r,t)-流体执行碰撞步后在时刻t，位置r时的分布函数；
f(r,t)-流体在时刻t，位置r时的分布函数；
Figure BDA0002195598030000039
-中间变量，具体定义为/>
Figure BDA00021955980300000310
M-转换矩阵，具体定义为
Figure BDA0002195598030000032
Figure BDA0002195598030000033
-松弛过程相关的主对角线矩阵，具体定义为/>
Figure BDA0002195598030000034
其中s8＝1/τ；
τ-松弛时间，
Figure BDA0002195598030000035
ρ-动力粘度；
F’-计算外力项，具体定义为
Figure BDA0002195598030000036
I-单位矩阵；
Figure BDA0002195598030000037
-外力项分项，对于牛顿流体/>
Figure BDA0002195598030000038
S2.4计算应变率张量和剪切率，涉及公式如下：
Figure BDA0002195598030000041
式中Sαβ-应变率张量；
eiα,eiβ-分别为x方向和y方向离散速度ei的分量；
DII-应变率张量第二不变量；
l-仿真维度，这里l＝2；
Figure BDA0002195598030000042
-剪切率；
S2.5计算迁移步，具体公式为f(r+eiδt,t+δt)＝f+(r,t)
式中f(r+eiδt,t+δt)-流体在时刻t+δt，位置r+eiδt时的分布函数。
S2.6采用非平衡反弹形式进行壁面及边界处理；
S2.7计算相对误差
Figure BDA0002195598030000043
若其小于10-4，则跳转至S2.9，否则执行S2.8；
S2.8跳至步骤S2.2，执行下一次循环计算；
S2.9计算所需的宏观量，所述宏观量至少包括速度和压力分布，并利用图形呈现具体的流线图、压力场和速度场分布图，完成模拟过程。
进一步地，所述内容(2)根据三种常见非牛顿流体的流变方程，将它们的非牛顿特征转换成特殊的外力项，代入(1)中计算，完成数值模拟过程，具体内容如下：
A幂律流体
对于幂律流体，其离散外力项可以表达为
Figure BDA0002195598030000051
式中μp0-幂律流体流变方程的粘度系数。
B宾汉流体
对于宾汉流体，其离散外力项可以表达为
Figure BDA0002195598030000052
式中μb0-宾汉流体流变方程的粘度系数；
τb0-宾汉流体流变方程的初始屈服应力；
m-一个与流体粘度、应力相关的参数；
Figure BDA0002195598030000053
-剪切率。
C Herschel-Bulkley流体
对于Herschel-Bulkley流体，其离散外力项可以表达为
Figure BDA0002195598030000054
式中μh0-Herschel-Bulkley流体流变方程的粘度系数；
τh0-Herschel-Bulkley流体流变方程的初始屈服应力；
m-一个与流体粘度、应力相关的参数；
Figure BDA0002195598030000055
-剪切率。
进一步地，本方法可通过编程实现，在利用MATLAB软件编程时，对于S2.2-S2.7步可利用Cell元胞数组计算代替for循环，可以有效提高运算效率。
与现有技术相比，本发明的有益效果在于：该方法可以有效改善应用格子玻尔兹曼方法模拟非牛顿流体时的稳定性和准确性，结合介观格子玻尔兹曼方法原理简单、计算方便、易于实现的特点，可指导常见非牛顿流体的流动仿真。
附图说明
图1为本发明的流程图。
具体实施方式
下面结合实施例及附图对本发明作进一步详细的描述，但本发明的的实施方式不限于此。
如图1所示，一种可用于常见非牛顿流体的数值模拟方法，其特征在于包括步骤：
(1)基于格子玻尔兹曼方法，建立含有外力项的多松弛参数格子玻尔兹曼方法；
(2)根据三种常见非牛顿流体的流变方程，将它们的非牛顿特征转换成离散外力项，代入(1)中计算，完成数值模拟过程。
所述步骤(1)包括如下分步骤：
S2.1定义各物理量初值，及划分网格区域；
S2.2计算平衡态分布函数
Figure BDA0002195598030000061
式中feq(r,t)-流体在时刻t，位置r时平衡态分布函数；
ωi-权重系数，具体定义为
Figure BDA0002195598030000062
ρ-流体的密度；
ei-离散速度，具体定义为
Figure BDA0002195598030000071
u-流体的流动速度；
cs-格子声速，具体定义为
Figure BDA0002195598030000072
c-格子速度，具体定义为c＝δx/δt；
δx-格子步长；
δt-时间步长；
S2.3计算碰撞步，具体公式为
Figure BDA0002195598030000073
式中f+(r,t)-流体执行碰撞步后在时刻t，位置r时的分布函数；
f(r,t)-流体在时刻t，位置r时的分布函数；
Figure BDA0002195598030000077
-中间变量，具体定义为/>
Figure BDA0002195598030000078
M-转换矩阵，具体定义为
Figure BDA0002195598030000074
Figure BDA0002195598030000075
-松弛过程相关的主对角线矩阵，具体定义为/>
Figure BDA0002195598030000076
其中s8＝1/τ；
τ-松弛时间，
Figure BDA0002195598030000081
ρ-动力粘度；
F’-计算外力项，具体定义为
Figure BDA0002195598030000082
I-单位矩阵；
Figure BDA0002195598030000083
-外力项分项，对于牛顿流体/>
Figure BDA0002195598030000084
/>
S2.4计算应变率张量和剪切率，涉及公式如下：
Figure BDA0002195598030000085
式中Sαβ-应变率张量；
eiα,eiβ-分别为x方向和y方向离散速度ei的分量；
DII-应变率张量第二不变量；
l-仿真维度，这里l＝2；
Figure BDA0002195598030000086
-剪切率；
S2.5计算迁移步，具体公式为f(r+eiδt,t+δt)＝f+(r,t)
式中f(r+eiδt,t+δt)-流体在时刻t+δt，位置r+eiδt时的分布函数。
S2.6采用非平衡反弹形式进行壁面及边界处理；
S2.7计算相对误差
Figure BDA0002195598030000087
若其小于10-4，则跳转至S2.9，否则执行S2.8；
S2.8跳至步骤S2.2，执行下一次循环计算；
S2.9计算所需的宏观量，所述宏观量至少包括速度和压力分布，并利用图形呈现具体的流线图、压力场和速度场分布图，完成模拟过程。
步骤(2)中针对幂律流体，其离散外力项表达为
Figure BDA0002195598030000091
式中μp0-幂律流体流变方程的粘度系数。
步骤(2)中针对宾汉流体，其离散外力项表达为
Figure BDA0002195598030000092
式中μb0-宾汉流体流变方程的粘度系数；
τb0-宾汉流体流变方程的初始屈服应力；
m-一个与流体粘度、应力相关的参数；
Figure BDA0002195598030000093
-剪切率。
步骤(2)中针对Herschel-Bulkley流体，其离散外力项表达为
Figure BDA0002195598030000094
式中μh0-Herschel-Bulkley流体流变方程的粘度系数；
τh0-Herschel-Bulkley流体流变方程的初始屈服应力；
m-一个与流体粘度、应力相关的参数；
Figure BDA0002195598030000095
-剪切率。
利用MATLAB软件编程，对于S2.2-S2.7步利用Cell元胞数组计算代替for循环，以提高运算效率。
本发明并不局限于上述实施例，在本发明公开的技术方案的基础上，本领域的技术人员根据所公开的技术内容，不需要创造性的劳动就可以对其中的一些技术特征作出一些替换和变形，这些替换和变形均在本发明的保护范围内。

Claims (5)

1.一种可用于常见非牛顿流体的数值模拟方法，其特征在于包括步骤：
(1)基于格子玻尔兹曼方法，建立含有外力项的多松弛参数格子玻尔兹曼方法；包括如下分步骤：
S1.1定义各物理量初值，及划分网格区域；
S1.2计算平衡态分布函数
Figure QLYQS_1
式中feq(r,t)-流体在时刻t，位置r时平衡态分布函数；
ωi-权重系数，具体定义为
Figure QLYQS_2
ρ-流体的密度；
ei-离散速度，具体定义为
Figure QLYQS_3
u-流体的流动速度；
cs-格子声速，具体定义为
Figure QLYQS_4
c-格子速度，具体定义为c＝δx/δt；
δx-格子步长；
δt-时间步长；
S1.3计算碰撞步，具体公式为
Figure QLYQS_5
式中f+(r,t)-流体执行碰撞步后在时刻t，位置r时的分布函数；
f(r,t)-流体在时刻t，位置r时的分布函数；
Figure QLYQS_6
-中间变量，具体定义为/>
Figure QLYQS_7
M-转换矩阵，具体定义为/>
Figure QLYQS_8
Figure QLYQS_9
-松弛过程相关的主对角线矩阵，具体定义为/>
Figure QLYQS_10
其中s8＝1/τ；
τ-松弛时间，
Figure QLYQS_11
ρ-动力粘度；
F’-计算外力项，具体定义为
Figure QLYQS_12
I-单位矩阵；
Figure QLYQS_13
-外力项分项，对于牛顿流体/>
Figure QLYQS_14
S1.4计算应变率张量和剪切率，涉及公式如下：
Figure QLYQS_15
式中Sαβ-应变率张量；
eiα,eiβ-分别为x方向和y方向离散速度ei的分量；
DII-应变率张量第二不变量；
l-仿真维度，这里l＝2；
Figure QLYQS_16
-剪切率；
S1.5计算迁移步，具体公式为f(r+eiδt,t+δt)＝f+(r,t)
式中f(r+eiδt,t+δt)-流体在时刻t+δt，位置r+eiδt时的分布函数；
S1.6采用非平衡反弹形式进行壁面及边界处理；
S1.7计算相对误差
Figure QLYQS_17
若其小于10-4，则跳转至S1.9，否则执行S1.8；
S1.8跳至步骤S1.2，执行下一次循环计算；
S1.9计算所需的宏观量，所述宏观量至少包括速度和压力分布，并利用图形呈现具体的流线图、压力场和速度场分布图，完成模拟过程；
(2)根据三种常见非牛顿流体的流变方程，将它们的非牛顿特征转换成离散外力项，执行步骤(1)中的各分步骤，完成数值模拟过程。
2.根据权利要求1所述的可用于常见非牛顿流体的数值模拟方法，其特征在于，所述步骤(2)中针对幂律流体，其离散外力项表达为
Figure QLYQS_18
式中μp0-幂律流体流变方程的粘度系数。
3.根据权利要求1所述的可用于常见非牛顿流体的数值模拟方法，其特征在于，所述步骤(2)中针对宾汉流体，其离散外力项表达为
Figure QLYQS_19
式中μb0-宾汉流体流变方程的粘度系数；
τb0-宾汉流体流变方程的初始屈服应力。
4.根据权利要求1所述的可用于常见非牛顿流体的数值模拟方法，其特征在于，所述步骤(2)中针对Herschel-Bulkley流体，其离散外力项表达为
Figure QLYQS_20
式中μh0-Herschel-Bulkley流体流变方程的粘度系数；
τh0-Herschel-Bulkley流体流变方程的初始屈服应力；
m-一个与流体粘度、应力相关的参数；
Figure QLYQS_21
-剪切率。
5.根据权利要求1-4任一项所述的可用于常见非牛顿流体的数值模拟方法，其特征在于，利用MATLAB软件编程，对于S2.2-S2.7步利用Cell元胞数组计算代替for循环，以提高运算效率。
CN201910847006.5A 2019-09-09 2019-09-09 可用于常见非牛顿流体的数值模拟方法 Active CN110532727B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910847006.5A CN110532727B (zh) 2019-09-09 2019-09-09 可用于常见非牛顿流体的数值模拟方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910847006.5A CN110532727B (zh) 2019-09-09 2019-09-09 可用于常见非牛顿流体的数值模拟方法

Publications (2)

Publication Number Publication Date
CN110532727A CN110532727A (zh) 2019-12-03
CN110532727B true CN110532727B (zh) 2023-05-23

Family

ID=68667660

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910847006.5A Active CN110532727B (zh) 2019-09-09 2019-09-09 可用于常见非牛顿流体的数值模拟方法

Country Status (1)

Country Link
CN (1) CN110532727B (zh)

Families Citing this family (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111174648B (zh) * 2020-01-13 2022-05-06 内蒙古工业大学 一种基于压板结构的降低药柱气孔率的浇注方法
CN112287622B (zh) * 2020-12-28 2021-03-09 中国人民解放军国防科技大学 基于链接方向人工压缩的快速湍流数值模拟方法和装置
CN112765841B (zh) * 2020-12-31 2024-05-07 西北工业大学深圳研究院 一种用以解决流体变物性计算的预处理格子玻尔兹曼方法
CN113051842A (zh) * 2021-03-05 2021-06-29 清华大学 一种非牛顿流体模拟方法及装置
CN113468829B (zh) * 2021-06-24 2024-04-26 西南石油大学 基于孔隙网络模型的非稳态非牛顿两相流体驱替模拟方法

Also Published As

Publication number Publication date
CN110532727A (zh) 2019-12-03

Similar Documents

Publication Publication Date Title
CN110532727B (zh) 可用于常见非牛顿流体的数值模拟方法
Zeng et al. Numerical simulation of proppant transport in hydraulic fracture with the upscaling CFD-DEM method
Wang et al. Lattice Boltzmann simulation of particle-laden turbulent channel flow
Ammar et al. A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids
Guo et al. Force imbalance in lattice Boltzmann equation for two-phase flows
Xu et al. A level-set continuum method for two-phase flows with insoluble surfactant
Misra et al. A vortex-based subgrid stress model for large-eddy simulation
Chen et al. An evaluation of the MPM for simulating dynamic failure with damage diffusion
Afra et al. An immersed boundary-lattice Boltzmann method combined with a robust lattice spring model for solving flow–structure interaction problems
US9283695B1 (en) Computer-implemented simulation method and non-transitory computer medium capable of predicting fiber orientation for use in a molding process
Wen et al. Lattice-type-dependent momentum-exchange method for moving boundaries
Dimakopoulos et al. The PAL (Penalized Augmented Lagrangian) method for computing viscoplastic flows: A new fast converging scheme
CN104281730B (zh) 一种大转动变形的板壳结构动响应的有限元分析方法
Li et al. A fractional-step lattice Boltzmann method for multiphase flows with complex interfacial behavior and large density contrast
CN104573269B (zh) 一种基于强耦合整体技术的索膜结构抗风设计方法
CN104239625A (zh) 一种基于修正流体运动方程线性迭代的稳态求解方法
Mompean et al. Numerical prediction of three-dimensional time-dependent viscoelastic extrudate swell using differential and algebraic models
CN107657075A (zh) 模拟地下水介质交界面处达西速度的区域分解有限元法
De Paulo et al. A marker-and-cell approach to viscoelastic free surface flows using the PTT model
Cai et al. Numerical study on bi-material interface crack using symplectic analytical singular element
Sun et al. Material point method and smoothed particle hydrodynamics simulations of fluid flow problems: a comparative study
Nguyen et al. Numerical simulation of annular bubble plume by vortex in cell method
Wang et al. The immersed boundary-lattice Boltzmann method for solving solid-fluid interaction problem with Navier-slip boundary condition
Zhou et al. Numerical investigation of non-Newtonian power law flows using B-spline material point method
Wang et al. Comparisons of two representative methods classified as immersed boundary and domain methods

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant