[发明专利]一种基于最小二乘法的Lagrange反解分析方法

摘要

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

主权项

1.一种基于最小二乘法的Lagrange反解分析方法，其特征在于包括以下步骤：1)构造路径线：将沿等时线的积分转换成沿径线和沿迹线的积分；路径线是人为构造的，在各个Lagrange位置处记录到波形的特征点以及特征点之间，将其按等时间划分并按照对应点相互连接起来的曲线；设某一力学量其中，h为Lagrange位置，t为时间，假设实验测量了n个Lagrange位置处的信号，则即有n条迹线数据；对应第i条迹线来说，其开始时间为t_oi，结束时间为t_i，因此有效时间长度为(t_i-t_oi)，i＝1，2，…，n；将迹线按等时间分割成N个点后，其每条迹线间隔为：Δti=ti-toiN-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)构造目标函数：利用最小二乘法将流场中的粒子速度、压力结合动量守恒方程得到一个整体函数；在忽略热传导、体积力、内部能源和能穴的假定条件下，一维平面波拉格朗日坐标下的守恒方程为：动量守恒：ρ0∂u∂t+∂σ∂h=0---(1)]]>质量守恒：∂ϵ∂t+∂u∂h=0---(2)]]>能量守恒：∂E∂t+σρ0∂u∂h=0---(3)]]>式中ε、u、σ分别表示应变、粒子速度和应力；E表示比内能，h和t分别表示拉格朗日坐标和时间，ρ₀表示试件的初始密度，其中应力应变以压为正，相速度以右行波为正；根据路径线法，将沿等时线的计算转化为沿径线和沿迹线的计算，即(∂σ∂h)t=(∂σ∂h)J-(∂σ∂t)h(∂t∂h)J---(4)]]>式中下标“t”表示沿时间的偏导数，下标“J”表示沿径线的偏导数，下标“h”表示沿迹线的偏导数；将公式(1)按照路径法即可表示为：ρ0∂u∂t+(∂σ∂h)J-(∂σ∂t)h(∂t∂h)J=0---(5)]]>对于传统的方法，在测量了粒子速度历程情况下将公式(5)写为如下差分形式：σj,k=σj,k-1+[-ρ0(∂uj,k-1∂t)h+12(∂σj,k-1∂t)h(∂tj,k-1∂h+∂tj,k∂h)J](hk-hk-1)---(6)]]>将其按差分形式展开得uj+1,k-uj,k=-12ρ0[(∂σj,k∂h)J+(∂σj+1,k∂h)J](tj+1,k-tj,k)---(7)]]>+12ρ0(σj+1,k-σj,k)[(∂tj,k∂h)J+(∂tj+1,k∂h)J]]]>按照公式(7)求解得到的应力场满足自洽检验法；假定沿径线的应力剖面σ为关于h的n-1次多项式函数，即应力的n阶导数σj,k=Σi=1nbijhki-1---(8)]]>则其沿经线的偏导数为：(∂σj,k∂h)J=Σi=1nbij(i-1)hki-2---(9)]]>令Δu_j，k＝u_j+1，k-u_j，k，Δt_j，k＝t_j+1，k-t_j，k，Δdthj,k=(∂tj,k∂h)J+(∂tj+1,k∂h)J]]>则公式(7)变为Δuj,k=-12ρ0[Σi=1nbij(i-1)hki-2+Σi=1nbij+1(i-1)hki-2]Δtj,k+12ρ0[Σi=1nbij+1hki-1-Σi=1nbijhki-1]Δdthj,k---(10)]]>令aj,k=-12ρ0Δtj,k,]]>cj,k=12ρ0Δdthj,k]]>则方程(10)可以写为：Δuj,k=[Σi=1nbij(i-1)hki-2aj,k-Σi=1nbijhki-1cj,k]+[Σi=1nbij+1(i-1)hki-2aj,k+Σi=1nbij+1hki-1cj,k]---(11)]]>当实验测量数据为粒子速度u时，根据公式(11)构造出的目标函数为：f=Σj=1LΣk=1M{[Σi=1nbij(i-1)hki-2aj,k-Σi=1nbijhki-1cj,k]---(12)]]>+[Σi=1nbij+1(i-1)hki-2aj,k+Σi=1nbij+1hki-1cj,k]-Δuj,k}2]]>式中L为迹线上的数据点数，M为迹线的条数；3)反解应力和应变：通过目标函数反解出应力和应变关系；根据目标函数(12)，利用最小二乘法求解其系数，当j＝1时其偏导数的矩阵形式为：Σk=1M[-cj,k][(i-1)hki-2a1,k-hki-1c1,k]...Σk=1M[(ii-1)hkii-2a1,k-hkii-1c1,k][(i-1)hki-2a1,k-hki-1c1,k]...Σk=1Mc1,k[(i-1)hki-2a1,k-hki-1c1,k]...Σk=1M[(ii-1)hkii-2a1,k+hkii-1c1,k][(i-1)hki-2a1,k-hki-1c1,k]...Tb1,1...bii,1...b1,2...bii,2...---(13)]]>=Σk=1MΔu1,k[(i-1)hki-2a1,k-hki-1c1,k]]]>当j＝2，3，…，L-1时其偏导数为：当j＝L时其偏导数满足的方程为：Σk=1M[-cL-1,k][(i-1)hki-2aL-1,k+hki-1cL-1,k]...Σk=1M[(ii-1)hkii-2aL-1,k-hkii-1cL-1,k][(i-1)hki-2aL-1,k+hki-1cL-1,k]...Σk=1M[cL-1,k][(i-1)hki-2aL-1,k+hki-1cL-1,k]...Σk=1M[(ii-1)hkii-2aL-1,k+hkii-1cL-1,k][(i-1)hki-2aL-1,k+hki-1cL-1,k]...Tb1,L-1...bii,L-1...b1,L...bii,L...]]>=Σk=1MΔuL-1,k[(i-1)hki-2aL-1,k+hki-1cL-1,k]---(15)]]>对于各个拉格朗日位置，由于在第一条径线上，其初始时刻的粒子速度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)，将其按差分形式展开得：ϵj+1,k-ϵj,k=-12ρ0[(∂uj,k∂h)J+(∂uj+1,k∂h)J](tj+1,k-tj,k)---(16)]]>+12ρ0(uj+1,k-uj,k)[(∂tj,k∂h)J+(∂tj+1,k∂h)J]]]>由于上式中t，h，u_j，k已知，根据计算得到偏导数则公式(16)可以直接求解出应变ε_j，k；同理，将公式(3)按照差分形式展开即可计算得到比内能E_j，k。