水-热-力耦合作用下寒区弧底梯形渠道结构优化设计 Structural optimization design for trapezoidal canal with arc-bottom during thermal-moisture-mechanical coupling in cold regions 1 2 3 corresp first-author 王正中，教授，博士生导师，从事旱区寒区水工程安全研究；E-mail：wangzz0910@163.com 王正中，教授，博士生导师，从事旱区寒区水工程安全研究；E-mail：wangzz0910@163.com wangzz0910@163.com 1 2 1 2 1 2 1 2 4 西北农林科技大学 旱区寒区水工程安全研究中心，陕西 杨凌，712100 西北农林科技大学 旱区寒区水工程安全研究中心，陕西 杨凌，712100 Cold and Arid Regions Water Engineering Safety Research Center, Northwest Agriculture & Forestry University, Yangling 712100, China Cold and Arid Regions Water Engineering Safety Research Center, Northwest Agriculture & Forestry University, Yangling 712100, China 西北农林科技大学 旱区农业水土工程教育部重点实验室，陕西 杨凌，712100 西北农林科技大学 旱区农业水土工程教育部重点实验室，陕西 杨凌，712100 Key Laboratory of Agricultural Soil and Water Engineering in Arid and Semiarid Areas of Ministry of Education, Northwest Agriculture & Forestry University, Yangling 712100, China Key Laboratory of Agricultural Soil and Water Engineering in Arid and Semiarid Areas of Ministry of Education, Northwest Agriculture & Forestry University, Yangling 712100, China 中国科学院西北生态环境资源研究院 冻土工程国家重点实验室，甘肃 兰州，730000 中国科学院西北生态环境资源研究院 冻土工程国家重点实验室，甘肃 兰州，730000 State Key Laboratory of Frozen Soil Engineering, Northwest Institute of Eco-Environment and Resources, Lanzhou 730000, China State Key Laboratory of Frozen Soil Engineering, Northwest Institute of Eco-Environment and Resources, Lanzhou 730000, China 宁夏回族自治区水利科学研究院，宁夏 银川，750021 宁夏回族自治区水利科学研究院，宁夏 银川，750021 Ningxia Institute of Hydraulic Research, Yinchuan 750021, China Ningxia Institute of Hydraulic Research, Yinchuan 750021, China 杨幼平 为合理优化设计寒区弧底梯形渠道的形体结构，基于冻土水热力三场耦合及冻土-衬砌相互作用的渠道冻胀数值模型，结合分层序列法将渠道断面水力参数最优解集作为变量空间，以衬砌结构的强度、刚度及几何构造为约束条件，以衬砌适应冻胀变形能力为目标函数，建立寒区弧底梯形渠道水力-抗冻胀双优数学模型。随后对寒区各类典型渠道工程进行双优结构标准化设计，得到弧底梯形渠道在不同负温、地下水埋深、土质、断面规模下的双优边坡系数、实佳比及衬砌厚度的标准参数，供工程设计参考。研究结果表明：双优结构尺寸参数均随地下水埋深的增大而减小、随基土冻胀率的增大而增大、随冬季负温的降低而增大、随断面规模的增大而增大。以新疆某渠道工程为例，优化后的断面实佳比为1.02，较原设计值增大1.8%；最大法向冻胀量为2.52 cm，较原设计值增大20.0%；最大拉应力为0.95 MPa，较原设计值减小60.4%；衬砌整体柔度为305.47 cm/MPa，较原设计值增大42.3%，水力性能及抗冻胀性能均达到最优。 In order to optimize the shape structure of trapezoidal canal with arc-bottom in cold region, a bi-objective optimization mathematical model was established for the water conveyance performance and anti-frost heave performance of the trapezoidal canal with arc-bottom in cold region. This mathematical method was based on the canal frost heave numerical model of thermal-moisture-mechanical coupling of frozen soil and the interaction between frozen soil and lining and the hierarchical sequence method. It was taken the optimal solution set of the hydraulic parameters of the canal's section as the variable space, was taken the strength, stiffness and geometric structure of the lining structure as the constraint conditions, and was taken the ability of the lining to adapt to the frost heave deformation as the objective function. Subsequently, standard designs of bi-objective optimal structures were carried out for various typical canal projects in cold regions, and a standard parameter table including the bi-objective optimal slope coefficient, the ratio of practical economic section to hydraulic optimum section and the lining thickness of the trapezoidal canal with arc-bottom was obtained at different temperatures, groundwater depths, soil quality and cross-section scale. The standard parameter table can provide reference for engineering design. The results show that the size parameters of the bi-objective optimal structure decrease with the increase of groundwater buried depth, increase with the increase of foundation soil frost heave rate, increase with the decrease of negative temperature in winter and the increase of section size. Taking a canal project in Xinjiang as an example, the ratio of practical economic section to hydraulic optimum section is 1.02, which is 1.8% higher than the original design value. The maximum normal frost heave is 2.52 cm, which is 20.0% higher than the original design value. The maximum tensile stress is 0.95 MPa, which is 60.4% lower than the original design value. The overall flexibility of the lining is 305.47 cm/MPa, which is 42.3% higher than the original design value. The water delivery performance and anti-frost heave performance of the optimized canal are both in the optimal state. 冻土 水-热-力耦合 冻土-衬砌相互作用 弧底梯形渠道 水力-抗冻胀双优 frozen soil thermal-moisture-mechanical coupling interaction between frozen soil and lining trapezoidal canal with arc-bottom hydraulic and anti-frost heave bi-objective optimization 渠道作为主要输水建筑物之一，因其造价低，运行管理方便，被广泛应用于大型长距离调水工程 [ 1 1 ] 。然而，我国北方季节性冻土区分布较广，以冻土为地基的混凝土衬砌渠道结构强度低、刚度较大，抗拉及适应冻胀变形能力差，衬砌结构极易发生冻胀、开裂、脱空甚至渠坡滑塌等工程问题，严重影响水资源利用效率及工程安全运行 [ 2 2 ] 。弧底梯形渠道是一种输水输沙性能优良的渠道型式，且因其弧底的“拱效应”，表现出结构复位能力强、冻胀适应性强等受力性能好的优点，被广泛应用于北方旱寒灌区及供水工程。但是，当前寒区弧底梯形渠道的相关设计规范 [ 3 3 - 4 4 3-4 ] 多依据工程经验，缺乏科学的结构计算模型，更无综合考虑渠道水沙与抗冻胀最优的定量优化设计方法，导致冻害现象频繁发生，影响工程的正常运行。如陕西泾惠渠和山西大禹渡扬水工程弧底梯形渠道在冻胀力作用下渠底中心处出现了裂缝 [ 5 5 - 6 6 5-6 ] 。 因此，亟待提出一种输水输沙能力与适应冻胀变形能力俱佳的弧底梯形渠道结构标准化设计方法。当前渠道断面优化研究基本可以分为两类：第一类为水力优化研究，其研究内容主要集中在渠道水力计算方面，旨在提高渠道输水效率，兼顾考虑渠道占地面积与工程造价等经济成本因素，并采用传统的优化算法对衬砌渠道进行断面优化设计 [ 7 7 - 10 10 7-10 ] ；在此类研究中基本不考虑冻胀或仅凭工程经验进行抗冻胀方面的设计。第二类为抗冻胀优化研究 [ 11 11 - 14 14 11-14 ] ，如：辛英华等 [ 11 11 ] 采用工程力学计算方法，对大U形渠道断面进行优化设计；刘东等 [ 12 12 ] 根据文献[ 15 15 ]提出的冻胀破坏工程力学模型，引入计算效率较高的粒子群算法，以经济成本为目标，对梯形渠道进行优化设计。安鹏等 [ 13 13 ] 依据半经验的冻胀量计算公式，对弧底梯形及弧形坡脚梯形渠道进行冻胀破坏分析，并以工程造价为目标函数对其进行优化设计。王玉宝等 [ 14 14 ] 建立了U形渠道冻胀破坏工程力学模型，并提出了一种小U型渠道断面设计优化方法。在此类研究中，渠道冻胀力多为经验或半经验公式，未考虑冻土水热力三场耦合及衬砌与冻土相互作用的冻胀破坏动态机制，也未兼顾渠道输水输沙与抗冻胀性能双优目标。 针对上述问题，本文作者基于冻土水热力三场耦合及衬砌-冻土相互作用的渠道冻胀数值模型，结合分层序列法，提出衬砌整体柔度指标以描述其适应冻胀变形的能力，建立寒区弧底梯形渠道水力-抗冻胀双优标准化设计方法；利用多物理场耦合软件COMSOL进一步对寒区不同冬季负温、地下水埋深、渠基土质及断面规模下的弧底梯形渠道进行双优结构设计，探究双优渠道的结构尺寸参数在上述因素下的变化规律。最后，以新疆某大型弧底梯形渠道为例进行水力抗冻胀双优设计，验证该方法的合理性与准确性，以期为寒区弧底梯形渠道工程的结构抗冻胀设计提供科学依据与理论参考。 1 冻土水热力三场耦合冻胀数学模型 1 1.1 　冻土水-热耦合控制方程 2 渠道表面与外界环境进行热量交换，利用Fourier公式计算土体内部温度场随时间的变化，并建立考虑水冰相变的二维非稳态热传导方程： <math display="block" id="M1"><mi>ρ</mi><msub><mrow><mi>C</mi></mrow><mrow><mi mathvariant="normal">e</mi><mi mathvariant="normal">q</mi></mrow></msub><mfrac><mrow><mo>∂</mo><mi>T</mi></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>=</mo><mfrac><mrow><mo>∂</mo></mrow><mrow><mo>∂</mo><mi>x</mi></mrow></mfrac><mo stretchy="false">(</mo><msub><mrow><mi>λ</mi></mrow><mrow><mi mathvariant="normal">e</mi><mi mathvariant="normal">q</mi></mrow></msub><mfrac><mrow><mo>∂</mo><mi>T</mi></mrow><mrow><mo>∂</mo><mi>x</mi></mrow></mfrac><mo stretchy="false">)</mo><mo>+</mo><mfrac><mrow><mo>∂</mo></mrow><mrow><mo>∂</mo><mi>y</mi></mrow></mfrac><mo stretchy="false">(</mo><msub><mrow><mi>λ</mi></mrow><mrow><mi mathvariant="normal">e</mi><mi mathvariant="normal">q</mi></mrow></msub><mfrac><mrow><mo>∂</mo><mi>T</mi></mrow><mrow><mo>∂</mo><mi>y</mi></mrow></mfrac><mo stretchy="false">)</mo><mo>+</mo><mi>L</mi><msub><mrow><mi>ρ</mi></mrow><mrow><mi mathvariant="normal">i</mi></mrow></msub><mfrac><mrow><mo>∂</mo><msub><mrow><mi>θ</mi></mrow><mrow><mi mathvariant="normal">i</mi><mn mathvariant="normal">2</mn></mrow></msub></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706456&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706455&type= 68.15666962 8.12800026 式中： ρ 和 ρ i 分别为基土及冰密度，kg/m 3 ； C eq 为渠基土的等效比定压热容，J/(kg·K)； λ eq 为土体等效导热系数，W/(m 2 ·K)； L 为水冰相变潜热，取 333 kJ/kg； T 为温度，℃； θ i2 为迁移水冻结成冰量，cm 3 /cm 3 。 考虑土颗粒、水及冰体积含量的变化及水冰相变产生的潜热对土体热参数的影响，依据混合律公式对渠基土的热物理参数进行简化计算 [ 16 16 ] ： <math display="block" id="M2"><msub><mrow><mi>λ</mi></mrow><mrow><mi mathvariant="normal">e</mi><mi mathvariant="normal">q</mi></mrow></msub><mo>=</mo><msub><mrow><mi>θ</mi></mrow><mrow><mi mathvariant="normal">s</mi></mrow></msub><msub><mrow><mi>λ</mi></mrow><mrow><mi mathvariant="normal">s</mi></mrow></msub><mo>+</mo><msub><mrow><mi>θ</mi></mrow><mrow><mi mathvariant="normal">w</mi></mrow></msub><msub><mrow><mi>λ</mi></mrow><mrow><mi mathvariant="normal">w</mi></mrow></msub><mo>+</mo><msub><mrow><mi>θ</mi></mrow><mrow><mi mathvariant="normal">i</mi><mn mathvariant="normal">1</mn></mrow></msub><msub><mrow><mi>λ</mi></mrow><mrow><mi mathvariant="normal">i</mi></mrow></msub><mo>+</mo><msub><mrow><mi>θ</mi></mrow><mrow><mi mathvariant="normal">i</mi><mn mathvariant="normal">1</mn></mrow></msub><msub><mrow><mi>λ</mi></mrow><mrow><mi mathvariant="normal">i</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706464&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706457&type= 43.77266693 3.72533321 <math display="block" id="M3"><msub><mrow><mi>C</mi></mrow><mrow><mi mathvariant="normal">e</mi><mi mathvariant="normal">q</mi></mrow></msub><mo>=</mo><mfrac><mrow><mn mathvariant="normal">1</mn></mrow><mrow><mi>ρ</mi></mrow></mfrac><mo stretchy="false">(</mo><msub><mrow><mi>θ</mi></mrow><mrow><mi mathvariant="normal">s</mi></mrow></msub><msub><mrow><mi>C</mi></mrow><mrow><mi>p</mi><mo>,</mo><mi mathvariant="normal">s</mi></mrow></msub><mo>+</mo><msub><mrow><mi>θ</mi></mrow><mrow><mi mathvariant="normal">w</mi></mrow></msub><msub><mrow><mi>C</mi></mrow><mrow><mi>p</mi><mo>,</mo><mi mathvariant="normal">w</mi></mrow></msub><mo>+</mo><msub><mrow><mi>θ</mi></mrow><mrow><mi mathvariant="normal">i</mi><mn mathvariant="normal">1</mn></mrow></msub><msub><mrow><mi>C</mi></mrow><mrow><mi>p</mi><mo>,</mo><mi mathvariant="normal">i</mi></mrow></msub><mo stretchy="false">)</mo><mo>-</mo><mi>L</mi><msub><mrow><mi>ρ</mi></mrow><mrow><mi mathvariant="normal">i</mi></mrow></msub><mfrac><mrow><mo>∂</mo><msub><mrow><mi>θ</mi></mrow><mrow><mi mathvariant="normal">i</mi><mn mathvariant="normal">1</mn></mrow></msub></mrow><mrow><mo>∂</mo><mi>T</mi></mrow></mfrac></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706459&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706443&type= 63.16133118 8.12800026 <math display="block" id="M4"><msub><mrow><mi>θ</mi></mrow><mrow><mi mathvariant="normal">s</mi></mrow></msub><mo>+</mo><msub><mrow><mi>θ</mi></mrow><mrow><mi mathvariant="normal">w</mi></mrow></msub><mo>+</mo><msub><mrow><mi>θ</mi></mrow><mrow><mi mathvariant="normal">i</mi><mn mathvariant="normal">1</mn></mrow></msub><mo>=</mo><mn mathvariant="normal">1</mn></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706460&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706466&type= 21.92866707 3.72533321 式中： λ s 、 λ i 和 λ w 分别为土颗粒、冰相及水相的导热系数，W/(m 2 ·K)； θ s 、 θ w 及 θ i1 分别为土颗粒、原位水及原位冰的体积分数，cm 3 /cm 3 ； C p ,s 、 C p ,w 和 C p ,i 分别为土颗粒、水及冰的比定压热容，J/(kg·K)。 在渠基土内，距离地下水位较近的土体含水量近似饱和，水分迁移可按饱和土渗流进行计算，且在冻结过程中的水分补给量较大，冻胀较为剧烈；而较远处土体含水量较低，水分迁移符合非饱和土渗流规律，水分迁移量小且冻胀量可以忽略。为此，本文不考虑冻结过程中距离地下水位大于毛细水上升高度的渠基土水分迁移，而小于或等于毛细水上升高度的渠基土内水分迁移按饱和土渗流进行计算。此种处理所得冻胀量计算结果符合《渠系工程抗冻胀设计规范》 [ 4 4 ] 中的现场实测规律。因此，结合达西定律与质量守恒原理建立渠基土迁移水运动方程，如 式(5) 式(5) 和(6)所示 [ 17 17 ] 。 <math display="block" id="M5"><mfrac><mrow><mo>∂</mo><msub><mrow><mi>θ</mi></mrow><mrow><mi mathvariant="normal">i</mi><mn mathvariant="normal">2</mn></mrow></msub></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>=</mo><mfrac><mrow><msub><mrow><mi>ρ</mi></mrow><mrow><mi mathvariant="normal">w</mi></mrow></msub></mrow><mrow><msub><mrow><mi>ρ</mi></mrow><mrow><mi mathvariant="normal">i</mi></mrow></msub></mrow></mfrac><mfrac><mrow><mo>∂</mo><msub><mrow><mi>θ</mi></mrow><mrow><mi mathvariant="normal">w</mi><mn mathvariant="normal">1</mn></mrow></msub></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>=</mo><mo>-</mo><mi>k</mi><mfrac><mrow><msub><mrow><mi>ρ</mi></mrow><mrow><mi mathvariant="normal">w</mi></mrow></msub></mrow><mrow><msub><mrow><mi>ρ</mi></mrow><mrow><mi mathvariant="normal">i</mi></mrow></msub></mrow></mfrac><mo stretchy="false">(</mo><mfrac><mrow><msup><mrow><mo>∂</mo></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msup><mi>h</mi></mrow><mrow><mo>∂</mo><msup><mrow><mi>x</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msup></mrow></mfrac><mo>+</mo><mfrac><mrow><msup><mrow><mo>∂</mo></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msup><mi>h</mi></mrow><mrow><mo>∂</mo><msup><mrow><mi>y</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msup></mrow></mfrac><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706468&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706446&type= 56.81133270 8.97466660 <math display="block" id="M6"><mi>k</mi><mo>=</mo><mfenced open="{" close="" separators="|"><mrow><mtable columnalign="left"><mtr><mtd><msub><mrow><mi>k</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub><mo stretchy="false">[</mo><mn mathvariant="normal">1</mn><mo>-</mo><mo stretchy="false">(</mo><mi>T</mi><mo>-</mo><msub><mrow><mi>T</mi></mrow><mrow><mi mathvariant="normal">f</mi></mrow></msub><msup><mrow><mo stretchy="false">)</mo><mo stretchy="false">]</mo></mrow><mrow><mi>β</mi></mrow></msup><mo>;</mo><mi>T</mi><mo>≤</mo><msub><mrow><mi>T</mi></mrow><mrow><mi mathvariant="normal">f</mi></mrow></msub><mo>,</mo><mi>y</mi><mo><</mo><msub><mrow><mi>y</mi></mrow><mrow><mi mathvariant="normal">s</mi><mi mathvariant="normal">e</mi><mi mathvariant="normal">p</mi></mrow></msub></mtd></mtr><mtr><mtd><msub><mrow><mi>k</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub><mo>;</mo><mi>T</mi><mo>></mo><msub><mrow><mi>T</mi></mrow><mrow><mi mathvariant="normal">f</mi></mrow></msub><mo>,</mo><mi>y</mi><mo><</mo><msub><mrow><mi>y</mi></mrow><mrow><mi mathvariant="normal">s</mi><mi mathvariant="normal">e</mi><mi mathvariant="normal">p</mi></mrow></msub></mtd></mtr><mtr><mtd><mn mathvariant="normal">0</mn><mo>,</mo><mi>y</mi><mo>≥</mo><msub><mrow><mi>y</mi></mrow><mrow><mi mathvariant="normal">s</mi><mi mathvariant="normal">e</mi><mi mathvariant="normal">p</mi></mrow></msub></mtd></mtr></mtable></mrow></mfenced></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706470&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706469&type= 51.47733307 16.76399994 式中： θ w1 为迁移水含量，cm 3 /cm 3 ； k 0 为未冻土的渗透系数，m/s； h 为压力水头，m； T f 为冰水相变温度，℃，本文设为0 ℃； y sep 为冰透镜体坐标位置，m； β 为试验参数。 研究表明 [ 18 18 ] ，冻土中土颗粒周围仍然存在部分未冻水膜，为下部未冻区的水分迁移提供通道。水分迁移驱动力则可由Clepyron方程推导为温度梯度的函数，故未冻土中的水迁移至冻结锋面处并相变成冰所产生的热量方程如下 [ 19 19 ] ： <math display="block" id="M7"><mi>L</mi><msub><mrow><mi>ρ</mi></mrow><mrow><mi>i</mi></mrow></msub><mfrac><mrow><mo>∂</mo><msub><mrow><mi>θ</mi></mrow><mrow><mi>i</mi><mn mathvariant="normal">2</mn></mrow></msub></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>=</mo><mo>-</mo><mi>k</mi><mfrac><mrow><msub><mrow><mi>ρ</mi></mrow><mrow><mi mathvariant="normal">w</mi></mrow></msub><msup><mrow><mi>L</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msup></mrow><mrow><mi>g</mi><msub><mrow><mi>T</mi></mrow><mrow><mi mathvariant="normal">f</mi></mrow></msub></mrow></mfrac><mo>·</mo><mo>∇</mo><mo stretchy="false">(</mo><mo>∇</mo><mi>T</mi><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706480&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706493&type= 39.28533173 8.72066593 由上述分析可知：方程(1)、(5)及(7)构成了渠基土水-热耦合控制方程，描述了冻胀过程中的热量传递与水分迁移的动态演变规律。 1.2 　冻土冻胀本构方程 2 土体在冻结过程中垂直于土体温度梯度方向形成冰透镜体，并随着温度的降低，冻结锋面持续向下推进，在冻结锋面内的任意方向处具有对称性及面内各向同性 [ 20 20 ] 。据此建立考虑初始应变的渠基冻土本构方程 [ 21 21 ] 。 <math display="block" id="M8"><mfenced open="[" close="]" separators="|"><mrow><mtable columnalign="left"><mtr><mtd><msub><mrow><mi>ε</mi></mrow><mrow><mi>x</mi></mrow></msub><mo>-</mo><msub><mrow><mi>ε</mi></mrow><mrow><mi>x</mi><mn mathvariant="normal">0</mn></mrow></msub></mtd></mtr><mtr><mtd><msub><mrow><mi>ε</mi></mrow><mrow><mi>y</mi></mrow></msub><mo>-</mo><msub><mrow><mi>ε</mi></mrow><mrow><mi>y</mi><mn mathvariant="normal">0</mn></mrow></msub></mtd></mtr><mtr><mtd><msub><mrow><mi>ε</mi></mrow><mrow><mi>x</mi><mi>y</mi></mrow></msub><mo>-</mo><msub><mrow><mi>ε</mi></mrow><mrow><mi>x</mi><mi>y</mi><mn mathvariant="normal">0</mn></mrow></msub></mtd></mtr></mtable></mrow></mfenced><mo>=</mo><mfenced open="[" close="]" separators="|"><mrow><mtable><mtr><mtd><mfrac><mrow><mn mathvariant="normal">1</mn></mrow><mrow><msub><mrow><mi>E</mi></mrow><mrow><mi mathvariant="normal">f</mi></mrow></msub></mrow></mfrac></mtd><mtd><mo>-</mo><mfrac><mrow><mi>ν</mi></mrow><mrow><msub><mrow><mi>E</mi></mrow><mrow><mi mathvariant="normal">f</mi></mrow></msub></mrow></mfrac></mtd><mtd><mn mathvariant="normal">0</mn></mtd></mtr><mtr><mtd><mo>-</mo><mfrac><mrow><mi>ν</mi></mrow><mrow><msub><mrow><mi>E</mi></mrow><mrow><mi mathvariant="normal">f</mi></mrow></msub></mrow></mfrac></mtd><mtd><mfrac><mrow><mn mathvariant="normal">1</mn></mrow><mrow><msub><mrow><mi>E</mi></mrow><mrow><mi mathvariant="normal">f</mi></mrow></msub></mrow></mfrac></mtd><mtd><mn mathvariant="normal">0</mn></mtd></mtr><mtr><mtd><mn mathvariant="normal">0</mn></mtd><mtd><mn mathvariant="normal">0</mn></mtd><mtd><mfrac><mrow><mo stretchy="false">(</mo><mn mathvariant="normal">1</mn><mo>+</mo><mi>ν</mi><mo stretchy="false">)</mo></mrow><mrow><msub><mrow><mi>E</mi></mrow><mrow><mi mathvariant="normal">f</mi></mrow></msub></mrow></mfrac></mtd></mtr></mtable></mrow></mfenced><mfenced open="[" close="]" separators="|"><mrow><mtable columnalign="left"><mtr><mtd><msub><mrow><mi>σ</mi></mrow><mrow><mi>x</mi></mrow></msub></mtd></mtr><mtr><mtd/></mtr><mtr><mtd><msub><mrow><mi>σ</mi></mrow><mrow><mi>y</mi></mrow></msub></mtd></mtr><mtr><mtd/></mtr><mtr><mtd><msub><mrow><mi>τ</mi></mrow><mrow><mi>x</mi><mi>y</mi></mrow></msub></mtd></mtr></mtable></mrow></mfenced></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706471&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706482&type= 59.68999863 29.37933159 式中： ε x 0 和 ε y 0 分别为 x 和 y 方向的初始应变； E f 为随温度呈非线性变化的冻土弹性模量，MPa； <math id="M9"><msub><mrow><mi>σ</mi></mrow><mrow><mi>x</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706485&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706484&type= 3.04800010 3.72533321 和 <math id="M10"><msub><mrow><mi>σ</mi></mrow><mrow><mi>y</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706486&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706473&type= 3.04800010 3.72533321 分别为 x 和 y 方向的正应力； <math id="M11"><msub><mrow><mi>τ</mi></mrow><mrow><mi>x</mi><mi>y</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706488&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706498&type= 3.47133350 3.72533321 为剪应力。 以冻胀率表示物理方程中的初始应变，定义垂直于温度梯度方向的初始应变，利用坐标转换得到 x 和 y 坐标平面下的初始应变，进而描述平面应变问题下的土体冻胀过程。方程如下 [ 22 22 ] ： <math display="block" id="M12"><mfenced open="{" close="" separators="|"><mrow><mtable columnalign="left"><mtr><mtd><msub><mrow><mi>ε</mi></mrow><mrow><mi>x</mi><mn mathvariant="normal">0</mn></mrow></msub><mo>=</mo><msub><mrow><mi>ε</mi></mrow><mrow><mi mathvariant="normal">g</mi><mi mathvariant="normal">r</mi><mi mathvariant="normal">a</mi><mi mathvariant="normal">d</mi><mi>T</mi></mrow></msub><mo>⋅</mo><msup><mrow><mi mathvariant="bold-italic">l</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msup></mtd></mtr><mtr><mtd><msub><mrow><mi>ε</mi></mrow><mrow><mi>y</mi><mn mathvariant="normal">0</mn></mrow></msub><mo>=</mo><msub><mrow><mi>ε</mi></mrow><mrow><mi mathvariant="normal">g</mi><mi mathvariant="normal">r</mi><mi mathvariant="normal">a</mi><mi mathvariant="normal">d</mi><mi>T</mi></mrow></msub><mo>⋅</mo><msup><mrow><mi mathvariant="bold-italic">n</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msup></mtd></mtr><mtr><mtd><msub><mrow><mi>ε</mi></mrow><mrow><mi>x</mi><mi>y</mi><mn mathvariant="normal">0</mn></mrow></msub><mo>=</mo><msub><mrow><mi>ε</mi></mrow><mrow><mi mathvariant="normal">g</mi><mi mathvariant="normal">r</mi><mi mathvariant="normal">a</mi><mi mathvariant="normal">d</mi><mi>T</mi></mrow></msub><mo>⋅</mo><mi mathvariant="bold-italic">l</mi><mo>⋅</mo><mi mathvariant="bold-italic">n</mi></mtd></mtr></mtable></mrow></mfenced></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706501&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706476&type= 23.96066666 17.27199936 式中： ε grad T 为渠基土的冻胀率， <math id="M13"><mi>T</mi><mo>≤</mo><msub><mrow><mi>T</mi></mrow><mrow><mi mathvariant="normal">f</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706491&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706509&type= 8.38199997 3.80999994 ； l 和 n 为温度梯度向量的方向余弦向量。 综上可知：方程(1)、(5)、(7)及(9)共同构成描述冻土冻胀的水-热-力三场耦合数学模型。LIU等 [ 23 23 ] 通过室内、现场试验结合有限元程序验证了此模型的正确性；王正中等 [ 19 19 ] 利用此模型提出了寒区渠道冻胀破坏的尺寸效应。本文将以此模型为基础结合实际边界条件及多物理耦合有限元软件对各类典型弧底梯形渠道进行水力-抗冻胀双优设计研究。 2 弧底梯形渠道水力-抗冻胀双优数学模型 1 渠道水力-抗冻胀双优设计在数学上属于典型的多目标优化问题，若直接对该问题进行优化求解，则计算参数范围广，计算成本较大；且计算结果精度较低，难以运用到实际工程设计中 [ 24 24 - 25 25 24-25 ] 。故本文采用分层序列法进行优化设计，在优化设计过程中，先对渠道进行水力优化设计，得到渠道水力最优断面的尺寸参数解集；再将其作为抗冻胀设优化计算的求解空间，从而得到满足水力与抗冻胀双目标优化的渠道衬砌结构形式。 2.1 　渠道水力最优结构设计 2 仅针对渠道横断面进行抗冻胀优化设计，忽略渠道纵向坡降与表面糙率的影响。根据文献[ 3 3 ]推导弧底梯形渠道的水力最优及经济实用断面表达式，如式( 10 10 )~( 13 13 )所示。 <math display="block" id="M14"><msub><mrow><mi>H</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub><mo>=</mo><mn mathvariant="normal">1.414</mn><msup><mrow><mo stretchy="false">(</mo><mfrac><mrow><msub><mrow><mi>ω</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub></mrow><mrow><mi>θ</mi><mo>+</mo><mn mathvariant="normal">2</mn><mi>m</mi></mrow></mfrac><mo stretchy="false">)</mo></mrow><mrow><mn mathvariant="normal">1</mn><mo>/</mo><mn mathvariant="normal">2</mn></mrow></msup></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706492&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706510&type= 32.00400162 7.36600018 <math display="block" id="M15"><mi>α</mi><mo>=</mo><mfrac><mrow><mi>ω</mi></mrow><mrow><msub><mrow><mi>ω</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub></mrow></mfrac></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706525&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706513&type= 10.66800022 8.12800026 <math display="block" id="M16"><msub><mrow><mi>K</mi></mrow><mrow><mi mathvariant="normal">r</mi></mrow></msub><mo>=</mo><mfrac><mrow><mi>r</mi></mrow><mrow><mi>H</mi></mrow></mfrac></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706515&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706538&type= 10.92199993 7.28133297 <math display="block" id="M17"><mi>H</mi><mo>=</mo><mfrac><mrow><mo stretchy="false">(</mo><mi>θ</mi><mo>+</mo><mn mathvariant="normal">2</mn><mi>m</mi><mo stretchy="false">)</mo><msup><mrow><mi>α</mi></mrow><mrow><mn mathvariant="normal">5</mn><mo>/</mo><mn mathvariant="normal">2</mn></mrow></msup></mrow><mrow><mo stretchy="false">(</mo><mn mathvariant="normal">2</mn><mi>m</mi><mo>-</mo><mn mathvariant="normal">2</mn><mroot><mrow><mn mathvariant="normal">1</mn><mo>+</mo><msup><mrow><mi>m</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msup><mo>+</mo><mi>θ</mi></mrow><mrow/></mroot><mo stretchy="false">)</mo><msub><mrow><mi>K</mi></mrow><mrow><mi mathvariant="normal">r</mi></mrow></msub><mo>+</mo><mn mathvariant="normal">2</mn><mroot><mrow><mn mathvariant="normal">1</mn><mo>+</mo><msup><mrow><mi>m</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msup></mrow><mrow/></mroot></mrow></mfrac><msub><mrow><mi>H</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706529&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706528&type= 64.76999664 10.41399956 式中： H 0 和 H 分别为水力最优及经济实用断面水深，m； ω 0 和 ω 分别为水力最优及经济实用断面的面积，m 2 ； m 为坡度系数； θ 为弧底段的圆心角， θ =2arctan(1/ m )，rad； α 为实际过水断面面积与水力最佳断面的面积比，简称实佳比； K r 为弧底半径与水深之比； r 为弧底半径，m。 由式( 10 10 )~( 13 13 )可以看出：当断面面积 ω 一定时，选定一组 m 、 α 便可确定水力最优及经济断面水深 H 0 、 H 。此外，将衬砌板厚度 d 也作为优化过程中的独立设计参数。因此，弧底梯形衬砌渠道断面的基本设计参数分别为 m 、 α 、 d 。依据相关规范及实际工程经验 [ 4 4 ] ，各尺寸参数的解集空间为： m =[1，3]， α =[1，1.04]， d =[6，20] cm。 2.2 　水力-抗冻胀双优数学模型建立及计算 2 寒区混凝土衬砌渠道的抗冻胀性能评价指标主要包括法向冻胀位移及截面正应力 [ 3 3 - 4 4 3-4 ] ，从“削减、适应冻胀”的抗冻胀思路以及“允许一定冻胀位移”的工程设计标准出发，认为法向冻胀位移在规范允许值以内时，截面正应力越小，衬砌结构越能适应冻胀变形。为此，本文综合上述两种指标，构建衬砌结构整体柔度指标，表征寒区弧底梯形衬砌渠道对基土冻胀变形的适应性，如 式(14) 式(14) 所示。 <math display="block" id="M18"><mi>C</mi><mo>=</mo><mfrac><mrow><msub><mo>∫</mo><mrow><msub><mrow><mi>S</mi></mrow><mrow><mi mathvariant="normal">u</mi></mrow></msub><mo>+</mo><msub><mrow><mi>S</mi></mrow><mrow><mi mathvariant="normal">d</mi></mrow></msub></mrow></msub><mrow><mi>w</mi><mi mathvariant="normal">d</mi><mi>S</mi></mrow></mrow><mrow><msub><mo>∫</mo><mrow><msub><mrow><mi>S</mi></mrow><mrow><mi mathvariant="normal">u</mi></mrow></msub><mo>+</mo><msub><mrow><mi>S</mi></mrow><mrow><mi mathvariant="normal">d</mi></mrow></msub></mrow></msub><mrow><mfenced open="|" close="|" separators="|"><mrow><msub><mrow><mi>σ</mi></mrow><mrow><mi mathvariant="normal">n</mi></mrow></msub></mrow></mfenced><mi mathvariant="normal">d</mi><mi>S</mi></mrow></mrow></mfrac></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706532&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706544&type= 25.23066711 17.10266685 式中： C 为衬砌结构整体柔度指标，表示衬砌结构整体单位应力下的法向位移，其值越大，表明衬砌结构适应冻胀变形的能力越强，反之，则适应能力越弱，cm/MPa； w 为衬砌板各点处的法向冻胀位移，cm； σ n 为截面正应力，正值表示拉应力，负值表示压应力，MPa； S u 和 S d 分别为衬砌结构上、下表面的展开长度，m。 综上所述，建立以柔度指标为目标函数、以法向冻胀位移及正应力为约束条件、以水力最优尺寸参数解集为求解空间的双优数学模型，如 式(15) 式(15) 所示。 <math display="block" id="M19"><mtable columnalign="left"><mtr><mtd><mi>F</mi><mo stretchy="false">(</mo><mi>m</mi><mo>,</mo><mi>α</mi><mo>,</mo><mi>d</mi><mo stretchy="false">)</mo><mo>=</mo><mi mathvariant="normal">m</mi><mi mathvariant="normal">a</mi><mi mathvariant="normal">x</mi><mi>C</mi><mo stretchy="false">(</mo><mi>m</mi><mo>,</mo><mi>α</mi><mo>,</mo><mi>d</mi><mo stretchy="false">)</mo></mtd></mtr><mtr><mtd><mi mathvariant="normal">s</mi><mo>.</mo><mi mathvariant="normal">t</mi><mo>.</mo><mfenced open="{" close="" separators="|"><mrow><mtable columnalign="left"><mtr><mtd><mi>m</mi><mo>=</mo><mo stretchy="false">[</mo><mn mathvariant="normal">1,3</mn><mo stretchy="false">]</mo></mtd></mtr><mtr><mtd><mi>α</mi><mo>=</mo><mo stretchy="false">[</mo><mn mathvariant="normal">1,1.04</mn><mo stretchy="false">]</mo></mtd></mtr><mtr><mtd><mi>d</mi><mo>=</mo><mo stretchy="false">[</mo><mn mathvariant="normal">6,20</mn><mo stretchy="false">]</mo><mtext> </mtext><mi mathvariant="normal">c</mi><mi mathvariant="normal">m</mi></mtd></mtr><mtr><mtd><mi>s</mi><mo>≤</mo><mn mathvariant="normal">3.0</mn><mtext> </mtext><mi mathvariant="normal">c</mi><mi mathvariant="normal">m</mi></mtd></mtr><mtr><mtd><mo>-</mo><mn mathvariant="normal">9.6</mn><mtext> </mtext><mi mathvariant="normal">M</mi><mi mathvariant="normal">P</mi><mi mathvariant="normal">a</mi><mo>≤</mo><msub><mrow><mi>σ</mi></mrow><mrow><mi mathvariant="normal">n</mi></mrow></msub><mo>≤</mo><mn mathvariant="normal">1.1</mn><mtext> </mtext><mi mathvariant="normal">M</mi><mi mathvariant="normal">P</mi><mi mathvariant="normal">a</mi></mtd></mtr></mtable></mrow></mfenced></mtd></mtr></mtable></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706535&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706522&type= 43.34933472 29.54866600 由 式(15) 式(15) 可知：衬砌渠道形体结构优化设计属于高度约束的非线性规划问题，拟采用COMSOL优化模块中的Nelder-Mead无梯度优化算法进行优化计算。其中，目标函数的相对收敛容差 δ 设为0.01；最大计算迭代次数设为1 000；约束处理方式采用罚函数法。 3 考虑冻土-衬砌相互作用的有限元模型 1 为充分考虑冻土与衬砌的接触效应以及冻胀过程中的相互作用机制，以新疆某大型弧底梯形衬砌渠道为工程背景，采用大型有限元软件COMSOL对上述方程进行联立求解，得到温度-水分-基土-衬砌结构相互耦合作用的渠道冻胀模型，并明确寒区弧底梯形渠道在不同形体结构下的冻胀破坏规律。 3.1 　考虑冻土-衬砌接触的界面本构模型 2 渠基土与衬砌结构在负温作用下会“黏结”在一起并具有一定的冻结强度，并同时受到基土冻胀力作用，当渠基冻土与衬砌作用力超过法向或切向冻结强度时将会发生脱离或滑移 [ 26 26 ] 。基于此，采用修正的无厚度弹性薄层单元对其进行模拟。弹性薄层单元的本质是在冻土与衬砌结构的接触面处设置具有一定刚度的法向与切向弹簧单元，弹簧刚度可随其拉伸量发生非线性变化，并且该模型能够根据相互接触结构的相对位移计算二者间的反作用力，再根据反作用力调整二者的相对位移。其本构方程如下： <math display="block" id="M20"><mfenced open="{" close="}" separators="|"><mrow><mtable columnalign="left"><mtr><mtd><msubsup><mrow><mi>σ</mi></mrow><mrow><mi mathvariant="normal">n</mi></mrow><mrow><mi mathvariant="normal">'</mi></mrow></msubsup></mtd></mtr><mtr><mtd><msub><mrow><mi>σ</mi></mrow><mrow><mi mathvariant="normal">t</mi></mrow></msub></mtd></mtr></mtable></mrow></mfenced><mo>=</mo><mfenced open="[" close="]" separators="|"><mrow><mtable><mtr><mtd><msub><mrow><mi>K</mi></mrow><mrow><mi mathvariant="normal">n</mi></mrow></msub></mtd><mtd><mn mathvariant="normal">0</mn></mtd></mtr><mtr><mtd><mn mathvariant="normal">0</mn></mtd><mtd><msub><mrow><mi>K</mi></mrow><mrow><mi mathvariant="normal">t</mi></mrow></msub></mtd></mtr></mtable></mrow></mfenced><mfenced open="{" close="}" separators="|"><mrow><mtable><mtr><mtd><msub><mrow><mi>u</mi></mrow><mrow><mi mathvariant="normal">n</mi><mi mathvariant="normal">l</mi></mrow></msub><mo>-</mo><msub><mrow><mi>u</mi></mrow><mrow><mi mathvariant="normal">n</mi><mi mathvariant="normal">s</mi></mrow></msub></mtd></mtr><mtr><mtd><msub><mrow><mi>u</mi></mrow><mrow><mi mathvariant="normal">t</mi><mi mathvariant="normal">l</mi></mrow></msub><mo>-</mo><msub><mrow><mi>u</mi></mrow><mrow><mi mathvariant="normal">t</mi><mi mathvariant="normal">s</mi></mrow></msub></mtd></mtr></mtable></mrow></mfenced></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706568&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706536&type= 40.21666718 9.56733322 式中， <math id="M21"><msubsup><mrow><mi>σ</mi></mrow><mrow><mi mathvariant="normal">n</mi></mrow><mrow><mi mathvariant="normal">'</mi></mrow></msubsup></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706571&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706557&type= 3.13266683 4.06400013 和 σ t 分别为法向和切向弹性薄层反力； K n 和 K t 分别为薄层单元法向和切向单位长度的刚度系数，kN/m 3 ； u nl 和 u tl 分别为衬砌板的法向与切向位移，m； u ns 和 u ts 分别为土体法向和切向位移，m。根据冻土-衬砌相互作用特点，对薄层单元法向和切向刚度系数进行修正，方程如下： <math display="block" id="M22"><msub><mrow><mi>K</mi></mrow><mrow><mi mathvariant="normal">n</mi></mrow></msub><mo>=</mo><mfenced open="{" close="" separators="|"><mrow><mtable columnalign="left"><mtr><mtd><msub><mrow><mi>E</mi></mrow><mrow><mi mathvariant="normal">f</mi></mrow></msub><mo>,</mo><msub><mrow><mi>σ</mi></mrow><mrow><mi mathvariant="normal">n</mi></mrow></msub><mo><</mo><msub><mrow><mi>τ</mi></mrow><mrow><mi mathvariant="normal">f</mi></mrow></msub></mtd></mtr><mtr><mtd><mn mathvariant="normal">0</mn><mo>,</mo><msub><mrow><mi>σ</mi></mrow><mrow><mi>n</mi></mrow></msub><mo>≥</mo><msub><mrow><mi>τ</mi></mrow><mrow><mi mathvariant="normal">f</mi></mrow></msub></mtd></mtr></mtable></mrow></mfenced></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706574&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706572&type= 23.36800003 9.65200043 <math display="block" id="M23"><msub><mrow><mi>K</mi></mrow><mrow><mi mathvariant="normal">t</mi></mrow></msub><mo>=</mo><mfenced open="{" close="" separators="|"><mrow><mtable columnalign="left"><mtr><mtd><msubsup><mrow><mi>K</mi></mrow><mrow><mi mathvariant="normal">t</mi></mrow><mrow><mi mathvariant="normal">'</mi></mrow></msubsup><mo>,</mo><msub><mrow><mi>σ</mi></mrow><mrow><mi mathvariant="normal">t</mi></mrow></msub><mo><</mo><msub><mrow><mi>τ</mi></mrow><mrow><mi mathvariant="normal">f</mi></mrow></msub></mtd></mtr><mtr><mtd><mi>κ</mi><mo>·</mo><mfrac><mrow><mfenced separators="|"><mrow><msub><mrow><mi>u</mi></mrow><mrow><mi mathvariant="normal">t</mi><mi mathvariant="normal">l</mi></mrow></msub><mo>-</mo><msub><mrow><mi>u</mi></mrow><mrow><mi mathvariant="normal">t</mi><mi mathvariant="normal">s</mi></mrow></msub></mrow></mfenced></mrow><mrow><mfenced open="|" close="|" separators="|"><mrow><msub><mrow><mi>u</mi></mrow><mrow><mi mathvariant="normal">t</mi><mi mathvariant="normal">l</mi></mrow></msub><mo>-</mo><msub><mrow><mi>u</mi></mrow><mrow><mi mathvariant="normal">t</mi><mi mathvariant="normal">s</mi></mrow></msub></mrow></mfenced><mo>+</mo><mi>ε</mi><mi>'</mi></mrow></mfrac><mo>,</mo><msub><mrow><mi>σ</mi></mrow><mrow><mi mathvariant="normal">t</mi></mrow></msub><mo>≥</mo><msub><mrow><mi>τ</mi></mrow><mrow><mi mathvariant="normal">f</mi></mrow></msub></mtd></mtr></mtable></mrow></mfenced></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706588&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706576&type= 41.90999985 16.76399994 式中： E f 为渠基冻土弹性模量，MPa； τ f 为冻土与衬砌的冻结强度，MPa； <math id="M24"><msubsup><mrow><mi>K</mi></mrow><mrow><mi mathvariant="normal">t</mi></mrow><mrow><mi mathvariant="normal">'</mi></mrow></msubsup></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706620&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706601&type= 3.64066648 4.06400013 为基土-衬砌间未达到冻结强度时的剪切刚度，MPa/m； κ 为冻土与衬砌在滑动状态下的摩擦因数； <math id="M25"><mi>ε</mi><mi>'</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706591&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706581&type= 2.28600001 2.62466669 为大于0的极小数。 基于上述界面本构模型，王羿等 [ 2 2 ] 针对寒区衬砌渠道在不同双膜防渗层布设形式下的冻胀受力与变位进行分析，得出寒区双膜防渗层的合理布设形式；江浩源等 [ 22 22 ] 提出了寒区大型衬砌渠道的“适缝”概念，并对其布设方式进行探究。故上述接触本构模型可应用于工程计算。 3.2 　基本参数选取及网格划分 2 新疆某大型弧底梯形渠道断面面积为 63.40 m 2 ，渠深为5.5 m，坡度系数为1.5，实佳比为1.002。衬砌材料采用C20现浇混凝土，衬砌板厚度为10 cm。具体尺寸及有限元网格划分如 图1 图1 所示。 数据单元：m。 依据土体恒温层厚度研究 [ 27 27 ] ，从渠底向下取10 m作为恒温层，渠顶各向两侧延伸4.5 m作为左右边界。根据现场实测资料，渠基土为低液限黏土，土体孔隙率为0.55。考虑渠道排水系统布置及灌溉水渗漏情况，地下水位距渠底以下1.3 m，土体初始含水率为20%，超过土体初始冻胀含水率，渗透系数取6×10 -7 m/s，试验参数 β 取-8 [ 19 19 ] 。冻土与衬砌板接触面的参数 <math id="M26"><msubsup><mrow><mi>K</mi></mrow><mrow><mi mathvariant="normal">t</mi></mrow><mrow><mi mathvariant="normal">'</mi></mrow></msubsup></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706620&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706601&type= 3.64066648 4.06400013 取120 MPa/m、 τ f 取 0.15 MPa、 f 取0.8 [ 22 22 ] 。将混凝土视作各向同性材料，弹性模量取2.1×10 4 MPa，泊松比取0.23。未冻土弹性模量取15 MPa，冻土弹性模量随温度变化规律如 表1 表1 所示。未冻土及冻土的热物理参数根据 式(2) 式(2) 及 式(3) 式(3) 进行计算，材料具体热物理参数见 表2 表2 [ 28 28 - 29 29 28-29 ] 。 温度/℃ 弹性模量/MPa 泊松比 0 15 0.33 -1 19 0.33 -2 26 0.33 -3 33 0.33 -5 46 0.33 材料 导热系数/ (W·m -2 ·K -1 ) 比定压热容/ (kJ·kg -1 ·K -1 ) 密度/(kg·m -3 ) 水 0.58 4.20 1 000 冰 2.20 2.10 916 土颗粒 1.54 0.92 2 700 混凝土 1.58 0.97 2 400 依据文献[ 4 4 ]，渠基土含水率与土质及地下水位相关，采用分段函数对不同计算点处的土体含水率进行拟合，如 式(19) 式(19) 所示；渠基土体冻胀率与土质、冻深及地下水位有关，依次对低液限黏土、粉土及含细粒砂土的冻胀率进行拟合，如式( 20 20 )~( 22 22 )所示。 <math display="block" id="M27"><msub><mrow><mi>θ</mi></mrow><mrow><mi mathvariant="normal">w</mi></mrow></msub><mo>=</mo><mfenced open="{" close="" separators="|"><mrow><mtable columnalign="left"><mtr><mtd><msub><mrow><mi>θ</mi></mrow><mrow><mi mathvariant="normal">s</mi></mrow></msub><mo>,</mo><mn mathvariant="normal">0</mn><mo>≤</mo><mi mathvariant="normal">Δ</mi><mi>h</mi><mo>≤</mo><msub><mrow><mi>h</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub></mtd></mtr><mtr><mtd><msub><mrow><mi>θ</mi></mrow><mrow><mi mathvariant="normal">s</mi></mrow></msub><mo>⋅</mo><msup><mrow><mi mathvariant="normal">e</mi></mrow><mrow><mo stretchy="false">[</mo><mo>-</mo><mn mathvariant="normal">0.08</mn><mo stretchy="false">(</mo><mi mathvariant="normal">Δ</mi><mi>h</mi><mo>-</mo><msub><mrow><mi>h</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub><mo stretchy="false">)</mo><mo stretchy="false">]</mo></mrow></msup><mo>,</mo><mi mathvariant="normal">Δ</mi><mi>h</mi><mo>></mo><msub><mrow><mi>h</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub></mtd></mtr></mtable></mrow></mfenced></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706633&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706625&type= 42.50266647 10.32933331 <math display="block" id="M28"><msub><mrow><mi>ε</mi></mrow><mrow><mi mathvariant="normal">g</mi><mi mathvariant="normal">r</mi><mi mathvariant="normal">a</mi><mi mathvariant="normal">n</mi><mi mathvariant="normal">d</mi></mrow></msub><mo>=</mo><mn mathvariant="normal">1.038</mn><mtext> </mtext><mn mathvariant="normal">6</mn><msup><mrow><mi mathvariant="normal">e</mi></mrow><mrow><mo>-</mo><mn mathvariant="normal">0.013</mn><msub><mrow><mi>Z</mi></mrow><mrow><mi mathvariant="normal">w</mi></mrow></msub></mrow></msup><msup><mrow><msub><mrow><mi>Z</mi></mrow><mrow><mi mathvariant="normal">d</mi></mrow></msub></mrow><mrow><mo>-</mo><mn mathvariant="normal">0.292</mn><mtext> </mtext><mn mathvariant="normal">5</mn></mrow></msup><mo>,</mo><msup><mrow><mi>R</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msup><mo>=</mo><mn mathvariant="normal">0.998</mn><mtext> </mtext><mn mathvariant="normal">9</mn></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706610&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706627&type= 62.73799896 5.24933338 <math display="block" id="M29"><msub><mrow><mi>ε</mi></mrow><mrow><mi mathvariant="normal">g</mi><mi mathvariant="normal">r</mi><mi mathvariant="normal">a</mi><mi mathvariant="normal">n</mi><mi mathvariant="normal">d</mi></mrow></msub><mo>=</mo><mn mathvariant="normal">1.855</mn><mtext> </mtext><mn mathvariant="normal">7</mn><msup><mrow><mi mathvariant="normal">e</mi></mrow><mrow><mo>-</mo><mn mathvariant="normal">0.014</mn><msub><mrow><mi>Z</mi></mrow><mrow><mi mathvariant="normal">w</mi></mrow></msub></mrow></msup><msup><mrow><msub><mrow><mi>Z</mi></mrow><mrow><mi mathvariant="normal">d</mi></mrow></msub></mrow><mrow><mo>-</mo><mn mathvariant="normal">0.435</mn><mtext> </mtext><mn mathvariant="normal">0</mn></mrow></msup><mo>,</mo><msup><mrow><mi>R</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msup><mo>=</mo><mn mathvariant="normal">0.998</mn><mtext> </mtext><mn mathvariant="normal">7</mn></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706645&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706637&type= 62.73799896 5.24933338 <math display="block" id="M30"><msub><mrow><mi>ε</mi></mrow><mrow><mi mathvariant="normal">g</mi><mi mathvariant="normal">r</mi><mi mathvariant="normal">a</mi><mi mathvariant="normal">n</mi><mi mathvariant="normal">d</mi></mrow></msub><mo>=</mo><mn mathvariant="normal">0.106</mn><mtext> </mtext><mn mathvariant="normal">3</mn><msup><mrow><mi mathvariant="normal">e</mi></mrow><mrow><mo>-</mo><mn mathvariant="normal">0.019</mn><msub><mrow><mi>Z</mi></mrow><mrow><mi mathvariant="normal">w</mi></mrow></msub></mrow></msup><msup><mrow><msub><mrow><mi>Z</mi></mrow><mrow><mi mathvariant="normal">d</mi></mrow></msub></mrow><mrow><mn mathvariant="normal">0.040</mn><mtext> </mtext><mn mathvariant="normal">7</mn></mrow></msup><mo>,</mo><msup><mrow><mi>R</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msup><mo>=</mo><mn mathvariant="normal">0.976</mn><mtext> </mtext><mn mathvariant="normal">3</mn></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706648&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706639&type= 61.55266953 5.24933338 式中： θ w 为不同点处含水率； θ s 为土体饱和含水率，与土体孔隙率有关；Δ h 为计算点到地下水位处线的距离，m； h 0 为影响土体冻结的临界地下水埋深，低液限黏土、粉土及含细粒砂土的 <math id="M31"><msub><mrow><mi>h</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706642&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706641&type= 3.21733332 3.72533321 分别取2.0、1.5及1.0 m； Z w 为地下水埋深，m； Z d 为冻深，m； <math id="M32"><msub><mrow><mi>ε</mi></mrow><mrow><mi mathvariant="normal">g</mi><mi mathvariant="normal">r</mi><mi mathvariant="normal">a</mi><mi mathvariant="normal">n</mi><mi mathvariant="normal">d</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706676&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706650&type= 6.68866634 3.72533321 为冻胀率； <math id="M33"><msup><mrow><mi>R</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msup></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706679&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706653&type= 3.55599999 2.96333337 为拟合方差。 3.3 　边界条件及初始条件 2 温度边界条件：渠道上边界与外界环境进行对流换热引起渠基土的持续降温，采用牛顿冷却定律进行描述 [ 25 25 ] ，如 式(23) 式(23) 所示。 <math display="block" id="M34"><mi mathvariant="bold-italic">n</mi><mo>⋅</mo><mo stretchy="false">(</mo><mi>λ</mi><mo>∇</mo><mi>T</mi><mo stretchy="false">)</mo><mo>=</mo><msub><mrow><mi>h</mi></mrow><mrow><mi mathvariant="normal">c</mi></mrow></msub><mo stretchy="false">(</mo><msub><mrow><mi>T</mi></mrow><mrow><mi mathvariant="normal">o</mi><mi mathvariant="normal">u</mi><mi mathvariant="normal">t</mi></mrow></msub><mo>-</mo><mi>T</mi><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706682&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706656&type= 32.51200104 4.48733330 式中： n 为模型边界法向向量； T out 和 T 分别为外部环境温度和模型上边界温度，℃； h c 为对流换热系数，W/(m 2 ·K)。 h c 与渠道衬砌表面的风速有关，关系式如下： <math display="block" id="M35"><msub><mrow><mi>h</mi></mrow><mrow><mi mathvariant="normal">c</mi></mrow></msub><mo>=</mo><mn mathvariant="normal">3.06</mn><mi>v</mi><mo>+</mo><mn mathvariant="normal">4.11</mn></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706669&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706690&type= 24.21466637 3.72533321 根据现场监测结果，渠顶处的平均风速为1.85 m/s，渠底平均风速为0.8 m/s，考虑渠道断面形式的特点，渠坡上各点处的风速以二次函数方式进行差值计算。 依据渠道当地气象站2019年10月1日至2020年9月30日的日平均温度监测值，采用正弦函数进行拟合 [ 23 23 ] ，计算2019年11月至2020年3月的日平均气温值，并作为外界环境温度 T out ，如 式(25) 式(25) 所示。 <math display="block" id="M36"><msub><mrow><mi>T</mi></mrow><mrow><mi mathvariant="normal">o</mi><mi mathvariant="normal">u</mi><mi mathvariant="normal">t</mi></mrow></msub><mo>=</mo><msub><mrow><mi>a</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub><mo>+</mo><msub><mrow><mi>b</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub><mi mathvariant="normal">s</mi><mi mathvariant="normal">i</mi><mi mathvariant="normal">n</mi><mo stretchy="false">(</mo><mfrac><mrow><mn mathvariant="normal">2</mn><mi mathvariant="normal">π</mi></mrow><mrow><mn mathvariant="normal">366</mn></mrow></mfrac><mo stretchy="false">(</mo><mi>t</mi><mo>+</mo><msub><mrow><mi>c</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub><mo stretchy="false">)</mo><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706693&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=44706670&type= 43.18000031 7.28133297 式中： a 1 、 b 1 和 c 1 为拟合参数，不同地区下的参数取值及拟合方差如 表3 表3 所示。 灌区 a 1 b 1 c 1 R 2 新疆奎屯灌区 8.13 20.10 174.0 0.941 3 内蒙古河套灌区 9.04 16.94 174.7 0.938 4 宁夏引黄灌区 11.25 14.24 177.0 0.935 8 模型左右边界按热绝缘边界处理，下边界恒温层的温度取当地多年平均气温值6.5 ℃ [ 22 22 ] ；渠基土内部各点的初始温度由模型下边界温度与模型上边界初始温度沿深度方向进行线性插值。 位移边界条件：模型上边界为自由边界；左右边界设为辊支撑；下边界设为固定边界。在冻胀变形计算中，对不同结构尺寸下的渠道模型计算地应力平衡状态为初始状态，初始地应力通过自重应力获得。 3.4 　计算思路与方案 2 计算方案分为三个部分：1) 以新疆某大型弧底梯形渠道为研究背景，利用上述三场耦合模型探究不同结构尺寸下的渠道温度、变形及截面正应力分布，明晰寒区弧底梯形渠道的冻胀破坏规律。2) 利用水力-抗冻胀双优设计方法研究典型地区、工程规模和不同水热力土参数下的双优设计成果及规律。3) 利用提出的双优设计方法对1)中的原型渠道进行优化计算实验，对比分析优化前后的变形及受力结果，验证其合理性。 4 结果分析与讨论 1 4.1 　寒区弧底梯形渠道冻胀破坏机理分析 2 4.1.1　不同渠槽形状下的渠道温度场分布规律 3 忽略衬砌板厚度微小变化对渠道稳定温度场的影响，仅探究不同断面参数组合的渠道温度场分布，如 图2 图2 所示。由 图2 图2 可知：等温线在不同渠道表面处的分布均较为密集，且随着渠基深度的增加分布逐渐稀疏，表明渠道表面温度低且梯度较大。同时，不同断面形式下的渠顶冻深均比渠底的大，这是由于渠顶处风速较大，对流换热系数较大，且该位置为双向冻结区域；同时，渠顶距离地下水位较远，由 式(2) 式(2) 及(3)确定的导热系数大且比热容较小，导致其在相同时间内的冻结速率较大，相反渠底处则导热系数小且比热容大，冻结速率小，冻深发展较慢。 (a) m =1， α =1.00；(b) m =2， α =1.04 进一步提取不同断面形式下渠顶、渠底处的冻深计算结果，如 图3 图3 所示。由 图3 图3 可知：随着 m 与 α 的逐渐增大，渠道断面由窄深式逐渐过渡至宽浅式，渠顶处冻深由1.52 m逐渐降低至1.08 m，渠底处冻深则由0.50 m不断增大至0.95 m，在二者的综合作用下，渠道整体冻深逐渐趋于一致。 4.1.2　不同渠槽形状下的衬砌变形及正应力分布规律 3 由前文可知，不同的渠槽形状将引起渠道温度场分布的改变，并直接影响衬砌结构的冻胀变形及受力分布特征，故提取衬砌厚度 d =10 cm、不同断面参数组合下的法向冻胀位移及截面正应力分布，分别如 图4 图4 和 图5 图5 所示。 (a) 上表面；(b) 下表面 注： 拉为正，压为负，虚线表示混凝土抗拉强度。 由 图4 图4 可知：不同断面形式下的渠道冻胀变形均表现为整体上抬、弧底向上隆起、两侧坡板向内凸起的变形特征。对于弧底段而言，其最大法向冻胀位移位于渠底中心处，而渠坡段则位于渠顶。当渠槽形状发生改变时，渠道表面的法向冻胀位移将产生显著变化。其中，当 m 保持不变、 α 逐渐增大时，因弧底段长度增加，弧底段特别是弧底中心处的法向冻胀显著增大；而当 α 保持不变、 m 逐渐增大时，因渠道边坡变缓且坡板长度增加，整体法向冻胀位移均呈明显增大的趋势。由此可知，随着断面参数 m 、 α 的增大，渠道断面将由窄深式转变为宽浅式，法向冻胀量均呈现出增大的趋势，冻胀变形释放程度较大。 结合截面正应力分布( 图5 图5 )可知，不同断面形式下的衬砌表面均存在两个拉应力极大值点。其中，最大极值点是由下部基土向上冻胀隆起而引起的，始终位于渠底中心上表面；次大极值点则由衬砌间相互挤压作用而引起的，位于衬砌下表面，具体位置则与坡板及弧底板长度有关。当 m 保持不变，实佳比 α 逐渐增大时，弧底段长度增大，“拱效应”增强，冻胀变形释放较多，拉应力极值均由应力限值外收缩至应力限值内。同时，次大极值点随着实佳比 α 的增大由坡板向弧底中心处转移。当实佳比 α 不变、边坡系数 m 不断增大时，衬砌对基土的冻胀约束减弱，冻胀变形增大，衬砌表面正应力减小，且随着 m 的增大次极值点由弧底板向坡板处逐渐转移。 由此可知，寒区宽浅式渠道更能够适应来自基土的冻胀变形作用，并通过释放冻胀变形以削减衬砌受力。然而，若一味将渠道断面宽浅化，会过度解除衬砌与基土间的冻结约束，导致弧底段“反拱效应”的丧失，从而产生过大的弯曲变形，降低衬砌结构的稳定性。因此，在进行弧底梯形渠道结构设计时，需充分考虑渠槽坡板与弧底板长度对结构受力及变形的综合影响，合理选择渠道断面形式。 4.1.3　不同衬砌厚度下的衬砌变形及正应力分布规律 3 以 m =1.5、 α =1.02的渠道断面为例，提取不同衬砌厚度下的法向冻胀位移及表面正应力分布，分别如 图6 图6 和 图7 图7 所示。 (a) 上表面；(b) 下表面 由 图6 图6 可知：不同衬砌厚度下的最大法向冻胀位移均位于渠底中心处，在衬砌厚度 d 由6 cm增加至20 cm的过程中，最大法向冻胀位移由 3.10 cm降至2.36 cm，减小23.9%；对于渠坡板而言，法向冻胀位移及分布基本不受到衬砌厚度变化的影响。结合截面正应力分布( 图7 图7 )可知：弧底中心处的上表面正应力由-1.11 MPa逐渐增大至2.96 MPa，超过混凝土极限抗拉强度，将会拉裂破坏。同时，距坡脚1/3的弧底板处挤压作用加剧，其下表面正应力由-0.31 MPa增加至1.54 MPa，亦会发生拉裂破坏。 由此可知，渠道衬砌厚度的增加将会强化坡板与底板间的约束，抑制衬砌板的不均匀冻胀变形。然而，过厚的衬砌板则限制了弧底板“拱效应”的发挥，导致弧底板因拉应力过大而发生拉裂破坏。故可在设计过程中适当加大衬砌厚度以抵抗来自基土的冻胀变形，但不能过度增加衬砌厚度，否则仍将导致渠道发生冻胀破坏。 4.2 　寒区弧底梯形渠道水力-抗冻胀双优标准化设计 2 分别对寒区不同温、水、土、断面规模条件下的典型渠道工程进行优化设计，得到双优结构标准尺寸参数表，并探究渠道结构尺寸在各因素下的变化规律，从而指导寒区弧底梯形渠道工程的结构抗冻胀设计。 在优化设计过程中，选取西北3处典型高寒地区，即宁夏引黄灌区、内蒙古河套灌区及新疆奎屯灌区，不同地区温度正弦函数的参数取值如 表3 表3 所示；依据断面面积分别取60、30及10 m 2 作为大、中、小型渠道面积；地下水埋深依据各灌区实际情况分别取距地表以下2.5、3.5、4.5及5.5 m四种工况；土质依据规范 [ 4 4 ] 取典型低液限黏土、粉土及含细粒砂土三种土质，冻胀率取值如式( 20 20 )~( 22 22 )所示；其余参数则参照前文及文献[ 16 16 - 19 19 16-19 ]。计算结果如 表4 表4 所示。 注： “—”表示该工况下不满足抗冻胀约束，无优化计算结果。 3 地区 3 断面面积/m 2 3 渠基土质 12 地下水埋深/m 3 2.5 3 3.5 3 4.5 3 5.5 m α d /cm m α d /cm m α d /cm m α d /cm 9 宁 夏 引 黄 灌 区 3 60 低液限黏土 — — — 2.15 1.031 12 1.76 1.024 10 1.37 1.019 8 粉土 — — — 2.08 1.025 10 1.65 1.020 8 1.29 1.016 6 含细粒砂土 2.16 1.027 10 1.71 1.023 8 1.39 1.013 6 1.20 1.010 6 3 30 低液限黏土 — — — 1.99 1.029 12 1.58 1.024 10 1.28 1.015 8 粉土 — — — 1.64 1.022 10 1.46 1.016 8 1.23 1.011 6 含细粒砂土 1.87 1.025 10 1.53 1.021 8 1.35 1.012 6 1.17 1.009 6 3 10 低液限黏土 — — — — — — 1.69 1.021 12 1.16 1.015 10 粉土 — — — 1.70 1.020 12 1.27 1.015 10 1.12 1.009 8 含细粒砂土 1.82 1.025 12 1.47 1.016 10 1.24 1.010 8 1.08 1.005 6 9 内 蒙 古 河 套 灌 区 3 60 低液限黏土 — — — — — — 2.31 1.026 12 1.45 1.019 10 粉土 — — — 2.12 1.027 12 1.76 1.023 10 1.37 1.016 6 含细粒砂土 2.39 1.033 12 1.85 1.025 10 1.53 1.017 6 1.32 1.013 6 3 30 低液限黏土 — — — — — — 2.14 1.025 14 1.36 1.017 12 粉土 — — — 1.98 1.027 12 1.66 1.019 10 1.26 1.013 8 含细粒砂土 2.23 1.029 12 1.69 1.021 10 1.42 1.018 8 1.18 1.012 6 3 10 低液限黏土 — — — — — — 1.81 1.024 16 1.23 1.017 14 粉土 — — — 1.78 1.024 14 1.37 1.018 10 1.18 1.011 8 含细粒砂土 1.53 1.019 12 1.53 1.019 12 1.27 1.012 8 1.15 1.007 6 9 新 疆 奎 屯 灌 区 3 60 低液限黏土 — — — — — — 2.54 1.029 14 1.63 1.019 10 粉土 — — — 2.38 1.031 12 1.98 1.024 10 1.48 1.017 8 含细粒砂土 2.58 1.034 14 1.99 1.026 12 1.65 1.018 10 1.36 1.016 8 3 30 低液限黏土 — — — — — — 2.37 1.026 16 1.55 1.017 12 粉土 — — — 2.25 1.027 14 1.81 1.023 10 1.34 1.013 8 含细粒砂土 2.31 1.033 14 1.87 1.025 12 1.58 1.021 8 1.27 1.013 8 3 10 低液限黏土 — — — — — — 1.88 1.025 16 1.32 1.016 14 粉土 — — — 1.94 1.027 14 1.48 1.021 12 1.23 1.015 8 含细粒砂土 — — — 1.83 1.021 14 1.45 1.019 10 1.21 1.010 8 4.2.1　地下水埋深及渠基土质影响分析 3 由 表4 表4 可以看出：对于同一环境温度、渠基土质及断面的弧底梯形渠道而言，当地下水埋深由2.5 m逐渐增大至5.5 m时，双优渠道的边坡系数 m 、实佳比 α 及衬砌厚度 d 均呈现出减小的趋势。这是由于当地下水埋深较浅时，基土在冻胀变形过程中水分补给量较大，导致渠道整体冻胀加剧，增大渠道边坡系数 m 及实佳比 α 能够使得渠道断面转变为宽浅式，整体冻胀变形及受力分布更加均匀。同时，合理增大衬砌厚度能够增强弧底段衬砌的抗冻胀能力，在二者综合作用下使得渠道整体抗冻胀效果最优。反之，当地下水埋深较深时，渠道整体冻胀程度较轻，设置较小尺寸参数的窄深式衬砌结构即可满足抗冻胀目标。 当地下水埋深过浅时，水分在基土冻胀过程中能够补给至渠道全断面，即形成开放系统 [ 30 30 ] ，渠道冻胀破坏严重，仅依靠改变渠道形体结构已无法满足抗冻胀约束条件。此时应增加额外的保温、排水、分缝及布设双膜垫层等工程措施 [ 1 1 ] ，削弱基土冻胀强度，再进行结构优化设计。由此可知，在进行寒区渠道设计时，应充分考虑渠基地下水埋深的影响，合理设置渠道结构尺寸参数。同时，本文在计算过程中因缺乏土体热参数实测值，故采用混合律公式对渠基土热物理参数进行简化计算。这将导致渠道冻胀变形及受力结果在数值上较实际情况偏大 [ 30 30 - 32 32 30-32 ] 。因而在优化计算过程中，优化结果是趋于安全的，能够应用于工程实际。 由此可知，在相同断面规模下，寒区宽浅式弧底梯形渠道的抗冻胀能力更强的原因如下：一是其自身能够较好地适应来自基土的冻胀变形作用，通过释放衬砌整体变形以降低其受力；二是相较于窄深式断面，其距地下水位更远，衬砌下部基土冻胀强度降低，使得衬砌整体冻胀破坏程度减轻。 另外，由于在设计过程中难以通过现场试验准确获取各地区不同渠基土质在冻结的过程中的冻胀率，故依据规范 [ 4 4 ] 对不同典型土质的冻胀率进行拟合(见式( 20 20 )~( 22 22 ))。由拟合公式可知，当地下水埋深及冻深条件相同时，冻胀率从大到小顺序为低液限黏土、粉土、含细粒砂土。这是由于冻胀率与土体中活性矿物含量及冻前初始含水等因素密切相关，例如，当低液限黏土中活性矿物成分的含量过高时，土体分散度增大，比表面积增大，在冻结过程中的水分迁移和冰分凝作用增强，导致其冻胀率将大于粉土冻胀率 [ 33 33 ] 。由 表4 表4 可知：对于同一环境温度、断面面积及地下水埋深下的弧底梯形渠道而言，不同土质渠道相应的结构尺寸参数按从大到小顺序为低液限黏土、粉土、含细粒砂土。因此，在寒区弧底梯形渠道设计中，应着重考虑渠基土质的影响，渠基土在冻结过程中的冻胀率越大，则结构设计尺寸参数也应越大。 4.2.2　外界环境温度及断面规模影响分析 3 由 表4 表4 还可以看出：对于同一断面规模、渠基土质、地下水埋深的弧底梯形渠道而言，不同地区对应的双优结构参数按从大到小顺序为新疆奎屯灌区、内蒙古河套灌区、宁夏引黄灌区。这表明冬季外界环境温度越低，基土冻深越大，渠道冻胀破坏越剧烈，应设置更为宽浅式的断面形式以适应冻胀变形，同时增大衬砌厚度以抵抗来自基土的冻胀作用，使得渠道整体柔度达到最大。 由 表4 表4 还可知：当环境温度、渠基土质、地下水埋深一致时，双优结构参数随断面规模的增大而逐渐增大。这表明大型渠道对基土的冻胀敏感性较强，更易发生冻胀破坏。因此，在寒区大型弧底梯形渠道设计中，宜采用宽浅式渠道断面，并合理增大衬砌厚度，以确保渠道安全稳定运行。本文仅给出断面面积为60、30及10 m 2 的渠道结构优化结果，若需获取其他断面面积下的最优结构参数，则根据 表4 表4 中数值进行线性插值即可。 4.2.3　工程案例计算 3 利用上述双优设计方法对新疆某弧底梯形渠道进行优化设计，提取优化前后的衬砌法向冻胀位移、截面正应力及结构尺寸参数，分别如 图8 图8 和 图9 图9 及 表5 表5 所示。 (a) 上表面；(b) 下表面 断面 边坡系数 m 实佳比 α 衬砌厚度 d /cm 弧底板长度/m 最大拉应力/MPa 最大位移/cm 湿周 χ /m 整体柔度指标值/ (cm∙MPa -1 ) 原设计断面 1.5 1.002 10 3.52 2.40 2.10 18.59 214.69 优化后断面 2.2 1.020 14 5.02 0.95 2.52 19.31 305.47 由 图8 图8 和 图9 图9 可知：该渠道在优化前的最大法向冻胀位移为2.10 cm，位于弧底中心处；弧底段衬砌上表面受拉、下表面受压，应力极值分别为2.40 MPa和-2.76 MPa；靠近坡脚的渠坡下半段因受到挤压作用而表现为上表面受压、下表面受拉，应力极值分别为1.25 MPa和-6.05 MPa。因此，弧底中心上表面与坡脚附近下表面将会发生拉裂破坏，计算结果与现场实际破坏情况基本一致 [ 22 22 ] 。 由 表5 表5 可知：对渠道进行双优设计后，坡度系数 m 、实佳比 α 、衬砌厚度 d 及断面湿周 χ 分别增大了46.7%、1.8%、40.0%及3.9%，输水性能仍保持在水力最优范围内。同时，优化后断面的最大法向冻胀位移仍位于渠底中心处，较优化前增大了20.0%，符合设计要求；衬砌上、下表面拉应力极大值分别较优化前降低了60.4%和52.0%，混凝土材料强度得到充分利用；衬砌整体柔度较优化前增大了42.3%，表明其变形与受力分布更加均匀，适应冻胀变形的能力得到显著提高。 5 结论 1 1) 衬砌下表面拉应力极大值点的位置与坡板及弧底板长度有关。当边坡系数 m 保持不变，下表面拉应力极大值点随着实佳比 α 的增大由坡板向弧底中心处转移；当实佳比 α 保持不变，下表面拉应力极大值点随着边坡系数 m 的增大由弧底板向坡板处逐渐转移。 2) 双优结构尺寸参数均随地下水埋深的增大而减小、随渠基土质冻胀率的增大而增大、随冬季负温的降低而增大、随断面规模的增大而增大。 3) 寒区高地下水位地区优选弧底梯形宽浅式渠道断面不仅可使地下水位距离渠底较远，基土冻胀强度降低，而且衬砌释放应力和适应冻胀变形的能力强，可使衬砌整体冻胀破坏程度减轻。 4) 新疆的工程算例结果表明：满足水力及抗冻胀双优的渠道与原设计值相比，衬砌上、下表面拉应力极大值分别降低了60.4%和52.0%，最大法向冻胀位移增大了20.0%，符合设计要求；衬砌整体柔度增大了42.3%，适应冻胀变形的能力得到显著增强。 参考文献 王正中 , 江浩源 , 王羿 , 等 . 旱寒区输水渠道防渗抗冻胀研究进展与前沿 [J]. 农业工程学报 , 2020 , 36 ( 22 ): 120 - 132 . 旱寒区输水渠道防渗抗冻胀研究进展与前沿 WANG Zhengzhong , JIANG Haoyuan , WANG Yi , et al . Research progresses and frontiers in anti-seepage and anti-frost heave of canals in cold-arid regions [J]. Transactions of the Chinese Society of Agricultural Engineering , 2020 , 36 ( 22 ): 120 - 132 . Research progresses and frontiers in anti-seepage and anti-frost heave of canals in cold-arid regions 王羿 , 王正中 , 刘铨鸿 , 等 . 基于弹性薄层接触模型研究衬砌渠道双膜防冻胀布设 [J]. 农业工程学报 , 2019 , 35 ( 12 ): 133 - 141 . 基于弹性薄层接触模型研究衬砌渠道双膜防冻胀布设 WANG Yi , WANG Zhengzhong , LIU Quanhong , et al . Double membranes cushion layout for preventing frost heave damage based on elastic thin contact model [J]. Transactions of the Chinese Society of Agricultural Engineering , 2019 , 35 ( 12 ): 133 - 141 . Double membranes cushion layout for preventing frost heave damage based on elastic thin contact model GB/T 50600—2010 . 渠道防渗工程技术规范 [S]. GB/T 50600—2010 . Technical specifications for canal seepage prevention engineering [S]. SL 23—2006 . 渠系工程抗冻胀设计规范 [S]. SL 23—2006 . Code for frost heave resistance design of canal works [S]. 程敏 . 泾惠渠灌区现代化改造和建设进程的分析与思考 [J]. 陕西水利 , 2020 ( 11 ): 139 - 140 . 泾惠渠灌区现代化改造和建设进程的分析与思考 CHENG Min . Analysis and thinking on modernization and construction process of Jinghuiqu irrigation area [J]. Shaanxi Water Resources , 2020 ( 11 ): 139 - 140 . Analysis and thinking on modernization and construction process of Jinghuiqu irrigation area 陈向阳 . 大禹渡弧底梯形渠道裂缝产生原因及处理措施 [J]. 山西水利 , 2004 , 20 ( 3 ): 64 - 66 . 大禹渡弧底梯形渠道裂缝产生原因及处理措施 CHEN Xiangyang . Causes and treatment measures of cracks in trapezoidal channel at the bottom of Dayudu arc [J]. Shanxi Water Resources , 2004 , 20 ( 3 ): 64 - 66 . Causes and treatment measures of cracks in trapezoidal channel at the bottom of Dayudu arc BAI Ruiqiang , LAI Yuanming , PEI Wansheng , et al . Investigation on frost heave of saturated-unsaturated soils [J]. Acta Geotechnica , 2020 , 15 ( 11 ): 3295 - 3306 . Investigation on frost heave of saturated-unsaturated soils MONADJEMI P . General formulation of best hydraulic channel section [J]. Journal of Irrigation and Drainage Engineering , 1994 , 120 ( 1 ): 27 - 35 . General formulation of best hydraulic channel section 余长洪 , 周明耀 , 姜健俊 , 等 . 灌区节水改造中防渗渠道断面的优化设计 [J]. 农业工程学报 , 2004 , 20 ( 1 ): 91 - 94 . 灌区节水改造中防渗渠道断面的优化设计 YU Changhong , ZHOU Mingyao , JIANG Jianjun , et al . Optimal design of anti-seepage channel section for water-saving transformation in irrigation district [J]. Transactions of the Chinese Society of Agricultural Engineering , 2004 , 20 ( 1 ): 91 - 94 . Optimal design of anti-seepage channel section for water-saving transformation in irrigation district BHATTACHARJYA R K , SATISH M G . Optimal design of a stable trapezoidal channel section using hybrid optimization techniques [J]. Journal of Irrigation and Drainage Engineering , 2007 , 133 ( 4 ): 323 - 329 . Optimal design of a stable trapezoidal channel section using hybrid optimization techniques 辛英华 , 王正中 . U形衬砌渠道结构及水力最佳断面的分析 [J]. 节水灌溉 , 2008 ( 2 ): 36 - 38 . U形衬砌渠道结构及水力最佳断面的分析 XIN Yinghua , WANG Zhengzhong . The structural and hydraulic analysis on frost-heaving concrete lining of the U-shape canal [J]. Water Saving Irrigation , 2008 ( 2 ): 36 - 38 . The structural and hydraulic analysis on frost-heaving concrete lining of the U-shape canal 刘东 , 胡宇祥 , 付强 , 等 . 北方灌区混凝土衬砌渠道断面优化及参数分析 [J]. 农业工程学报 , 2015 , 31 ( 20 ): 107 - 114 . 北方灌区混凝土衬砌渠道断面优化及参数分析 LIU Dong , HU Yuxiang , FU Qiang , et al . Optimization and parameter analysis for channel cross section with concrete lining in northern irrigation district [J]. Transactions of the Chinese Society of Agricultural Engineering , 2015 , 31 ( 20 ): 107 - 114 . Optimization and parameter analysis for channel cross section with concrete lining in northern irrigation district 安鹏 , 邢义川 , 张爱军 , 等 . 弧底梯形渠道抗冻胀结构优化与数值模拟 [J]. 灌溉排水学报 , 2017 , 36 ( 11 ): 56 - 62 . 弧底梯形渠道抗冻胀结构优化与数值模拟 AN Peng , XING Yichuan , ZHANG Aijun , et al . Optimal analysis of anti-frost heaping canal with cross section of trapezoidal slopes and curved bed [J]. Journal of Irrigation and Drainage , 2017 , 36 ( 11 ): 56 - 62 . Optimal analysis of anti-frost heaping canal with cross section of trapezoidal slopes and curved bed 王玉宝 , 王亮 , 程森浩 , 等 . 整体式U形渠道不同断面结构抗冻胀性研究 [J]. 东北大学学报(自然科学版) , 2020 , 41 ( 9 ): 1341 - 1347 . 整体式U形渠道不同断面结构抗冻胀性研究 WANG Yubao , WANG Liang , CHENG Senhao , et al . Study on frost heave resistance of different section structures of U-shaped canal [J]. Journal of Northeastern University(Natural Science) , 2020 , 41 ( 9 ): 1341 - 1347 . Study on frost heave resistance of different section structures of U-shaped canal 王正中 . 梯形渠道砼衬砌冻胀破坏的力学模型研究 [J]. 农业工程学报 , 2004 , 20 ( 3 ): 24 - 29 . 梯形渠道砼衬砌冻胀破坏的力学模型研究 WANG Zhengzhong . Establishment and application of mechanics models of frost heaving damage of concrete lining trapezoidal open canal [J]. Transactions of the Chinese Society of Agricultural Engineering , 2004 , 20 ( 3 ): 24 - 29 . Establishment and application of mechanics models of frost heaving damage of concrete lining trapezoidal open canal LI Zhuo , LIU Sihong , FENG Youting , et al . Numerical study on the effect of frost heave prevention with different canal lining structures in seasonally frozen ground regions [J]. Cold Regions Science and Technology , 2013 , 85 : 242 - 249 . Numerical study on the effect of frost heave prevention with different canal lining structures in seasonally frozen ground regions LI Shuangyang , LAI Yuanming , ZHANG Mingyi , et al . Centrifuge and numerical modeling of the frost heave mechanism of a cold-region canal [J]. Acta Geotechnica , 2019 , 14 ( 4 ): 1113 - 1128 . Centrifuge and numerical modeling of the frost heave mechanism of a cold-region canal LIU Zhen , YU Xiong . Coupled thermo-hydro-mechanical model for porous materials under frost action: theory and implementation [J]. Acta Geotechnica , 2011 , 6 ( 2 ): 51 - 65 . Coupled thermo-hydro-mechanical model for porous materials under frost action: theory and implementation 王正中 , 刘少军 , 王羿 , 等 . 寒区弧底梯形衬砌渠道冻胀破坏的尺寸效应研究 [J]. 水利学报 , 2018 , 49 ( 7 ): 803 - 813 . 寒区弧底梯形衬砌渠道冻胀破坏的尺寸效应研究 WANG Zhengzhong , LIU Shaojun , WANG Yi , et al . Size effect on frost heave damage for lining trapezoidal canal with arc-bottom in cold regions [J]. Journal of Hydraulic Engineering , 2018 , 49 ( 7 ): 803 - 813 . Size effect on frost heave damage for lining trapezoidal canal with arc-bottom in cold regions TABER S . The mechanics of frost heaving [J]. The Journal of Geology , 1930 , 38 ( 4 ): 303 - 317 . The mechanics of frost heaving 王正中 , 袁驷 , 陈涛 . 冻土横观各向同性非线性本构模型的实验研究 [J]. 岩土工程学报 , 2007 , 29 ( 8 ): 1215 - 1218 . 冻土横观各向同性非线性本构模型的实验研究 WANG Zhengzhong , YUAN Si , CHEN Tao . Study on the constitutive model of transversely isotropic frozen soil [J]. Chinese Journal of Geotechnical Engineering , 2007 , 29 ( 8 ): 1215 - 1218 . Study on the constitutive model of transversely isotropic frozen soil 江浩源 , 王正中 , 王羿 , 等 . 大型弧底梯形渠道“适缝”防冻胀机理及应用研究 [J]. 水利学报 , 2019 , 50 ( 8 ): 947 - 959 . 大型弧底梯形渠道“适缝”防冻胀机理及应用研究 JIANG Haoyuan , WANG Zhengzhong , WANG Yi , et al . Study on anti-frost heave mechanism and application of “reasonable joints” for large-size trapezoidal canal with arc-bottom [J]. Journal of Hydraulic Engineering , 2019 , 50 ( 8 ): 947 - 959 . Study on anti-frost heave mechanism and application of “reasonable joints” for large-size trapezoidal canal with arc-bottom LIU Quanhong , WANG Zhengzhong , LI Zhanchao , et al . Transversely isotropic frost heave modeling with heat-moisture-deformation coupling [J]. Acta Geotechnica , 2020 , 15 ( 5 ): 1273 - 1287 . Transversely isotropic frost heave modeling with heat-moisture-deformation coupling SARUYA T , REMPEL A W , KURITA K . Hydrodynamic transitions with changing particle size that control ice lens growth [J]. The Journal of Physical Chemistry B , 2014 , 118 ( 47 ): 13420 - 13426 . Hydrodynamic transitions with changing particle size that control ice lens growth 王羿 , 刘瑾程 , 刘铨鸿 , 等 . 温-水-土-结构耦合作用下寒区梯形衬砌渠道结构形体优化 [J]. 清华大学学报(自然科学版) , 2019 , 59 ( 8 ): 645 - 654 . 温-水-土-结构耦合作用下寒区梯形衬砌渠道结构形体优化 WANG Yi , LIU Jincheng , LIU Quanhong , et al . Shape optimization of a trapezoidal canal structure for coupled temperature-water-soil conditions in cold regions [J]. Journal of Tsinghua University(Science and Technology) , 2019 , 59 ( 8 ): 645 - 654 . Shape optimization of a trapezoidal canal structure for coupled temperature-water-soil conditions in cold regions LIU Jiankun , WANG Tengfei , TAI Bowen , et al . A method for frost jacking prediction of single pile in permafrost [J]. Acta Geotechnica , 2020 , 15 ( 2 ): 455 - 470 . A method for frost jacking prediction of single pile in permafrost ZHOU Jiazuo , LI Dongqing . Numerical analysis of coupled water, heat and stress in saturated freezing soil [J]. Cold Regions Science and Technology , 2012 , 72 : 43 - 49 . Numerical analysis of coupled water, heat and stress in saturated freezing soil 王正中 , 沙际德 , 蒋允静 , 等 . 正交各向异性冻土与建筑物相互作用的非线性有限元分析 [J]. 土木工程学报 , 1999 , 32 ( 3 ): 55 - 60 . 正交各向异性冻土与建筑物相互作用的非线性有限元分析 WANG Zhengzhong , SHA Jide , JIANG Yunjing , et al . Nonlinear finite element analysis of interaction of orthotropic frozen ground and construction [J]. China Civil Engineering Journal , 1999 , 32 ( 3 ): 55 - 60 . Nonlinear finite element analysis of interaction of orthotropic frozen ground and construction ZHANG Mingyi , PEI Wansheng , LI Shuangyang , et al . Experimental and numerical analyses of the thermo-mechanical stability of an embankment with shady and sunny slopes in a permafrost region [J]. Applied Thermal Engineering , 2017 , 127 : 1478 - 1487 . Experimental and numerical analyses of the thermo-mechanical stability of an embankment with shady and sunny slopes in a permafrost region 石梁宏 , 李双洋 , 王冲 , 等 . 冻土热物理性质的统计特征及分布规律 [J]. 中南大学学报(自然科学版) , 2018 , 49 ( 12 ): 3060 - 3067 . 冻土热物理性质的统计特征及分布规律 SHI Lianghong , LI Shuangyang , WANG Chong , et al . Statistical characters and distribution laws of thermophysical properties of frozen soil [J]. Journal of Central South University(Science and Technology) , 2018 , 49 ( 12 ): 3060 - 3067 . Statistical characters and distribution laws of thermophysical properties of frozen soil FANG Hongyuan , ZHAO Peng , ZHANG Chao , et al . A cleaner polyurethane elastomer grouting material with high hardening strain for the fundamental rehabilitation: the comprehensive mechanical properties study [J]. Construction and Building Materials , 2022 , 318 : 125951 . A cleaner polyurethane elastomer grouting material with high hardening strain for the fundamental rehabilitation: the comprehensive mechanical properties study 刘为民 , 何平 , 张钊 . 土体导热系数的评价与计算 [J]. 冰川冻土 , 2002 , 24 ( 6 ): 770 - 773 . 土体导热系数的评价与计算 LIU Weimin , HE Ping , ZHANG Zhao . A calculation method of thermal conductivity of soils [J]. Journal of Glaciolgy and Geocryology , 2002 , 24 ( 6 ): 770 - 773 . A calculation method of thermal conductivity of soils 徐学祖 . 冻土物理学 [M]. 北京 : 科学出版社 , 2001 : 125 - 136 . 冻土物理学 XU Xuezu . Frozen soil physics [M]. Beijing : Science Press , 2001 : 125 - 136 . Frozen soil physics 10.11817/j.issn.1672-7207.2023.05.025 TV68 S277 A 1672-7207（2023）05-1916-14 王正中, 孙涛, 郑艾磊, 等. 水-热-力耦合作用下寒区弧底梯形渠道结构优化设计[J]. 中南大学学报(自然科学版), 2023, 54(5): 1916-1929. 王正中, 孙涛, 郑艾磊, 等. 水-热-力耦合作用下寒区弧底梯形渠道结构优化设计[J]. 中南大学学报(自然科学版), 2023, 54(5): 1916-1929. WANG Zhengzhong, SUN Tao, ZHENG Ailei, et al. Structural optimization design for trapezoidal canal with arc-bottom during thermal-moisture-mechanical coupling in cold regions[J]. Journal of Central South University(Science and Technology), 2023, 54(5): 1916-1929. WANG Zhengzhong, SUN Tao, ZHENG Ailei, et al. Structural optimization design for trapezoidal canal with arc-bottom during thermal-moisture-mechanical coupling in cold regions[J]. Journal of Central South University(Science and Technology), 2023, 54(5): 1916-1929. 国家重点研发计划项目(2017YFC0405102)；国家自然科学基金资助项目(U2003108，51279168)；冻土工程国家重点实验室开放基金资助项目(SKLFSE201801) Project(2017YFC0405102) supported by the National Key Research and Development Program of China; Projects(U2003108, 51279168) supported by the National Natural Science Foundation of China; Project(SKLFSE201801) supported by the Open Research Fund of State Key Laboratory of Frozen Soil Engineering 2022-08-14 2022-10-24 2023-05-26 2023-06-20 38910721 中南大学学报(自然科学版) 5 54