一种关于数值解场变量的高精度求解方法与流程

1.本发明属于场量应用分析、数据处理技术领域，具体涉及一种关于数值解场变量的高精度求解方法。


背景技术：

2.自然界中，几乎所有问题都可以转化微分方程的求解问题，如物体的运动问题、热传导问题、人口增长问题等等，所以微分方程的求解对于促使我们认识事物发展变化的规律、提高相关学科认识，具有关键作用。对于一些简单的微分方程，可以通过一些数学技巧解出，如齐次方程、一阶线性微分方程等，但能得到解析解的微分方程毕竟只是少数，对于从实际问题中提取的大量微分方程，往往无法得到方程的解析解。随着计算机技术的发展，数值方法的应用逐渐广泛，成为数学与计算机之间的桥梁，为各种数学问题的求解提供可行手段。
3.经过几十年的发展，为解决从流体流动、固体变形、人口变化等实际问题提取出的数学问题，研究学者们提出的数值分析方法包括：有限单元法、边界元方法、控制容积有限元方法、光滑粒子流体动力学方法、无网格伽辽金方法等，利用这些方法所得到的直接结果是一系列单元结点或离散点上对应的场量值。然而，在部分方法(如有限单元法)的求解过程中或对结果的后处理中，还需要高精度获得非初始单元结点或离散点位置处的场变量及其导数的数值，如当利用有限单元法进行多场耦合分析，而所采用的单元网格精度不一致，需要相互传递场的数据时；当要求后处理绘制等值线时；当需要根据位移场数据绘制应变场时；等等。无疑，过程中，非单元结点或离散点位置处的场变量值及其导数值的计算精度将对最终结果产生关键影响。
4.为获得非初始单元结点或离散点处的场变量值，常规做法是通过形函数插值或移动最小二乘方法，但结果精度往往不高；为获得单元结点或离散点处的导数值，其精度则取决于形函数精度，如有限单元法中经常使用的常应变单元，在计算导数值分布上的效果差强人意。如果能够高精度获得任意点处的数值解场变量，保证从实际问题提取出的数学问题的求解过程中或后处理过程中的结果精度，对于我们正确应对实际问题、认识其相关规律，有重大意义。值得注意的是，这里提到的任意点处的数值解场变量，区别于前述的通过有限单元法、边界元方法等数值方法直接得到的一系列单元结点或离散点上对应的场量值，而是以他们为已知量，采用数值分析的方式求解得到的任意点处的函数场量值及其导数值的统称。


技术实现要素：

5.针对现有问题，本发明提供一种关于数值解场变量的高精度求解方法，思想是联合泰勒级数展开技术与移动最小二乘技术，建立并求解关于已知局部离散点坐标及场量值的线性方程组，以高效高精度地得到某目标点处的数值解场变量，并服务于非线性计算的中间步骤，或满足高精度后处理的需求。针对所求解问题的不同来源，即不同实际问题，所
述数值解场变量有不同的实际物理意义，如应力、应变、水位、温度、人口密度等。
6.本发明是通过下述技术方案实现的：
7.一种关于数值解场变量的高精度求解方法，联合泰勒级数展开技术与移动最小二乘技术，获得任一点处场变量关于周围局部离散点场量值的离散表达式，将局部离散点场量值代入离散表达式，计算求解任意目标点处的数值解场变量。
8.一种关于数值解场变量的高精度求解方法，具体包括如下步骤：
9.步骤一，确定目标点的坐标信息，及研究域内所有单元结点或离散点的坐标信息及其场量值、边界信息；
10.步骤二，根据步骤一中目标点的坐标信息和边界信息，结合用户目的，确定局部离散点集合的选取规则、最少点要求和加权函数类型；
11.步骤三，在目标点处建立以目标点为原点的局部坐标系，根据步骤二确定的局部离散点集合的选取规则和最少点要求确定目标点的局部离散点集合；
12.步骤四，根据泰勒展开式、步骤二确定的加权函数类型、步骤三建立的局部坐标系、步骤三确定的目标点的局部离散点集合，建立目标点处的数值解场变量的离散表达式；
13.步骤五，将步骤三局部离散点上的场量值代入步骤四的离散表达式，求解目标点处的数值解场变量；
14.步骤六，结果可视化：将计算结果以数值形式、绘制彩图或绘制等高线图进行显示。
15.所述步骤一中研究域内所有单元结点或离散点的边界信息包括：对研究域边界的描述，及已知的边界条件。
16.所述步骤二中局部离散点集合选取规则的确定，包括如下内容：
17.当研究域内有限单元中的网格信息为已知时，采用邻接单元覆盖方法：在数量上满足最少点要求的情况下，根据所需选取与目标点直接通过单元进行连接的点、目标点所属单元的单元结点即一次点、与一次点通过单元进行连接的点即二次点、与二次点通过单元进行连接的点即三次点或者与三次点通过单元进行连接的点即四次点中的一种或几种作为局部离散点集合；
18.当研究域内网格信息不是已知时，采用相似域方法来确定局部离散点的集合，即当研究域为长条形时，采用长方形或椭圆形作为局部域形状；当研究域为圆形时，采用正方形、圆形或椭圆形为局部域形状；在保证选取局部域内的局部离散点数量满足最少点要求的前提下，选取局部域内所需数量的任意离散点作为局部离散点集合。
19.所述步骤二中局部离散点集合的最少点要求为：用户根据预期精度与待求解数值解场变量类型，在泰勒展开式中确定截断阶数，泰勒展开式中保留项数即为局部离散点集合最少点要求的个数。
20.所述步骤二中确定的加权函数类型包括：不加权、三次样条加权、四次样条加权或负指数函数加权。
21.所述步骤四建立的目标点处数值解场变量的离散表达式如下：
22.g＝d
t
·
s
t
·
w
·
s
·
d
23.h＝d
t
·
s
t
·
w
24.25.其中：d＝diag{1,|δx|
max
,|δy|
max
,|δx|
max2
,|δy|
max2
,|δx|
max
|δy|
max
}， |δx|
max
＝max{|δx1|,|δx2|,...,|δx
n
|}，|δy|
max
＝max{|δy1|,|δy2|,...,|δy
n
|}，t为转置符号，
26.w＝diag([w
1 w2... w
n
])，
[0027]
δx
i
＝x
i
‑
x0，
[0028]
δy
i
＝y
i
‑
y0，(x
i
,y
i
)是第i个局部离散点的坐标，i＝1,2,
…
,n，u0、、为目标点处的数值解场变量，为第i个局部离散点处场量的真值，w
i
为第i个局部离散点的不为零的权函数。
[0029]
本发明的有益效果：
[0030]
1.本发明原理简单、精度高，可适用于多维度、多阶次数学问题求解，在非线性场和高阶场中性能卓越，有利于正确认识实际问题，总结客观规律；
[0031]
2.本发明所述离散表达式具有统一解析表达，可在保证计算精度的同时，缩短计算用时，兼备计算效率高的优势；
[0032]
3.本发明可以任意数值方法计算得到的系列单元结点或离散点上对应的场量值为基础，高效高精度获得任意点处的数值解场变量的值，实现从不同实际问题中提取的数学问题求解过程中的数据传递及后处理目的；
[0033]
4.本发明还可联系任意数值方法，服务于实际问题分析求解过程中非线性计算中间步骤，为数值分析提供新思路，发展前景广阔。
附图说明
[0034]
图1为本发明研究域网格信息已知时局部离散点集合的定义示意图。
[0035]
图2为本发明研究域为长方形时局部离散点集合的定义示意图。
[0036]
图3为本发明研究域为圆形时局部离散点集合的定义示意图。
[0037]
图4为本发明实施例1的单元示意图。
[0038]
图5为本发明实施例2隧道分析粗网格系统示意图。
[0039]
图6为本发明实施例2隧道分析细网格系统示意图。
[0040]
图7为本发明实施例2细网格隧道系统x方向位移分布解析解。
[0041]
图8为本发明实施例2细网格隧道系统y方向位移分布解析解。
[0042]
图9为本发明在实施例2隧道x方向位移后处理中的结果应用。
[0043]
图10为本发明在实施例2隧道y方向位移后处理中的结果应用。
具体实施方式
[0044]
一种关于数值解场变量的高精度求解方法，联合泰勒级数展开技术与移动最小二乘技术，获得任一点处场变量关于周围局部离散点场量值的离散表达式，将局部离散点场量值代入离散表达式，计算得到任意目标点处的数值解场变量，包括以下步骤：
[0045]
步骤一，确定目标点的坐标信息，及研究域内所有单元结点或离散点的坐标信息及其场量值、边界信息；
[0046]
步骤二，根据步骤一中目标点的坐标信息和边界信息，结合用户目的，确定局部离散点集合的选取规则、最少点要求和加权函数类型；
[0047]
步骤三，在目标点处建立以目标点为原点的局部坐标系，根据步骤二确定的局部离散点集合的选取规则和最少点要求确定目标点的局部离散点集合；
[0048]
步骤四，根据泰勒展开式、步骤二确定的加权函数类型、步骤三建立的局部坐标系、步骤三确定的目标点的局部离散点集合，建立目标点处的数值解场变量的离散表达式(见公式 (8))；
[0049]
步骤五，将步骤三局部离散点上的场量值代入步骤四的离散表达式，求解目标点处的数值解场变量；
[0050]
步骤六，结果可视化：将计算结果以数值形式、绘制彩图或绘制等高线图进行显示。
[0051]
所述步骤一中研究域内所有单元结点或离散点的边界信息包括：对研究域边界的描述，及已知的边界条件。
[0052]
所述步骤二中局部离散点集合选取规则的确定，包括如下内容：如附图1
‑
3所示，
[0053]
如图1所示，当研究域内有限单元中的网格信息为已知时，采用邻接单元覆盖方法：在数量上满足最少点要求的情况下，选取与目标点直接通过单元进行连接的点、目标点所属单元的单元结点即一次点、与一次点通过单元进行连接的点即二次点或与二次点通过单元进行连接的点即三次点甚至与三次点通过单元进行连接的点即四次点作为局部离散点集合；注意这里的选择应呈递进关系，不能跨过低次点而只选择高次点作为局部离散点，如不能跨过一次点而只选择二次点或三次点作为离散点，且不是每种类型点都要选取的；
[0054]
如图2和图3所示，当研究域内网格信息不是已知时，采用相似域方法来确定局部离散点的集合，如当研究域为长条形时，采用长方形或椭圆形作为局部域形状；当研究域为圆形时，采用正方形、圆形或椭圆形为局部域形状；局部域的形状参数由用户指定，在保证选取局部域内的局部离散点数量满足最少点要求的前提下，选取局部域内所需数量的任意离散点作为局部离散点集合。
[0055]
所述步骤二中局部离散点集合的最少点要求为：用户根据预期精度(具有一阶导数精度、具有二阶导数精度等)与待求解数值解场变量类型(函数场变量值、一阶导数值、二阶导数值等)，在泰勒展开式中确定截断阶数，泰勒展开式中保留项数即为局部离散点集合最少点要求的个数。
[0056]
所述步骤二中确定的加权函数类型包括：不加权、三次样条加权、四次样条加权或负指数函数加权。
[0057]
所述步骤四的具体内容如下：
[0058]
设目标点0周围有n个离散点，构成目标0点的局部离散点集合，有泰勒完全展开
式：
[0059][0060]
其中，δx
i
＝x
i
‑
x0，δy
i
＝y
i
‑
y0，(x
i
,y
i
)是第i个(i＝1,2,
…
,n)局部离散点的坐标，(δx
i
,δy
i
) 构成第i个局部离散点在以目标点0(x0,y0)为原点的局部坐标系中的局部坐标，u0、、为目标点0处的数值解场变量，分别简写为u0、u
,x
|0、u
,y
|0、u
,xx
|0、 u
,yy
|0、u
,xy
|0；
[0061]
式(1)中未知量为目标点0的数值解场变量(函数场变量u0及其各阶导数)。
[0062]
假设式(1)保留了前k项，为第i个局部离散点处场量的真值，而为通过截断k项后的式(1)计算得到的第i个局部离散点处场量的近似值，为使总体偏差最小，采用最小二乘法求解，并引入加权函数来描述不同局部离散点的重要程度，则要求：
[0063][0064]
其中，w
i
为第i个局部离散点的不为零的权函数，通过用户确定的加权函数类型计算得到，进一步地，式(2)可写成式(3)的矩阵形式；
[0065][0066]
其中，
[0067][0068]
w＝diag([w
1 w2... w
n
])
ꢀꢀꢀ
(5)
[0069][0070][0071]
其中t为转置符号，令g＝s
t
·
w
·
s，h＝s
t
·
w，代入式(3)，有
[0072][0073]
为避免矩阵g的病态与奇异，引入比例矩阵
[0074]
d＝diag{1,|δx|
max
,|δy|
max
,|δx|
max2
,|δy|
max2
,|δx|
max
|δy|
max
}
ꢀꢀꢀ
(9)
[0075]
其中，|δx|
max
＝max{|δx1|,|δx2|,...,|δx
n
|}，|δy|
max
＝max{|δy1|,|δy2|,...,|δy
n
|}，比例矩阵修正后，有如下表达式：
[0076][0077]
从而可以确定目标点0处的数值解场变量(目标点0处的场变量u0及其各阶导数)。
[0078]
本发明的特点是：
[0079]
1.本发明联合泰勒级数展开技术与移动最小二乘技术，可高效高精度获得目标点的数值解场变量的值，满足数据传递及高精度后处理方面的需求，基础功能明确；
[0080]
2.本发明为线性方程组求解运算，计算原理简单易懂，所述离散表达式可以具有统一解析表达，针对不同目标点，通过代入局部坐标系与具体场量值，获得该点处的唯一表达式，具有方便快捷简单的特点，计算效率高；
[0081]
3.本发明所述求解任意点处的数值解场变量的方法与所应对的实际问题无关，在步骤1
‑
4 中均未涉及到场量具体数值及实际物理意义，即针对具有相同网格分布或离散点分布的数学问题，在某点上的离散表达式具有一致性，通过改变代入离散表达式中的已知场量的数值，来实现在不同实际问题中的应用。
[0082]
为使本发明要解决的技术问题、技术方案和优点更加清楚，下面结合附图对本发明做进一步描述。
[0083]
实施例1：三结点三角形单元，有限单元法，温度场
[0084][0085]
根据附图4，目标点坐标为(7,5)。选择邻接单元覆盖算法、四次样条加权曲线(即式(11)，δd
max
＝α
·
max{δd1,δd2,...,δd8}，本实施例中取α＝1.2，为加权函数、二阶精度保证。根据覆盖算法，确定目标点的局部离散点集合为(4,5)、(8,
7.2)、(9,4)、(1.2,5)、(1,1)、(5,3)、(6,1)、(10,1.5)、(12,4)、(11.5,7.25)、(9,10.5)、(4.75,9)，共12个。将局部离散点的场量值u
i
(i＝1,2,
…
,12)在目标点处展开，根据二阶精度，在三阶导数处截断，故有泰勒展开式
[0086][0087]
即k＝6，将局部离散点和目标点的坐标信息代入式(4)、式(11)、式(5)、式(7)、式(9)中，分别得到：
[0088][0089]
w＝diag{[0.9863 0.9941 0.9957 0.8425 0.6777 0.9891 0.9541 0.9306 0.9000 0.9047 0.8376 0.9317]}
[0090]
d＝diag{[1 0.1667 0.1818 0.0278 0.0331 0.0303]}
[0091]
代入式(10)中，可得到g、h的数学表达式。假设研究域内温度场分布形式为u＝x2 y(温度单位：摄氏度)，即
[0092]
u＝[21 71.2 85 6.44 2 28 37 101.5 148 139.5 91.5 31.5625]
t
[0093]
通过式(8)，可计算得到目标点处的函数场变量为54.0000，即：目标点处的温度为 54.0000℃，与解析解54℃的绝对误差为1.4211e
‑
14，相对误差为2.6316e
‑
16，远小于千分之一，准确率很高。
[0094]
实施例2：应变场的后处理应用
[0095]
本实施例中包括含整个研究区域的粗网格和针对开挖扰动区划分的细网格两套网格系统，如附图5、附图6所示，其中附图5中虚线部分为细化网格边界，即附图6的边界。步骤如下：第一步，根据粗网格模型和地应力场进行有限元分析得到整个区域的应力场和位移场；第二步，根据得到的位移场，采用本文所提方案求得细化网格系统中单元结点的位移场；第三步，根据细化位移场求出地应力分布，便于以后的开挖模拟。因此，在这里有限元计算是在粗网格上进行的，再划分细网格的目的是利用细化分布的单元结点来作为目标点，求解目标点上的位移场、应变场、应力场等，以提高应力集中区域的结果精度，正确认识位
移、应力等在近开挖区域部分的分布特征。
[0096]
为说明本发明的有效性，假设附图5左右边界均匀受压，而上下边界为自由边界，选用此边界条件，主要是因为该问题的近似解形式被广泛研究。附图7
‑
10分别展示了实施例 2的x方向位移、y方向位移的解析解和利用本发明计算得到的结果，利用本发明进行计算时，使用了相似域方法构建局部离散点集合，选择了四次样条曲线作为加权函数。附图7
‑
10 表明：利用本发明，针对粗网格下的位移结果进行单元内部加密后处理，可高精度得到开挖扰动区内的位移分布情况，因此本发明具有很好的应用价值。
[0097]
值得指出的是，上述虽然结合附图对本发明的具体实施方式进行了描述，但并非对本发明保护范围的限制。利用本发明所述的技术方案，或本领域的技术人员在本发明技术方案的启发下，设计出类似的技术方案，而达到上述技术效果的，均是落入本发明的保护范围内。