科技论文报告会范文分析解析

更新时间:2023-03-08 04:53:24 阅读量: 人文社科 文档下载

说明:文章内容仅供预览,部分内容可能不全。下载后的文档,内容与下面显示的完全一致。下载之前请确认下面内容是否您想要的,是否完整无缺。

工程学院科技论文报告会参会论文

紫红色泥岩蠕变损伤本构关系研究

李男,徐辉,胡新丽,马越野

(中国地质大学(武汉)工程学院,湖北 武汉430074)

摘 要:岩体的蠕变性是影响边坡变形与长期稳定性的重要因素之一。以湖北巴东地区的五里堆2号滑坡为工程背景,对滑坡内典型的软弱夹层紫红色泥岩,开展室内岩石剪切蠕变试验,研究其蠕变变形特性,进行了蠕变本构模型的辨识,最终选定西原模型并通过计算得到其蠕变参数,结果显示西原模型能够较好地模拟泥岩的初始和稳定蠕变阶段,但无法模拟加速蠕变阶段。于是,在西原模型的基础上,引入损伤演化方程,从而建立了可以描述该泥岩全过程蠕变变形的损伤模型,并具体介绍了该模型参数的求解方法,拟合结果证明了新模型的正确和合理性。 关键词:岩体力学;泥岩;蠕变;损伤本构模型;参数求解方法

STUDY ON CREEP DAMAGE CONSTITUTIVE RELATION OF

FUCHSIA MUDSTONE

LI Nan, XU Hui, HU Xin-li

(College of Engineering, China University of Geosciences, Wuhan, Hubei, 430074, China)

Abstract: The time-dependent behavior of rock mass is one of the important factors that influence the deformation and long term stability of landslide. As an engineering example, Wu-li-dui landslide is studied,which is located in Ba-dong city, Hubei Province, with typical soft interlayer- mudstone. The creep deformation characteristics of rock are studied through laboratory shear creep tests. Then, the Xiyuan creep constitutive model is chosen from many ones, and the parameters are identified accordingly, the results show that Xiyuan model can simulate initialized and stable creep deformation, but cannot simulate accelerated creep deformation. Therefore, a new damage developing function is tumbled in Xiyuan creep model to found the damage creep model which can describe whole-process creep deformation characteristics, and the methods to calculate the model parameters are also given in detail. The results indicate that the damage creep model is very correct and reasonable.

Key words: rock mechanics; mudstone; creep; damage constitutive model; parameter calculating method

1.引言

工程实践与研究表明,岩石流变是岩石地下工程、构筑物基础、边坡及滑坡产生大变形乃至失稳的重要原因之一。目前,国内外因对岩土体流变特性研究不够而导致的延误施工甚至工程失败的先例不胜枚举,如1963年意大利瓦依昂(Vajont)库岸的蠕滑破坏,巴拿马运河滑坡区的片状粘土在运河边坡完工几年后突然崩塌等等。在第1届国际岩石力学会议上,Zischinsky用流变学模型描述了高边坡的变形,并指出岩土体的蠕变在高边坡变形中起重要作用,自此,国内外对边坡岩体流变特性的研究越来越多,成果也十分显著[1]~[12]。

五里堆滑坡区位于湖北省恩施市巴东县官渡口

基金项目:国家自然科学基金项目(40872175)

镇五里堆村,处于巫峡与西陵峡之间的长江北岸,与巴东新城区隔江相望,滑坡区总面积达0.476平方公里,经勘查资料证实其为多次滑移后形成的一个滑坡群。其下伏基岩为紫红色粉砂岩夹泥岩,软硬相间,单层厚度大于0.3m。

由于该滑坡滑床为软硬互层且软弱结构面发育,属于典型的“红层岩体”。相关研究表明,软弱岩层是红层岩体的薄弱部位,且具有明显的蠕变特性,其发育特征和空间形态对红层岩体的稳定性有重要的控制作用[13~14],根据对滑带土矿物成分的研究分析,可以认为该滑坡基本上是沿软弱泥岩夹层发生蠕滑变形的,因此很有必要对该软弱岩层进行蠕变试验并开展蠕变特性及其本构关系的研究,以为五里堆2号滑坡的形成机制和蠕滑机理的探索以

作者简介:李男,女,1989年生,中国地质大学(武汉)工程学院工程地质专业052073班本科生。 指导老师:胡新丽,中国地质大学(武汉)工程学院教授。

工程学院科技论文报告会参会论文

及治理工作的开展提供重要的理论和数据支持。 连续加载,当连续24h总变形量小于0.001mm时,开始加下一级荷载增量。试样取于五里堆2号滑坡的典型部位,岩性为紫红色泥岩,致密块状结构,共分为5组,其垂向应力分别为表盘上读数2,4,6,8,10MPa。取其中一组试验数据进行分析研究,如图1所示:

2.泥岩蠕变试验及其结果分析

2.1蠕变试验概况和蠕变特性分析

本蠕变试验仪器采用中国科学院岩土研究所研制的JQ系列剪切流变仪,精度0.001mm。试验操作严格按照规程[15]中的要求进行,加载方法采用分级

图1 0.617MPa法向应力下蠕变试验结果

Fig.1 Results of creep test for 0.617MPa norm stress

从图1可以看出,每一级应力施加瞬间,泥岩试件都会出现瞬时弹性应变,其量值随着剪应力水平的增大而增大。当剪应力较小时,试样主要以瞬时弹性应变为主,一般在1~2d内就能够达到稳定蠕变阶段。随着剪应力的增加,初始蠕变阶段所经历的时间逐渐加大,蠕变特征逐渐明显,蠕变总应变量也逐渐增大。当剪应力等于3.580MPa时,试样反常地经历了大约15个小时的非衰减蠕变过程,然后趋于稳定,其原因可能是泥岩试件发生了微小结构性破坏,随后经过重新调整后恢复稳定状态。当剪应力等于7.486MPa时,泥岩试件经历了初始蠕变阶段、稳定蠕变阶段和加速蠕变阶段后最终破坏。对于泥岩试件,其加速蠕变阶段经历的时间较为短暂,破坏前出现明显的征兆,表现为延性破坏形式。 2.2 长期剪切强度的确定

岩石的蠕变性与其所承受的应力水平密切相

关,当应力水平低于某一界限值时,蠕变过程中蠕变速率持续衰减直至为零,而当应力水平达到或高于该应力界限值时,就会出现稳态蠕变(??const?0)或加速蠕变现象,这一应力水平界限值称为岩石的长期强度[16]。目前,确定长期强度主要有以下两种方法[9]:一种是将梯级加载的试验结果在对数坐标下绘制应力-变形等时曲线,将对数坐标中曲线上明显拐点对应的应力作为长期强度,这是因为在非衰减蠕变荷载阶段变形值会急剧增长;另一种方法是根据流动速度来判断,绘制应力与应变速率的关系曲线,曲线和应力轴的交点即为长期强度。这里采用第一种方法来确定泥岩的长期剪切强度值,绘制的等时曲线如图2所示,由此确定0.617MPa垂直压力作用下泥岩的长期剪切强度值τs为6.30MPa。

第工程学院科技论文报告会参会论文

图2 应变-应力等时曲线

Fig.2 Isochronal curves of strain-stress

3 泥岩蠕变本构模型

3.1衰减蠕变阶段本构模型的辨识与拟合

在黏弹性体假设的前提下,利用Boltzmann迭加

原理对实验数据进行处理,可以得到各剪切力水平下的剪应变-时间蠕变曲线,如图3所示。

图3 不同剪应力下蠕变曲线及拟合结果

Fig.3 Creep curves of different shear stresses and their fitting results

在实际应用中,需要根据岩石或岩体的真实流变特征及变形性态来选择模型进行实际工程问题的分析。故根据室内或现场试验所得的数据进行蠕变模型的辨识就显得尤为重要。目前,蠕变模型的辨识主要有直接筛选法、后验排除法或者两种方法的综合运用[16]。本文将采用直接筛选法,即根据应变

-时间曲线特征并参照表1直接进行模型的辨识。

表1 常用蠕变模型的蠕变特性

Table1 creep properties of in common used models 模型 Maxwell模型 Kelvin模型 弹性应变 有 无 应变速率 不变 递减 黏性流动 有 无 工程学院科技论文报告会参会论文

H-K模型 有 递减 无 H|M模型 有 递减 无 Burgers模型 有 递减 有 黏塑性模型 无 不变 有 宾汉姆模型 有 不变 有 西原模型 有 递减 有 从图2可以看出:在任何一级剪切荷载作用下,都存在着瞬时弹性应变现象;当剪应力小于长期强度时,随着时间的增长,应变速率逐渐降低,近似呈负指数的形式趋于某一渐近线;当剪应力τ0=6.835MPa>τs时,曲线中出现了略微的黏性流动现象。根据上述分析,并参照表1,只有Burgers模型和西原模型符合泥岩的蠕变特性,又由于应力小于长期强度时Burgers模型有4个参数,而西原模型退化为三参量H-K模型,参数的增多会增加计算分析和实际应用的难度,因此,本文选用西原模型来模拟泥岩的蠕变过程,其拟合结果见图2。从图2中可以看出,西原模型能够较好地模拟泥岩的初始蠕变阶段和衰减或非衰减的稳定蠕变阶段,但不能模拟复杂的非线性加速蠕变阶段,那是因为剪应力大于长期强度时西原模型只是增加了一个简单的线性黏塑性元件,故只能模拟线性黏性流动现象。 3.2全过程损伤模型 3.2.1 损伤模型的建立

为了更好地模拟泥岩蠕变变形的全过程性态,本文将提出一种全过程蠕变损伤模型。

基于黏弹塑性介质模型理论,我们可以认为黏弹性部分与黏塑性部分为串联的关系,由此得到如下的关系式:

???e??p (1) ???e??p (2)

式中,σ为总应力,σe为黏弹性部分的应力,σp

为黏塑性部分的应力;ε为总应变,εe为黏弹性部分的应变,εp为黏塑性部分的应变。

将式(2)对时间求导后,可得:

???e??p (3) 对于黏弹性流变变形部分,其算子形式的通用表达式为:

?e?P(D)Q(D)? (4)

式中,P(D),Q(D)均为黏弹性模型算子函数,nP(D)??p?kk,nQ(D)?k?0?t?q?kk,D??为kkk?0?t?t对时间的微分算子,σs为岩石材料的屈服应力。 将式(4)对时间求导后,可得: ?e?P(D)Q(D)????s (5)

对于具有应变软化特性的岩石的黏塑性部分,则有

?p?1?(???s)??s (6)

2(1?w(t))?式中,w(t)为损伤演化函数,用以表征岩石损伤的演变规律,考虑到应变量的大小往往能够表征岩体的变形特征、微裂隙的发展以及判断岩石破坏与否,因此本文由文献[17]将损伤演化方程引申为:

w(t)?k?1(1?w(t))?(?) (7) c其中

??(?k2??(???1c?)????)cc?? (8) ??0??1c式中,εc为稳定蠕变与加速蠕变阶段的分界点

所对应的应变量。

从而可以解得:

?w(t)???1?exp?????k??k1??(?)2?1?????t?tc (9)

??c????0t?tc式中,tc为稳定蠕变与加速蠕变阶段的分界点的时刻。

对于西原模型,其流变算子函数分别为:

工程学院科技论文报告会参会论文

P(D)?1??1?DE0?E1?Q(D)?E0E1?E0?1????s?D(E?0?E1)E0?E1?? (10)

P(D)?1?(?2??1??2)???12DED?0E1E0E1???s?Q(D)?????122??ED1??式中,E0表示瞬时弹性模量,E1表示黏弹性模量,η1表示黏弹性系数。

由式(3)、(5)、(6)、(9)可推算得到一维流变全过程的微分本构方程:

E0?1??E0E1?E??1???????s1?E0E1?E0E1?E0???1?????1??(1??1??2????(11)

s?)?????s?E1E0E1E0E12?2(1?w(t))??3.2.2 损伤模型参数的求解

模型参数的求解主要有直接解析法和数值法。 i.

当σ<σs时,对式(11)进行Laplace变换和

逆变换可求解得到:

?(t)??E??(1?exp(?E1?t))???s (12)

0E11此时,西原模型退化为三参量H-K模型,其参数E0、E1和η1的求解可以用解析法完成。

由式(12)可得到

?(0)???E0?? (13)

?(?)?????E?0E1??从而可以解得:

E??0??(0)?E???? (14)

1??(?)????(?)??(0)?E?0?将式(12)变换为

E1?(?)??(t)???(?)??(0)?e??t1 (15)

对式(15)两边取对数,可得

In??(?)??(t)??In??(?)??(0)??E1?t (16)

1令

Yi?In??(?)??(ti)??a??E1???? (17)

b?In??(?)?1?(0)????则式(16)变为

Yi?ati?b (18)

对式(18)进行最小二乘线性回归法即可得到参

数a,从而解得

?1??E1a (19) ii. 当σ≥σs且t

?(t)??(E??E(1?exp(?E1t))????s)t (20) 01?1?2当t=0时,?(t)??,从而可以确定参数EE0。

0?E1t当t→∞时,e?1?0,所以应变-时间曲线趋

向于渐近线?(t)?????(???s)t,将试验数据进

E0E1?2行线性回归分析后得:

Y?a?bt (21)

式中,a??E??E,b?(???s),从而01?2可以确定参数E1和η2。

将式(21)减去式(20),可得

Y??(t)??Eexp(?E1t)) (22)

1?1对式(22)两边取对数,可得

In[Y??(t)]?In?E1E?t (23)

1?1对式(23)按照Qi?In[Yi??(ti)]?A?Bti进行

工程学院科技论文报告会参会论文

线性回归分析后可以得到

??A?In?E1? (24) ??B??E1??1?从而可以确定参数η1,至此,所有参数均已求得。 iii.

当σ≥σs且t≥tc时,式(11)变为

h??????zn?zn?1?n?1n?2?????k2????k1??n??1?E??????E1?1?? (27) s????c??e?zn??????1h12?zn?1?zn???????k2??2??k1??n?1??1?E????????E??1?s1??c???e?zn?1??????1?12??为了便于求解,对式(27)作如下变换

E1????s?E1??e ???1??12????k2?k1????1????c???? (25)

式(25)的解析解不易求出,因此本文采用数值方法来求解其参数k1,k2。

这里采用均值欧拉(Euler)方法。首先将试验数据按线性内插法变换成等步长的数据,即相邻两数据点间的距离xn+1-xn=h=const。为了求解式(25)这个二阶微分方程,先把它化为一阶微分方程组的形式,

d?然后进行求解。令z(t)?,则式(25)变为

dth??????zn?zn?1?n?1n?2??????k2???(28) ?k1??n?1??1?E????????E?s1??c???zn?1?zn?h?1e?zn?1??????112??????由于ε0,ε1,ε2……εn+1均可直接从曲线上得到,

从而可以根据式(28)上半部分求得z0,z1,z2……zn+1。

对式(28)下半部分移项后两边取对数,得到

k2E1????s????n?1???zn?1?znE1?In??zn?1??In?k1????1?(29)

?1??????h12??c???d???z(t)?(tc)??c?dt?????k2?(26) ?k1????1????s?dzE1????s?????c???E1?e?zz(t)?c?dt???1?2?12将向前欧拉法和向后欧拉法叠加求算术平均值后即为均值欧拉方法,应用该方法将式(26)变为

最后,对式(29)做非线性幂回归分析即可求得参数k1、k2,至此,全过程蠕变损伤本构模型的所有参数都已得到。

用上述方法对建立的蠕变全过程损伤本构模型进行参数的计算,然后对0.617MPa垂直压力作用下泥岩的蠕变曲线进行全过程模拟,其结果如表2和图4所示。从结果可以看出,本文所建立的蠕变全过程损伤本构模型与试验结果的吻合情况相当理想,这也彰显了该模型的合理性和正确性。

表2 蠕变全过程损伤本构模型计算参数

Table2 Calculating parameters on creep damage constitutive model τ0/MPa 0.325 0.976 1.627 2.278 2.929 3.580 4.231 4.882 5.533 6.184 6.835 7.486 E0/MPa 63.9253 126.85803 187.44033 247.06254 302.07126 351.6409 393.15413 435.54052 468.90919 508.34023 513.22537 537.54888 E1/MPa 337.3952 31812.38 17785.78 59568.61 47014.4 12634.23 44216.37 17106.8 60675.68 6586.47 31581.8 223256.8 η1/MPa.h 18.87374 7454.558 761752.2 465756.7 560473.2 126291.2 611635.5 771222.4 1488657 15217.07 25047.54 102236.7 η2/MPa.h 1974395.31 246637.793 k1 -2.150438 k2 -41.8005 R2 0.9634 0.9218 0.8125 0.7780 0.8900 0.9812 0.9702 0.9919 0.9010 0.9612 0.9841 0.9916

工程学院科技论文报告会参会论文

图4 全过程蠕变损伤本构模型拟合结果

Fig.4 Fitting results of mudstone for 0.617MPa norm stress on the whole process creep damage constitutive mode

4 结论

本文对五里堆2号滑坡中的软弱夹层紫红色泥岩进行了剪切蠕变试验,并对其蠕变特性和本构关系开展了研究,主要结论如下:

a) 选取一组泥岩剪切蠕变试验曲线作为本文研究分析的基础,首先对其剪切蠕变特性进行了分析,岩样破坏前出现明显的加速蠕变阶段,有显著的破坏征兆,表现为延性破坏形式。

b) 运用等时曲线确定在0.617MPa法向应力作用下该泥岩的长期剪切强度值为6.30MPa。

c) 采用直接筛选法进行蠕变模型的辨识,最终选用西原模型来描述泥岩的蠕变变形行为,结果显示:西原模型能够较好地模拟泥岩的初始蠕变和稳定蠕变或出现线性黏性流动的稳态蠕变现象,但不能模拟复杂的非线性加速蠕变现象。

d) 基于黏弹塑性介质模型理论,建立了一个能够模拟泥岩全过程蠕变变形行为的损伤模型,并具体介绍了模型参数的求解方法。将其运用到泥岩剪切蠕变试验曲线中,结果显示该模型能够很好地模拟泥岩的全过程蠕变变形行为,同时彰显了本文所提出的蠕变全过程损伤模型的正确性和合理性。

e) 获得了能够充分描述五里堆滑坡软弱夹层蠕变变形特征的本构模型参数值,为该滑坡的形成机制和启滑过程的数值模拟奠定了重要基础。

参 考 文 献

[1]

[5] [4] [3] [2]

Rock M ass of TGP’S Lock[J]. Journal of Yangtze River Scientific Research Institute, 1999, 16(2): 31-34. (in Chinese)

夏熙伦,徐平,丁秀丽.岩石流变特性及高边坡稳定性流变分析[J].岩石力学与工程学报,1996,15(4): 312- 322.

XIA Xi-lun, XU Ping, DING Xiu-li. Rheological characteristics of rock and stability rheological analysis for high slope[J].Chinese Journal of Rock Mechanics and Engineering,1996, 15(4):312-322. (in Chinese)

张友谊,刘俊新,冯涛. 峡口滑坡不同降雨条件下蠕变行为分析[J].铁道工程学报,2009,第1期:23-26. ZHANG You-yi, LIU Jun-xin, FENG Tao. Numerical Simulation of Creep Behavior of Xiakou Landslide Under Different Rainfall Conditions[J]. Journal of railway engineering society, 2009, NO.1: 23- 26. (in Chinese) 丁秀丽,付敬,刘建等.软硬互层边坡岩体的蠕变特性研究及稳定性分析[J].岩石力学与工程学报, 2005,24(19):3410-3417.

DING Xiu-li, FU Jing, LIU Jian, et, al. Study on creep behavior of alternatively distributed soft and hard rock layers and slope stability analysis[J].Chinese Journal of Rock Mechanics and Engineering,2005,24(19):3410- 3417. (in Chinese)

胡斌,冯夏庭,等.龙滩水电站左岸高边坡泥板岩体蠕变参数的智能反演[J].岩石力学与工程学报,2005,24(17):3064-3070.

HU Bin, FENG Xia-tin, et al. Intelligent inversion of creeping parameters of left bank high slope shale rock masses at long-tan hydropower station[J].Chinese Journal of Rock Mechanics and Engineering,2005, 24(17):3064- 3070. (in Chinese)

[6]

徐平,甘军,丁秀丽.三峡工程船闸高边坡岩体长期变形及稳定性有限元分析[J].长江科学院院报,1999,16(2):31-34.

XU Ping, GAN Jun, DING Xiu-li. FEM Analysis on Long-term Deformation and Stability for High Slope

赵法锁,张伯友,彭建兵,等.仁义河特大桥南桥台

工程学院科技论文报告会参会论文

边坡软岩流变性研究[J].岩石力学与工程学报,2002,21(10):1 527-1 532.

Zhao Fa-suo, Zhang Bo-you, Peng Jian-bing, et a1.Rheological study on soft rocks in the south abutment slope of Renyi River great bridge[J]. Chinese Journal of Rock Mechanics and Engineering,2002,21(10):1527-1532.(in Chinese)

[7]

韩爱果,聂德新,任光明,韩坤立.大型滑坡滑带土剪切流变特性研究[J].工程地质学报, 2001,9 (4):345-348.

HAN Ai-guo, NIE De-xin, REN Guang-ming, HAN Kun-li. Study on shear rheological behaviors of soil in slip zone of a large-scale landslide[J]. Journal of Engineering Geology, 2001, 9 (4):345-348.(in Chinese)

[8]

张强勇,向文,杨文东,张建国.坝区岩体蠕变参数反演与边坡开挖流变计算分析[J]. 武汉大学学报(工学版),2008,41(5):72-76.

ZHANG Qiang-yong, XIANG Wen, YANG Wen-dong, ZHANG Jian-guo. Creep parameters inversion for dam zone rockmass and rheological computational analysis of slope excavation[J]. Engineering Journal of Wuhan University, 2008, 41(5): 72-76. (in Chinese)

[9]

严绍军,项伟,唐辉明,满作武,徐瑞春.大岩淌滑坡滑带土蠕变性质研究[J].岩土力学,2008, 29(1): 58- 62.

YAN Shao-jun, XIANG Wei, TANG Hui-ming, et, al. Research on creep behavior of slip band soil of Dayantang landslide[J]. Rock and Soil Mechanics, 2008,29(1):58-62. (in Chinese)

[10] 许 强,曾裕平.具有蠕变特点滑坡的加速度变化特征

及临滑预警指标研究[J].岩石力学与工程学报,2009,28(6):1099-1106.

XU Qiang,ZENG Yu-ping.Research on acceleration variation characteristics of creep landslide and early-warning prediction indicator of critical sliding[J]. Chinese Journal of Rock Mechanics and Engineering,2009,28(6):1099- 1106.(in Chinese)

[11] Cheng-Wei Chen, Hani Salim, John J. Bowders et al.

Creep Behavior of Recycled Plastic Lumber in Slope Stabilization Applications[J]. Journal of materials in civil engineering, 2007,19(2): 130-138.

[12] Jin Feng, Zhang Chuhan, M.ASCE, et al. Creep Modeling

in Excavation Analysis of a High Rock Slope[J]. Journal of materials in civil engineering, 2003,129(9):849-857.

[13] 程强,寇小兵,黄绍槟.中国红层的分布及地质环境特

征[J] .工程地质学报,2004,12(1):34-40.

CHENG Qiang, KOU Xiao-bing, HUANG Shao-bin. The distributions and geologic environment characteristics of red beds in China[J]. Journal of Engineering Geology, 2004, 12(1):34-40.(in Chinese)

[14] 程强,周德培,寇小兵.红层软岩地区公路建设中的主

要岩土工程问题[C].第8次全国岩石力学与工程学术大会论文集.北京:科学出版社,2004: 128-131. CHENG qiang, ZHOU De-pei, KOU Xiao-bing. The main geotechnical engineering problems of highway construction in red beds[C]. Proceedings of the 8th National Conference of Rock Mechanics and Engineering. Beijing: Science Press, 2004: 128-131. (in Chinese)

[15] 中华人民共和国国家标准编写组. SL264–2001 水利

水电工程岩石试验规程[S]. 北京:中国水利水电出版社,2001.

The National Standards Compilation Group of People?s Republic of China. SL264–2001 Water conservancy and hydropower project order on rock test[S]. Beijing:China Water Power Press,2001.(in Chinese)

[16] 王芝银,李云鹏. 岩体流变理论及其数值模拟[M].北

京:科学出版社,2008.

WANG Zhi-yin, LI Yun-peng. Rock rheological theory and numerical simulation[M]. Beijing:Science Press, 2008. (in Chinese)

[17] 徐海滨,朱维申,白世伟.岩体黏弹塑-损伤本构模型及

其有限元分析[J].岩土力学,1992,13(1):11-20. XU Hai-bin,ZHU Wei-shen,BAI Shi- wei.Visco-elastic- plastic constitutive model and finite-element analysis of rock[J]. Rock and Soil Mechanics, 1992, 13(1):11-20. (in Chinese)

本文来源:https://www.bwwdw.com/article/6hg.html

Top