CN101853334A - 化学反应扩散中的免方程多尺度模拟方法 - Google Patents

Abstract

本发明公开一种化学反应扩散中的免方程多尺度模拟方法，从介观尺度出发建立能真实逼近反应扩散系统的细小尺度模型，即格子玻尔兹曼模型；分别对Schlogl和Selkov反应扩散系统构造粗粒化时间步进；采用传统数值技术中的欧拉法与Schlogl和Selkov反应扩散系统的粗粒化时间步进耦合成Schlogl和Selkov反应扩散系统的粗粒化投影积分，从而实现不同投影时长和演化时长的浓度常稳态多尺度模拟，并与格子玻尔兹曼直接模拟进行对比；其中，采用最小二乘法对宏观系统浓度值进行拟合。本发明成功地回避了从细小尺度模型中获得反应扩散模型方程这一过程，并提高了格子玻尔兹曼直接模拟的运算效率，为在宏观尺度上对方程解析困难、在细小尺度上直接模拟计算量巨大的反应扩散仿真研究提供了一种新的思路。

Description

化学反应扩散中的免方程多尺度模拟方法

技术领域

本发明涉及一种新型的复杂系统多尺度模拟框架——免方程(Equation-free)方法，尤其涉及化学反应扩散问题中的免方程多尺度模拟方法。属于复杂系统多尺度模拟技术领域。

背景技术

反应扩散系统是自然界中典型的复杂系统，它的应用范围覆盖了许多学科，如生态系统中捕食者猎物模型，物理系统中气体放电模型，舆论传播模型等。化学反应就是一种典型的反应扩散系统，取不同的控制参数，系统能呈现出不同的时空行为，如双稳性、化学波、非均匀定态、化学振荡等复杂现象，所表现出来的规律与自然界中广泛存在的生物体有序结构、自组织现象的规律有共同之处，可以作为自然界中反应扩散系统的典型代表。同时化学反应过程中的控制参数方便调控，便于探索其变化规律，因而化学反应扩散的研究愈来愈受关注，不少研究者在这方面做了大量工作。因此，对化学反应扩散问题的研究无疑极大地有助于研究化学反应系统中的复杂现象，如极限环、图灵斑图、螺旋波的产生与控制等。

大多数化学反应扩散系统都具有大量交互成分，其内部关联复杂、不确定、总体行为具有时空多尺度特性，即不能通过系统的局部特性，形式地或者抽象地描述整个系统特性。研究这类具有时-空多尺度特性的系统，多尺度模拟是一个重要的、甚至是唯一的研究手段。反应扩散问题的传统模拟方法主要依赖于单一尺度的模型方程，如宏观尺度的反应扩散动力学方程。这方面研究从多，并取得了很多研究成果。然而，传统模拟方法存在以下问题：一、所得结果不仅会掩盖宏观尺度下的结构效应，还会抹平更大尺度上的结构效应，从而造成显著的误差。二、值得注意的是，由于这些描述反应扩散的演化方程都是在人为假设和理想约束的条件下通过反应扩散的细小尺度模型(fine scale models，如微观模型、介观模型和随机模型等)演绎推导出来的，因此它们表达或刻画出来的理想系统不能较有效地逼近真实系统，并且单从方程本身来考虑根本无法透析其系统的内在机理，从而给化学反应扩散系统的进一步研究分析带来了无法跨越的鸿沟。三、反应扩散问题的宏观模型方程多数属于非线性偏微分方程组，这种方程的数值解析往往非常困难甚至不可解。因此，对于这些反应扩散系统，传统单尺度方法已很难适用，需要寻求新的以多尺度模拟为基础的研究策略。

近几年来，化学反应扩散系统在介观尺度上的研究(即格子玻尔兹曼模型)发展迅速。对于化学反应，用格子玻尔兹曼方法可以实现模型的构造，且具有精度高、程序代码简单、稳定性好等特点。特别是对二维和三维问题，虽然空间复杂性大大增加，但算法的复杂性并没有增加。另外，该方法有很强的物理背景，对于无对流的反应系统，表现为粒子的动量是对称的。而对不同的组分，如需要的话，可以有不同的能量和分布。格子玻尔兹曼从介观尺度的机理出发逐步揭示出宏观尺度上反应扩散现象的机理，将介观尺度上流体分子微团的各种特性逐步反映到宏观尺度的现象中，从而对实现多尺度模拟起到了非常重要的作用。自从20世纪90年代初由Qian等人和Chen等人提出单弛豫时间近似代替碰撞项以来获得了广泛的应用，在很多领域取得了成功。目前，较成功地用于反应扩散问题研究的格子玻尔兹曼方法是1993年S.P.Dawson等人提出的，其方法是在流行的格子玻尔兹曼的碰撞项中加入反应碰撞项，其中反应碰撞项以参加反应的粒子的密度作为自变量Rs(ρ)。利用这种格子玻尔兹曼，S.P.Dawson等人从理论上导出了宏观的反应扩散方程，并且所用的反应项直观简捷。总的来说，相对于反应扩散动力学方程，反应扩散系统的格子玻尔兹曼模型能够较简便地获取和更有效地表达。此外，与实验研究相比它们还具有能随意初始化的优点，这样就为反应扩散系统的多尺度研究开辟了一条新的道路。

由于反应扩散系统的绝大多数相关研究工作如获取系统宏观特性、预测系统趋势、系统优化等都是在宏观水平上进行的，这样就需要把介观尺度上的信息跨尺度地传输到宏观尺度水平。目前，虽然格子玻尔兹曼的传统直接模拟已可在较简单的水平上实现跨尺度传输，并且可以在模拟过程上直接地观察到系统变化，但由于缺少运用相关的数值技术，对于系统特性的获取往往表现得很迟钝，并且它们很难完成甚至不可能实现一些要求较高的研究工作如系统分岔分析、系统优化、控制系统设计等。另外更关键的是，设想在处理拥有巨大运算量、更复杂的新型多组分反应扩散系统时，它们的运算效率必然很低，甚至有时会因运算量过大而导致无法运行，从而使得它们根本不能有效地达到系统研究的预期目的。本发明将考虑采用以细小尺度模型(即反应扩散系统的格子玻尔兹曼模型)为基础的新型多尺度模拟框架——免方程方法，并将其分别应用在Schlogl和Selkov两类反应扩散问题的模拟研究中。

发明内容

本发明的目的是，克服传统直接模拟方法在获取复杂系统特性方面的本质缺陷，首次应用免方程方法中的粗粒化投影积分(coarse projective integration)分别对Schlogl和Selkov两类反应扩散问题中的浓度常稳态进行多尺度模拟，为研究拥有巨大运算量、更复杂的新型多组分反应扩散系统打下坚实的基础。

本发明的技术方案为：反应扩散问题的免方程多尺度模拟方法，其特征在于对Schlogl和Selkov反应扩散问题中的浓度常稳态采用免方程方法中的粗粒化投影积分进行多尺度模拟。其模拟方法包括如下步骤：

(一)结合本研究关于Schlogl、Selkov反应扩散系统特性的要求，对Schlogl、Selkov反应扩散系统进行免方程分析，并对反应扩散系统进行尺度划分，即宏观尺度和介观尺度。

(二)根据Schlogl、Selkov反应扩散系统在宏观尺度和介观尺度上的时空特性，筛选出宏观系统量(即反应物浓度)和介观量离散分布(即粒子密度分布函数)并从中分析得出具体的提升表达式和约束表达式。从系统的底层(介观尺度)出发建立能真实逼近Schlogl和Selkov反应扩散系统的细小尺度模型，即格子玻尔兹曼模型。

(三)构造Schlogl、Selkov反应扩散系统的粗粒化时间步进(coarse timestepper)。该粗粒化时间步进是由Schlogl、Selkov反应扩散系统的格子玻尔兹曼模型中的演化、提升和约束构成。

(四)采用传统数值技术中的欧拉法与Schlogl和Selkov反应扩散系统的粗粒化时间步进耦合成Schlogl和Selkov反应扩散系统的粗粒化投影积分，从而实现不同投影时长和演化时长的浓度常稳态多尺度模拟，并与格子玻尔兹曼直接模拟进行对比。其中，采用最小二乘法对宏观系统浓度值进行拟合。

本发明所提供的免方程多尺度模拟方法，具有如下有益效果：

结果分析表明，本发明能够快而准地获得Schlogl、Selkov反应扩散的浓度常稳态并且能优化格子玻尔兹曼直接模拟的运算效率。

1)由于本发明采用的是粗粒化投影积分，在模拟过程中成功地回避了从细小尺度模型中获得反应扩散模型方程这一过程，为复杂系统多尺度模拟提供了一种新的思路。

2)由于本发明的应用研究着眼于从介观尺度上分析反应扩散系统，使得模拟结果能够更真实且有效地反映该复杂系统的本质特性，即浓度常稳态。

3)又由于粗粒化投影积分中能融入了欧拉法等外推数值技术，从而捕获反应扩散系统的浓度常稳态表现得更快捷，更重要的是在优化格子玻尔兹曼直接模拟的运算效率方面它是一种极为有效的方法。

本发明受到国家自然科学基金项目“复杂能源系统多尺度优化方法研究(资助号：50876117)”的资助。

附图说明

图1为传统直接模拟流程图；

图2为本发明免方程多尺度模拟的基本框架流程图；

图3为粗粒化投影积分的基本原理图；

图4为Schlogl反应扩散问题的第一项实验结果图；

图5为Schlogl反应扩散问题的第二项实验结果图；

图6为Selkov反应扩散问题的第一项实验结果图；

图7为Selkov反应扩散问题的第二项实验结果图；

表1为Schlogl反应扩散问题的各实验运算时间；

表2为Selkov反应扩散问题的各实验运算时间；

具体实施方式

下面结合附图和具体实施方式对本发明作进一步说明。

如图2所示，一种反应扩散问题的免方程多尺度模拟方法，对Schlogl和Selkov反应扩散中的浓度常稳态采用免方程方法中的粗粒化投影积分进行多尺度模拟。

其模拟方法包括如下步骤：

(一)以免方程(Equation-free)为基本模拟框架，针对化学反应扩散系统的特性，将系统空间尺度划分为宏观尺度和介观尺度，以此为基础建立免方程多尺度模拟方法，实现通过细小尺度的模拟得出系统宏观性能。

(二)根据Schlogl、Selkov反应扩散系统在宏观尺度和介观尺度上的时空特性，筛选出宏观系统量(即反应物浓度)和介观量离散分布(即粒子密度分布函数)并从中分析得出具体的提升表达式和约束表达式。从系统的底层(介观尺度)出发建立能真实逼近Schlogl和Selkov反应扩散系统的细小尺度模型，即格子玻尔兹曼模型。

(三)分别对Schlogl和Selkov反应扩散系统构造粗粒化时间步进(coarse timestepper)。该粗粒化时间步进是由Schlogl和Selkov反应扩散系统的格子玻尔兹曼模型中的演化、提升和约束构成。

(四)采用传统数值技术中的欧拉法与Schlogl和Selkov反应扩散系统的粗粒化时间步进耦合成Schlogl和Selkov反应扩散系统的粗粒化投影积分，从而实现不同投影时长和演化时长的浓度常稳态多尺度模拟，并与格子玻尔兹曼直接模拟进行对比。其中，采用最小二乘法对宏观系统浓度值进行拟合。

进一步，在步骤(一)中，免方程的基本框架由能跨越时间尺度差的粗粒化时间步进和与之相应的数值技术构成，其中，粗粒化时间步进相当于整个框架中的“内核”作用，类似于系统辨识中的“黑箱”，而数值技术的作用则可看似粗粒化时间步进的“外壳”。此处粗粒化时间步进是一个由提升、演化和约束组成的跨尺度组合，它解决了细小尺度模型和宏观尺度之间的尺度限制，框架中的数值技术要根据系统所要求的研究任务进行筛选；基本框架的一般实现过程首先需要做两个准备：

a)根据具体的研究任务对系统进行尺度划分，即宏观尺度和细小尺度，在相应的细小尺度上建立模型并根据不同尺度的时空特性筛选出具体的物理量和粒子离散分布分别作为宏观系统量u和细小量离散分布U。

b)根据宏观系统量u和细小量离散分布U之间的映射关系分析得出具体的提升表达式，则提升μ可表达为：U＝μu；再由已得到的提升μ推导出合适的约束表达式，约束M可表达为：u＝MU。其中提升μ与约束M互逆：μM＝I。

完成以上准备工作后就可按如下描述进行：

1)对宏观系统量设定一个基本符合研究要求的初值u(t₀)，t₀为系统的初始时刻。

2)先通过提升μ把初值u(t₀)离散转化为细小量离散分布，具体表达式：U(t₀)＝μu(t₀)。

3)再将细小量离散分布U(t₀)投入细小尺度模型中，运行演化时步(即步长)δt得到新的细小量离散分布，其过程可表达为：

其中为演化算子。

4)最后，把新得到的细小量离散分布函数通过约束M得到宏观系统量的新值，其过程具体表达为：u(t₀+δt)＝MU(t₀+δt)。把刚获得的新值返回2步中继续运行。

以上的2)、3)、4)步就是粗粒化时间步进的三个基本步骤。把这三步的表达形式结合起来，我们可以推导出粗粒化时间步进的一般表达式为：

接下来，再把从连续几次粗粒化时间步进中得到的一系列宏观系统值抽离出来，用所选的数值技术如Euler法、Runge-Kutta法、Newton-Raphson法等对它们进行相应的数值处理，然后把所得到的处理值重新输入到步骤2中运行。整个模拟过程在不停地循环运行中达到对系统研究分析的目的。另外，此基本框架能够通过适当的修正和补充扩展到空间尺度，甚至是时空尺度。

对化学反应扩散系统进行尺度划分，即宏观尺度和介观尺度。在介观尺度上建立能真实逼近反应扩散系统的细小尺度模型，再以此为基础试构造反应扩散系统的粗粒化时间步进。根据反应扩散系统特性的研究要求(即反应扩散浓度常稳态问题)选用外推数值方法与反应扩散系统的粗粒化时间步进耦合，从而实现不同投影时长和演化时长的粗粒化投影积分多尺度模拟。

进一步，在步骤(二)中，选用反应扩散格子玻尔兹曼模型分别作为反应扩散系统的细小尺度模型，对于有反应源项的扩散方程，基本思想是在碰撞项中加入反应项具体形式如下：

$f_{s,i}(x+c_i\mathrm{\Δt},t+\mathrm{\Δt})-f_{s,i}(x,t)=-1/{\mathrm{\τ}}_s[f_{s,i}(x,t)-f_{s,i}^{\mathrm{eq}}(x,t)]+\mathrm{\Δt}{R}_s/b$

上式可作为Schlogl反应扩散的细小尺度模型；对上式采用多尺度方法得到反应扩散宏观方程：

$\frac{{\∂\mathrm{\ρ}}_s}{\∂t}-D_s{\▿}^2{\mathrm{\ρ}}_s=R_s$

其中f_s，i(x，t)为s组分时刻t位于x处具有速度c_i的粒子密度分布函数；

为s组分的局部平衡态分布函数；ρ_s为反应物浓度；D_s为扩散系数；τ_s为松弛时间；b为格子模型中的速度数；R_s为反应扩散方程里的反应项，由具体化学反应决定。

本发明采用格子玻尔兹曼D2Q5模型对Schlogl反应扩散系统建立细小尺度模型，格子玻尔兹曼D1Q3模型对Selkov反应扩散系统建立细小尺度模型，模拟结果均在处理器为AMDAthlon(tm)，主频为2.01GHz，操作环境是Windows XP的计算机上应用MATLAB编程运行得到。

进一步，在步骤(三)中，粗粒化时间步进是免方程方法基本框架的核心部分，类似于系统辨识中的“黑箱”。在实施免方程多尺度模拟时，构造合理有效且能真实逼近复杂系统的粗粒化时间步进是一个重要的前提。化学反应扩散的粗粒化时间步进可分别由在反应扩散细小尺度模型中的演化过程和已提炼的提升、约束构成。本发明中反应扩散的粗粒化时间步进以格子玻尔兹曼模型为仿真器，结合介观尺度信息—宏观尺度的交互耦合机制，从而实现了对于反应扩散系统的跨尺度组态，它担负着透析反应扩散系统和映射系统底层至宏观特性的任务。

进一步，在步骤(四)中，粗粒化投影积分(可简称CPI)是免方程方法中一种具有代表性的基础方法。简单地讲就是将基本框架中的某种选定的数值方法进行外推时所呈现出来的一种具体的免方程方法，其基本思想是把粗粒化时间步进获取的宏观系统值外推投影从而实现捕获系统特性和预测系统的目的。其中较常见的是用Euler法进行外推处理，其过程可表示为：

$u^{0,N+1}=u^{k,N}+(\mathrm{\ΔT}-t_k)\tilde{F}\left(u^{k,N}\right)$

其中，

可通过已获得的宏观系统值得到：

$\tilde{F}\left(u^{k,N}\right)=\frac{u^{k+1,N}-u^{k,N}}{\mathrm{\δt}}\≈\frac{\mathrm{du}(\mathrm{N\ΔT}+t_k)}{\mathrm{dt}}=F\left(u(\mathrm{N\ΔT}+t_k)\right)\≈F\left(u^{k,N}\right)$

在具体应用中可以根据研究目的地要求灵活地处理

在本应用研究中采用最小二乘拟合。最后，再把投影获得的u^0，N+1输入到下一个模拟阶段中继续按照此模式运行直到能够满足研究要求。此方法称为投影向前欧拉方法，它是粗粒化投影积分中的一种，能够很清晰地表达粗粒化投影积分的基本原理。其中，某个阶段起始点的宏观系统量初值为u^0，N，在这它表示经过N个模拟阶段时的模拟值，且每个阶段的宏观时间长为ΔT，因此宏观系统量初值还可写成u(NΔT+0)。t_k＝kδt，相当于k个演化时步，显然ΔT＞＞t_k。另外，为了方便理解我们把包含外推过程的模拟阶段时间长ΔT称为投影时长，而在投影时长中需要在细小尺度模型中演化运算的时间长t_k称为演化时长。

本发明实施反应扩散浓度常稳态的粗粒化投影积分模拟，在化学反应扩散系统的粗粒化时间步进上耦合宏观尺度的传统外推数值方法，达到对化学反应扩散系统特性快速而准确地模拟，同时大大提高了格子玻尔兹曼直接模拟的运算效率。

下面以对Schlogl、Selkov反应扩散中的浓度常稳态进行多尺度模拟为例，并与格子玻尔兹曼直接模拟进行对比，进一步说明该模拟方法包括反应扩散系统的免方程分析、建立反应扩散的细小尺度模型、构造反应扩散的粗粒化时间步进和实施粗粒化投影积分多尺度模拟。

(一)反应扩散系统的免方程分析

1.免方程方法的基本框架

近几年国外提出了一种计算机辅助分析复杂系统的新型多尺度模拟框架——免方程方法。它能在细小尺度模型上运用以往基于宏观模型方程的传统数值技术直接达到对复杂系统的研究分析。其中的关键思路是如何绕过从细小尺度模型中获得宏观模型方程这一传统过程，为复杂系统多尺度模拟提供了一种新的思路。免方程多尺度模拟中主要针对这样一类情形，即所研究系统的机理可以在细小尺度水平模拟，然而其宏观模拟方程则难以从细小水平规则中完整地获得。这种情形普遍存在从材料科学到社会生态学，从工程科学到计算化学的众多模拟问题中，在这类问题中，机理在细小尺度是已知的，而要将这些机理映射到粗粒化的宏观水平则无法实现。人们有时甚至不知道在什么样的尺度上能够有效地预测其机理。而试图进行直接计算机仿真，在细小尺度和宏观系统尺度之间的巨大鸿沟之间建立联系则受到计算机能力的极大限制，这种限制是当今复杂系统仿真的一个主要障碍。显然若能通过以细小尺度模拟来执行宏观系统尺度的研究分析工作，即通过细小的时间和空间尺度上适当初始化仿真器并重复调用该仿真器以获取在宏观时间和空间尺度上系统的信息，它是多尺度模拟和优化研究中的前沿性问题。此多尺度模拟框架能够避免复杂的详细宏观模型方程的推导(这些方程在概念上存在，但无法描述或不能得到封闭的方程组)。其基本思想在于利用对物理问题的细小水平仿真(如格子玻尔兹曼方法、动力学蒙特卡洛方法、以及分子动力学模拟方法等)完成宏观尺度的研究分析工作实际上。这种多尺度模拟框架的本质相当于系统辨析，通过所谓约束(restriction，通过平均、筛选、光滑处理从细小尺度模型获得粗糙宏观状态集)和提升过程(lifting，借助于适当的分布函数从粗糙集中获得细小量粒子分布)建立起细小尺度模型和传统连续科学计算、数值分析的“桥梁”。目前免方程的具体方法中有以下四种：

(1)粗粒化分岔分析(coarse bifurcation analysis)，Shroff-Keller回归投影方法的延伸，传统数值分岔分析和粗粒化时间步进(coarse timestepper)的结合；

(2)粗粒化投影积分，传统数值投影积分和粗粒化时间步进(coarse timestepper)的结合；

(3)区域化动态学(patch dynamics)，存在有时空区域的齿隙(gap-tooth)格式，细小尺度模型在很小的区域(the teeth)中求解，这些区域是由较大的“隙”(gap)分隔而成的，粗粒化投影积分和齿隙格式的结合；

(4)区域化分岔分析(patch bifurcation analysis)。粗粒化分岔分析和齿隙格式的结合。

目前，免方程方法常见的应用是时间尺度上的多尺度模拟，因此一般可认为在时间尺度上的模拟框架为其基本框架。简单地讲该基本框架由能跨越时间尺度差的粗粒化时间步进和与之相应的数值技术构成，其中，粗粒化时间步进相当于整个框架中的“内核”作用，类似于系统辨识中的“黑箱”，而数值技术的作用则可看似粗粒化时间步进的“外壳”。此处粗粒化时间步进是一个由提升(lifting)、演化(evolution)和约束(restriction)组成的跨尺度组合，它解决了细小尺度模型和宏观尺度之间的尺度限制。其中，提升是该组合中首先需实现的过程，它是一种跨尺度过程，能够把宏观系统量离散转化成与之相匹配的细小量离散分布，其精度直接决定着粗粒化时间步进、甚至整个框架的模拟效果，例如在气体动态模拟中，用提升可把某气体的浓度离散成大量气体粒子的随机离散分布，如果提升的精度不适合研究要求势必将影响整个动态研究。接下来要进行的就是演化过程，在粗粒化时间步进中只有演化不需要跨越尺度差，这是由于它的整个过程都是在细小尺度上进行的。它把从提升中得到的细小量离散分布投入到细小尺度模型中进行运算，经过一些短的演化时间后从中能得到新的细小量离散分布。细小尺度模型在这就相当于一个在细小水平上能进行自行演化计算的仿真器。最后，要执行的就是约束，它是粗粒化时间步进中的最后一个过程，其要达到的目的恰好与提升相反，即把新得到的细小量离散分布重新反馈到宏观水平，该过程在整个粗粒化时间步进中同样具有关键作用。

对于框架中的数值技术，我们要根据系统所要求的研究任务进行筛选。如果粗粒化时间步进能够配合到相应有效的数值技术，那么就可使该基本框架的模拟能够很巧妙地完成以前那些直接模拟不能有效达到甚至不可能实现的系统研究目的。另外值得强调的是，免方程方法并不是一种具体的多尺度模拟方法，而是一种多尺度模拟框架。其中具有代表性的免方程方法有粗粒化投影积分和粗粒化分岔分析(coarse bifurcation analysis)。再则提升和约束这两个过程的价值主要体现在它们能够将信息跨尺度传输，因此它们相当于在数值技术和细小尺度模型之间的架起了一座跨尺度“桥梁”。以下为基本框架的一般实现过程：

首先需要做两个准备：

a)先结合具体的研究任务对系统进行尺度划分，即宏观尺度和细小尺度，再在合适的细小尺度上建立细小尺度模型并根据不同尺度的时空特性筛选出具体的物理量和粒子离散分布分别作为宏观系统量u和细小量离散分布U。

b)根据宏观系统量u和细小量离散分布U之间的映射关系分析得出具体的提升表达式，则提升μ可表达为：U＝μu。再由已得到的提升μ推导出合适的约束表达式，约束M可表达为：u＝MU。其中提升μ与约束M互逆：μM＝I。

完成以上准备工作后就可按如下描述进行：

1)宏观系统量设定一个基本符合研究要求的初值u(t₀)，t₀为系统的初始时刻。

2)先通过提升μ把初值u(t₀)离散转化为细小量离散分布，具体表达式：U(t₀)＝μu(t₀)。

3)再将细小量离散分布U(t₀)投入细小尺度模型中，运行演化时步(即步长)δt得到新的细小量离散分布，其过程可表达为：

其中

为演化算子。

4)最后，把新得到的细小量离散分布函数通过约束M得到宏观系统量的新值，其过程具体表达为：u(t₀+δt)＝MU(t₀+δt)。把刚获得的新值返回2步中继续运行。

以上第2)、3)、4)步就是粗粒化时间步进的三个基本步骤。把这三步的表达形式结合起来，我们可以推导出粗粒化时间步进的一般表达式为：

接下来，再把从连续几次粗粒化时间步进中得到的一系列宏观系统值抽离出来，用所选的数值技术如Euler法、Runge-Kutta法、Newton-Raphson法等对它们进行相应的数值处理，然后把所得到的处理值重新输入到第2)步中运行，最后，整个模拟过程就这样在不停地循环运行中达到对系统研究分析的目的，如图2所示(图中的数值技术假设只需一次粗粒化时间步进结果即可)。另外，此基本框架能够通过适当的修正和补充扩展到空间尺度，甚至是时空尺度。本发明将运用粗粒化投影积分对Schlogl、Selkov反应扩散浓度常稳态进行数值模拟(只考虑时间尺度)。

2.反应扩散系统

[1]Schlog反应扩散

Schlogl反应扩散模型为：

式(1)

式(2)

其中，设定式(1)、(2)的正反应速率常数分别为k₁和k₂，逆反应速率常数分别为k_-1和k_-2。用ρ_a和ρ_b分别表示反应物A和B的浓度，ρ₁则表示X的浓度。该模型有无量纲化的模型方程：

$\frac{{\∂\mathrm{\ρ}}_1}{\∂t}-D_1{\▿}^2{\mathrm{\ρ}}_1=R_1=k_2{\mathrm{\ρ}}_a{\mathrm{\ρ}}_1^2-k_{-1}{\mathrm{\ρ}}_1+k_1{\mathrm{\ρ}}_b-k_{-2}{\mathrm{\ρ}}_1^3,$ 式(3)

式中D₁为X的扩散系数。当保持各边界上ρ₁不变，模型方程参数取：

[k₂ρ_a，k_-1，k₁ρ_b，k_-2]＝[0.015，0.006875，0.0009375，0.01]，式(4)可解析得到常稳态值：

ρ₁(1)＝0.25，ρ₁(2)＝0.75。              式(5)

[2]Selkov反应扩散

Selkov反应扩散模型为：

式(6)

式(7)

式(8)

同样，设定式(6)、(7)、(8)的正反应速率常数分别为k₃、k₄和k₅，逆反应速率常数分别为k_-3、k_-4和k_-5。反应物C和D的浓度用ρ_c和ρ_d表示，M和N的浓度用ρ₂和ρ₃表示。因M和N都具有不同的扩散速率，该模型有无量纲化的模型方程：

$\frac{{\∂\mathrm{\ρ}}_2}{\∂t}-D_2{\▿}^2{\mathrm{\ρ}}_2=R_2=k_3{\mathrm{\ρ}}_c-k_{-3}{\mathrm{\ρ}}_2-k_4{\mathrm{\ρ}}_2{\mathrm{\ρ}}_3^2+k_{-4}{\mathrm{\ρ}}_3^3,$ 式(9)

$\frac{{\∂\mathrm{\ρ}}_3}{\∂t}-D_3{\▿}^2{\mathrm{\ρ}}_3=R_3=k_{-5}{\mathrm{\ρ}}_d-k_5{\mathrm{\ρ}}_3+k_4{\mathrm{\ρ}}_2{\mathrm{\ρ}}_3^2+k_{-4}{\mathrm{\ρ}}_3^3,$ 式(10)

式中D₂和D₃分别为M和N的扩散系数。当保持各边界上ρ₂和ρ₃不变，模型方程参数取：

[k₃ρ_c，k_-3，k₄，k_-4]＝[0.002656673，0.000665，0.015，0.015]，式(11)

[k_-5ρ_d，k₅，k₄，k_-4]＝[0.00053134，0.00665，0.015，0.015]，式(12)可求解得到常稳态值：

${\mathrm{\ρ}}_2^{*}=1.331141,$ ${\mathrm{\ρ}}_3^{*}=0.346286.$ 式(13)

对化学反应扩散系统进行尺度划分，即宏观尺度和介观尺度。在介观尺度上建立能真实逼近反应扩散系统的细小尺度模型，再以此为基础试构造反应扩散系统的粗粒化时间步进。根据反应扩散系统特性的研究要求(即反应扩散浓度常稳态问题)采用外推数值方法与反应扩散系统的粗粒化时间步进耦合，从而实现不同投影时长和演化时长的粗粒化投影积分多尺度模拟。

(二)建立反应扩散的细小尺度模型

本发明选用Schlogl、Selkov反应扩散的格子玻尔兹曼模型分别作为它们的细小尺度模型。基本描述如下：

格子玻尔兹曼方法是一种基于分子运动论和统计力学理论的流体计算方法。与以宏观连续方程为基础的传统计算流体力学方法不同，格子玻尔兹曼是基于流体介观模型的方法。与传统的计算流体力学方法相比，格子玻尔兹曼具有许多独特的优势，如边界条件容易实现、具有完全并行性等。目前，除了在一般的流体力学问题中得到了成功的应用外，格子玻尔兹曼已经在多相流、多孔介质流、悬浮粒子流、反应流、磁流体力学和生物力学等领域取得了很大的成功。格子玻尔兹曼方法是由McNamara和Zanetti于1988年提出的，这种模型中模拟的对象不是数目庞大的流体分子个体，而是数目大大减少的流体粒子，即微观充分大、宏观充分小的流体分子微团。格子玻尔兹曼从描述流体分子微团出发，有效的模拟了多个单粒子分布函数的流动和碰撞过程，并通过提升和约束将介观尺度的单粒子分布函数与宏观尺度的密度、速度、内能有机的结合起来，建立了不同层次之间参数的联系机制，实现了连接介观尺度和宏观尺度的作用，解决了不同层次模型耦合的问题，实现了从介观尺度向宏观尺度的跨越。格子玻尔兹曼从介观尺度的机理出发逐步揭示出宏观尺度上各种反应扩散现象的本质，将介观尺度上流体分子微团的各种特性逐步反映到宏观尺度的现象中，从而对实现多尺度模拟起到了非常重要的作用。

1)Schlogl反应扩散的细小尺度模型

采用D2Q5模型，有五个粒子离散速度：

$c_i=\left\{\begin{array}{c}\left(\mathrm{0,0}\right),i=0\\ (\cos (i-1)\mathrm{\π}/2,\sin (i-1)\mathrm{\π}/2)c,i=\mathrm{1,2,3,4}\end{array}\right.$ 式(14)

其中，可取c＝Δx/Δt＝1，i表示粒子离散速度方向。设f_s，i(x，t)为s组分时刻t位于x处具有速度c_i的粒子密度分布函数。在Schlogl反应扩散系统中组分序号s由“1”表示，在此s组分即为反应物X。

此时能给出扩散过程的演化方程为：

$f_{s,i}(x+c_i\mathrm{\Δt},t+\mathrm{\Δt})-f_{s,i}(x,t)=-1/{\mathrm{\τ}}_s[f_{s,i}(x,t)-f_{s,i}^{\mathrm{eq}}(x,t)].$ 式(15)

其中，方程式右边为碰撞项，它表示发生碰撞后f_s，i(x，t)的表化率。

为s组分的局部平衡态分布函数，τ_s为s组分的松弛时间。

为了正确反映扩散方程，必须适当选取局部平衡态分布函数，且使其满足质量守恒条件，s组分的浓度为：

${\mathrm{\ρ}}_s=\underset{i}{\mathrm{\Σ}}{f}_{s,i}=\underset{i}{\mathrm{\Σ}}{f}_{s,i}^{\mathrm{eq}}.$ 式(16)

假定反应扩散中无流动，可简单取s组分的局部平衡分布函数的形式为：

$f_{s,i}^{\mathrm{eq}}={\mathrm{\ρ}}_s/b.$ 式(17)

其中b表示粒子离散速度数，此细小尺度模型中可取5。

在此需指出，(16)式即为粗粒化时间步进的约束表达式，而(17)式则为提升表达式。然后采用多尺度方法可得到扩散过程的宏观方程：

$\frac{{\∂\mathrm{\ρ}}_s}{\∂t}-D_s{\▿}^2{\mathrm{\ρ}}_s=0,$ 式(18)

扩散系数计算式：

$D_s=c_s^2({\mathrm{\τ}}_s-0.5).$ 式(19)

其中

为模型参数，在此细小尺度模型取0.4。对于有反应源的扩散过程，基本思想是在碰撞项中加入反应项

具体形式如下：

$f_{s,i}(x+c_i\mathrm{\Δt},t+\mathrm{\Δt})-f_{s,i}(x,t)=-1/{\mathrm{\τ}}_s[f_{s,i}(x,t)-f_{s,i}^{\mathrm{eq}}(x,t)]+\mathrm{\Δt}{R}_s/b,$ 式(20)

上式可作为Schlogl反应扩散的细小尺度模型。同样对(20)采用多尺度方法得到反应扩散宏观方程：

$\frac{{\∂\mathrm{\ρ}}_s}{\∂t}-D_s{\▿}^2{\mathrm{\ρ}}_s=R_s.$ 式(21)

2)Selkov反应扩散的细小尺度模型

与上述Schlogl反应扩散的细小尺度模型相比。Selkov反应扩散的细小尺度模型除了有如下几个关键性的变化外其他部分都类似：

a)Selkov反应扩散采用的是格子玻尔兹曼D1Q3模型，有三个离散粒子速度：0，-1，1，b＝3。

b)Selkov反应扩散中有两个组分M和N，s表示成“2”和“3”。扩散系数计算式推导

为0.67。

c)Selkov反应扩散细小尺度模型中的两个子模型通过M和N两个组分反应项之间的耦合而联接在一起。

(三)构造反应扩散的粗粒化时间步进

粗粒化时间步进是免方程方法基本框架的核心部分，类似于系统辨识中的“黑箱”。在实施免方程多尺度模拟时，构造合理有效且能真实逼近复杂系统的粗粒化时间步进是一个重要的前提。它的准确程度直接决定着整个免方程多尺度模拟的研究效果。具体应用时应注意两点问题：第一，作为仿真器的细小尺度模型，不一定要是诸如蒙特卡洛模型、分子动力学模型、格子玻尔兹曼模型等常见模型，只要在细小尺度水平上能准确反映系统、满足研究目标要求的模型均可；第二，提升和约束与细小尺度模型的耦合匹配十分关键，它们并不是越精细越好用，关键要根据具体的研究目标而定，此外，信息跨尺度传输不会是一对一的映射，而是一对多或者多对一的映射。

本发明通过以上Schlogl、Selkov反应扩散系统各两个细小尺度模型的建立，可以很清晰地构造出Schlogl、Selkov反应扩散的粗粒化时间步进，它们可分别由在Schlogl、Selkov反应扩散细小尺度模型中的演化过程和上面已给出的提升、约束构成。但需注意，在本免方程应用研究中Schlogl、Selkov反应扩散的粗粒化时间步进是彼此独立的。

本发明中反应扩散的粗粒化时间步进以格子玻尔兹曼模型为仿真器，结合介观尺度信息—宏观尺度的交互耦合机制，从而实现了对于反应扩散系统的跨尺度组态，它担负着透析反应扩散系统和映射系统底层至宏观特性的任务。由于有在介观层面上研究反应扩散系统的独特背景，为反应扩散系统的粗粒化时间步进的可靠性打下了坚实的基础。

(四)实施粗粒化投影积分的多尺度模拟

粗粒化投影积分(可简称CPI)是免方程方法中一种具有代表性的基础方法。简单地讲就是将基本框架中的某种选定的数值方法进行外推时所呈现出来的一种具体的免方程方法，其基本思想是把粗粒化时间步进获取的宏观系统值外推投影从而实现捕获系统特性和预测系统的目的。其基本原理如下：

传统的直接模拟完全依赖于细小尺度模型，它们需依次完成整个模拟过程中的每个演化时步运算。其大致过程描述如下(如图1)：先假设把直接模拟过程按时间均分成若干个模拟阶段，其中某个阶段起始点的宏观系统量初值为u^0，N，在这它表示经过N个模拟阶段时的模拟值，且每个阶段的宏观时间长为ΔT，因此宏观系统量初值还可写成u(NΔT+0)。接下来依照基本模拟步骤：u^0，N通过提升可以得到：U^0，N＝μu^0，N，U^0，N通过该阶段的第一个演化时步δt可以得到：

U^1，N通过约束可以得到：u^1，N＝MU^1，N。再把u^1，N通过提升得到：U^1，N＝μu^1，N......U^k-1，N通过第k个演化时步δt得到：

U^k，N通过约束得到：u^k，N＝MU^k，N......就这样一直迭代运行直到模拟结束，其中t_k＝kδt，相当于k个演化时步。上述的直接模拟虽然可以在模拟过程中直接地观察到系统变化，但是对于系统特性的获取往往表现得很迟钝，并且它们很难完成甚至不可能实现一些要求较高的研究工作如系统分岔分析、系统优化、控制系统设计等。另外，更关键的是它的运算效率一般都较低，甚至有时会因运算量过大而导致无法运行。这是由于直接模拟的每个演化时步都必须真实地进行，所以即使时长为t_k的运算量不是很大，但在ΔT＞＞t_k的情况下，要想通过直接模拟获得u^0，N+1就会很困难。因此，人们考虑运用外推宏观系统值的方法来解决上述问题。其中较常见的是用Euler法进行外推处理，其过程可表示为：

$u^{0,N+1}=u^{k,N}+(\mathrm{\ΔT}-t_k)\tilde{F}\left(u^{k,N}\right).$ 式(22)

其中，

可通过已获得的宏观系统值得到：

$\tilde{F}\left(u^{k,N}\right)=\frac{u^{k+1,N}-u^{k,N}}{\mathrm{\δt}}\≈\frac{\mathrm{du}(\mathrm{N\ΔT}+t_k)}{\mathrm{dt}}=F\left(u(\mathrm{N\ΔT}+t_k)\right)\≈F\left(u^{k,N}\right),$ 式(23)

在具体应用中可以根据研究目的地要求灵活地处理

在本应用研究中采用最小二乘拟合。最后，再把投影获得的u^0，N+1输入到下一个模拟阶段中继续按照此模式运行直到能够满足研究要求。此方法称为投影向前欧拉方法，它是粗粒化投影积分中的一种，能够很清晰地表达粗粒化投影积分的基本原理(见图3)。另外，为了方便理解我们把包含外推过程的模拟阶段时间长ΔT称为投影时长，而在投影时长中需要在细小尺度模型中演化运算的时间长t_k称为演化时长。采用粗粒化投影积分方法对Schlogl和Selkov反应扩散问题中的浓度常稳态进行多尺度模拟，并与格子玻尔兹曼直接模拟进行对比分析。

本发明实施反应扩散浓度常稳态的粗粒化投影积分模拟，在化学反应扩散系统的粗粒化时间步进上耦合宏观尺度的传统外推数值方法，达到对化学反应扩散系统特性快速而准确地捕获，并提高格子玻尔兹曼直接模拟的运算效率。为研究拥有巨大运算量、更复杂的新型多组分反应扩散系统打下坚实的基础。

本发明的研究结果如下：

1.Schlogl反应扩散问题的研究结果

在Schlogl反应扩散的免方程模拟应用中，本发明进行了两项数值实验，第一项实验是用粗粒化投影积分模拟捕获X单点浓度(指系统中一个格点的X浓度，取(2，8)点)的常稳态并与直接模拟结果做对比。基本参数：系统大小取100*100，演化时步δt取0.001s，演化时长取20步，浓度初始分布的系综平均为0.2，浓度允许误差为±5×10-4，如图4所示。

从图4可以看出，除了投影时长为0.4s的粗粒化投影积分模拟达到稳态需要的模拟时间稍长外，其他粗粒化投影积分模拟都比直接模拟经历的模拟时间要短，并且都能准确地获取X单点浓度的常稳态。这说明粗粒化投影积分模拟不仅能准确地得到宏观系统特性，而且当投影时长选取适当时还能比直接模拟表现得更快捷。

我们再把系统扩大到128*128，浓度初始分布的系综平均增大到0.8，以X系统浓度(指系统中所有X浓度的均值)为研究对象进行第二项实验，如图5所示。

由图5可知，直接模拟需要经过近4.5s才能进入稳态区域(0.7495～0.7505)，而当投影时长为0.1s、0.2s时粗粒化投影积分模拟分别只需近2s、3.5s就可实现。再次说明了粗粒化投影积分对宏观系统特性具有快且准地捕获能力。接下来我们还对这两项实验中各次模拟的运算效率分别进行了对比，见表1(本研究的模拟结果均在处理器为AMD Athlon(tm)，主频为2.01GHz，操作环境是Windows XP的计算机上应用MATLAB编程运行得到)。

由表1可知，无论在第一项还是第二项中采用粗粒化投影积分模拟所耗费的运算时间(指在计算机上投入运行的时间)都远小于直接模拟，而且随着投影时长的逐渐加大，所需运算时间也随之缩短。说明粗粒化投影积分这种免方程方法能提高运算效率。因此对于那些传统直接模拟中运算量过大、收敛困难等问题，它是一种极为有效的处理方法。

最后我们再选取投影时长为0.2s的粗粒化投影积分模拟与直接模拟做下对比。内容为：当它们首次达到常稳态时各自耗费的运算时间，在第一项中分别约为7s和88s，在第二项中分别约为17s和221s。很明显粗粒化投影积分模拟节省了近12倍的运算资源。需指出，本研究中的运算时间不同于模拟时间。

2.在Selkov反应扩散中的研究结果

在此应用研究中做了四项数值实验(因对于M和N实验效果类似，故在这只展示有关M的实验)，第一项实验是以M单点浓度(取(5)点)为对象对比粗粒化投影积分模拟和直接模拟分别获取其常稳态的效果。基本参数为：系统为一维，格子大小取1000，演化时步δt取0.001s，演化时长取20步，浓度初始分布的系综平均分别为1.34、0.34，浓度允许误差为±5×10-6，如图6所示。

由图6可知，当投影时长为0.1s、0.2s时粗粒化投影积分模拟能够很快地达到常稳态，明显要优于直接模拟，表明粗粒化投影积分模拟对于具有多组分的反应扩散同样适用，显示了其一般性。再则由表2能看出直接模拟所耗费的运算时间将近是投影时长0.3s的粗粒化投影积分模拟的16倍，但是它们的实验效果几乎一致。这样再结合Schlogl反应扩散研究中的结果分析，从而能更加肯定粗粒化投影积分在优化直接模拟运算效率方面的优势。另外，在实际应用中需注意对投影时长的选取。当投影时长过大时将影响模拟效果，而过小则运算效率又得不到明显地改善。

本研究为了进一步分析粗粒化投影积分技术特点以掌握其应用，本文还考察了演化时长对CPI模拟的影响。在第二项实验中保持外推时长(即投影时长与演化时长的差值)不变，参数为：(1)投影时长取0.09s时，演化时长取10步；(2)投影时长取0.1s时，演化时长取20步；(3)投影时长取0.12s时，演化时长取40步；(4)投影时长取0.14s时；演化时长取60步。对象为M系统浓度，其它参数与第一项相同，如图7所示。

从图7可以看出，当演化时长由10步增大到20步时粗粒化投影积分获取稳态所需的模拟时间明显减少，但是当演化时长继续加大时效果却不再明显，几乎同时达到稳态。另外，由表2可得知演化时长逐渐增大时运算时间随之增加。由此可推断出演化时长能直接影响粗粒化投影积分模拟的效果，演化时长越大实验效果越好但耗费的运算时间也随之越长。其中，还值得注意的是演化时长具有临界值，超过这个临界值将削弱演化时长对粗粒化投影积分模拟的影响。因此在粗粒化投影积分应用研究中，演化时长需尽量接近其临界值。

上述两个实施例表明，本发明成功地将免方程方法中的粗粒化投影积分应用于反应扩散的系统特性研究中，首次对Schlogl、Selkov反应扩散问题中的浓度常稳态进行免方程多尺度模拟。本发明验证了该框架的可行性，并直接在细小尺度模型上进行了宏观水平的研究工作，从而实现了免方程的基本属性——无需获得宏观模型方程。更重要的是，相比传统的直接模拟，它能够更快速而准确地获得Schlogl、Selkov反应扩散的浓度常稳态并且能优化格子玻尔兹曼直接模拟的运算效率。本发明进一步完善免方程方法的基本理论并为它在国内各研究领域中的探索应用提供参考。

Claims (5)

1.化学反应扩散中的免方程多尺度模拟方法，其特征在于，包括如下步骤：

(一)以免方程为基本模拟框架，针对化学反应扩散系统的特性，将系统空间尺度划分为宏观尺度和介观尺度，以此为基础建立免方程多尺度模拟方法，实现通过细小尺度的模拟得出系统宏观性能；

(二)根据Schlogl、Selkov反应扩散系统在宏观尺度和介观尺度上的时空特性，筛选出宏观系统反应物浓度和介观粒子密度分布函数，并从中分析得出具体的提升表达式和约束表达式；从介观尺度出发建立能真实逼近Schlogl和Selkov反应扩散系统的细小尺度模型，即格子玻尔兹曼模型；

(三)分别对Schlogl和Selkov反应扩散系统构造粗粒化时间步进；该粗粒化时间步进是由Schlogl和Selkov反应扩散系统的格子玻尔兹曼模型中的演化、提升和约束构成；

(四)采用传统数值技术中的欧拉法与Schlogl和Selkov反应扩散系统的粗粒化时间步进耦合成Schlogl和Selkov反应扩散系统的粗粒化投影积分，从而实现不同投影时长和演化时长的浓度常稳态多尺度模拟，并与格子玻尔兹曼直接模拟进行对比；其中，采用最小二乘法对宏观系统浓度值进行拟合。

2.根据权利要求1所述的化学反应扩散中的免方程多尺度模拟方法，其特征在于，在步骤(一)中，免方程的基本框架由能跨越时间尺度差的粗粒化时间步进和与之相应的数值技术构成；其中，粗粒化时间步进相当于整个框架中的“内核”作用，类似于系统辨识中的“黑箱”，而数值技术的作用则可看似粗粒化时间步进的“外壳”；此处粗粒化时间步进是一个由提升、演化和约束组成的跨尺度组合；基本框架的一般实现过程：

首先需要做两个准备：

a)根据具体的研究任务对系统进行尺度划分，即宏观尺度和细小尺度，在相应的细小尺度上建立模型并根据不同尺度的时空特性筛选出具体的物理量和粒子离散分布分别作为宏观系统量u和细小量离散分布U；

b)根据宏观系统量u和细小量离散分布U之间的映射关系分析得出具体的提升表达式，则提升μ可表达为：U＝μu；再由已得到的提升μ推导出合适的约束表达式，约束M可表达为：u＝MU；其中提升μ与约束M互逆：μM＝I；

完成以上准备工作后就可按如下描述进行：

1)对宏观系统量设定一个基本符合研究要求的初值u(t₀)，t₀为系统的初始时刻；

2)先通过提升μ把初值u(t₀)离散转化为细小量离散分布，具体表达式：U(t₀)＝μu(t₀)；

3)再将细小量离散分布U(t₀)投入细小尺度模型中，运行演化时步δt得到新的细小量离散分布，其过程可表达为：

其中

为演化算子；

4)最后，把新得到的细小量离散分布函数通过约束M得到宏观系统量的新值，其过程

具体表达为：u(t₀+δt)＝MU(t₀+δt)；把刚获得的新值返回第2)步中继续运行；以上的第2)、3)、4)步就是粗粒化时间步进的三个基本步骤；把这三步的表达形式结合起来，推导出粗粒化时间步进的一般表达式为：

3.根据权利要求1所述的化学反应扩散中的免方程多尺度模拟方法，其特征在于，在步骤(二)中，选用反应扩散格子玻尔兹曼模型分别作为反应扩散系统的细小尺度模型，对于有反应源项的扩散方程，是在碰撞项中加入反应项具体形式如下：

$f_{s,i}(x+c_i\mathrm{\Δt},t+\mathrm{\Δt})-f_{s,i}(x,t)=-1/{\mathrm{\τ}}_s[f_{s,i}(x,t)-f_{s,i}^{\mathrm{eq}}(x,t)]+\mathrm{\Δt}{R}_s/b$

上式可作为Schlogl和Selkov反应扩散的细小尺度模型；

对上式采用多尺度方法得到反应扩散宏观方程：

$\frac{{\∂\mathrm{\ρ}}_s}{\∂t}-D_s{\▿}^2{\mathrm{\ρ}}_s=R_s$

其中f_s，i(x，t)为s组分时刻t位于x处具有速度c_i的粒子密度分布函数；

为s组分的局部平衡态分布函数；ρ_s为反应物浓度；D_s为扩散系数；τ_s为松弛时间；b为格子模型中的速度数；R_s为反应扩散方程里的反应项。

4.根据权利要求1所述的化学反应扩散中的免方程多尺度模拟方法，其特征在于，在步骤(三)中，粗粒化时间步进以格子玻尔兹曼模型为仿真器，结合介观尺度信息—宏观尺度的交互耦合机制，实现扩散系统的跨尺度组态。

5.根据权利要求1所述的化学反应扩散中的免方程多尺度模拟方法，其特征在于，在步骤(四)中，粗粒化投影积分是免方程方法中一种具有代表性的基础方法，是用Euler法进行外推处理，其过程可表示为：

$u^{0,N+1}=u^{k,N}+(\mathrm{\ΔT}-t_k)\tilde{F}\left(u^{k,N}\right)$

其中，

通过已获得的宏观系统值得到：

$\tilde{F}\left(u^{k,N}\right)=\frac{u^{k+1,N}-u^{k,N}}{\mathrm{\δt}}\≈\frac{\mathrm{du}(\mathrm{N\ΔT}+t_k)}{\mathrm{dt}}=F\left(u(\mathrm{N\ΔT}+t_k)\right)\≈F\left(u^{k,N}\right)$

采用最小二乘法拟合；最后，再把投影获得的u^0，N+1输入到下一个模拟阶段中继续按照此模式运行直到能够满足研究要求。

Images