分数阶蠕变本构方法及模型

1.本发明涉及分数阶蠕变本构方法及模型。


背景技术：

2.沥青混合料在长期的循环荷载和复杂的服役环境中，沥青混合料的力学性能较初始状态逐渐发生劣化，为后期的工程结构维护以及结构性能的演化发展预测带来难题。特别像高速铁路无砟轨道结构中ca砂浆层，因其工程的特殊性以及处于多层结构的之间的隐蔽性，更加剧了其性能发展的预测的挑战性。从材料性能发展角度考虑，建立沥青混合料精准的本构模型是解决其材料性能演化以及结构性能预测的途径之一，是构建现代化基础设施精准维护体系的基础性工作。
3.目前研究途径包括：(1)直接对沥青混合料的试验数据进行线性或非线性函数拟合，构建沥青混合料的应力应变映射关系。其忽略了沥青混合料黏弹性变形过程的时间相关性，导致后期的结构性能计算和发展预测出现较大误差；(2)基于不可逆热力学理论和内变量理论的本构模型构建方法，其若设置内变量过多则会导致计算量大幅增加，同时各内变量的演化方程均是基于耗散势函数对内变量的微分形式，而各内变量之间的相互作用机制尚不能清晰，导致耗散势函数建立时加入太多的主观性因素，不利于推广应用；(3)基于不同力学性质元件建立的概念性本构模型，难以表征沥青混合料变形过程中的非线性现象。对此有学者开始将非线性力学元件(如非线性弹簧元件)和理论引入本构模型的构建中，以此理论为基础利用分数阶黏弹性元件(如abel黏壶)改进原有模型或建立新形式的沥青混合料模型已经有丰富的研究成果。在表征沥青混合料蠕变的完整过程中，目前所建立的模型多是基于蠕变过程的三阶段(即衰减蠕变阶段、稳定蠕变阶段以及加速蠕变阶段)展开的，将前两阶段用黏弹性变形表示，第三阶段用黏塑性变形表示。这种分段表示的研究方法虽也能将整个过程拟合的精度非常高，但在物理意义上缺少对沥青混合料蠕变过程中黏弹与黏塑性变形以及两者之间的耦合作用机制的解释，无法做到模型的普适性。


技术实现要素：

4.本发明所要解决的技术问题总的来说是提供一种分数阶蠕变本构方法及模型。本发明提出了一组由两个不同分数阶参数的黏弹性元件并联组成的新型黏弹性体，并参考burgers模型和西原模型的结构形式，考虑黏弹塑变形阶段的黏弹-黏塑性耦合作用机制，构建了沥青混合料黏弹塑性力学行为特征的分数阶导数本构模型。所建本构模型对沥青混合料蠕变过程各阶段力学行为的性质做了重新解释，描述了不同条件下沥青混合料的蠕变全过程，完善了各阶段的物理意义。通过示例试验并与已有模型进行比较，验证了模型的适用性，对其力学行为提供帮助。
5.为解决上述问题，本发明所采取的技术方案是：
6.方法包括以下步骤；
7.s1，设立分数阶黏弹性体；
8.s2，黏弹塑性本构模型；
9.s3操作步骤；
10.s4，本构模型验证。
11.本发明提出的考虑黏弹性-黏塑性耦合作用机制的分数阶蠕变本构模型具有普适性，不仅可适用于沥青类混合料的蠕变本构方程的拟合，而且也适合蠕变变形曲线类似的其他材料的本构建模。本发明提出的模型建立时以不同阶段的应变性质进行区分，各阶段物理意义明确，可为其他材料本构建模提供思路参考。
附图说明
12.图1是本发明的分数阶黏弹性体示意图。
13.图2是本发明的黏弹塑性本构模型图。
14.图3是本发明的黏弹塑性材料蠕变过程图。
15.图4是本发明的ⅱ部分与iii部分分界点确定示意图。
16.图5是本发明的ⅱ部分与iii部分分界点确定方法示意图。
17.图6是本发明的沥青混合料本构模型拟合结果示意图。
18.图7是本发明的沥青混合料本构模型拟合结果示意图。
19.图8是本发明的试验样机参考示意图。
20.图9是本发明的试件结果参考示意图。
具体实施方式
21.如图1-9所示，本实施例的分数阶蠕变本构方法及模型，
22.s1，设立分数阶黏弹性体；
23.基于沥青混合料在黏弹性变形阶段表现出明显的非线性特征，为此，本发明基于分数阶理论，提出了一组由两个不同参数的分数阶黏壶并联组成的新黏弹性体，如图1所示。
24.分数阶黏弹性元件的本构方程为：
25.σ＝ηd
γ
ε/dt
γ
，0≤γ≤1
ꢀꢀ
(1)；
26.式中，σ、ε分别为应力和应变，η为材料黏弹性系数，t为时间，γ为求导阶数；当γ＝0时，σ＝ηε，即为hooke弹簧元件；
27.当γ＝1时，即为牛顿流体元件；
28.通过调整参数γ黏弹性模型较线性组合黏弹性体和仅一个分数阶黏弹性元件与线性元件组合的黏弹性体一般化，普适性更强。
29.将公式(1)涉及的两并联元件关系可得，
30.σ＝σ1+σ2ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
式(2)；
31.ε＝ε1(t)＝ε2(t)
ꢀꢀ
式(3)；
32.σ、ε分别为黏弹性体总应力、总应变；σ1、ε1分别为黏弹性元件1应力、应变；σ2、ε2分别为黏弹性元件2应力、应变；
33.对公式(2)、(3)采用riemann-liouville型分数阶导数算子定义模型中分数阶黏弹性元件，则：
[0034][0035][0036]
η1、η2分别为两个黏弹性元件的材料黏弹性系数；γ1、γ2为求导阶数。
[0037]
联立式(1)-(5)得：
[0038][0039]
对式(6)两边进行laplace变换，
[0040][0041]
式中s为laplace变换域自变量；η1、η2分别为两个元件的材料黏弹性系数材料黏弹性系数，γ1、γ2分别为两个元件的求导阶数。
[0042]
考虑初始状态时，ε1(0)＝0，且当施加应力值σ为常数σ0时，则式(7)变为：
[0043][0044][0045]
为应变ε(t)经laplace变换后的值
[0046]
利用mittag-leffler函数的laplace变换式，求式(9)的逆变换：
[0047][0048]
式中：p为laplace变换域自变量；k为求导阶数；a，b为数学公式中常系数。
[0049]
式中是mittag-leffler函数的k阶导数，mittag-leffler函数的定义为：
[0050][0051]
z为自变量，a,b为常系数
[0052]
式中γ()为gamma函数，定义为：
[0053][0054]
由式(9)、(10)、(11)，推得的laplace逆变换式，并由式(3)得的分数阶粘弹性体的本构方程：
[0055][0056]
取mittag-leffler函数展开式的前两项进行近似表达，则式(13)简化为：
＞0则表示求导运算，γ5＜0则表示积分运算。
[0079]
排除材料的屈服硬化效应，设材料的屈服应力强度为σs，则模型iii部分黏塑性本构方程为：
[0080][0081]
对式(18)进行laplace变换，再进行laplace逆变换可得黏塑性体的本构方程。
[0082][0083]
黏弹-黏塑性耦合效应应变取与(ⅱ)黏弹性本构方程不同参数的分数阶黏弹性体应变表示，本构方程为：
[0084][0085]
式中ε
′
(t)为模型iii部分黏弹-黏塑性耦合效应应变；η3、η4分别为两个黏弹性元件的材料参数；γ3、γ4为导数的阶数。
[0086]
联立式(19)、式(20)，则模型iii部分的应变为：
[0087][0088]
联立式(15)、式(16)、式(21)得到材料蠕变过程的本构方程，分两种情况，即：
[0089][0090]
s3操作步骤，在图3中。
[0091]
s3.1，第一阶段，进行蠕变试验
[0092]
s3.1.1，制作试件。可以采用室内制作或者现场取样的方式获取试验用试件。采用室内制作的方式，需要依据行业标准制作试件(目前主要为圆柱体)，并在一定条件下进行养护；采用现场取样的方式，则直接利用钻芯取样的方法获取试件，但钻芯取样虽然是对现场结构具有一定的破坏性，但其获取的试验数据与实际更为吻合；
[0093]
s3.1.2.测量尺寸。对获取的试件进行尺寸测量，选择是否符合相关规定对试件尺寸误差的规定；
[0094]
s3.1.3.确定加载设备。选择精度能够满足要求的包含温控箱的万能试验机或者其他设备(具备荷载、试件应变及应变发生时间同步测量的功能，且可以控制加载速率)；
[0095]
s3.1.4.安装试件。将试件对中放置于加载试验台上，并控制温控箱内温度(放置时间与试件大小相关)，待试件温度与温箱设置温度一致后记录温度；
[0096]
s3.1.5.进行初始加载及测试。首先，依据拟合曲线使用目的，确定初始加载的大小，其中，依据施加的恒定荷载应力σ0与屈服强度σs的大小关系，可以分为两种。当σ0＜σs时，仅可以得到i、ii部分；当σ0≥σs时，可以得到i、ii、iii部分；
[0097]
开始加载，并进行荷载、应变、时间参数的测量，测量方式可以采用离散测量也可以采用连续测量的方式，完成加载的时间为t0，对应的应变为ε0，初始加载速率应稍高，缩短初始加载时间，减小蠕变应变的发生；
[0098]
s3.1.6.蠕变测试。初始荷载加载完成后，保持荷载稳定不变，持续记录荷载稳定时间和对应的应变，获取持续的时间-应变(t
c-εc)关系曲线；
[0099]
由于蠕变试验持续时间较长，为避免采样数据量过大可适当减小采样频率，同时蠕变试验应采用应力控制下的应变数据采集。
[0100]
s3.2，第二阶段，数据处理
[0101]
s3.2.1，时间修正。将蠕变发生的时间设为零点，即应变(ε)-时间(t)关系曲线向右侧偏移t0，得到修正后的时间-应变(t
x-ε
x
)关系曲线；
[0102]
s3.2.2.数据拟合。包括σ0＜σs和σ0≥σs两种工况。
[0103]
(1)当σ0＜σs时，基于t
x-ε
x
曲线，利用式(22)的中σ0＜σs情况对应的材料本构方程，利用l-m算法进行参数拟合，也可以转化为目标优化的方式进行拟合，将得到对应的参数e、η1、η2、γ1、γ2，分别为弹性模量、材料黏弹性系数；
[0104]
(2)当σ0≥σs时，分三类分别进行。
[0105]
①
屈服应力σs与其对应的屈服应变εs均已知：
[0106]
如图4所示，依据屈服应变εs均在t
x-ε
x
曲线中确定ii部分与iii部分的分界点对应的时间为t1，利用0～t1间的t
x-ε
x
曲线数据并按照式(22)的中σ0＜σs情况进行拟合，确定参数e、η1、η2、γ1、γ2取值。
[0107]
将参数e、η1、η2、γ1、γ2值代入式(22)的中σ0≥σs情况对应的材料本构方程，并利用图4中t1～td间的t
x-ε
x
曲线数据进行拟合，确定参数η3、η4、η5、γ3、γ4、γ5的取值。
[0108]
②
屈服应力σs与其对应的屈服应变εs均未知：
[0109]
将试件在同等温度条件下，按照设定的加载速率开展单轴压缩试验，至试件发生破坏；加载过程中连续记录应力-应变测试数据；并根据应力-应变曲线确定材料的屈服应力σs与其对应的屈服应变εs；
[0110]
此时，按照
“①
屈服应力σs与其对应的屈服应变εs均已知”情况进行材料本构方程中待定常数的确定，从而确定材料的本构方程。
[0111]
③
屈服应力σs已知，屈服应变εs均未知：
[0112]
方法一：是按照
“②
屈服应力σs与其对应的屈服应变εs均未知”进行分析，重新开展单轴压缩试验，确定屈服应力与其对应的屈服应变。
[0113]
以图5为例。方法二：利用已有测试数据确定屈服应变。
[0114]
设整个曲线测点个数为n，第n个测点对应的时刻为tn，则利用拟合结果与测试结果误差最小为目标，确定ii部分与iii部分的分界点，建立的优化模型：
[0115][0116]
在图5所示测试曲线中，ii部分与iii部分的分界点应在变化比较平缓的部分，因此为简化式(23)优化模型的计算方法，可以在时刻0～tn的范围：
[0117]
1)计算应变的前向差分，并确定差分最小值对应的时刻tm及对应的应变εm；
[0118]
2)确定屈服应变可能存在范围
[0119]
在图5中，式中ε
max
为所测试的最大应变；k为经验系数，其值越小所确定的范围越大，优化耗时越长；依据实际拟合经验，建议取值为3-5即可满足要求，从而确定屈服应变所对应时间范围[t
x
，ts]，即为tn的取值范围。
[0120]
式(23)的求解方法可以采用遗传算法或蚁群算法，其得到的优化结果不仅包含了屈服应变对应的时刻，也包含了具体的拟合参数；
[0121]
s4，本构模型验证
[0122]
(1)对沥青混合料进行拟合，拟合结果如图6所示。
[0123]
具体工况为：试验施加应力值为0.537mpa，屈服应力为0.358mpa，最终拟合曲线表达式，
[0124][0125]
屈服前黏弹性阶段与屈服后黏弹塑性两段数据拟合结果的表征拟合效果的参数r
12
，r
12
越接近1越好，分别为0.999、0.997。
[0126]
(2)沥青砂拟合，如图7所示。具体工况为：试验施加应力值为0.2mpa，屈服应力为0.05mpa，最终拟合曲线表达式；
[0127]
ε＝0.0573+0.2113t
0.5264-0.1864t
0.5262
+4535.013t
13.0474-0.0015t
0.1462
+0.1723t
0.1477
(29)
[0128]
屈服前黏弹性阶段与屈服后黏弹塑性两段数据拟合结果的表征拟合效果的参数
越接近1越好，分别为0.989、0.998。
[0129]
本发明充分描述是为了更加清楚的公开，而对于现有技术就不再一一列举。
[0130]
最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；作为本领域技术人员对本发明的多个技术方案进行组合是显而易见的。而这些修改或者替换，并不使相应技术方案的本质脱离本发明实施例技术方案的精神和范围。

技术特征：
1.一种分数阶蠕变本构方法，其特征在于：方法包括以下步骤；s1，设立分数阶黏弹性体；s2，黏弹塑性本构模型；s3操作步骤；s4，本构模型验证。2.根据权利要求1所述的分数阶蠕变本构方法，其特征在于：在步骤s1中，首先，基于沥青混合料在黏弹性变形阶段表现出明显的非线性特征，设立一组由两个不同参数的分数阶黏壶并联组成的新黏弹性体，并建立分数阶黏弹性元件的本构方程为：σ＝ηd
γ
ε/dt
γ
，0≤γ≤1
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
式(1)；式中，σ、ε分别为应力和应变，η为材料黏弹性系数，t为时间，γ为求导阶数；当γ＝0时，σ＝ηε，即为hooke弹簧元件；当γ＝1时，即为牛顿流体元件；然后，利用公式(1)将涉及的两并联元件关系得，σ＝σ1+σ2ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
式(2)；ε＝ε1(t)＝ε2(t)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
式(3)；σ、ε分别为黏弹性体总应力、总应变；σ1、ε1分别为黏弹性元件1应力、应变；σ2、ε2分别为黏弹性元件2应力、应变；其次，对公式(2)、(3)采用riemann-liouville型分数阶导数算子定义模型中分数阶黏弹性元件，则：弹性元件，则：η1、η2分别为两个黏弹性元件的材料参数；γ1、γ2为求导阶数；再次，联立式(1)—(5)得：对式(6)两边进行laplace变换，式中s为laplace变换域自变量；η1、η2分别为两个元件的材料黏弹性系数，γ1、γ2分别为两个元件的求导阶数；再次，在考虑初始状态时，ε1(0)＝0，且当施加应力值σ为常数σ0时，则式(7)变为：时，则式(7)变为：时，则式(7)变为：为应变ε(t)经laplace变换后的值
之后，利用mittag-leffler函数的laplace变换式，即式(10)，求式(9)的逆变换：式中：p为laplace变换域自变量；k为求导阶数；a,b为常系数；在后，式中是mittag-leffler函数的k阶导数，mittag-leffler函数的定义为:z为自变量，a,b为常系数式中γ()为gamma函数，定义为：γ(z)＝∫
0∞
e-t
t
z
dt，re(z)＞0
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
式(12)；由式(9)、(10)、(11)，推得的laplace逆变换式，并由式(3)得的分数阶粘弹性体的本构方程：取mittag-leffler函数展开式的前两项进行近似表达，则式(13)简化为：3.根据权利要求1所述的分数阶蠕变本构方法，其特征在于：在步骤s2中，s2.1，鉴于沥青相基质的存在，设定黏弹塑性材料在设定小应力作用下只发生黏弹性的蠕变变形；在超过材料屈服应力强度的蠕变过程中，总的变形由瞬时弹性进入黏弹性过渡段，直至应变超过相应的应变屈服值时，塑性变形开始出现，即进入黏塑性阶段；s2.2，在进入屈服状态后，将s1的分数阶黏弹性体引入以代表该非线性耦合机制；然后，参考burgers模型和西原模型的结构形式，构建了沥青混合料的本构模型；其次，以材料应变属性将整个过程划分为瞬时弹性应变、黏弹性应变以及黏弹塑性应变；其中，模型i的瞬时弹性应变符合hooke定律，为模型ii的黏弹性应变由新型分数阶黏弹性体表征，为模型iii的黏弹塑性应变由一个新型分数阶黏弹性体和一个由分数阶黏弹性元件并联一个塑性元件的黏塑性体串联进行表征；s2.3，在本构方程建立过程中，对模型施加恒定应力σ0，其中，针对瞬时弹性应变；由广义hooke定律得：式中ε
i
为模型i部分的应变；e为材料修正的弹性模量；针对黏弹性应变：由步骤s1中的本构方程的式(1)得：
式中ε
ii
为模型ii部分的应变；η1、η2分别为两个黏弹性元件的黏弹性系数；γ1、γ2为导数的阶数；针对黏弹塑性应变模型iii部分应变分为表征黏弹-黏塑性耦合作用产生的应变ε
′
(t)和黏塑性体应变ε5(t)；黏塑性体中黏弹性元件本构方程表达式为：式中σ5、ε5(t)为黏弹性元件所受应力、应变；η5为黏弹性元件的黏弹性系数；γ5为微积分阶数，γ5＞0则表示求导运算，γ5＜0则表示积分运算；排除材料的屈服硬化效应，设材料的屈服应力强度为σ
s
，则模型iii部分黏塑性本构方程为：式中<>为选择函数，当σ
0-σ
s
小于零时结果为零，当σ
0-σ
s
大于零时为σ
0-σ
s
；对式(18)进行laplace变换，再进行laplace逆变换可得黏塑性体的本构方程；黏弹-黏塑性耦合效应应变取与(ⅱ)黏弹性本构方程不同参数的分数阶黏弹性体应变表示，本构方程为：式中ε
′
(t)为模型iii部分黏弹-黏塑性耦合效应应变；η3、η4分别为两个黏弹性元件的黏弹性系数；γ3、γ4为导数的阶数；联立式(19)、式(20)，则模型iii部分的应变为：联立式(15)、式(16)、式(21)得到材料蠕变过程的本构方程，分两种情况，即：
4.根据权利要求1所述的分数阶蠕变本构方法，其特征在于：在步骤3中；s3.1，第一阶段，进行蠕变试验；s3.1.1，制作试件；s3.1.2，测量尺寸；s3.1.3，确定加载设备；s3.1.4，安装试件；s3.1.5，进行初始加载及测试；首先，依据拟合曲线使用目的，确定初始加载的大小，其中，依据施加的恒定荷载应力σ0与屈服强度σ
s
的大小关系，分为两种；当σ0<σ0时，得到i、ii部分；当σ0≥σ0时，得到i、ii、iii部分；然后，开始加载，并进行荷载、应变、时间参数的测量，完成加载的时间为t0，对应的应变为ε0；s3.1.6，蠕变测试；初始荷载加载完成后，保持荷载稳定不变，持续记录荷载稳定时间和对应的应变，获取持续的时间-应变(t
c-ε
c
)关系曲线。5.根据权利要求4所述的分数阶蠕变本构方法，其特征在于：s3.2，第二阶段，数据处理；s3.2.1,时间修正；将蠕变发生的时间设为零点，即应变(ε)-时间(t)关系曲线向右侧偏移t0，得到修正后的时间-应变(t
x-ε
x
)关系曲线；s3.2.2.数据拟合；包括σ0<σ
s
和σ0≥σ
s
两种工况；当σ0<σ
s
时，基于(t
x-ε
x
)关系曲线，利用式(22)的中σ0＜σ
s
情况对应的材料本构方程，利用l-m算法进行参数拟合或转化为目标优化的方式进行拟合，将得到对应的材料参数e、η1、η2、γ1、γ2，分别为弹性模量、材料黏弹性系数，求导阶数；当σ0≥σ
s
时，分三类分别进行；当屈服应力σ
s
与其对应的屈服应变ε
s
均已知：依据屈服应变ε
s
均在(t
x-ε
x
)关系曲线中确定ⅱ部分与iii部分的分界点对应的时间为t1，利用0～t1间的(t
x-ε
x
)关系曲线数据并按照式(22)的中σ0＜σ
s
情况进行拟合，确定参数e、η1、η2、γ1、γ2取值；将参数e、η1、η2、γ1、γ2值代入式(22)的中σ0≥σ
s
情况对应的材料本构方程，并利用(t
x-ε
x
)关系曲线t1～t
d
间的数据进行拟合，确定参数η3、η4、η5、γ3、γ4、γ5的取值；当屈服应力σ
s
与其对应的屈服应变ε
s
均未知：将试件在同等温度条件下，按照设定的加载速率开展单轴压缩试验，至试件发生破坏；加载过程中连续记录应力-应变测试数据；并根据应力-应变曲线确定材料的屈服应力σ
s
与其对应的屈服应变ε
s
；此时，按照屈服应力σ
s
与其对应的屈服应变ε
s
均已知的情况进行材料本构方程中待定常数的确定，从而确定材料的本构方程；当屈服应力σ
s
已知，屈服应变ε
s
均未知：策略一：是按照当屈服应力σ
s
与其对应的屈服应变ε
s
均未知的情况进行分析，重新开展单轴压缩试验，确定屈服应力与其对应的屈服应变；策略二：利用已有测试数据确定屈服应变；设整个曲线测点个数为n，第n个测点对应的时刻为t
n
，则利用拟合结果与测试结果误差
最小为目标，确定ⅱ部分与iii部分的分界点，建立的优化模型：简化式(23)优化模型的计算方法，在时刻0～t
n
的范围：计算应变的前向差分，并确定差分最小值对应的时刻t
m
及对应的应变ε
m
；确定屈服应变可能存在范围式中ε
max
为所测试的最大应变；k为经验系数，经验系数k取值为3-5，从而保证屈服应变所对应时间范围[t
x
,t
s
]，即为t
n
的取值范围；式(23)的求解方法采用遗传算法或蚁群算法，其屈服应变对应的时刻与变的前向差分。6.根据权利要求3所述的分数阶蠕变本构方法，其特征在于：在s4中，s4.1，拟合沥青混合料；其中，试验施加应力值为0.537mpa，屈服应力为0.358mpa，最终拟合曲线表达式，s4.2，拟合沥青砂，其中，试验施加应力值为0.2mpa，屈服应力为0.05mpa，最终拟合曲线表达式；
7.一种分数阶蠕变本构模型，其特征在于：借助于权利要求1-6任一项所述的方法进行搭建；所述的模型包括黏弹塑性本构模型。

技术总结
本发明涉及分数阶蠕变本构方法及模型，其方法包括以下步骤；S1，设立分数阶黏弹性体；S2，黏弹塑性本构模型；S3操作步骤；S4，本构模型验证。本发明设计合理、结构紧凑且使用方便。结构紧凑且使用方便。结构紧凑且使用方便。


技术研发人员：谢铠泽 赵佳 赵维刚 杨勇 田秀淑
受保护的技术使用者：石家庄铁道大学
技术研发日：2023.06.09
技术公布日：2023/9/23