CN102663257A - 一种基于最小二乘法的Lagrange反解分析方法 - Google Patents

Abstract

本发明提供一种基于最小二乘法的Lagrange反解分析方法，属于冲击加载下计算材料本构关系技术领域。包括以下步骤：1)构造路径线、即将沿等时线的积分转换成沿径线和沿迹线的积分；2)构造目标函数、即利用最小二乘法将流场中的粒子速度、压力结合动量守恒方程得到的一个整体函数；3)反解应力和应变、即通过目标函数反解出应力和应变关系。该方法能够对已知粒子速度条件下，无需已知第一条应力边界即可以求解出应力和应变时程关系，减少了实验测量的数据量，且测量粒子速度较应力更方便、精度更高。

Description

一种基于最小二乘法的Lagrange反解分析方法

技术领域

本发明涉及一种冲击加载下计算材料本构关系的计算方法，具体来说，是涉及一种基于最小二乘法的Lagrange反解分析方法。

背景技术

近年来，随着人们安全意识的逐渐提高，人们对固体材料在高应变率下的本构模型开展了越来越多的研究。Lagrange分析方法主要就是在事先不作任何本构假定的情况下，通过测量材料不同位置处的一些力学信息(如应力、应变或粒子速度等)的变化，由冲击动力学的三个守恒方程经过适当的数学推演及数值计算，求得其他未知的力学量。该方法主要是用来研究分析材料对冲击加载的动态力学响应形态。

Lagrange分析方法中需要已知各个拉格朗日位置处的试验数据，而该试验数据可以在试验中通过嵌入在材料内部不同拉格朗日位置处的传感器测量试件的应力时程、应变时程或粒子速度时程获得。因此，通过Lagrange分析方法来研究材料率型本构关系的研究，具体可以有三种求解途径：(1)当测量的物理量为应力时程时，根据Lagrange分析方法来求解可以得到粒子速度，应变等时程关系，该方法不存在数学上的计算障碍。(2)当测量的物理量为粒子速度时程时，根据Lagrange分析方法来求解得到应力、应变时程关系，在求解应力时程关系时还需要已知初始位置处的一条应力时程作为边界条件。然而在同一个位置难以同时测量两个物理量。(3)当测量的物理量为应变时程时，根据Lagrange分析方法来求解粒子速度、应力时程关系，不仅需要已知初始位置处的一条粒子速度时程作为边界条件来求解粒子速度时程关系，还需要已知初始位置处的一条应力时程作为边界条件来求解应力时程关系。因此，对于Lagrange分析的三种方法，第一种方法最易于实现；第二种方法需要已知一个应力边界条件；第三种方法最为困难，需要已知一个应力边界条件、一个粒子速度边界条件。

关于Lagrange分析方法的研究，最早是由Fowles R和Williams R.F.在70年代初期提出的，但是由于难以实现而未能得到实际应用。Cowperthwaite和Williams将这一方法的应用范围作了进一步推广。1973年Grady在Fowles的基础上提出了路径线法，用来计算Lagrange分析中所需要的导数。在此基础上，Seaman等人提出了曲面拟合法，即根据实验测量波形的波形特征进行分区曲面拟合，Seaman假定所拟合的曲面是单调光滑的，且沿径线的三阶导数

后来，学者们提出了反解分析方法，反解法避开了Seaman假定而求得应力流场，但由于其结果依赖于假定的未知应力函数形式而可能产生非唯一解，因而它也曾遭到Seaman的否定。GuptaYM提出自洽检验法，由于速度波形推应力波形时需要作一定的假设，因此该方法通过应力场再回过来反推速度场，从而验证结果的可靠性。1989年C.A.Forest提出了冲量时间积分函数法，并将其应用于数据处理和方差估计中，该方法在测量应力来求解是非常有效的，但是当测量的量完全不涉及应力时，要建立精确的函数形式仍有困难。

基于上述原因，传统的方法均是采用第一种方法，即通过锰铜传感器来测量应力并求解出其他物理量(粒子速度、应变等)。但是，粒子速度计的时间响应比测量应力的锰铜传感器好，有效记录时间也长一些，因此，粒子速度计更适合用于测量冲击波后流场的信息。那么，提出一种以粒子速度为已知条件的计算方法成为关键。

发明内容

针对已知粒子速度作为已知条件时传统Lagrange分析方法的不足，我们提供一种基于最小二乘法的Lagrange反解分析方法，能够在已知粒子速度条件下，无需再增加应力边界条件，降低了实验测量的难度，且计算精度也比传统方法的更高。

本发明通过以下技术手段实现：一种基于最小二乘法的Lagrange反解分析方法，包括以下步骤：

1)构造路径线：将沿等时线的积分转换成沿径线和沿迹线的积分；

路径线是人为构造的，在各个Lagrange位置处记录到波形的特征点以及特征点之间，将其按等时间划分并按照对应点相互连接起来的曲线；设某一力学量

其中，h为Lagrange位置，t为时间，其在空间的关系如图1所示；假设实验测量了n个Lagrange位置处的信号，则即有n条迹线数据；对应第i条迹线来说，其开始时间为t_oi，结束时间为t_i，因此有效时间长度为(t_i-t_oi)，i＝1，2，…，n；将迹线按等时间分割成N个点后，其每条迹线间隔为：

$\mathrm{\Δ}{t}_i=\frac{t_i-t_{\mathrm{oi}}}{N-1}$ (i＝1，2，…，n)

在(h，t)平面上同时得到n×N个离散点：(h_i，t_oi+(j-1)Δt_i)，简记为(i，j)；对于每一个固定的j，可得到n个上离散点，因此用最小二乘法拟合一条曲线，这条曲线就叫做第j条(h，t)径线，见图2所示；因为在所研究的流场区域内，前导冲击波波速变化不太大，只要适当选取各条迹线的有效时间长度，用二次多项式来构造(h，t)径线已能满足要求，即：

t_j＝b_1jh²+b₂jh+b₃j  (i＝1，2，…，n)

对某一确定的路径线j(h，t)，当h一定时，j和t有确定的对应关系，故可以用j代替t，而把变量表示为

这样即可以求解出偏导数

即：

2)构造目标函数：利用最小二乘法将流场中的粒子速度、压力结合动量守恒方程得到一个整体函数；

有关拉氏分析的基本原理，Fowles于1973年已经做过较全面的阐述，在忽略热传导、体积力、内部能源和能穴的假定条件下，一维平面波拉格朗日坐标下的守恒方程为：

动量守恒： ${\mathrm{\ρ}}_0\frac{\∂u}{\∂t}+\frac{\∂\mathrm{\σ}}{\∂h}=0---\left(1\right)$

质量守恒： $\frac{\∂\mathrm{\ϵ}}{\∂t}+\frac{\∂u}{\∂h}=0---\left(2\right)$

能量守恒： $\frac{\∂E}{\∂t}+\frac{\mathrm{\σ}}{{\mathrm{\ρ}}_0}\frac{\∂u}{\∂h}=0---\left(3\right)$

式中ε、u、σ分别表示应变、粒子速度和应力；E表示比内能，h和t分别表示拉格朗日坐标和时间，ρ₀表示试件的初始密度，其中应力应变以压为正，相速度以右行波为正；

根据路径线法，将沿等时线的计算转化为沿径线和沿迹线的计算，即

${\left(\frac{\∂\mathrm{\σ}}{\∂h}\right)}_t={\left(\frac{\∂\mathrm{\σ}}{\∂h}\right)}_J-{\left(\frac{\∂\mathrm{\σ}}{\∂t}\right)}_h{\left(\frac{\∂t}{\∂h}\right)}_J---\left(4\right)$

式中下标“t”表示沿时间的偏导数，下标“J”表示沿径线的偏导数，下标“h”表示沿迹线的偏导数；

将公式(1)按照路径法即可表示为：

${\mathrm{\ρ}}_0\frac{\∂u}{\∂t}+{\left(\frac{\∂\mathrm{\σ}}{\∂h}\right)}_J-{\left(\frac{\∂\mathrm{\σ}}{\∂t}\right)}_h{\left(\frac{\∂t}{\∂h}\right)}_J=0---\left(5\right)$

对于传统的方法，在测量了粒子速度历程情况下将公式(5)写为如下差分形式：

${\mathrm{\σ}}_{j,k}={\mathrm{\σ}}_{j,k-1}+[-{\mathrm{\ρ}}_0{\left(\frac{{\∂u}_{j,k-1}}{\∂t}\right)}_h+\frac{1}{2}{\left(\frac{\∂{\mathrm{\σ}}_{j,k-1}}{\∂t}\right)}_h{(\frac{\∂t_{j,k-1}}{\∂h}+\frac{{\∂t}_{j,k}}{\∂h})}_J](h_k-h_{k-1})---\left(6\right)$

在已知粒子速度情况下，按照公式(6)求解应力时还需要一条应力边界条件；GuptaYM提出自洽检验法，由于速度波形推应力波形时需要作一定的假设，而根据应力场求解其他力学量是无需假设，因此该方法还需要通过应力场再回过来反推速度场，检验计算结果的正确性；为了能够满足自洽检验法，这里仍旧从应力场来求解粒子速度的方程出发，将其应力场按差分形式展开得

$u_{j+1,k}-u_{j,k}=-\frac{1}{2{\mathrm{\ρ}}_0}[{\left(\frac{\∂{\mathrm{\σ}}_{j,k}}{\∂h}\right)}_J+{\left(\frac{\∂{\mathrm{\σ}}_{j+1,k}}{\∂h}\right)}_J](t_{j+1,k}-t_{j,k})---\left(7\right)$

$+\frac{1}{2{\mathrm{\ρ}}_0}({\mathrm{\σ}}_{j+1,k}-{\mathrm{\σ}}_{j,k})[{\left(\frac{\∂t_{j,k}}{\∂h}\right)}_J+{\left(\frac{\∂t_{j+1,k}}{\∂h}\right)}_J]$

按照公式(7)求解得到的应力场自然满足自洽检验法；

假定沿径线的应力剖面σ为关于h的n-1次多项式函数，即应力的n阶导数

${\mathrm{\σ}}_{j,k}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{b}_{\mathrm{ij}}{h}_k^{i-1}---\left(8\right)$

则其沿经线的偏导数为：

${\left(\frac{\∂{\mathrm{\σ}}_{j,k}}{\∂h}\right)}_J=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{b}_{\mathrm{ij}}(i-1)h_k^{i-2}---\left(9\right)$

令Δu_j，k＝u_j+1，k-u_j，k，Δt_j，k＝t_j+1，k-t_j，k， $\mathrm{\Δ}{\mathrm{dth}}_{j,k}={\left(\frac{\∂t_{j,k}}{\∂h}\right)}_J+{\left(\frac{\∂t_{j+1,k}}{\∂h}\right)}_J$

则公式(7)变为

$\mathrm{\Δ}{u}_{j,k}=-\frac{1}{2{\mathrm{\ρ}}_0}[\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{b}_{\mathrm{ij}}(i-1)h_k^{i-2}+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{b}_{\mathrm{ij}+1}(i-1)h_k^{i-2}]\mathrm{\Δ}{t}_{j,k}+\frac{1}{2{\mathrm{\ρ}}_0}[\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{b}_{\mathrm{ij}+1}{h}_k^{i-1}-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{b}_{\mathrm{ij}}{h}_k^{i-1}]\mathrm{\Δ}{\mathrm{dth}}_{j,k}---\left(10\right)$

令 $a_{j,k}=-\frac{1}{2{\mathrm{\ρ}}_0}\mathrm{\Δ}{t}_{j,k},$ $c_{j,k}=\frac{1}{2{\mathrm{\ρ}}_0}\mathrm{\Δ}{\mathrm{dth}}_{j,k}$ 则方程(10)可以写为：

$\mathrm{\Δ}{u}_{j,k}=[\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{b}_{\mathrm{ij}}(i-1)h_k^{i-2}{a}_{j,k}-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{b}_{\mathrm{ij}}{h}_k^{i-1}{c}_{j,k}]+[\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{b}_{\mathrm{ij}+1}(i-1)h_k^{i-2}{a}_{j,k}+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{b}_{\mathrm{ij}+1}{h}_k^{i-1}{c}_{j,k}]---\left(11\right)$

当实验测量数据为粒子速度u时，根据公式(11)构造出的目标函数为：

$f=\underset{j=1}{\overset{L}{\mathrm{\Σ}}}\underset{k=1}{\overset{M}{\mathrm{\Σ}}}\left\{\right[\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{b}_{\mathrm{ij}}(i-1)h_k^{i-2}{a}_{j,k}-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{b}_{\mathrm{ij}}{h}_k^{i-1}{c}_{j,k}]---\left(12\right)$

$+[\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{b}_{\mathrm{ij}+1}(i-1)h_k^{i-2}{a}_{j,k}+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{b}_{\mathrm{ij}+1}{h}_k^{i-1}{c}_{j,k}]-{\mathrm{\Δu}}_{j,k}{\}}^2$

式中L为迹线上的数据点数，M为迹线的条数；

3)反解应力和应变：通过目标函数反解出应力和应变关系；

根据目标函数(12)，利用最小二乘法求解其系数，当j＝1时其偏导数的矩阵形式为：

${\left[\begin{array}{c}\underset{k=1}{\overset{M}{\mathrm{\Σ}}}[-c_{j,k}][(i-1)h_k^{i-2}{a}_{1,k}-h_k^{i-1}{c}_{1,k}]\\ .\\ .\\ .\\ \underset{k=1}{\overset{M}{\mathrm{\Σ}}}[(\mathrm{ii}-1)h_k^{\mathrm{ii}-2}{a}_{1,k}-h_k^{\mathrm{ii}-1}{c}_{1,k}][(i-1)h_k^{i-2}{a}_{1,k}-h_k^{i-1}{c}_{1,k}]\\ .\\ .\\ .\\ \underset{k=1}{\overset{M}{\mathrm{\Σ}}}{c}_{1,k}[(i-1)h_k^{i-2}{a}_{1,k}-h_k^{i-1}{c}_{1,k}]\\ .\\ .\\ .\\ \underset{k=1}{\overset{M}{\mathrm{\Σ}}}[(\mathrm{ii}-1)h_k^{\mathrm{ii}-2}{a}_{1,k}+h_k^{\mathrm{ii}-1}{c}_{1,k}][(i-1)h_k^{i-2}{a}_{1,k}-h_k^{i-1}{c}_{1,k}]\\ .\\ .\\ \end{array}\right]}^T\left[\begin{array}{c}{b}_{\mathrm{1,1}}\\ .\\ .\\ .\\ b_{\mathrm{ii},1}\\ .\\ .\\ .\\ b_{\mathrm{1,2}}\\ .\\ .\\ .\\ b_{\mathrm{ii},2}\\ .\\ .\\ .\end{array}\right]---\left(13\right)$

$=\underset{k=1}{\overset{M}{\mathrm{\Σ}}}\mathrm{\Δ}{u}_{1,k}[(i-1)h_k^{i-2}{a}_{1,k}-h_k^{i-1}{c}_{1,k}]$

当j＝2，3，…，L-1时其偏导数为：

当j＝L时其偏导数满足的方程为：

${\left[\begin{array}{c}\underset{k=1}{\overset{M}{\mathrm{\Σ}}}[-c_{L-1,k}][(i-1)h_k^{i-2}{a}_{L-1,k}+h_k^{i-1}{c}_{L-1,k}]\\ .\\ .\\ .\\ \underset{k=1}{\overset{M}{\mathrm{\Σ}}}[(\mathrm{ii}-1)h_k^{\mathrm{ii}-2}{a}_{L-1,k}-h_k^{\mathrm{ii}-1}{c}_{L-1,k}][(i-1)h_k^{i-2}{a}_{L-1,k}+h_k^{i-1}{c}_{L-1,k}]\\ .\\ .\\ .\\ \underset{k=1}{\overset{M}{\mathrm{\Σ}}}\left[c_{L-1,k}\right][(i-1)h_k^{i-2}{a}_{L-1,k}+h_k^{i-1}{c}_{L-1,k}]\\ .\\ .\\ .\\ \underset{k=1}{\overset{M}{\mathrm{\Σ}}}[(\mathrm{ii}-1)h_k^{\mathrm{ii}-2}{a}_{L-1,k}+h_k^{\mathrm{ii}-1}{c}_{L-1,k}][(i-1)h_k^{i-2}{a}_{L-1,k}+h_k^{i-1}{c}_{L-1,k}]\\ .\\ .\\ .\end{array}\right]}^T\left[\begin{array}{c}{b}_{1,L-1}\\ .\\ .\\ .\\ b_{\mathrm{ii},L-1}\\ .\\ .\\ .\\ b_{1,L}\\ .\\ .\\ .\\ b_{\mathrm{ii},L}\\ .\\ .\\ .\end{array}\right]$

$=\underset{k=1}{\overset{M}{\mathrm{\Σ}}}\mathrm{\Δ}{u}_{L-1,k}[(i-1)h_k^{i-2}{a}_{L-1,k}+h_k^{i-1}{c}_{L-1,k}]---\left(15\right)$

对于各个拉格朗日位置，由于在第一条径线上，其初始时刻的粒子速度u_1，k、应力σ_1，k均等于0(k＝1，2，…，M)，则b_jj，1＝0，jj＝1，2，…，n-1；

根据实验测量了粒子速度u，同时t，h，u_j，k已知，通过公式计算可以得到Δu_j，k，Δt_j，k，

Δdth_j，k，a_j，k，c_j，k；根据公式(13)、(14)和(15)，方程中只包含了b_j，k是未知的，因此通过联立求解即可得到应力的n-1次多项式函数的系数，但其函数的有效次项与迹线数M直接相关，即可以精确实现应力的M次多项式函数；

由于实验中测量到了粒子速度u，根据公式(2)，将其按差分形式展开得：

${\mathrm{\ϵ}}_{j+1,k}-{\mathrm{\ϵ}}_{j,k}=-\frac{1}{2{\mathrm{\ρ}}_0}[{\left(\frac{\∂u_{j,k}}{\∂h}\right)}_J+{\left(\frac{\∂u_{j+1,k}}{\∂h}\right)}_J](t_{j+1,k}-t_{j,k})---\left(16\right)$

$+\frac{1}{2{\mathrm{\ρ}}_0}(u_{j+1,k}-u_{j,k})[{\left(\frac{\∂t_{j,k}}{\∂h}\right)}_J+{\left(\frac{\∂t_{j+1,k}}{\∂h}\right)}_J]$

由于上式中t，h，u_j，k已知，根据计算得到偏导数

则公式(16)可以直接求解出应变ε_j，k；同理，将公式(3)按照差分形式展开即可计算得到比内能E_j，k；

与现有技术相比，本发明具有的有益效果为：

1)拉格朗日分析方法是一种研究材料的动态力学性能，特别是研究材料率型本构关系的常用方法。本发明的分析方法利用了拉格朗日分析方法，基于最小二乘法的Lagrange方法在已知粒子速度求解应力时，避免了同一个拉格朗日位置测量两个物理量的问题。

2)本发明的分析方法结合了反解法和Gupta Y M提出的自洽检验法，算法的理论精度可以达到沿径线的M阶导数

(M为迹线数)，且计算结果也是自洽的。

3)根据气体炮冲击加载实验，将混凝土实验测量的粒子速度作为输入量，计算得到了流场中其他物理量历程(应力、应变等)。对应本次实验只测量了4个Lagrange位置处的粒子速度时程曲线，因此，其计算的应力精度为4阶导数

因此，该方法比Seaman法的3阶导数为零具有更高的计算精度。

4)本发明的分析方法利用动量守恒方程通过上述计算方法得到的应力函数表达式为n-1次多项式，算法理论实现了应力函数表达式的迹线数M阶导数为零而Lagrange分析方法一般不得少于4条试验数据，因此，根据需求设计4条以上的多条实验数据，可使该方法具有更高的计算精度。

附图说明

图1为Lagrange实验波形及路径线的示意图(图中只示出4个Lagrange位置迹线)；

图2为(h，t)平面示意图(图中只示出4个Lagrange位置迹线)；

图3为轻气炮实验装置及测试系统示意图；

图4为粒子速度历程示意图；

图5为粒子应力历程示意图；

图6为粒子应变历程示意图；

图7为混凝土材料在冲击加载、卸载过程的应力-应变曲线；

图中：1.第i条路径线；1.高压气室；2.弹丸；3.飞片；4.靶架；5.靶板试件；6.永久磁铁；7.回收室；

具体实施方式

以下结合附图和实施例对本发明做进一步说明，但并不对本发明造成任何限制。

为了进一步说明该方法的有效性，本实施例对具体的一组实验数据进行分析。

实验数据通过直径为100mm的一级气体炮实验得到。实验装置如图3所示，图中飞片3与靶板试件5均为混凝土材料，飞片3安装在弹丸2的端部，飞片3直径为98.0mm，厚度为6.0mm，靶板试件5直径为98.0mm，厚度为12.0mm。

按照下述步骤进行反解分析：

1)构造路径线：将沿等时线的积分转换成沿径线和沿迹线的积分；

在靶板试件5中，每间隔为3.2mm位置处设置一粒子速度计来测量该位置处的速度波形，试验中一共设置测量了4个位置的粒子速度波形。

将记录到波形的特征点以及特征点之间，将其按等时间划分并按照对应点相互连接起来形成曲线；如图4所示，为通过一组粒子速度计测量到的各个Lagrange位置(分别为0.0mm，3.2mm，6.4mm，9.6mm)处的粒子速度历程，在图中分别用g1，g2，g3，g4表示。

设某一力学量

其中，h为Lagrange位置，t为时间，其在空间的关系如图1所示。实验测量了4个Lagrange位置处的信号，则即有4条迹线数据。对应第i条迹线来说，其开始时间为t_oi，结束时间为t_i，因此有效时间长度为(t_i-t_oi)，i＝1，2，…，n。将迹线按等时间分割成N个点后，其每条迹线间隔为：

$\mathrm{\Δ}{t}_i=\frac{t_i-t_{\mathrm{oi}}}{N-1}$ (i＝1，2，…，n)

在(h，t)平面上同时得到n×N个离散点：(h_i，t_oi+(j-1)Δt_i)，简记为(i，j)。对于每一个固定的j，可得到n个上离散点，因此用最小二乘法拟合一条曲线，这条曲线就叫做第j条(h，t)径线，见图2所示。因为在所研究的流场区域内，前导冲击波波速变化不太大，只要适当选取各条迹线的有效时间长度，用二次多项式来构造(h，t)径线已能满足要求，即：

t_j＝b_1jh²+b₂jh+b₃j  (i＝1，2，…，n)

对某一确定的路径线j(h，t)，当h一定时，j和t有确定的对应关系，故可以用j代替t，而把变量

表示为

这样即可以求解出偏导数

即：

2)构造目标函数：利用最小二乘法将流场中的粒子速度、压力结合动量守恒方程得到一个整体函数；

有关拉氏分析的基本原理，Fowles于1973年已经做过较全面的阐述，在忽略热传导、体积力、内部能源和能穴的假定条件下，一维平面波拉格朗日坐标下的守恒方程为：

动量守恒： ${\mathrm{\ρ}}_0\frac{\∂u}{\∂t}+\frac{\∂\mathrm{\σ}}{\∂h}=0---\left(1\right)$

质量守恒： $\frac{\∂\mathrm{\ϵ}}{\∂t}+\frac{\∂u}{\∂h}=0---\left(2\right)$

能量守恒： $\frac{\∂E}{\∂t}+\frac{\mathrm{\σ}}{{\mathrm{\ρ}}_0}\frac{\∂u}{\∂h}=0---\left(3\right)$

式中ε、u、σ分别表示应变、粒子速度和应力；E表示比内能，h和t分别表示拉格朗日坐标和时间，ρ₀表示试件的初始密度，其中应力应变以压为正，相速度以右行波为正；

根据路径线法，将沿等时线的计算转化为沿径线和沿迹线的计算，即

${\left(\frac{\∂\mathrm{\σ}}{\∂h}\right)}_t={\left(\frac{\∂\mathrm{\σ}}{\∂h}\right)}_J-{\left(\frac{\∂\mathrm{\σ}}{\∂t}\right)}_h{\left(\frac{\∂t}{\∂h}\right)}_J---\left(4\right)$

式中下标“t”表示沿时间的偏导数，下标“J”表示沿径线的偏导数，下标“h”表示沿迹线的偏导数；

将公式(1)按照路径法即可表示为：

${\mathrm{\ρ}}_0\frac{\∂u}{\∂t}+{\left(\frac{\∂\mathrm{\σ}}{\∂h}\right)}_J-{\left(\frac{\∂\mathrm{\σ}}{\∂t}\right)}_h{\left(\frac{\∂t}{\∂h}\right)}_J=0---\left(5\right)$

对于传统的方法，在测量了粒子速度历程情况下将公式(5)写为如下差分形式：

${\mathrm{\σ}}_{j,k}={\mathrm{\σ}}_{j,k-1}+[-{\mathrm{\ρ}}_0{\left(\frac{{\∂u}_{j,k-1}}{\∂t}\right)}_h+\frac{1}{2}{\left(\frac{\∂{\mathrm{\σ}}_{j,k-1}}{\∂t}\right)}_h{(\frac{\∂t_{j,k-1}}{\∂h}+\frac{{\∂t}_{j,k}}{\∂h})}_J](h_k-h_{k-1})---\left(6\right)$

在已知粒子速度情况下，按照公式(6)求解应力时还需要一条应力边界条件；GuptaYM提出自洽检验法，由于速度波形推应力波形时需要作一定的假设，而根据应力场求解其他力学量是无需假设，因此该方法还需要通过应力场再回过来反推速度场，检验计算结果的正确性；为了能够满足自洽检验法，这里仍旧从应力场来求解粒子速度的方程出发，将其应力场按差分形式展开得

$u_{j+1,k}-u_{j,k}=-\frac{1}{2{\mathrm{\ρ}}_0}[{\left(\frac{\∂{\mathrm{\σ}}_{j,k}}{\∂h}\right)}_J+{\left(\frac{\∂{\mathrm{\σ}}_{j+1,k}}{\∂h}\right)}_J](t_{j+1,k}-t_{j,k})---\left(7\right)$

$+\frac{1}{2{\mathrm{\ρ}}_0}({\mathrm{\σ}}_{j+1,k}-{\mathrm{\σ}}_{j,k})[{\left(\frac{\∂t_{j,k}}{\∂h}\right)}_J+{\left(\frac{\∂t_{j+1,k}}{\∂h}\right)}_J]$

按照公式(7)求解得到的应力场自然满足自洽检验法；

假定沿径线的应力剖面σ为关于h的n-1次多项式函数，即应力的n阶导数

${\mathrm{\σ}}_{j,k}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{b}_{\mathrm{ij}}{h}_k^{i-1}---\left(8\right)$

则其沿经线的偏导数为：

${\left(\frac{\∂{\mathrm{\σ}}_{j,k}}{\∂h}\right)}_J=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{b}_{\mathrm{ij}}(i-1)h_k^{i-2}---\left(9\right)$

令Δu_j，k＝u_j+1，k-u_j，k，Δt_j，k＝t_j+1，k-t_j，k， $\mathrm{\Δ}{\mathrm{dth}}_{j,k}={\left(\frac{{\∂t}_{j,k}}{\∂h}\right)}_J+{\left(\frac{{\∂t}_{j+1,k}}{\∂h}\right)}_J$

则公式(7)变为

$\mathrm{\Δ}{u}_{j,k}=-\frac{1}{2{\mathrm{\ρ}}_0}[\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{b}_{\mathrm{ij}}(i-1)h_k^{i-2}+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{b}_{\mathrm{ij}+1}(i-1)h_k^{i-2}]\mathrm{\Δ}{t}_{j,k}+\frac{1}{2{\mathrm{\ρ}}_0}[\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{b}_{\mathrm{ij}+1}{h}_k^{i-1}-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{b}_{\mathrm{ij}}{h}_k^{i-1}]\mathrm{\Δ}{\mathrm{dth}}_{j,k}---\left(10\right)$

令 $a_{j,k}=-\frac{1}{2{\mathrm{\ρ}}_0}\mathrm{\Δ}{t}_{j,k},$ $c_{j,k}=\frac{1}{2{\mathrm{\ρ}}_0}\mathrm{\Δ}{\mathrm{dth}}_{j,k}$ 则方程(10)可以写为：

$\mathrm{\Δ}{u}_{j,k}=[\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{b}_{\mathrm{ij}}(i-1)h_k^{i-2}{a}_{j,k}-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{b}_{\mathrm{ij}}{h}_k^{i-1}{c}_{j,k}]+[\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{b}_{\mathrm{ij}+1}(i-1)h_k^{i-2}{a}_{j,k}+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{b}_{\mathrm{ij}+1}{h}_k^{i-1}{c}_{j,k}]---\left(11\right)$

当实验测量数据为粒子速度u时，根据公式(11)构造出的目标函数为：

$f=\underset{j=1}{\overset{L}{\mathrm{\Σ}}}\underset{k=1}{\overset{M}{\mathrm{\Σ}}}\left\{\right[\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{b}_{\mathrm{ij}}(i-1)h_k^{i-2}{a}_{j,k}-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{b}_{\mathrm{ij}}{h}_k^{i-1}{c}_{j,k}]---\left(12\right)$

$+[\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{b}_{\mathrm{ij}+1}(i-1)h_k^{i-2}{a}_{j,k}+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{b}_{\mathrm{ij}+1}{h}_k^{i-1}{c}_{j,k}]-{\mathrm{\Δu}}_{j,k}{\}}^2$

式中L为迹线上的数据点数，M为迹线的条数；

3)反解应力和应变：通过目标函数反解出应力和应变关系；

根据目标函数(12)，利用最小二乘法求解其系数，当j＝1时其偏导数的矩阵形式为：

${\left[\begin{array}{c}\underset{k=1}{\overset{M}{\mathrm{\Σ}}}[-c_{j,k}][(i-1)h_k^{i-2}{a}_{1,k}-h_k^{i-1}{c}_{1,k}]\\ .\\ .\\ .\\ \underset{k=1}{\overset{M}{\mathrm{\Σ}}}[(\mathrm{ii}-1)h_k^{\mathrm{ii}-2}{a}_{1,k}-h_k^{\mathrm{ii}-1}{c}_{1,k}][(i-1)h_k^{i-2}{a}_{1,k}-h_k^{i-1}{c}_{1,k}]\\ .\\ .\\ .\\ \underset{k=1}{\overset{M}{\mathrm{\Σ}}}{c}_{1,k}[(i-1)h_k^{i-2}{a}_{1,k}-h_k^{i-1}{c}_{1,k}]\\ .\\ .\\ .\\ \underset{k=1}{\overset{M}{\mathrm{\Σ}}}[(\mathrm{ii}-1)h_k^{\mathrm{ii}-2}{a}_{1,k}+h_k^{\mathrm{ii}-1}{c}_{1,k}][(i-1)h_k^{i-2}{a}_{1,k}-h_k^{i-1}{c}_{1,k}]\\ .\\ .\\ \end{array}\right]}^T\left[\begin{array}{c}{b}_{\mathrm{1,1}}\\ .\\ .\\ .\\ b_{\mathrm{ii},1}\\ .\\ .\\ .\\ b_{\mathrm{1,2}}\\ .\\ .\\ .\\ b_{\mathrm{ii},2}\\ .\\ .\\ .\end{array}\right]---\left(13\right)$

$=\underset{k=1}{\overset{M}{\mathrm{\Σ}}}\mathrm{\Δ}{u}_{1,k}[(i-1)h_k^{i-2}{a}_{1,k}-h_k^{i-1}{c}_{1,k}]$

当j＝2，3，…，L-1时其偏导数为：

当j＝L时其偏导数满足的方程为：

${\left[\begin{array}{c}\underset{k=1}{\overset{M}{\mathrm{\Σ}}}[-c_{L-1,k}][(i-1)h_k^{i-2}{a}_{L-1,k}+h_k^{i-1}{c}_{L-1,k}]\\ .\\ .\\ .\\ \underset{k=1}{\overset{M}{\mathrm{\Σ}}}[(\mathrm{ii}-1)h_k^{\mathrm{ii}-2}{a}_{L-1,k}-h_k^{\mathrm{ii}-1}{c}_{L-1,k}][(i-1)h_k^{i-2}{a}_{L-1,k}+h_k^{i-1}{c}_{L-1,k}]\\ .\\ .\\ .\\ \underset{k=1}{\overset{M}{\mathrm{\Σ}}}\left[c_{L-1,k}\right][(i-1)h_k^{i-2}{a}_{L-1,k}+h_k^{i-1}{c}_{L-1,k}]\\ .\\ .\\ .\\ \underset{k=1}{\overset{M}{\mathrm{\Σ}}}[(\mathrm{ii}-1)h_k^{\mathrm{ii}-2}{a}_{L-1,k}+h_k^{\mathrm{ii}-1}{c}_{L-1,k}][(i-1)h_k^{i-2}{a}_{L-1,k}+h_k^{i-1}{c}_{L-1,k}]\\ .\\ .\\ .\end{array}\right]}^T\left[\begin{array}{c}{b}_{1,L-1}\\ .\\ .\\ .\\ b_{\mathrm{ii},L-1}\\ .\\ .\\ .\\ b_{1,L}\\ .\\ .\\ .\\ b_{\mathrm{ii},L}\\ .\\ .\\ .\end{array}\right]$

$=\underset{k=1}{\overset{M}{\mathrm{\Σ}}}\mathrm{\Δ}{u}_{L-1,k}[(i-1)h_k^{i-2}{a}_{L-1,k}+h_k^{i-1}{c}_{L-1,k}]---\left(15\right)$

对于各个拉格朗日位置，由于在第一条径线上，其初始时刻的粒子速度u_1，k应力σ_1，k均等于0(k＝1，2，…，M)，则b_jj，1＝0，jj＝1，2，…，n-1；

根据实验测量了粒子速度u，同时t，h，u_j，k已知，通过公式计算可以得到Δu_j，k，Δt_j，k，Δdth_j，k，a_j，k，c_j，k；根据公式(13)、(14)和(15)，方程中只包含了b_j，k是未知的，因此通过联立求解即可得到应力的n-1次多项式函数的系数，但其函数的有效次项与迹线数M直接相关，即可以精确实现应力的M次多项式函数；

由于实验中测量到了粒子速度u，根据公式(2)，将其按差分形式展开得：

${\mathrm{\ϵ}}_{j+1,k}-{\mathrm{\ϵ}}_{j,k}=-\frac{1}{2{\mathrm{\ρ}}_0}[{\left(\frac{\∂u_{j,k}}{\∂h}\right)}_J+{\left(\frac{\∂u_{j+1,k}}{\∂h}\right)}_J](t_{j+1,k}-t_{j,k})---\left(16\right)$

$+\frac{1}{2{\mathrm{\ρ}}_0}(u_{j+1,k}-u_{j,k})[{\left(\frac{\∂t_{j,k}}{\∂h}\right)}_J+{\left(\frac{\∂t_{j+1,k}}{\∂h}\right)}_J]$

由于上式中t，h，u_j，k已知，根据计算得到偏导数

则公式(16)可以直接求解出应变ε_j，k；同理，将公式(3)按照差分形式展开即可计算得到比内能E_j，k；

将试验得到的粒子速度作为初始输入量，利用上述方法计算分别得到应力历程、应变历程，分别如图5、6所示。

由于求解应力时所采用的方程与已知应力求解粒子速度所采用的的方程是相同的，因此，以应力来求解得到的结果与粒子速度的结果相同，保证了自洽法。实验一共测量了四个拉格朗日位置的粒子速度时程关系，因此，所得到的结果沿径线能够达到4阶导数为零，即

如果需要更高阶精度，还需要增加测量数据。

根据应力历程、应变历程即可以得到混凝土材料在冲击加载、卸载过程的应力-应变曲线，如图7所示，应力-应变曲线具有明显的速率相关效应。从图中可以看到冲击波卸载后混凝土中存在残余应变，曲线呈现滞回型。

以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简单推演或替换，都应当视为属于本发明的保护范围。

Claims (1)

1.一种基于最小二乘法的Lagrange反解分析方法，其特征在于包括以下步骤：

1)构造路径线：将沿等时线的积分转换成沿径线和沿迹线的积分；

路径线是人为构造的，在各个Lagrange位置处记录到波形的特征点以及特征点之间，将其按等时间划分并按照对应点相互连接起来的曲线；设某一力学量

其中，h为Lagrange位置，t为时间，假设实验测量了n个Lagrange位置处的信号，则即有n条迹线数据；对应第i条迹线来说，其开始时间为t_oi，结束时间为t_i，因此有效时间长度为(t_i-t_oi)，i＝1，2，…，n；将迹线按等时间分割成N个点后，其每条迹线间隔为：

$\mathrm{\Δ}{t}_i=\frac{t_i-t_{\mathrm{oi}}}{N-1}$ (i＝1，2，…，n)

在(h，t)平面上同时得到n×N个离散点：(h_i，t_oi+(j-1)Δt_i)，简记为(i，j)；对于每一个固定的j，可得到n个上离散点，因此用最小二乘法拟合一条曲线，这条曲线就叫做第j条(h，t)径线；因为在所研究的流场区域内，前导冲击波波速变化不太大，只要适当选取各条迹线的有效时间长度，用二次多项式来构造(h，t)径线已能满足要求，即：

t_j＝b_1jh²+b₂jh+b₃j  (i＝1，2，…，n)

对某一确定的路径线j(h，t)，当h一定时，j和t有确定的对应关系，故可以用j代替t，而把变量表示为

这样即可以求解出偏导数

即：

2)构造目标函数：利用最小二乘法将流场中的粒子速度、压力结合动量守恒方程得到一个整体函数；

在忽略热传导、体积力、内部能源和能穴的假定条件下，一维平面波拉格朗日坐标下的守恒方程为：

动量守恒： ${\mathrm{\ρ}}_0\frac{\∂u}{\∂t}+\frac{\∂\mathrm{\σ}}{\∂h}=0---\left(1\right)$

质量守恒： $\frac{\∂\mathrm{\ϵ}}{\∂t}+\frac{\∂u}{\∂h}=0---\left(2\right)$

能量守恒： $\frac{\∂E}{\∂t}+\frac{\mathrm{\σ}}{{\mathrm{\ρ}}_0}\frac{\∂u}{\∂h}=0---\left(3\right)$

式中ε、u、σ分别表示应变、粒子速度和应力；E表示比内能，h和t分别表示拉格朗日坐标和时间，ρ₀表示试件的初始密度，其中应力应变以压为正，相速度以右行波为正；

根据路径线法，将沿等时线的计算转化为沿径线和沿迹线的计算，即

${\left(\frac{\∂\mathrm{\σ}}{\∂h}\right)}_t={\left(\frac{\∂\mathrm{\σ}}{\∂h}\right)}_J-{\left(\frac{\∂\mathrm{\σ}}{\∂t}\right)}_h{\left(\frac{\∂t}{\∂h}\right)}_J---\left(4\right)$

式中下标“t”表示沿时间的偏导数，下标“J”表示沿径线的偏导数，下标“h”表示沿迹线的偏导数；

将公式(1)按照路径法即可表示为：

${\mathrm{\ρ}}_0\frac{\∂u}{\∂t}+{\left(\frac{\∂\mathrm{\σ}}{\∂h}\right)}_J-{\left(\frac{\∂\mathrm{\σ}}{\∂t}\right)}_h{\left(\frac{\∂t}{\∂h}\right)}_J=0---\left(5\right)$

对于传统的方法，在测量了粒子速度历程情况下将公式(5)写为如下差分形式：

${\mathrm{\σ}}_{j,k}={\mathrm{\σ}}_{j,k-1}+[-{\mathrm{\ρ}}_0{\left(\frac{{\∂u}_{j,k-1}}{\∂t}\right)}_h+\frac{1}{2}{\left(\frac{\∂{\mathrm{\σ}}_{j,k-1}}{\∂t}\right)}_h{(\frac{\∂t_{j,k-1}}{\∂h}+\frac{{\∂t}_{j,k}}{\∂h})}_J](h_k-h_{k-1})---\left(6\right)$

将其按差分形式展开得

$u_{j+1,k}-u_{j,k}=-\frac{1}{2{\mathrm{\ρ}}_0}[{\left(\frac{\∂{\mathrm{\σ}}_{j,k}}{\∂h}\right)}_J+{\left(\frac{\∂{\mathrm{\σ}}_{j+1,k}}{\∂h}\right)}_J](t_{j+1,k}-t_{j,k})---\left(7\right)$

$+\frac{1}{2{\mathrm{\ρ}}_0}({\mathrm{\σ}}_{j+1,k}-{\mathrm{\σ}}_{j,k})[{\left(\frac{\∂t_{j,k}}{\∂h}\right)}_J+{\left(\frac{\∂t_{j+1,k}}{\∂h}\right)}_J]$

按照公式(7)求解得到的应力场满足自洽检验法；

假定沿径线的应力剖面σ为关于h的n-1次多项式函数，即应力的n阶导数

${\mathrm{\σ}}_{j,k}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{b}_{\mathrm{ij}}{h}_k^{i-1}---\left(8\right)$

则其沿经线的偏导数为：

${\left(\frac{\∂{\mathrm{\σ}}_{j,k}}{\∂h}\right)}_J=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{b}_{\mathrm{ij}}(i-1)h_k^{i-2}---\left(9\right)$

令Δu_j，k＝u_j+1，k-u_j，k，Δt_j，k＝t_j+1，k-t_j，k， $\mathrm{\Δ}{\mathrm{dth}}_{j,k}={\left(\frac{\∂t_{j,k}}{\∂h}\right)}_J+{\left(\frac{\∂t_{j+1,k}}{\∂h}\right)}_J$

则公式(7)变为

$\mathrm{\Δ}{u}_{j,k}=-\frac{1}{2{\mathrm{\ρ}}_0}[\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{b}_{\mathrm{ij}}(i-1)h_k^{i-2}+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{b}_{\mathrm{ij}+1}(i-1)h_k^{i-2}]\mathrm{\Δ}{t}_{j,k}+\frac{1}{2{\mathrm{\ρ}}_0}[\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{b}_{\mathrm{ij}+1}{h}_k^{i-1}-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{b}_{\mathrm{ij}}{h}_k^{i-1}]\mathrm{\Δ}{\mathrm{dth}}_{j,k}---\left(10\right)$

令 $a_{j,k}=-\frac{1}{2{\mathrm{\ρ}}_0}\mathrm{\Δ}{t}_{j,k},$ $c_{j,k}=\frac{1}{2{\mathrm{\ρ}}_0}\mathrm{\Δ}{\mathrm{dth}}_{j,k}$ 则方程(10)可以写为：

$\mathrm{\Δ}{u}_{j,k}=[\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{b}_{\mathrm{ij}}(i-1)h_k^{i-2}{a}_{j,k}-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{b}_{\mathrm{ij}}{h}_k^{i-1}{c}_{j,k}]+[\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{b}_{\mathrm{ij}+1}(i-1)h_k^{i-2}{a}_{j,k}+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{b}_{\mathrm{ij}+1}{h}_k^{i-1}{c}_{j,k}]---\left(11\right)$

当实验测量数据为粒子速度u时，根据公式(11)构造出的目标函数为：

$f=\underset{j=1}{\overset{L}{\mathrm{\Σ}}}\underset{k=1}{\overset{M}{\mathrm{\Σ}}}\left\{\right[\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{b}_{\mathrm{ij}}(i-1)h_k^{i-2}{a}_{j,k}-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{b}_{\mathrm{ij}}{h}_k^{i-1}{c}_{j,k}]---\left(12\right)$

$+[\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{b}_{\mathrm{ij}+1}(i-1)h_k^{i-2}{a}_{j,k}+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{b}_{\mathrm{ij}+1}{h}_k^{i-1}{c}_{j,k}]-{\mathrm{\Δu}}_{j,k}{\}}^2$

式中L为迹线上的数据点数，M为迹线的条数；

3)反解应力和应变：通过目标函数反解出应力和应变关系；

根据目标函数(12)，利用最小二乘法求解其系数，当j＝1时其偏导数的矩阵形式为：

${\left[\begin{array}{c}\underset{k=1}{\overset{M}{\mathrm{\Σ}}}[-c_{j,k}][(i-1)h_k^{i-2}{a}_{1,k}-h_k^{i-1}{c}_{1,k}]\\ .\\ .\\ .\\ \underset{k=1}{\overset{M}{\mathrm{\Σ}}}[(\mathrm{ii}-1)h_k^{\mathrm{ii}-2}{a}_{1,k}-h_k^{\mathrm{ii}-1}{c}_{1,k}][(i-1)h_k^{i-2}{a}_{1,k}-h_k^{i-1}{c}_{1,k}]\\ .\\ .\\ .\\ \underset{k=1}{\overset{M}{\mathrm{\Σ}}}{c}_{1,k}[(i-1)h_k^{i-2}{a}_{1,k}-h_k^{i-1}{c}_{1,k}]\\ .\\ .\\ .\\ \underset{k=1}{\overset{M}{\mathrm{\Σ}}}[(\mathrm{ii}-1)h_k^{\mathrm{ii}-2}{a}_{1,k}+h_k^{\mathrm{ii}-1}{c}_{1,k}][(i-1)h_k^{i-2}{a}_{1,k}-h_k^{i-1}{c}_{1,k}]\\ .\\ .\\ .\end{array}\right]}^T\left[\begin{array}{c}{b}_{\mathrm{1,1}}\\ .\\ .\\ .\\ b_{\mathrm{ii},1}\\ .\\ .\\ .\\ b_{\mathrm{1,2}}\\ .\\ .\\ .\\ b_{\mathrm{ii},2}\\ .\\ .\\ .\end{array}\right]---\left(13\right)$

$=\underset{k=1}{\overset{M}{\mathrm{\Σ}}}\mathrm{\Δ}{u}_{1,k}[(i-1)h_k^{i-2}{a}_{1,k}-h_k^{i-1}{c}_{1,k}]$

当j＝2，3，…，L-1时其偏导数为：

当j＝L时其偏导数满足的方程为：

${\left[\begin{array}{c}\underset{k=1}{\overset{M}{\mathrm{\Σ}}}[-c_{L-1,k}][(i-1)h_k^{i-2}{a}_{L-1,k}+h_k^{i-1}{c}_{L-1,k}]\\ .\\ .\\ .\\ \underset{k=1}{\overset{M}{\mathrm{\Σ}}}[(\mathrm{ii}-1)h_k^{\mathrm{ii}-2}{a}_{L-1,k}-h_k^{\mathrm{ii}-1}{c}_{L-1,k}][(i-1)h_k^{i-2}{a}_{L-1,k}+h_k^{i-1}{c}_{L-1,k}]\\ .\\ .\\ .\\ \underset{k=1}{\overset{M}{\mathrm{\Σ}}}\left[c_{L-1,k}\right][(i-1)h_k^{i-2}{a}_{L-1,k}+h_k^{i-1}{c}_{L-1,k}]\\ .\\ .\\ .\\ \underset{k=1}{\overset{M}{\mathrm{\Σ}}}[(\mathrm{ii}-1)h_k^{\mathrm{ii}-2}{a}_{L-1,k}+h_k^{\mathrm{ii}-1}{c}_{L-1,k}][(i-1)h_k^{i-2}{a}_{L-1,k}+h_k^{i-1}{c}_{L-1,k}]\\ .\\ .\\ .\end{array}\right]}^T\left[\begin{array}{c}{b}_{1,L-1}\\ .\\ .\\ .\\ b_{\mathrm{ii},L-1}\\ .\\ .\\ .\\ b_{1,L}\\ .\\ .\\ .\\ b_{\mathrm{ii},L}\\ .\\ .\\ .\end{array}\right]$

$=\underset{k=1}{\overset{M}{\mathrm{\Σ}}}\mathrm{\Δ}{u}_{L-1,k}[(i-1)h_k^{i-2}{a}_{L-1,k}+h_k^{i-1}{c}_{L-1,k}]---\left(15\right)$

对于各个拉格朗日位置，由于在第一条径线上，其初始时刻的粒子速度u_1，k应力σ_1，k均等于0(k＝1，2，…，M)，则b_jj，1＝0，jj＝1，2，…，n-1；

根据实验测量了粒子速度u，同时t，h，u_j，k已知，通过公式计算可以得到Δu_j，k，Δt_j，k，Δdth_j，k，a_j，k，c_j，k；根据公式(13)、(14)和(15)，方程中只包含了b_j，k是未知的，因此通过联立求解即可得到应力的n-1次多项式函数的系数，但其函数的有效次项与迹线数M直接相关，即可以精确实现应力的M次多项式函数；

由于实验中测量到了粒子速度u，根据公式(2)，将其按差分形式展开得：

${\mathrm{\ϵ}}_{j+1,k}-{\mathrm{\ϵ}}_{j,k}=-\frac{1}{2{\mathrm{\ρ}}_0}[{\left(\frac{\∂u_{j,k}}{\∂h}\right)}_J+{\left(\frac{\∂u_{j+1,k}}{\∂h}\right)}_J](t_{j+1,k}-t_{j,k})---\left(16\right)$

$+\frac{1}{2{\mathrm{\ρ}}_0}(u_{j+1,k}-u_{j,k})[{\left(\frac{\∂t_{j,k}}{\∂h}\right)}_J+{\left(\frac{\∂t_{j+1,k}}{\∂h}\right)}_J]$

由于上式中t，h，u_j，k已知，根据计算得到偏导数

则公式(16)

可以直接求解出应变ε_j，k；同理，将公式(3)按照差分形式展开即可计算得到比内能E_j，k。

Images