CN110287590B - 基于算子分裂及改进半拉格朗日求解污染物传播的方法 - Google Patents

基于算子分裂及改进半拉格朗日求解污染物传播的方法 Download PDF

Info

Publication number
CN110287590B
CN110287590B CN201910551021.5A CN201910551021A CN110287590B CN 110287590 B CN110287590 B CN 110287590B CN 201910551021 A CN201910551021 A CN 201910551021A CN 110287590 B CN110287590 B CN 110287590B
Authority
CN
China
Prior art keywords
particle
velocity
equation
time
concentration
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
CN201910551021.5A
Other languages
English (en)
Other versions
CN110287590A (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.)
Tianjin University
Original Assignee
Tianjin 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 Tianjin University filed Critical Tianjin University
Priority to CN201910551021.5A priority Critical patent/CN110287590B/zh
Publication of CN110287590A publication Critical patent/CN110287590A/zh
Application granted granted Critical
Publication of CN110287590B publication Critical patent/CN110287590B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明涉及流体力学技术领域，为实现在满足数值精度的前提下使用改进半拉格朗日方法更简洁快速地模拟流体中的污染物传播问题，本发明，基于算子分裂及改进半拉格朗日求解污染物传播的方法，步骤如下：步骤一，初始化系统的相关变量和运行参数；步骤二，生成粒子信息；步骤三，列出求解方程并迭代计算：根据算子分裂思想将对流、扩散和反应三个过程分别采用不同的数值方法计算；步骤四，输出结果：每一个时间步结束，保存该时间步的计算结果并更新；结束时间循环，输出最终结果。本发明主要应用于模拟流体中的污染物传播场合。

Description

基于算子分裂及改进半拉格朗日求解污染物传播的方法
技术领域
本发明涉及流体力学技术领域，具体是涉及一种基于算子分裂求解污染物传播问题的改进半拉格朗日方法。
背景技术
传统半拉格朗日方法由于插值计算物质浓度导致缺乏守恒性，并且经过长时间的数值积分会带来误差累积。另外在处理非连续分布问题上，对流过程的插值计算会在间断点处发生数值振荡。通过速度插值的改进半拉格朗日方法能够克服间断点处发生的非物理振荡问题，计算过程更加简化，可以提高计算效率，基于算子分裂的改进半拉格朗日方法可处理非连续分布、非匀速流场等问题。
发明内容
为克服现有技术的不足，本发明旨在提出一种基于算子分裂求解污染物传播问题的改进半拉格朗日方法，该方法求解欧拉-拉格朗日系统下的对流扩散方程，充分考虑了网格方法以及传统半拉格朗日带来各种数值问题，在满足数值精度的前提下使用改进半拉格朗日方法更简洁快速地模拟流体中的污染物传播问题。为此，本发明采取的技术方案是，基于算子分裂及改进半拉格朗日求解污染物传播的方法，步骤如下：
步骤一，初始化系统的相关变量和运行参数；
步骤二，生成粒子信息；
步骤三，列出求解方程并迭代计算：
根据算子分裂思想将对流、扩散和反应三个过程分别采用不同的数值方法计算；
步骤四，输出结果：
①每一个时间步结束，保存该时间步的计算结果并更新；
②结束时间循环，输出最终结果。
步骤三具体地：
求解的一维对流扩散反应方程为：
其中，C为物质浓度，x和t分别为空间和时间坐标，u为流体流动速度，ε为物质扩散系数，G(C)为反应源项，代表单位时间单位体积空间内因化学反应生成的组分量，当G(C)＝0时，方程(1)即为对流扩散方程：
其中，等号右侧为扩散项，一维纯对流方程的欧拉形式为：
根据运动流体的物质导数将方程(3)转化为拉格朗日形式：
从方程(4)和(5)，在拉格朗日体系下，对流过程中物质不会随时间变化而销毁或创造，仅仅是随流体发生运动，解得速度之后，求得浓度在时间和空间上的变化规律；
改进拉格朗日方法针对速度场采用合适的数值插值方法插值粒子速度，插值表达式为：
Vh＝a1Vj-1+a2Vj+a3VXj-1+a4VXj (6)
其中，Vh为当前时间步粒子速度，系数a1＝α2(3-2α)，a2＝1-a1，a3＝α2(1-α)△x，a4＝-α(1-α)2△x，α为柯朗数的小数部分，Vj-1和Vj分别为所求粒子的相邻速度网格点的速度，VXj-1和VXj为相邻速度网格点速度导数；
一维纯扩散方程为：
将方程(7)中的时间离散通过向前差分，空间离散使用二阶中心差分得到：
其中n表示当前时刻，n-1表示上一时刻，j-1，j，j+1表示空间网格位置，为对流过程结果，得扩散结果：
最后选用欧拉法Euler、龙格-库塔二阶法RK2或龙格-库塔四阶法RK4求解反应项；
通过方程(6)、(9)和(10)分别计算出每个时间步下每个流体粒子的速度和浓度，进而求解出粒子新的位置。
用RK4进行时间积分，离散形式如下：
通过方程(6)、(9)和(10)分别计算出每个时间步下每个流体粒子的速度和浓度，进而求解出粒子新的位置，具体计算过程为：
①循环每个时间步；
②初始化计算域中的粒子后，根据已知速度场，首先通过方程(6)计算目标粒子的速度，再根据速度计算粒子位置的变化，最后更新该粒子的速度和位置；
③将上一步对流结果代入方程(9)计算出粒子扩散过程浓度变化，并更新粒子浓度；
④通过方程(10)计算反应过程粒子浓度变化，并更新粒子浓度；
⑤若有粒子运动至计算域外则删除该粒子，并在入口边界补充新的粒子；
循环计算域内的每一个流体粒子，重复做②③④⑤中的操作。
相关变量信息和运行参数以无量纲化形式设置：计算域是长度为1的一维管道,管道阀门位于0处，管道内有浓度分别为1，0，0的三种反应物u1，u2，u3，此时阀门处于关闭状态，管道内混合物处于稳态,管道外u1，u2，u3浓度分别为0，1，0，求解当阀门打开后，流体以速度1发生流动，计算域内反应物u2浓度变化：
步骤一，初始化。初始化系统的相关变量和运行参数，具体包括：
模拟的计算域为阀门右侧x＝[0,1]的一维矩形空间，内有三种处于稳定状态的反应物，三种物质粒子间距△x＝0.0005，计算时间步长为0.0002，计算时间为0.8，阀门左侧三种反应物浓度分别为u1＝0，u2＝1，u3＝0，打开阀门后，流体以速度1开始流动，三种反应物开始相互反应；
步骤二，生成粒子信息，具体包括：
初始化粒子步骤中，三种反应物分别生成流体粒子2000个，初始浓度分别为u1＝1，u2＝0，u3＝0；
步骤三，列出求解方程并迭代计算；
求解的一维对流反应方程拉格朗日形式为：
其中，U1，U2，U3为反应物浓度，x和t分别为空间和时间坐标，V为流体流动速度，方程(2)(3)(4)等号右侧均为反应项，代表单位时间单位体积空间内因化学反应生成的组分量，一维纯对流方程的欧拉形式为：
根据运动流体的物质导数可以将方程(5)转化为拉格朗日形式：
从方程(6)和(7)看出，解得速度之后，求得浓度在时间和空间上的变化规律；
改进拉格朗日方法针对速度场采用合适的数值插值方法插值粒子速度，插值表达式为：
Vh＝a1Vj-1+a2Vj+a3VXj-1+a4VXj (8)
其中，Vh为当前时间步粒子速度，系数a1＝α2(3-2α)，a2＝1-a1，a3＝α2(1-α)△x，a4＝-α(1-α)2△x，α为柯朗数的小数部分，Vj-1和Vj分别为所求粒子的相邻速度网格点的速度，VXj-1和VXj为相邻速度网格点速度导数；
采用RK4求解反应项，离散形式如下：
通过方程(1)、(2)、(3)、(8)和(9)分别计算出每个时间步下每个流体粒子的速度和浓度，进而求解出粒子新的位置。
本发明的特点及有益效果是：
本发明在求解欧拉-拉格朗日系统下的对流扩散方程过程中，充分考虑了网格方法以及传统半拉格朗日带来各种数值问题，在满足数值精度的前提下使用改进半拉格朗日方法更简洁快速地模拟流体中的污染物传播问题。
附图说明：
图1程序流程图。
图2一维对流反应问题物理模型。
图3CASE1中在不同时刻反应物u1(a)，反应物u2(b)，反应物u3(c)随位移的变化。
图4CASE1中在位置0.5处反应物u1(a)，反应物u2(b)，反应物u3(c)随时间的变化。
图5一维对流扩散问题物理模型。
图6 CASE2速度场。
图7 CASE2中在不同时刻物质浓度随位移的变化。
具体实施方式
本发明解决的技术问题是提供一种基于算子分裂求解污染物传播问题的改进半拉格朗日方法，该方法求解欧拉-拉格朗日系统下的对流扩散方程，充分考虑了网格方法以及传统半拉格朗日带来各种数值问题，在满足数值精度的前提下使用改进半拉格朗日方法更简洁快速地模拟流体中的污染物传播问题。
为了解决上述技术问题，本发明的技术方案是：
一种基于算子分裂求解污染物传播问题的改进半拉格朗日方法，包括以下步骤：
步骤一，初始化系统的相关变量和运行参数；
步骤二，生成粒子信息；
步骤三，列出求解方程并迭代计算：
根据算子分裂思想将对流、扩散和反应三个过程分别采用不同的数值方法计算，保证分裂后的方程更易求解且格式灵活，具有较好的稳定性。
求解的一维对流扩散反应方程为：
其中，C为物质浓度，x和t分别为空间和时间坐标，u为流体流动速度，ε为物质扩散系数，G(C)为反应源项，代表单位时间单位体积空间内因化学反应生成的组分量。当G(C)＝0时，方程(1)即为对流扩散方程：
其中，等号右侧为扩散项。
一维纯对流方程的欧拉形式为：
根据运动流体的物质导数可以将方程(3)转化为拉格朗日形式：
从方程(4)和(5)可以看出，在拉格朗日体系下，对流过程中物质不会随时间变化而销毁或创造，仅仅是随流体发生运动，解得速度之后，可以求得浓度在时间和空间上的变化规律。
改进拉格朗日方法针对速度场采用合适的数值插值方法插值粒子速度，插值表达式为：
Vh＝a1Vj-1+a2Vj+a3VXj-1+a4VXj (6)
其中，Vh为当前时间步粒子速度，系数a1＝α2(3-2α)，a2＝1-a1，a3＝α2(1-α)△x，a4＝-α(1-α)2△x，α为柯朗数的小数部分，Vj-1和Vj分别为所求粒子的相邻速度网格点的速度，VXj-1和VXj为相邻速度网格点速度导数。
一维纯扩散方程为：
将方程(7)中的时间离散通过向前差分，空间离散使用二阶中心差分得到：
其中n表示当前时刻，n-1表示上一时刻，j-1，j，j+1表示空间网格位置，为对流过程结果，所以可得扩散结果：
最后求解反应项，可用Euler、RK2、RK4等方法进行时间积分，以精度最高的RK4为例，离散形式如下：
通过方程(6)、(9)和(10)可以分别计算出每个时间步下每个流体粒子的速度和浓度，进而求解出粒子新的位置。具体计算过程为：
⑥循环每个时间步；
⑦初始化计算域中的粒子后，根据已知速度场，首先通过方程(6)计算目标粒子的速度，再根据速度计算粒子位置的变化，最后更新该粒子的速度和位置。
⑧将上一步对流结果代入方程(9)计算出粒子扩散过程浓度变化，并更新粒子浓度。
⑨通过方程(10)计算反应过程粒子浓度变化，并更新粒子浓度。
⑩若有粒子运动至计算域外则删除该粒子，并在入口边界补充新的粒子。
11循环计算域内的每一个流体粒子，重复做②③④⑤中的操作。
步骤四，输出结果：
③每一个时间步结束，保存该时间步的计算结果并更新。
④结束时间循环，输出最终结果。
进一步地，在上述方案中，步骤一所述初始化系统的相关变量和步骤二所述生成粒子信息具体包括：
①初始化与问题相关的变量信息和运行参数；
②生成流体粒子信息，在计算域初始化粒子分布，并添加初始信息；
③生成速度场信息，在有粒子移动至计算域外时删除该粒子并在入口边界补充粒子。
进一步地，在上述方案中，所述初始化变量信息和运行参数具体设置如下：
CASE 1：一维对流反应问题
本实验模拟问题的物理模型见图2，相关变量信息和运行参数以无量纲化形式设置：计算域是长度为1的一维管道。管道阀门位于0处，管道内有浓度分别为1，0，0的三种反应物u1，u2，u3，此时阀门处于关闭状态，管道内混合物处于稳态。管道外u1，u2，u3浓度分别为0，1，0，求解当阀门打开后，流体以速度1发生流动，计算域内反应物u2浓度变化。
下面结合附图对本发明做进一步详细的描述。
基于算子分裂求解污染物传播问题的改进半拉格朗日方法为：
步骤一，初始化。初始化系统的相关变量和运行参数，具体包括：
如图2所示，为本发明的实验模拟的计算域为阀门右侧x＝[0,1]的一维矩形空间，内有三种处于稳定状态的反应物，三种物质粒子间距△x＝0.0005，计算时间步长为0.0002，计算时间为0.8。阀门左侧三种反应物浓度分别为u1＝0，u2＝1，u3＝0，打开阀门后，流体以速度1开始流动，三种反应物开始相互反应。
步骤二，生成粒子信息，具体包括：
初始化粒子步骤中，三种反应物分别生成流体粒子2000个，初始浓度分别为u1＝1，u2＝0，u3＝0。
步骤三，列出求解方程并迭代计算。
求解的一维对流反应方程拉格朗日形式为：
其中，U1，U2，U3为反应物浓度，x和t分别为空间和时间坐标，V为流体流动速度，方程(2)(3)(4)等号右侧均为反应项，代表单位时间单位体积空间内因化学反应生成的组分量。一维纯对流方程的欧拉形式为：
根据运动流体的物质导数可以将方程(5)转化为拉格朗日形式：
从方程(6)和(7)可以看出，在拉格朗日体系下，对流过程中物质不会随时间变化而销毁或创造，仅仅是随流体发生运动，解得速度之后，可以求得浓度在时间和空间上的变化规律。
改进拉格朗日方法针对速度场采用合适的数值插值方法插值粒子速度，插值表达式为：
Vh＝a1Vj-1+a2Vj+a3VXj-1+a4VXj (8)
其中，Vh为当前时间步粒子速度，系数a1＝α2(3-2α)，a2＝1-a1，a3＝α2(1-α)△x，a4＝-α(1-α)2△x，α为柯朗数的小数部分，Vj-1和Vj分别为所求粒子的相邻速度网格点的速度，VXj-1和VXj为相邻速度网格点速度导数。
然后求解反应项，可用Euler、RK2、RK4等方法进行时间积分，以精度最高的RK4为例，离散形式如下：
通过方程(1)、(2)、(3)、(8)和(9)可以分别计算出每个时间步下每个流体粒子的速度和浓度，进而求解出粒子新的位置。具体计算过程为：
①循环每个时间步；
②初始化计算域中的粒子后，根据已知速度场，首先通过方程(8)计算目标粒子的速度，再根据速度计算粒子位置的变化，最后更新该粒子的速度和位置。
③将上一步对流结果代入方程(9)计算出粒子反应过程浓度变化，并更新粒子浓度。
④若有粒子运动至计算域外则删除该粒子，并在入口边界补充新的粒子。
⑤循环计算域内的每一个流体粒子，重复做②③④中的操作。
步骤四，输出结果：
①每完成一个时间步的计算就更新其结果；
②完成时间步的循环，输出最终结果。
CASE 2：一维变速对流扩散问题
本实验模拟问题的物理模型见图5，计算域是长度为1m的一维空间。污染物满足高斯分布，中心初始浓度为1，位于计算域的0.25m处，随流体流动并慢慢扩散，流体流速逐渐加快。求解污染物浓度随时间变化情况。
下面结合附图对本发明做进一步详细地描述。
基于算子分裂求解污染物传播问题的改进半拉格朗日方法为：
步骤一，初始化。初始化系统的相关变量和运行参数，具体包括：
如图五所示，实验模拟的是计算域长度为1m的一维空间，满足高斯分布的污染物中心位于计算域0.25m处，粒子间距△x＝0.01m。污染物扩散系数D＝0.00005，速度场v＝0.01(1+x)，速度网格间距△v＝0.01。计算时间步长为0.1s，计算时间为10s。
步骤二，生成粒子信息，具体包括：
初始分布流体粒子100个粒子，满足中心浓度为1边界浓度为0的高斯分布。
步骤三，列出求解方程并迭代计算。
求解的一维对流扩散方程为：
其中，C为物质浓度，x和t分别为空间和时间坐标，u为流体流动速度，ε为物质扩散系数。
一维纯对流方程的欧拉形式为：
根据运动流体的物质导数可以将方程(3)转化为拉格朗日形式：
从方程(3)(4)可以看出，在拉格朗日体系下，对流过程中物质不会随时间变化而销毁或创造，仅仅是随流体发生运动，解得速度之后，可以求得浓度在时间和空间上的变化规律。
改进拉格朗日方法针对速度场采用合适的数值插值方法插值粒子速度，插值表达式为：
Vh＝a1Vj-1+a2Vj+a3VXj-1+a4VXj (5)
其中，Vh为当前时间步粒子速度，系数a1＝α2(3-2α)，a2＝1-a1，a3＝α2(1-α)△x，a4＝-α(1-α)2△x，α为柯朗数的小数部分，Vj-1和Vj分别为所求粒子的相邻速度网格点的速度，VXj-1和VXj为相邻速度网格点速度导数。
一维纯扩散方程为：
将方程(7)中的时间离散通过向前差分，空间离散使用二阶中心差分得到：
其中n表示当前时刻，n-1表示上一时刻，j-1，j，j+1表示空间网格位置，为对流过程结果，所以可得扩散结果：
通过方程(5)和(8)可以计算出每个时间步下每个流体粒子的速度和浓度，进而求解出粒子新的位置。具体计算过程为：
①循环每个时间步；
②初始化计算域中的粒子后，根据已知速度场，首先通过方程(5)计算目标粒子的速度，再根据速度计算粒子位置的变化，最后更新该粒子的速度和位置。
③将上一步对流结果代入方程(8)计算出粒子扩散过程浓度变化，并更新粒子浓度。
④若有粒子运动至计算域外则删除该粒子，并在入口边界补充新的粒子。
⑤循环计算域内的每一个流体粒子，重复做②③④中的操作。
步骤四，输出结果：
①每完成一个时间步的计算就更新其结果。
②完成时间步的循环，输出最终结果。
尽管上面结合图对本发明进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨的情况下，还可以做出很多变形，这些均属于本发明的保护之内。

Claims (3)

1.一种基于算子分裂及改进半拉格朗日求解污染物传播的方法，其特征是，步骤如下：
步骤一，初始化系统的相关变量和运行参数；
步骤二，生成粒子信息；
步骤三，列出求解方程并迭代计算：
根据算子分裂思想将对流、扩散和反应三个过程分别采用不同的数值方法计算；
步骤四，输出结果：
①每一个时间步结束，保存该时间步的计算结果并更新；
②结束时间循环，输出最终结果；
步骤三包括如下具体步骤：
求解的一维对流扩散反应方程为：
其中，C为物质浓度，x和t分别为空间和时间坐标，u为流体流动速度，ε为物质扩散系数，G(C)为反应源项，代表单位时间单位体积空间内因化学反应生成的组分量，当G(C)＝0时，方程(1)即为对流扩散方程：
其中，等号右侧为扩散项，一维纯对流方程的欧拉形式为：
根据运动流体的物质导数将方程(3)转化为拉格朗日形式：
从方程(4)和(5)，在拉格朗日体系下，对流过程中物质不会随时间变化而销毁或创造，仅仅是随流体发生运动，解得速度之后，求得浓度在时间和空间上的变化规律；
改进拉格朗日方法针对速度场采用合适的数值插值方法插值粒子速度，插值表达式为：
Vh＝a1Vj-1+a2Vj+a3VXj-1+a4VXj (6)
其中，Vh为当前时间步粒子速度，系数a1＝α2(3-2α)，a2＝1-a1，a3＝α2(1-α)Δx，a4＝-α(1-α)2Δx，α为柯朗数的小数部分，Vj-1和Vj分别为所求粒子的相邻速度网格点的速度，VXj-1和VXj为相邻速度网格点速度导数；
一维纯扩散方程为：
将方程(7)中的时间离散通过向前差分，空间离散使用二阶中心差分得到：
其中n表示当前时刻，n-1表示上一时刻，j-1，j，j+1表示空间网格位置，为对流过程结果，得扩散结果：
最后选用欧拉法Euler、龙格-库塔二阶法RK2或龙格-库塔四阶法RK4求解反应项；
用RK4进行时间积分，离散形式如下：
通过方程(6)、(9)和(10)分别计算出每个时间步下每个流体粒子的速度和浓度，进而求解出粒子新的位置。
2.如权利要求1所述的基于算子分裂及改进半拉格朗日求解污染物传播的方法，其特征是，通过方程(6)、(9)和(10)分别计算出每个时间步下每个流体粒子的速度和浓度，进而求解出粒子新的位置，具体计算过程为：
①循环每个时间步；
②初始化计算域中的粒子后，根据已知速度场，首先通过方程(6)计算目标粒子的速度，再根据速度计算粒子位置的变化，最后更新该粒子的速度和位置；
③将上一步对流结果代入方程(9)计算出粒子扩散过程浓度变化，并更新粒子浓度；
④通过方程(10)计算反应过程粒子浓度变化，并更新粒子浓度；
⑤若有粒子运动至计算域外则删除该粒子，并在入口边界补充新的粒子；
循环计算域内的每一个流体粒子，重复做②③④⑤中的操作。
3.如权利要求2所述的基于算子分裂及改进半拉格朗日求解污染物传播的方法，其特征是，对于计算域是长度为1的一维管道，管道阀门位于0处，管道内有浓度分别为1，0，0的三种反应物u1，u2，u3，此时阀门处于关闭状态，管道内混合物处于稳态,管道外u1，u2，u3浓度分别为0，1，0，求解当阀门打开后，流体以速度1发生流动，计算域内反应物u2浓度变化：
步骤一，初始化，初始化系统的相关变量和运行参数，具体包括：
模拟的计算域为阀门右侧x＝[0,1]的一维矩形空间，内有三种处于稳定状态的反应物，三种物质粒子间距Δx＝0.0005，计算时间步长为0.0002，计算时间为0.8，阀门左侧三种反应物浓度分别为u1＝0，u2＝1，u3＝0，打开阀门后，流体以速度1开始流动，三种反应物开始相互反应；
步骤二，生成粒子信息，具体包括：
初始化粒子步骤中，三种反应物分别生成流体粒子2000个，初始浓度分别为u1＝1，u2＝0，u3＝0；
步骤三，列出求解方程并迭代计算；
求解的一维对流反应方程拉格朗日形式为：
其中，U1，U2，U3为反应物浓度，x和t分别为空间和时间坐标，V为流体流动速度，方程(12)(13)(14)等号右侧均为反应项，代表单位时间单位体积空间内因化学反应生成的组分量，一维纯对流方程的欧拉形式为：
根据运动流体的物质导数可以将方程(3)转化为拉格朗日形式：
从方程(4)和(5)看出，解得速度之后，求得浓度在时间和空间上的变化规律；改进拉格朗日方法针对速度场采用合适的数值插值方法插值粒子速度，插值表达式为：
Vh＝a1Vj-1+a2Vj+a3VXj-1+a4VXj (6)
其中，Vh为当前时间步粒子速度，系数a1＝α2(3-2α)，a2＝1-a1，a3＝α2(1-α)Δx，a4＝-α(1-α)2Δx，α为柯朗数的小数部分，Vj-1和Vj分别为所求粒子的相邻速度网格点的速度，VXj-1和VXj为相邻速度网格点速度导数；
采用RK4求解反应项，离散形式如下：
通过方程(11)、(12)、(13)、(6)和(10)分别计算出每个时间步下每个流体粒子的速度和浓度，进而求解出粒子新的位置。
CN201910551021.5A 2019-06-24 2019-06-24 基于算子分裂及改进半拉格朗日求解污染物传播的方法 Active CN110287590B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910551021.5A CN110287590B (zh) 2019-06-24 2019-06-24 基于算子分裂及改进半拉格朗日求解污染物传播的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910551021.5A CN110287590B (zh) 2019-06-24 2019-06-24 基于算子分裂及改进半拉格朗日求解污染物传播的方法

Publications (2)

Publication Number Publication Date
CN110287590A CN110287590A (zh) 2019-09-27
CN110287590B true CN110287590B (zh) 2023-08-18

Family

ID=68005454

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910551021.5A Active CN110287590B (zh) 2019-06-24 2019-06-24 基于算子分裂及改进半拉格朗日求解污染物传播的方法

Country Status (1)

Country Link
CN (1) CN110287590B (zh)

Families Citing this family (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110705185A (zh) * 2019-09-29 2020-01-17 天津大学 预测管道气锤的方法
CN111241742B (zh) * 2019-12-27 2021-11-19 西安交通大学 一种多相流计算方法
CN111191372B (zh) * 2020-01-01 2021-12-21 浙江大学 一种基于二维水动力学的大规模洪水场景模拟预警交互方法和系统
CN111650672B (zh) * 2020-05-22 2022-04-01 西北核技术研究院 采用时间堆栈实现空气污染物大气扩散快速预测的方法
CN112580770A (zh) * 2020-12-21 2021-03-30 西安交通大学 一种基于粒子法的变粒径分裂方法
CN117275756B (zh) * 2023-08-25 2024-06-04 中国科学院地理科学与资源研究所 基于hasm的传染病时空扩散模拟方法和系统
CN117331361B (zh) * 2023-12-01 2024-03-08 江苏迈鼎科技(集团)有限公司 一种基于物联网技术的沥青生产管理系统
CN118278257A (zh) * 2024-03-07 2024-07-02 南方海洋科学与工程广东省实验室(广州) 一种全局模拟水体污染物滞留特征的sph并行计算方法

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101754845A (zh) * 2007-07-02 2010-06-23 马格马铸造工艺有限公司 用于在模具填充过程的模拟中描述粒子统计取向分布的方法和装置
JP2015141122A (ja) * 2014-01-29 2015-08-03 株式会社奥村組 汚染物質の移流拡散解析方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7930155B2 (en) * 2008-04-22 2011-04-19 Seiko Epson Corporation Mass conserving algorithm for solving a solute advection diffusion equation inside an evaporating droplet
KR101453774B1 (ko) * 2013-06-12 2014-10-22 한국지질자원연구원 개선된 라그랑지안-율러리안 기법을 이용한 총유량 경계조건을 가지는 오염물 거동의 분석방법

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101754845A (zh) * 2007-07-02 2010-06-23 马格马铸造工艺有限公司 用于在模具填充过程的模拟中描述粒子统计取向分布的方法和装置
JP2015141122A (ja) * 2014-01-29 2015-08-03 株式会社奥村組 汚染物質の移流拡散解析方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
"基于 N-S 方程的烟雾模拟实时性改进";李晓艳;《中国优秀硕士学位论文全文数据库 信息科技辑》;20100815;第2章-第4章 *

Also Published As

Publication number Publication date
CN110287590A (zh) 2019-09-27

Similar Documents

Publication Publication Date Title
CN110287590B (zh) 基于算子分裂及改进半拉格朗日求解污染物传播的方法
Stam Stable fluids
Thürey et al. Detail-preserving fluid control
Treuille et al. Keyframe control of smoke simulations
US7565276B2 (en) Method of simulating detailed movements of fluids using derivative particles
Nielsen et al. A two-continua approach to Eulerian simulation of water spray
Coros et al. Synthesis of constrained walking skills
US8296112B2 (en) Flow simulation method, flow simulation system, and computer program product
US20080319722A1 (en) Water Particle Manipulation
EP2551825A2 (en) Fluid dynamics framework for animated special effects
US10147220B2 (en) Precomputing data for an interactive system having discrete control inputs
KR102026154B1 (ko) 천수흐름에서 천해파의 수치모의 방법
Kim et al. Fire sprite animation using fire-flake texture and artificial motion blur
US7460984B1 (en) Compensating for delay in modeling environments
CN105183965B (zh) 用于预测雾化过程的大涡模拟方法
Wendt et al. Finite volume flow simulations on arbitrary domains
Makungu et al. A generalized 1-dimensional particle transport method for convection diffusion reaction model
US20240126955A1 (en) Physics-Informed Machine Learning Model-Based Corrector for Deformation-Based Fluid Control
Stuyck et al. Model predictive control for robust art-directable fluids
KR101538141B1 (ko) 기체 밀도의 보정을 수행하는 유체 시뮬레이션 장치 및 방법
KR100705417B1 (ko) 물체 표면의 젖음과 마름 표현 장치 및 방법
Bonilla et al. Control methods for fluid-based image warping
CN113673140B (zh) 一种气压作用下的流体粒子实时交互视觉仿真方法
Huang et al. A local adaptive Catmull–Rom to reduce numerical dissipation of semi‐Lagrangian advection
Gowthaman et al. Time Reversal and Simulation Merging for Target-Driven Fluid Animation.

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