模拟复杂流体运动的求解方法

摘要

本公开提供了一种模拟复杂流体运动的求解方法，针对流体运动方程，先计算高精度的空间导数，进而计算出高精度的时间导数，再使用Taylor格式求出下一步计算解，并积分得到指定目标时刻的数值解。本公开的方案是适用于大气科学数值模式的时/空精度匹配的高阶算法，使数值模拟中的总误差得到控制，减小了天气和气候数值模式中由于计算误差导致的不确定性。

权利要求

1.一种模拟复杂流体运动的求解方法，针对流体运动方程 F为关于时间和空间的变量，L为包含F以及F关于空间变量导数的算子，先计算高精度的空间导数，进而计算出高精度的时间导数，再使用Taylor格式求出下一步计算解；包括：步骤A，给定初始条件，旋转合适的空间步长h和时间步长τ，设定积分的终值时刻，最高时间积分精度为M阶；

步骤B，在选定的M阶时间精度下，采用n阶高阶精度空间差分格式进行空间差分计算，n与最高时间积分精度M的值相匹配；

步骤C，设置时间积分精度为k阶，k＝1，2，…M，利用高精度的空间导数，进而计算出高精度的k阶时间导数；

步骤D，增加时间积分阶数k，重复步骤C，直到k＝M循环结束，采用Taylor级数法进行时间积分计算一步时间积分；以及步骤E，重复上述步骤B～D，积分得到指定目标时刻表示流体速度方向的数值解；

所述步骤B中，空间差分计算采用递归微分的方式计算，空间差分方法采用以下公式：其中u是一个有关x，t的函数u＝u(x ,t) ，“ux(1)(xi)”即为在xi格点上的函数u对x的一阶导数所述步骤D中，时间积分计算采用以下公式：

其中k＝1,2,…M，

对于线性方程描述的复杂流体运动 采用以下公式求得时间导数其中 k＝1,2,…M， 为空间导数；

对于非线性方程描述的复杂流体运动 其中A为非线性算子，采用以下公式求得时间导数：其中 k＝1,2,…M， 为空间导数。

2.根据权利要求1所述的求解方法，所述步骤A中，最高时间积分精度M为大于等于3的整数；空间差分阶数n为大于等于6的整数。

3.根据权利要求2所述的求解方法，所述步骤A中，最高时间积分精度M为5，10或20；空间差分阶数n为30，50，100或500。

4.根据权利要求2所述的求解方法，所述步骤A中，给定的所述初始条件为连续光滑的，周期性的初始条件。

5.根据权利要求2所述的求解方法，所述步骤B，当最高时间积分精度M越大时，设置匹配的空间差分精度n阶越高。

6.根据权利要求2所述的求解方法，使用了Multiple precision库，采用1024二进制位精度。

7.根据权利要求2至6中任一项所述的求解方法，所述步骤B中，采用n阶高阶精度致空间差分格式，通过递归微分方式进行空间差分计算。

说明书

模拟复杂流体运动的求解方法

技术领域

[0001] 本公开涉及大气科学领域，尤其涉及一种模拟复杂流体运动的高精度求解方法。

背景技术

[0002] 数值模式和数值模拟是定量研究复杂流体运动如天气和气候变化的主要工具，但是由于数学模型中各变量之间复杂的非线性作用关系，导致模拟结果中存在着不确定性。其中计算误差对数值模拟结果的影响可以从复杂的大气环流模式、耦合模式运行结果看出，也可以从简单的混沌动力系统、准地转模式的数值试验得到验证。因此，如何采取有效的方法来控制计算误差的增长，对天气过程和复杂流体运动的长时间的数值计算和精确的数值模拟至关重要。

[0003] 许多流动运动(从复杂的N-S方程到简单的平流方程)都可以写为算子形式：为了使用计算机求解此类方程，目前已有一些研究引进了高精度算法来对其进行数值模拟。但这些研究在空间差分上使用高精度格式的较多，而时间积分方案多为较低精度的计算格式(习惯上将阶数3阶或以下的算法称为低精度或低阶算法；高精度一般为4-
9阶算法，也称高阶算法；10阶以上称为超高阶或超高精度算法)。这就导致了空间、时间计算精度有时不匹配，会影响模式长时间计算的准确性。只有整体提高算法的精度(即同时提高时间和空间精度)，才可使数值模拟中的总误差得到控制，即需要研究适用于大气科学数值模式的时/空精度匹配的高阶算法，以减小天气和气候数值模式中由于计算误差导致的不确定性。

[0004] 公开内容

[0005] (一)要解决的技术问题

[0006] 本公开提供了一种模拟复杂流体运动的高精度Taylor-Li求解方法，以至少部分解决以上所提出的技术问题。

[0007] (二)技术方案

[0008] 根据本公开的一个方面，提供了一种模拟复杂流体运动的高精度求解方法，针对流体运动方程 F为关于时间和空间的变量，L为包含F以及F关于空间变量导数的算子，先计算高精度的空间导数，进而计算出高精度的时间导数，再使用Taylor格式求出下一步计算解；包括：

[0009] 步骤A，给定初始条件，旋转合适的空间步长h和时间步长τ，设定积分的终值时刻，最高时间积分精度为M阶；

[0010] 步骤B，在选定的M阶时间精度下，采用n阶高阶精度空间差分格式进行空间差分计算，n与最高时间积分精度M的值相匹配；

[0011] 步骤C，设置时间积分精度为k阶，k＝1，2，…M，利用高精度的空间导数，进而计算出高精度的k阶时间导数；

[0012] 步骤D，增加时间积分阶数k，重复步骤C，直到k＝M循环结束，采用Taylor级数法进行时间积分计算一步时间积分；

[0013] 步骤E，重复上述步骤B～D，积分得到指定目标时刻表示流体速度方向的数值解。

[0014] 在本公开的一些实施例中，所述步骤A中，最高时间积分精度M为大于等于3的整数；空间差分阶数n为大于等于6的整数，例如最高时间积分精度M为5，10或20；空间差分阶数n为30，50，100或500。当最高时间积分精度M越大时，设置匹配的空间差分精度n阶越高。

[0015] 在本公开的一些实施例中，给定的所述初始条件为连续光滑的，周期性的初始条件。

[0016] 在本公开的一些实施例中，所述高精度求解方法使用了Multiple precision(MP)库，采用1024二进制位精度。

[0017] 在本公开的一些实施例中，所述步骤B中，空间差分计算采用递归微分的方式计算，空间差分方法采用以下公式：

[0018]

[0019] 其中u是一个有关x，t的函数u＝u(x,t)，“ux(1)(xi)”即为在xi格点上的函数u对x的一阶导数
或者空间差分方法采用紧致空间差分方案。

[0020] 在本公开的一些实施例中，所述步骤D中，时间积分计算采用以下公式：

[0021]

[0022] 其中k＝1,2,…M，

[0023] 在本公开的一些实施例中，所述步骤C中，对于线性方程描述的复杂流体运动采用以下公式求得时间导数

[0024]

[0025] 其中 k＝1,2,…M， 为空间导数。

[0026] 在本公开的一些实施例中，对于非线性方程描述的复杂流体运动 其中A为非线性算子，采用以下公式求得时间导数：

[0027]

[0028] 其中 k＝1,2,…M， 为空间导数。

[0029] (三)有益效果

[0030] 从上述技术方案可以看出，本公开的模拟复杂流体运动的高精度求解方法至少具有以下有益效果其中之一：

[0031] (1)通过灵活地调整计算精度阶数，相比常规方法使用的固定低阶时间积分，可以实现5阶、10阶甚至20阶的时间积分算法，从而相应的空间差分精度也能从6阶拓展到30阶、50阶、100阶甚至500阶，因此该时/空精度匹配的高阶算法使得算法的精度得到整体提高，数值模拟的总误差得到控制；

[0032] (2)通过采用连续光滑、周期性较好的初值，可以获得很好的计算效果，精度阶数越高，误差越小。对于线性和非线性问题的计算复杂流体模拟结果表明，空间差分的精度可以远超6阶；

[0033] (3)使用了递归微分来提高计算速度，从而所实现的快速高精度差分格式比直接计算高阶空间微分项的方案速度有显著提升。这种计算方法不仅内存占用小，计算速度快，而且超过20阶的空间差分精度，对于常规方案已经无法在常规时间内完成的计算，本公开的方法仍然有效，当精度提高到100阶时仍能够在常规时将内完成计算；

[0034] (4)由于求解方法使用了Multiple precision(MP)库，采用1024二进制位精度，相当于200位以上的十进制有效数字，因此足以区分绝对误差小到10-200的数值解。

附图说明

[0035] 图1为本公开第一实施例对流体运动方程的计算流程图。

[0036] 图2为本公开第二实施例对线性方程的计算流程图。

[0037] 图3为本公开第二实施例线性情形下计算误差随空间精度阶数的变化。

[0038] 图4为本公开第三实施例对非线性方程的计算流程图。

[0039] 图5为本公开第三实施例非线性情形下计算误差随空间精度阶数的变化。

具体实施方式

[0040] 为使本公开的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本公开进一步详细说明。

[0041] 本公开某些实施例于后方将参照所附附图做更全面性地描述，其中一些但并非全部的实施例将被示出。实际上，本公开的各种实施例可以许多不同形式实现，而不应被解释为限于此数所阐述的实施例；相对地，提供这些实施例使得本公开满足适用的法律要求。

[0042] 在本公开的第一个示例性实施例中，提供了一种模拟复杂流体运动的高精度求解方法，图1为本公开第一实施例流体方程计算方法流程图。如图1所示：针对流体运动方程先计算高精度的空间导数，进而计算出高精度的时间导数，再使用Taylor格式求出下一步计算解。具体步骤如下：

[0043] 步骤A.给定初始条件，旋转合适的空间步长h和时间步长τ，设定积分的终值时刻，最高时间积分精度为M阶；

[0044] 其中：

[0045] 最高时间积分精度M为大于等于3的整数，例如，M＝5,10或20；

[0046] 给定的初始条件为连续光滑的，周期性好的初始条件；

[0047] 步骤B.在选定的M阶时间精度下，采用n阶高阶精度空间差分格式进行空间差分计算，n与最高时间积分精度M的值相匹配；

[0048] 其中：n为大于等于6的整数，例如n＝30，50，100或500，当最高时间积分精度M越大时，设置空间差分精度n阶越高；

[0049] 空间差分计算采用递归微分的方式计算，空间差分方法采用以下公式：

[0050]

[0051] 其中 (i＝j)，“ux(1)(xi)”即为在xi格点上的函数u对x的一阶导数

[0052] 本公开在计算空间高阶导数uy(m)(yi)时尽量避免使用系数 来直接计算，而是使用 来重复计算m次，此即递归微分的方式。这样做的好处如下：仅使用了两个控制参数m和i，因此比使用4个参数n，m，i，j的循环要高效。当m和n较大时，若在求解程序中计算是较费时间的，特别是当m大于10阶以上时更为明显。而且所需的内存为m×(n+1)×(n+1)×(n+1)，也较大，还会影响CPU的cache使用率，从而降低计算速度。此外， 在实际计算时公式简单，速度快，而且存储一次即可反复使用，内存占用比 小1-2个量级。

[0053] 步骤C.设置时间积分精度为k阶，k＝1，2，…M，利用高精度的空间导数，进而计算出高精度的k阶时间导数；

[0054] 步骤D.增加时间积分阶数k，重复步骤C，直到k＝M循环结束，采用Taylor级数法进行时间积分计算一步时间积分；

[0055] 步骤E.重复上述步骤B～D，积分得到指定目标时刻的数值解。

[0056] 为解决高精度求解方法在计算时还会碰到的舍入误差问题，如果仅采用双精度计算，会发现高精度算法得到的绝对误差会徘徊在10-15～10-16，这恰好是双精度计算时的相对误差极限。所以为了得到超高精度方案的计算效果，采用多精度计算是必要的。本公开的求解方法使用了Multipleprecision(MP)库，采用1024二进制位精度，相当于200位以上的十进制有效数字，足以区分绝对误差小到10-200的数值解。

[0057] 至此，本公开第一实施例线性方程情况下的模拟复杂流体运动的高精度求解方法介绍完毕。

[0058] 在本公开的第二个示例性实施例中，还提供了一种模拟复杂流体运动的高精度求解方法。

[0059] 为了便于理解，现以一个复杂流体运动的典型方程(Burgers方程)作为算例来进行说明。本实施例中，以线性为例进行说明。

[0060] 图2为本公开第二实施例线性方程计算方法流程图。如图2所示：

[0061] 第一步，对线性方程(1)进行空间差分计算，

[0062]

[0063] 公式中u是一个有关x，t的函数u＝u(x,t)。

[0064] 在本方案中，空间差分方法选用任意阶精度的有限差分公式(2)。其优点是可以实现高阶精度(甚至超高阶精度)的空间差分格式，从而方便寻找与高阶时间积分匹配的精度更高的差分精度。

[0065]

[0066] 其中

[0067] 第二步是计算时间积分，本求解方法采用的是Taylor级数法。Taylor级数法的核心是计算导数 因此可以利用线性方程 即 对时间导数 进行变换，得到公式(3)：

[0068]

[0069] 其中 k为时间积分的精度，若我们的计算中最高为M阶时间积分精度，那么k＝1,2,…M。计算方程(3)的关键是计算导数项 这可使用公式(2)来计算，“ux(1)(xi)”即为在xi格点上的函数u对x的一阶导数

[0070] 计算时k＝1,2,…M进行循环，当k＝M时，循环结束。此时，可以利用公式(4)求得下一时刻的数值解。

[0071]

[0072] 其中

[0073] 这就在格点位置为xi完成了时间步长为τ的M阶的Taylor格式的一步时间积分。重复上述过程，即可积分得到指定目标时刻的数值解。

[0074] 为验证结果中误差随时间积分精度的变化可以采用以下步骤：

[0075] a.给定初始条件，旋转合适的空间步长和时间步长，设定积分的终值时刻；

[0076] b.选择时将积分精度为3阶，进行计算。

[0077] c.在选定的时间精度下，改变不同的空间差分精度，由3阶到30阶，进行计算。

[0078] d.分析结果中误差的随空间精度的变化情况。

[0079] e.改变时间积分阶数，重复步骤b～d，得到结果中误差随时间积分精度的变化。

[0080] 图3所示为第二实施例线性情形下计算误差随空间精度阶数的变化，横坐标为空间精度阶数，纵坐标为误差取对数(log10)，图中，“□”、 “●”、“*”连线分别代表时间精度为3、4、5、6阶。

[0081] 图3所示求解平流方程(1)时，使用Gauss波初始条件： 计算时x方向的空间步长为h＝1/200，时间步长为τ＝1/400，计算区域为[0,1]，计算400步，从图
5可见，若M＝3阶，n达到6阶时，计算误差与n大于6阶时的计算误差相差不多，这与Feng07的结果一致；但M＝4阶时，n达到8阶左右计算误差接近平缓变化；M＝5阶时，n达到10阶左右计算误差接近平缓变化；M＝6阶时，n达到12阶左右计算误差接近平缓变化；这些结果说明，试验中总计算误差受时间方向差分精度的影响较大，当时间方向精度足够高时，空间方向的格式精度可以远超过6阶。

[0082] 至此，本公开第二实施例线性方程情况下的模拟复杂流体运动的高精度求解方法介绍完毕。

[0083] 在本公开的第三个示例性实施例中，提供了一种非线性方程情况下的模拟复杂流体运动的高精度求解方法。

[0084] 为了便于理解，现以一个复杂流体运动的典型方程(Burgers方程)作为算例来进行说明。本实施例中，以非线性为例进行说明。

[0085] 图4为本公开第三实施例线性方程计算方法流程图。如图4所示，本实施例中，对于一个描述复杂流体的非线性方程(5)，其空间差分方法也选为任意阶精度的有限差分公式(2)。

[0086]

[0087] 其中A为非线性算子，计算时间积分时，非线性方程与线性方程有所不同，使用的是推导得到的公式(6)。

[0088]

[0089] 其中 对非线性方程，仍然用空间差分的公式(2)来计算方程(6)中的 使k＝1,2,…M循环计算，直到k＝M时，循环结束，即可完成对非线性方程中的一步计算。重复上述过程，即可积分得到指定目标时刻的数值解。

[0090] 为验证结果中误差随时间积分精度的变化可以采用以下步骤：

[0091] a.给定初始条件，旋转合适的空间步长和时间步长，设定积分的终值时刻；

[0092] b.选择时将积分精度为3阶，进行计算。

[0093] c.在选定的时间精度下，改变不同的空间差分精度，由3阶到30阶，进行计算。

[0094] d.分析结果中误差的随空间精度的变化情况。

[0095] e.改变时间积分阶数，重复步骤b～d，得到结果中误差随时间积分精度的变化。

[0096] 图5所示为非线性情形下计算误差随空间精度阶数的变化，横坐标为空间精度阶数，纵坐标为误差取对数(log10)，图中，“□”、 “●”、“*”连线分别代表时间精度为3、4、5、6阶。

[0097] 图5所示求解无粘性Burgers方程(7)时，初始条件为：u(x,0)＝-sin(x)。试验时，计算区域选为[-π,π]，网格数为N＝800，时间步长为τ＝0.001，侧边界条件为周期边界条件。计算时进行检验计算误差的时刻选为：t＝0.8，共计算800步。从图5可以看出，当时间精度为3阶时，计算误差在空间精度达到6阶后就基本保持不变；当时间精度为4、5、6阶时，有效的空间差分精度阶数分别可达11、18、33。这个试验说明了对非线性Burgers方程，只要时间积分精度阶数足够高，空间方向的差分精度阶数可以超过6阶。

[0098] 为了达到简要说明的目的，上述第一实施例中任何可作相同应用的技术特征叙述皆并于此，无需再重复相同叙述。

[0099] 至此，本公开第三实施例非线性方程情况下的模拟复杂流体运动的高精度求解方法介绍完毕。

[0100] 在本公开的第四个示例性实施例中，提供了一种模拟复杂流体运动的高精度求解方法，先计算高精度的空间导数，进而计算出高精度的时间导数，再使用Taylor格式求出下一步计算解。具体步骤如下：

[0101] 步骤A.给定初始条件，旋转合适的空间步长和时间步长，设定积分的终值时刻，最高时间积分精度为M阶；

[0102] 步骤B.在选定的M阶时间精度下，采用n阶高阶精度致空间差分格式，通过递归微分进行空间差分计算，n与最高时间积分精度M的值相匹配；

[0103] 步骤C.设置时间积分精度为k阶，k＝1，2，…M，利用计算的高精度的空间导数，进而计算出高精度的时间导数；

[0104] 步骤D.增加时间积分阶数k，重复步骤C，直到k＝M循环结束，采用Taylor级数法进行时间积分计算一步时间积分；

[0105] 步骤E.重复上述步骤B～D，积分得到指定目标时刻表示流体速度方向的数值解。

[0106] 为了达到简要说明的目的，上述第一实施例中任何可作相同应用的技术特征叙述皆并于此，无需再重复相同叙述。

[0107] 至此，本公开第四实施例的模拟复杂流体运动的高精度求解方法介绍完毕。

[0108] 表1.采用不同算法实现的求解程序(线性情形)，计算时所需要的墙钟时间，其中DP含义为double precision，MP为Multiple precision，试验平台为使用Intel E5-2640 2.6GHz CPU的Linux系统，时间单位为：秒。

[0109]

[0110] 除非有所知名为相反之意，本说明书及所附权利要求中的数值参数是近似值，能够根据通过本公开的内容所得的所需特性改变。具体而言，所有使用于说明书及权利要求中表示组成的含量、反应条件等等的数字，应理解为在所有情况中是受到「约」的用语所修饰。一般情况下，其表达的含义是指包含由特定数量在一些实施例中±10％的变化、在一些实施例中±5％的变化、在一些实施例中±1％的变化、在一些实施例中±0.5％的变化。

[0111] 再者，单词“包含”不排除存在未列在权利要求中的元件或步骤。位于元件之前的单词“一”或“一个”不排除存在多个这样的元件。

[0112] 此外，除非特别描述或必须依序发生的步骤，上述步骤的顺序并无限制于以上所列，且可根据所需设计而变化或重新安排。并且上述实施例可基于设计及可靠度的考虑，彼此混合搭配使用或与其他实施例混合搭配使用，即不同实施例中的技术特征可以自由组合形成更多的实施例。

[0113] 在此提供的算法和显示不与任何特定计算机、虚拟系统或者其它设备固有相关。各种通用系统也可以与基于在此的启示一起使用。根据上面的描述，构造这类系统所要求的结构是显而易见的。此外，本公开也不针对任何特定编程语言。应当明白，可以利用各种编程语言实现在此描述的本公开的内容，并且上面对特定语言所做的描述是为了披露本公开的最佳实施方式。

[0114] 本公开可以借助于包括有若干不同元件的硬件以及借助于适当编程的计算机来实现。本公开的各个部件实施例可以以硬件实现，或者以在一个或者多个处理器上运行的软件模块实现，或者以它们的组合实现。本领域的技术人员应当理解，可以在实践中使用微处理器或者数字信号处理器(DSP)来实现根据本公开实施例的相关设备中的一些或者全部部件的一些或者全部功能。本公开还可以实现为用于执行这里所描述的方法的一部分或者全部的设备或者装置程序(例如，计算机程序和计算机程序产品)。这样的实现本公开的程序可以存储在计算机可读介质上，或者可以具有一个或者多个信号的形式。这样的信号可以从因特网网站上下载得到，或者在载体信号上提供，或者以任何其他形式提供。

[0115] 本领域那些技术人员可以理解，可以对实施例中的设备中的模块进行自适应性地改变并且把它们设置在与该实施例不同的一个或多个设备中。可以把实施例中的模块或单元或组件组合成一个模块或单元或组件，以及此外可以把它们分成多个子模块或子单元或子组件。除了这样的特征和/或过程或者单元中的至少一些是相互排斥之外，可以采用任何组合对本说明书(包括伴随的权利要求、摘要和附图)中公开的所有特征以及如此公开的任何方法或者设备的所有过程或单元进行组合。除非另外明确陈述，本说明书(包括伴随的权利要求、摘要和附图)中公开的每个特征可以由提供相同、等同或相似目的的替代特征来代替。并且，在列举了若干装置的单元权利要求中，这些装置中的若干个可以是通过同一个硬件项来具体体现。

[0116] 类似地，应当理解，为了精简本公开并帮助理解各个公开方面中的一个或多个，在上面对本公开的示例性实施例的描述中，本公开的各个特征有时被一起分组到单个实施例、图、或者对其的描述中。然而，并不应将该公开的方法解释成反映如下意图：即所要求保护的本公开要求比在每个权利要求中所明确记载的特征更多的特征。更确切地说，如下面的权利要求书所反映的那样，公开方面在于少于前面公开的单个实施例的所有特征。因此，遵循具体实施方式的权利要求书由此明确地并入该具体实施方式，其中每个权利要求本身都作为本公开的单独实施例。

[0117] 以上所述的具体实施例，对本公开的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本公开的具体实施例而已，并不用于限制本公开，凡在本公开的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本公开的保护范围之内。