西北农林科技大学 旱区寒区水工程安全研究中心, 旱区农业水土工程教育部重点实验室, 杨凌 712100
收稿日期：2018-10-31
基金项目：国家重点研发计划“水资源高效开发利用”重点专项（2017YFC0405100）；国家自然科学基金项目（51279168）；国家十二五科技支撑计划（2012BAD10B02）；中央高校基本科研业务科技创新专项（Z102021848，Z109021807）
作者简介：王羿(1987-), 男, 博士研究生
通信作者：王正中, 教授, E-mail:wangzz0910@163.com
摘要：为提高寒区大型梯形渠道结构输水及抗冻胀破坏能力，采用分层序列法构建渠道断面水力与抗冻胀双目标优化方法。为准确反映温、水、土等外环境与衬砌结构耦合作用所导致的渠道冻胀变形过程，基于冻土水-热-力耦合冻胀理论建立了衬砌渠道冻胀分析模型。以水力最优为第1层次优化，以所得形体几何参数作为抗冻胀优化的解集空间；在第2层次的抗冻胀优化过程中提出了衬砌结构整体刚度指标，并以整体刚度指标最小作为优化目标，以衬砌允许最大法向位移和拉应力作为约束条件；通过COMSOL软件二次开发对冻胀模型控制方程进行有限元求解。工程算例表明：优化结果体现了对渠道行水和衬砌抗冻胀的改进；与原设计相比，优化后渠道结构整体刚度系数减小30%~48%，提高了对冻胀变形的适应能力。该方法可为类似工程设计提供参考。
关键词：衬砌渠道温-水-土-结构耦合形体优化分层序列法水热力耦合理论
Shape optimization of a trapezoidal canal structure for coupled temperature-water-soil conditions in cold regions
WANG Yi, LIU Jincheng, LIU Quanhong, WANG Zhengzhong
Cold and Arid Regions Water Engineering Safety Research Center, Key Laboratory of Agricultural Soil and Water Engineering in Arid and Semiarid Areas of Ministry of Education, Northwest A & F University, Yangling 712100, China
Abstract: The shape of a large trapezoidal canal in a cold region should be optimized to improve the water delivery and minimize frost damage. This study used a double objective analysis to optimize the hydraulics and the frost damage prevention for a canal using the hierarchical sequence method. A frost heave model was developed for a concrete lining canal based on thermo-hydro-mechanics frost heave theory for freezing soil to predict the canal deformation for various temperature, water, and soil conditions. The optimal hydraulic design was predicted first with the canal section parameters then used as the solution set space for the frost damage prevention optimization. The overall stiffness index of the lining was set as the optimization object with the allowable maximum normal displacement and stress in the lining set as the optimization constraints. The model was solved using the secondary hierarchical sequence method in COMSOL. The results show that this method can consider both the water delivery and the frost damage effect on the canal section. The practical case results indicate that the overall stiffness index of the lining is reduced by 30%-48% from that of the original design using this method to get an optimized canal section, which improves the adaptability to frost damage of the lining structure. This method can also be used for similar analyses involving hydraulics and frost damage effects.
Key words: lined canalcoupling of temperature-water-soil and structureshape optimizationhierarchical sequence methodthermo-hydro-mechanics theory
随着西部大开发战略“十三五”规划的制定与实施，为了优化配置西北旱寒区水资源、提高水安全保障能力及加快灌区节水改造，一批重大长距离引调水工程、水源工程及大型灌区工程将开工建设。渠道作为主要输水结构，相较于管道与隧洞，具有施工造价低、日常管理与维护方便、输水流量大、水头损失小等优点，在条件允许情况下成为首选形式。但是，大型长距离宽浅式输水渠道虽然抗冻胀能力强，但开挖面积大，占地多，输水阻力大；而窄深式渠道虽然接近水力最优断面，但因不便施工、渠坡结构稳定性及抗冻胀能力较差，直接影响工程的安全运行。《渠系工程抗冻胀规范》[1]中没有综合考虑以上多方面因素的设计方法，工程上多凭经验来协同各影响因素。因此，有必要充分考虑渠道开挖面积、输水条件及衬砌抗冻胀破坏等因素，提出开挖量少、输水效率高、更适应冻胀变形的渠道断面形式，为旱寒区输水工程设计提供理论基础。
当前衬砌渠道断面优化多以水力优化为主，兼顾开挖面积和工程造价等经济性因素[2-5]，而对断面抗冻胀设计多根据《渠系工程抗冻胀规范》[1]中规定的以冻深、土质和地下水位经验性确定的冻胀力。部分研究者基于冻胀力假设建立衬砌结构工程力学模型和有限元模型，结合水力计算对渠道断面优化进行初步探索。辛英华等[6]通过对“U”形渠道衬砌的破坏特征及机理的分析，采用王正中等[7-8]的工程力学方法分析衬砌冻胀力，在水力最佳断面基础上对衬砌结构进行优化，使其冻胀力和结构应力分布较为均匀。刘东等[9]将王正中等[7-8]提出的冻胀模型仅用于确定衬砌板厚度参数，随后采用粒子群算法对考虑水力条件、开挖和衬砌量等因素的非线性优化模型进行求解，得到经济性最佳的衬砌板厚度及混凝土衬砌梯形渠道断面参数，但未考虑渠道断面的抗冻胀特性。王正中等[10]、刘旭东等[11]、李爽等[12]、郝晋彩等[13]等采用有限元软件计算渠道温度场，将冻土视为负温膨胀材料，并采用《渠系工程抗冻胀规范》[1]中的冻胀率作为材料的负温膨胀系数，建立渠基冻胀计算有限元模型。但是，这些研究通过有限元计算的衬砌应力场结果对渠道断面进行优化设计过程中没有考虑渠道水力特性，不符合实际工程设计需要。安鹏等[14]基于Stephen公式计算冻深，并建立冻深与冻胀量的经验公式，针对渠基封闭、半封闭和开放系统分别进行冻胀模拟。以冻胀位移和结构应力为约束，以渠道经济规模为优化目标，对弧底梯形渠道进行了优化设计，但所得结果并非抗冻胀最优的渠道断面。以上研究或受限于冻胀量取值的经验性，没有考虑实际渠道所处气候、水文地质条件，其优化结果难以与实际对接；或受限于优化目标单一，未能得到同时满足水力特性和抗冻胀特性最优要求。渠道冻胀量的发生发展受多种因素影响。冻胀本质在于水分迁移聚集于冻结缘所形成的分凝冰层的增长，分凝冰层增长量受制于土体含水量、渗透性、孔隙率和导热性等因素，水分迁移驱动力主要来源于冻结缘处温度梯度，而迁移水结冰所带来的潜热又反作用于温度场。此外，对于大型渠道，断面形式影响其表面热对流特性，导致渠基不同部位降温速率不同，最终造成冻胀量发展的不均匀性。因此，应考虑渠道所处环境下气温、水分、土质等因素，研究这些因素与渠道断面的耦合作用，准确计算冻胀量；再从水力学要素和冻胀量两方面逆推，对渠道断面尺寸进行优化，提高其输水性能和运行寿命。
本文基于冻土冻胀理论建立描述渠基土不均匀冻胀的水-热-力耦合数学模型，充分考虑气温、热流边界、渠基水分场、地下含水边界和土体渗流特性与渠槽形状之间的耦合作用。采用COMSOL Multiphysics软件自定义偏微分方程求解功能进行二次开发，对数学模型进行有限元求解。采用层次序列优化法，以渠道水力条件为第1层次优化，以所得结果作为解集空间。通过衬砌法向冻胀位移和应力构建衬砌整体刚度指标，表征结构对冻胀变形的适应能力，并以此指标为第2层次优化目标，采用COMSOL软件参数化求解方法进行优化求解，从而获得满足“水力+抗冻胀”双优的渠道断面尺寸参数。
1 基于冻土水-热-力耦合的衬砌渠道冻胀模型准确计算渠基土体冻胀量及衬砌结构的冻胀力学响应，是进行渠道抗冻胀断面优化的重要前提。为此，需要从冻土冻胀的基本理论出发，建立考虑温度、土质、含水量和冻结过程等因素的土体冻胀模型以及合理反映冻土与衬砌结构相互作用的接触模型，两者共同组成渠道衬砌冻胀破坏模型。渠道是典型的平面应变问题，可以用x-y平面上的Descartes坐标来描述。
1.1 渠基土温度场与水分场耦合控制方程仅考虑冻土中水分迁移和冰水相变过程，则渠基内热量输送方程为
| $\begin{array}{*{20}{c}}{\rho {C_p}\frac{{\partial T}}{{\partial t}} = \frac{\partial }{{\partial x}}\left( { - {\lambda _{{\rm{eq}}}}\frac{{\partial T}}{{\partial x}}} \right) + }\\{\frac{\partial }{{\partial y}}\left( { - {\lambda _{{\rm{eq}}}}\frac{{\partial T}}{{\partial y}}} \right) + {L_{\rm{f}}}{\rho _{\rm{i}}}\frac{{{\rm{d}}{\theta _{\rm{i}}}}}{{{\rm{d}}t}}.}\end{array}$ | (1) |
忽略空气相，仅考虑冻土与融土中土颗粒、水相与冰相所占体积及原位水冰相变对土体热参数的影响，式(1)中导热系数与等效定压热容可以表示为[15]：
| ${\lambda _{{\rm{eq}}}} = {\theta _{\rm{s}}}{\lambda _{\rm{s}}} + {\theta _{{\rm{i0}}}}{\lambda _{\rm{i}}} + {\theta _{\rm{w}}}{\lambda _{\rm{w}}}.$ | (2) |
| ${C_p} = \frac{1}{\rho }\left( {{\theta _{\rm{s}}}{C_{p,{\rm{s}}}} + {\theta _{{\rm{i0}}}}{C_{p,{\rm{i}}}} + {\theta _{\rm{w}}}{C_{p,{\rm{w}}}}} \right) + {L_{\rm{f}}}\frac{{{\rm{d}}{\theta _{{\rm{i0}}}}}}{{{\rm{d}}T}}.$ | (3) |
| ${\theta _{\rm{w}}} + {\theta _{{\rm{i0}}}} = 1 - {\theta _{\rm{s}}} = {\theta _{\rm{v}}}.$ | (4) |
式(1)仅考虑正冻土迁移水的相变潜热对温度场的影响，将其作为源项加入到方程(1)右项，而忽略该部分水分所带来的显热。根据文[17]，水分迁移驱动力由Clepyron方程推导为温度梯度函数，水分由未冻区迁移至冻结缘区后开始冻结并形成冰透镜体；同时，水分迁移满足Darcy定律，可以得到在冻结缘区水分迁移相变所产生的热量。
| $\begin{array}{*{20}{c}}{{L_{\rm{f}}}{\rho _{\rm{i}}}\frac{{{\rm{d}}{\theta _{\rm{i}}}}}{{{\rm{d}}t}} = \frac{{{\rho _{\rm{w}}}L_{\rm{f}}^2}}{{g{T_{\rm{f}}}}}\nabla \cdot \left( { - k\nabla T} \right)}\\{\left( { - 0.5{\;^ \circ }{\rm{C}} \le T \le 0{\;^ \circ }{\rm{C}}} \right).}\end{array}$ | (5) |
1.2 冻土冻胀本构关系渠基土因原位水及迁移水冻结所产生的冻胀变形可以用以下方程描述：
| $\left\{ \begin{array}{l}\frac{{\partial {\sigma _x}}}{{\partial x}} + \frac{{\partial {\tau _{xy}}}}{{\partial y}} = {f_x},\\\frac{{\partial {\sigma _y}}}{{\partial y}} + \frac{{\partial {\tau _{xy}}}}{{\partial x}} = {f_y}.\end{array} \right.$ | (6) |
| $\left\{ \begin{array}{l}{\sigma _x} = E\left( {{\varepsilon _x} - {\varepsilon _{x0}}} \right) - \nu \left( {{\varepsilon _y} - {\varepsilon _{y0}}} \right),\\{\sigma _y} = E\left( {{\varepsilon _y} - {\varepsilon _{y0}}} \right) - \nu \left( {{\varepsilon _x} - {\varepsilon _{x0}}} \right),\\{\tau _{xy}} = \frac{{E\left( {{\varepsilon _{xy}} - {\varepsilon _{xy0}}} \right)}}{{1 + \nu }}.\end{array} \right.$ | (7) |
| $\left\{ \begin{array}{l}{\varepsilon _{x0}} = {\varepsilon _0}/2 + {\varepsilon _{{\rm{grand}}T}} \cdot {l^2},\\{\varepsilon _{y0}} = {\varepsilon _0}/2 + {\varepsilon _{{\rm{grand}}T}} \cdot {n^2},\\{\varepsilon _{xy0}} = {\varepsilon _{yx0}} = {\varepsilon _{{\rm{grand}}T}} \cdot l \cdot n,\end{array} \right.T < 0{\;^{\rm{o}}}{\rm{C}}{\rm{.}}$ | (8) |
| ${\varepsilon _0} = {\theta _{\rm{w}}} \cdot 0.09.$ | (9) |
| $\begin{array}{*{20}{c}}{{\varepsilon _{{\rm{grand}}T}} = \int {\frac{{\partial {\theta _i}}}{{\partial t}}{\rm{d}}t} = }\\{\int {\left( {\frac{{{\rho _{\rm{w}}}{L_{\rm{f}}}}}{{{\rho _{\rm{i}}}g{T_{\rm{f}}}}}\nabla \cdot \left( { - k\nabla T} \right)} \right){\rm{d}}t} }\\{\left( { - 0.5\;{}^ \circ {\rm{C}} \le T \le 0\;{}^ \circ {\rm{C}}} \right).}\end{array}$ | (10) |
| $\begin{array}{*{20}{c}}{\frac{{{\rm{d}}{\varepsilon _{{\rm{grant}}T}}}}{{{\rm{d}}t}} = \frac{{{\rho _{\rm{w}}}{L_{\rm{f}}}}}{{{\rho _{\rm{i}}}g{L_{\rm{f}}}}}\nabla \cdot \left( { - k\nabla T} \right)}\\{\left( { - 0.5{\;^{\rm{o}}}{\rm{C}} \le T \le 0{\;^{\rm{o}}}{\rm{C}}} \right).}\end{array}$ | (11) |
1.3 COMSOL Multiphysics有限元求解方案采用多物理场分析软件COMSOL Multiphysics求解方程组，对已知环境温度、含水量和土质条件下渠道的冻胀过程进行数值计算。COMSOL Multiphysics软件是源于MATLAB微分方程求解模块开发的多物理场耦合有限元求解软件，可以将本文推导的控制微分方程组自定义编写入软件物理场模块，同时以函数和表达式的形式定义模型水、热参数和本构方程。最后，通过软件中参数化扫描求解功能，对解集空间内的渠道形式进行冻胀分析，得到各几何参数组合下渠道衬砌结构的冻胀位移和结构应力，最终得到满足优化目标和约束条件的几何参数组合。
2 水力与抗冻胀双目标断面优化数学模型2.1 水力最优设计参数及其解集空间采用分层序列法[18-19]对双目标优化问题进行求解，将水力最优目标作为第1层序列，在其解集空间内求解抗冻胀最优目标解集。梯形渠道最佳水力断面及满足工程需要的实用经济断面各尺寸参数满足以下关系[20]：
| ${\beta _{\rm{m}}} = \frac{{{b_{\rm{m}}}}}{{{h_{\rm{m}}}}} = 2\left[ {\sqrt {1 + {m^2}} - m} \right],$ | (12) |
| $\alpha = \frac{A}{{{A_{\rm{m}}}}} = \frac{{\left( {\beta + m} \right){h^2}}}{{\left( {{\beta _m} + m} \right)h_m^2}},$ | (13) |
| $\frac{h}{{{h_m}}} = {\alpha ^{5/2}}\left[ {1 - \sqrt {1 - {\alpha ^{ - 4}}} } \right],$ | (14) |
| $\beta = {\left( {\frac{{{h_m}}}{h}} \right)^2}\alpha \left( {2\sqrt {1 + {m^2}} - m} \right) - m,$ | (15) |
| $h = {\left( {\frac{{cQ}}{{\sqrt i }}} \right)^{3/8}}\frac{{{{\left( {\beta + 2\sqrt {1 + {m^2}} } \right)}^{1/4}}}}{{{{\left( {\beta + m} \right)}^{5/8}}}},b = \beta h.$ | (16) |
由式(12)—(16)可知，在已知Q、i、c的情况下，选定一组m和α即可唯一确定h，因此梯形渠道断面尺寸优化的独立参数为m、α。为减少参数建模和寻优计算量，根据大型渠道工程建设实际, 将独立参数的解集空间限定为m=[1, 3]，α=[1, 1.04]。
2.2 优化函数及约束条件确定评价渠道衬砌抗冻胀能力的主要指标为衬砌法向冻胀位移和应力，具体到混凝土板状结构的衬砌时以拉应力作为控制指标。由于相同冻胀荷载下，衬砌结构冻胀位移越大，应力越小，反之冻胀位移越小，应力越大。从“削减、适应冻胀”的原则出发，混凝土衬砌结构在满足位移要求的前提下，结构内拉应力越小，对冻胀变形的适应性越好，为此本文构造衬砌整体刚度指标表征其结构的冻胀适应性，
| $K = \frac{{\int_{{L_{\rm{u}}} + {L_{\rm{d}}}} {{\sigma _{\rm{n}}}{\rm{d}}L} }}{{\int_{{L_{\rm{u}}} + {L_{\rm{d}}}} {s{\rm{d}}L} }}.$ | (17) |
 |
| 图 1 衬砌结构正应力及表面长度示意图 |
| 图选项 |
在水力最佳及经济实用断面所确定的优化解集空间内，为确定使渠道满足抗冻胀最优要求的断面尺寸，结合文[1]中对梯形渠道衬砌冻胀位移材料强度的约束要求，建立以下优化函数：
| $\begin{array}{*{20}{c}}{F\left( {\alpha ,m} \right) = \min K\left( {\alpha ,m} \right),}\\{{\rm{s}}.\;{\rm{t}}.\left\{ \begin{array}{l}\beta = {\left( {\frac{{{h_{\rm{m}}}}}{h}} \right)^2}\alpha \left( {2\sqrt {1 + {m^2}} - m} \right) - m,\\m = \left[ {1,3} \right],\\\alpha = \left[ {1,1.04} \right],\\s \le 1\;{\rm{cm}},\\ - {f_{\rm{c}}} \le {\sigma _{\rm{n}}} \le {f_{\rm{t}}}.\end{array} \right.}\end{array}$ | (18) |
3 渠道形体优化的参数化数值建模及求解3.1 几何模型及网格划分以甘肃省景电工程干渠兰化段梯形混凝土衬砌渠道为例。原渠道设计流量为65 m3/s，糙率c为0.02，比降i为1/1 500，边坡系数m为1.5，底宽3 m，宽深比1.0，正常水深3.6 m，C20素混凝土衬砌板厚度8 cm，α为1.002。渠基土质为粉质黏土，因有灌溉回归水及渗漏，采取渠基排水措施后地下水位距离渠道底板1.3 m深。渠道断面设计参数在水力最优的解集空间内，因此可以直接进行抗冻胀断面优化。
保持流量不变，以参数m和α表示渠道其他断面尺寸，建立参数化的渠道有限元模型。模型建立流程如图 2a所示。采用实体单元对所建渠基与衬砌结构进行网格划分，如图 2b所示。
 |
| 图 2 渠道断面参数优化建模流程及几何模型 |
| 图选项 |
3.2 计算参数首先设置模型温度边界。依据文[21]，底部边界取渠顶以下10 m作为恒温层，温度为8 ℃，上表面采用热流边界，在软件中输入式(19)，
| $ - {q_{{\rm{up}}}} = {h_{\rm{ \mathsf{ α} }}} \cdot \left( {{T_{{\rm{ext}}}} - T} \right).$ | (19) |
 |
| 图 3 当地月平均气温 |
| 图选项 |
| ${h_{\rm{ \mathsf{ α} }}} = 3.06v + 4.11,$ | (20) |
用于土体冻胀计算的水、热、力物理场参数取值如表 1所示。
表 1 材料计算参数
| 材料 | 导热系数λ/ (W·(m·K)-1) | 比热容Cp/ (kJ·(kg·K)-1) | 弹性模量E/ MPa | Poisson比 | 渗透系数k/ (m·s-1) | 孔隙体积比 | 体积含水量 |
| 混凝土 | 1.58 | 0.97 | 2.1×104 | 0.23 | — | — | — |
| 已冻土 | 式(2) | 式(3) | 46 | 0.35 | 1×10-9 | 0.55 | 0.05 |
| 未冻土 | 式(2) | 式(3) | 15 | 0.35 | k(Δh) | 0.55 | f(Δh) |
| 冰 | 2.2 | 2.1 | — | — | — | — | — |
| 土颗粒 | 1.5 | 0.92 | — | — | — | — | — |
| 未冻水 | 0.6 | 4.2 | — | — | — | — | — |
表选项
表 1中渠基未冻土体积含水量与距离地下水位线的距离Δh有关，渠基为黏土时，地下水对冻结层无显著影响的临界值为2 m[1]。本文采用分段函数对渠基含水量进行定义，
| $f\left( {\Delta h} \right) = \left\{ {\begin{array}{*{20}{c}}\begin{array}{l}0.55,\\0.55 \cdot {{\rm{e}}^{\left[ { - 0.08\left( {\Delta h - 2} \right)} \right]}},\end{array}&\begin{array}{l}0 < \Delta h \le 2\;{\rm{m}}\;;\\\Delta h > 2\;{\rm{m}}\;.\end{array}\end{array}} \right.$ | (21) |
| $k\left( {\Delta h} \right) = 6.02 \times {10^{ - 5}} \cdot f{\left( {\Delta h} \right)^{10.61}}.$ | (22) |
图 4为不同边坡系数m和实佳比α时渠顶、渠底冻深计算值。图中白线为0 ℃线。图 5为渠道温度场分布特征。如图 5所示，由于渠顶距离地下水位线较远，土体含水量较小，由式(2)、(3)所确定的导热系数大且热容小，因此冻深最大；渠底距离地下水位线近且对流交换系数比渠顶小，因此冻深最小。但是，随着边坡系数m和实佳比α的增大，渠槽变得宽浅，渠顶逐渐靠近地下水位线，土体含水量增加，冻深有所减小；同时，开阔的渠底导致对流交换系数增大，热量散失加快，冻深加大；在两者综合作用下渠顶渠底的冻深逐渐趋近相同(图 4)。
 |
| 图 4 不同边坡系数m和实佳比α时渠顶、渠底冻深 |
| 图选项 |
 |
| 图 5 (网络版彩图)不同边坡系数m和实佳比α时渠道温度场分布特征 |
| 图选项 |
3.3.2 断面形体尺寸对衬砌冻胀位移及应力的影响根据3.3.1节冻深分析结果可知，窄深渠道渠底与渠顶冻深差异大，而宽浅渠道冻深差异小，并且随着断面变得宽浅，冻深总体变小。但是，由于渠基水分增加，迁移所形成的分凝冰层增厚导致冻胀量增大，加之断面约束减弱，因此宽浅渠道衬砌法向位移会有所增加。以往断面优化研究[6, 9, 11-13]忽略了冻胀量因断面尺寸变化导致的差异，在相同冻胀量前提下，仅考虑断面尺寸影响衬砌结构的受力性能，必然出现渠道越宽浅抗冻胀性能越好的优化结果。
根据文[1]，衬砌结构冻胀安全的评价指标包括最大法向冻胀位移和最大结构拉应力，不同设计参数下衬砌最大法向冻胀位移和拉应力如图 6a和6b所示。可以看出，当实佳比α=1时，设计断面为窄深式水力最优断面，随着边坡系数m的增加，渠坡减缓同时渠底宽度迅速减小(图 4)，衬砌最大位移和最大拉应力皆减小，且最大位移和最大应力始终满足约束条件。当α＞1时，设计断面为宽浅式实用经济断面，随着α和m的增加，衬砌最大位移增大并超出位移约束范围；最大应力减小并从应力约束范围外进入约束范围内。这说明依据工程经验一味增加底宽与边坡系数不能起到防冻胀效果，而是需要在合理的α和m的取值范围内进行甄选。
 |
| 图 6 不同设计参数下衬砌最大法向冻胀位移和截面最大拉应力 |
| 图选项 |
为综合考虑α和m对渠道冻胀适应性的影响，提取有限元计算的衬砌结构位移和应力数据，并依据式(17)进行计算，得到不同设计参数下结构上下表面的整体刚度指标K。如图 7a、7b所示，法向冻胀位移导致衬砌承受上表面拉应力和下表面压应力，由于衬砌结构安全受制于拉应力，因此以上表面整体刚度指标K作为优化目标。结合图 6a、6b和图 7a可以看出，在位移和应力约束范围内，使K取得最小值的断面设计参数取值如表 2所示。优化前后，各断面形状的湿周均在最佳水力断面所对应的最小湿周值附近微小变化，最大仅增大1.1%，然而整体刚度指标降低30%~48%，最大拉应力减小了36.4%~52.7%，并均小于C20素混凝土抗拉强度设计值1.1 MPa，可见优化后的渠道断面适应冻胀破坏能力得到了提高。
 |
| 图 7 衬砌上、下表面整体刚度指标 |
| 图选项 |
表 2 设计参数优化结果
| 实佳比α | 边坡系数m | 湿周χ | 最大位移/cm | 最大拉应力/MPa | 整体刚度指标/(MPa·cm-1) | |
| 原设计 | 1.002 | 1.5 | 17.78 | 0.70 | 1.10 | 37 |
| 最优组合1 | 1.00 | 2.0 | 17.77 | 0.50 | 0.70 | 25 |
| 最优组合2 | 1.01 | 2.5 | 17.84 | 0.60 | 0.65 | 26 |
| 最优组合3 | 1.02 | 2.5 | 17.91 | 0.95 | 0.55 | 19 |
| 最优组合4 | 1.03 | 2.0 | 17.98 | 0.95 | 0.52 | 19 |
表选项
4 温-水-土-结构耦合作用对衬砌结构优化的影响本文提出的寒区渠道形体优化设计方法，首次加入了冻土水-热-力耦合作用模型。与已往依靠经验或规范确定冻胀位移或冻胀力的研究方法相比[1, 6, 11-13]，本文方法综合考虑了渠道所处区域的气温、地下水条件和土质特征等因素与渠道衬砌结构的耦合作用，基于本文方法对冻胀过程和断面优化计算的结果更加合理。
首先，外界气温通过表面对流换热对渠基温度产生作用，渠道不同位置空气流速不同导致温度传导差异。在模型中，渠槽不同位置对流换热系数的设置(式(20))体现了渠槽形状对表面热通量传导的影响，并影响了最终的温度场计算结果。计算结果表明，窄深断面渠底保温较好，相同含水量条件下冻深较宽浅式断面浅。计算结果与实际现场地温监测情况(图 8)相符。
 |
| 图 8 景电工程渠道基土温度监测结果(2017年1月) |
| 图选项 |
其次，由于渠基土水分来源与渠道渗漏及周围农田灌溉回归水有关，冬季渠道停水时，水分在重力作用下汇聚于渠底以下的含水层内，导致渠底附近土体含水量大，而远离渠底的渠坡和渠顶含水量逐渐降低。本文模型通过式(21)对渠基含水量进行赋值来体现以上区别，一方面通过式(2)、(3)建立含水量与热参数间联系，另一方面通过式(9)建立初始含水量与原位冻胀率的关系，从而使土温与冻胀量计算结果考虑了不同位置因水分含量不同而导致的差异。
最后，在迁移冰透镜体所导致的冻胀量计算中，式(10)、(22)综合考虑了土质、地下水埋深、温度梯度及冻结时间对迁移水量的影响。图 9为不同设计参数下渠道衬砌冻胀位移分布。结果表明，渠基土冻胀量分布呈现渠底大、渠顶小的特点。随着断面实佳比α和边坡系数m的增加，一方面渠槽变得宽浅而对冻胀变形的适应性提高；另一方面由于渠底对流作用增强导致的冻深增大，加之渠顶含水量增加导致的冻胀率增大，使得衬砌法向冻胀位移增大，超过了可以恢复的最大位移值，从而不满足设计要求。这一结果在不考虑水-热-力3场耦合冻胀理论的模型中是无法得到的。因此，本文方法在渠基冻胀量的计算过程中考虑因素更多，更加符合实际。
 |
| 图 9 不同设计参数下渠道冻胀分布 |
| 图选项 |
5 结论通过建立的双目标优化模型和水-热-力3场耦合冻土冻胀模型，采用分层序列法结合多物理场耦合有限元软件对模型进行求解，进行了寒区衬砌渠道断面形体优化设计。
1) 应用冻土3场耦合冻胀理论模型，综合反映了实际渠道气温、含水量、地下水埋深和土质等因素耦合作用对渠基土冻胀变形及衬砌结构受力的影响，提高了优化结果的准确性。
2) 构造结构整体刚度指标作为衬砌冻胀适应性指标，并以整体刚度指标最小为优化目标，所得优化结果满足结构位移与应力的协调。
3) 将渠道水力特性与冻胀适应性共同作为优化目标，开发了相应大型渠道衬砌结构水力、抗冻胀双目标优化设计程序。通过工程算例计算得到了在经济实用、最大冻胀位移和最大结构拉应力约束条件下，满足最佳过水条件和最佳冻胀适用性双目标的最优设计参数组合。与原设计相比，最优设计保证湿周位于最小范围内，衬砌最大拉应力减小36.4%~52.7%，结构整体刚度指标减小30%~48%，提高了渠道适应冻胀破坏的能力。本文提出的方法可为渠道断面抗冻胀优化设计提供参考。
参考文献
| [1] | 中华人民共和国水利部.渠系工程抗冻胀设计规范: SL 23-2006[S].北京: 中国水利水电出版社, 2006. Ministry of Water Resources of the People's Republic of China.Design code for anti-frost-heave of canal and its structure: SL 23-2006[S].Beijing: China Water & Power Press, 2006.(in Chinese) |
| [2] | SWAMEE P K. Optimal irrigation canal sections[J]. Journal of Irrigation and Drainage Engineering, 1995, 121(6): 467-469. DOI:10.1061/(ASCE)0733-9437(1995)121:6(467) |
| [3] | SWAMEE P, MISHRA G, CHAHAR B. Comprehensive design of minimum cost irrigation canal sections[J]. Journal of Irrigation and Drainage Engineering, 2000, 126(5): 322-327. DOI:10.1061/(ASCE)0733-9437(2000)126:5(322) |
| [4] | 余长洪, 周明耀, 姜健俊, 等. 灌区节水改造中防渗渠道断面的优化设计[J]. 农业工程学报, 2004, 20(1): 91-94. YU C H, ZHOU M Y, JIANG J J, 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. (in Chinese) |
| [5] | 脱云飞, 王克勤, 宋维峰, 等. 灌溉渠道优化设计方法研究[J]. 中国农村水利水电, 2012(4): 58-60. TUO Y F, WANG K Q, SONG W F, et al. Research on optimal design and method of irrigation channel[J]. China Rural Water and Hydropower, 2012(4): 58-60. (in Chinese) |
| [6] | 辛英华, 王正中. U形衬砌渠道结构及水力最佳断面的分析[J]. 节水灌溉, 2008(2): 36-38, 45. XIN Y H, WANG Z Z. The structural and hydraulic analysis on frost-heaving concrete lining of the U-shape canal[J]. Water Saving Irrigation, 2008(2): 36-38, 45. (in Chinese) |
| [7] | 王正中, 李甲林, 陈涛, 等. 弧底梯形渠道砼衬砌冻胀破坏的力学模型研究[J]. 农业工程学报, 2008, 24(1): 18-23. WANG Z Z, LI J L, CHEN T, et al. Mechanics models of frost-heaving damage of concrete lining trapezoidal canal with arc-bottom[J]. Transactions of the Chinese Society of Agricultural Engineering, 2008, 24(1): 18-23. DOI:10.3321/j.issn:1002-6819.2008.01.004 (in Chinese) |
| [8] | 王正中. 梯形渠道砼衬砌冻胀破坏的力学模型研究[J]. 农业工程学报, 2004, 20(3): 24-29. WANG Z Z. 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. DOI:10.3321/j.issn:1002-6819.2004.03.006 (in Chinese) |
| [9] | 刘东, 胡宇祥, 付强, 等. 北方灌区混凝土衬砌渠道断面优化及参数分析[J]. 农业工程学报, 2015, 31(20): 107-114. LIU D, HU Y X, FU Q, 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. DOI:10.11975/j.issn.1002-6819.2015.20.016 (in Chinese) |
| [10] | 王正中, 沙际德, 蒋允静, 等. 正交各向异性冻土与建筑物相互作用的非线性有限元分析[J]. 土木工程学报, 1999, 32(3): 55-60. WANG Z Z, SHA J D, JIANG Y J, et al. Nonlinear finite element analysis of interaction of orthotropic frozen ground and construction[J]. China Civil Engineering Journal, 1999, 32(3): 55-60. DOI:10.3321/j.issn:1000-131X.1999.03.009 (in Chinese) |
| [11] | 刘旭东, 王正中, 闫长城, 等. 基于数值模拟的双层薄膜防渗衬砌渠道抗冻胀机理探讨[J]. 农业工程学报, 2011, 27(1): 29-35. LIU X D, WANG Z Z, YAN C C, et al. Exploration on anti-frost heave mechanism of lining canal with double films based on computer simulation[J]. Transactions of the Chinese Society of Agricultural Engineering, 2011, 27(1): 29-35. DOI:10.3969/j.issn.1002-6819.2011.01.005 (in Chinese) |
| [12] | 李爽, 王正中, 高兰兰, 等. 考虑混凝土衬砌板与冻土接触非线性的渠道冻胀数值模拟[J]. 水利学报, 2014, 45(4): 497-503. LI S, WANG Z Z, GAO L L, et al. Numerical simulation of canal frost heaving considering nonlinear contact between concrete lining board and soil[J]. Journal of Hydraulic Engineering, 2014, 45(4): 497-503. (in Chinese) |
| [13] | 郝晋彩, 娄宗科, 高凤. 梯形渠道冻胀数值模拟与边坡系数优化[J]. 灌溉排水学报, 2017, 36(12): 81-86. HAO J C, LOU Z K, GAO F. Numerical simulation and slope coefficient optimization of frost heave of trapezoidal channel[J]. Journal of Irrigation and Drainage, 2017, 36(12): 81-86. (in Chinese) |
| [14] | 安鹏, 邢义川, 张爱军, 等. 弧底梯形渠道抗冻胀结构优化与数值模拟[J]. 灌溉排水学报, 2017, 36(11): 56-62. AN P, XING Y C, ZHANG A J, 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. (in Chinese) |
| [15] | LIU Z, YU X. Coupled thermo-hydro-mechanical model for porous materials under frost action:Theory and implementation[J]. Acta Geotechnica, 2011, 6(2): 51-65. DOI:10.1007/s11440-011-0135-6 |
| [16] | KURYLYK B L, WATANABE K. The mathematical representation of freezing and thawing processes in variably-saturated, non-deformable soils[J]. Advances in Water Resources, 2013, 60: 160-177. DOI:10.1016/j.advwatres.2013.07.016 |
| [17] | 刘月, 王正中, 王羿, 等. 考虑水分迁移及相变对温度场影响的渠道冻胀模型[J]. 农业工程学报, 2016, 32(17): 83-88. LIU Y, WANG Z Z, WANG Y, et al. Frost heave model of canal considering influence of moisture migration and phase transformation on temperature field[J]. Transactions of the Chinese Society of Agricultural Engineering, 2016, 32(17): 83-88. (in Chinese) |
| [18] | 代鸿, 曹波. 多目标模糊优化在船舶结构设计中的应用[J]. 舰船科学技术, 2017, 39(7A): 4-6. DAI H, CAO B. Application of multi objective fuzzy optimization method in ship structure design[J]. Ship Science and Technology, 2017, 39(7A): 4-6. (in Chinese) |
| [19] | 张秀利, 梁迎春, 董申. 多目标结构模糊优化的分层序列法[J]. 机械设计与制造, 1999(5): 44-45. ZHANG X L, LIANG Y C, DONG S. Lexicograhic method for fuzzy multiobjective optimization of structures[J]. Machinery Design & Manufacture, 1999(5): 44-45. DOI:10.3969/j.issn.1001-3997.1999.05.019 (in Chinese) |
| [20] | 李炜. 水力计算手册[M]. 2版. 北京: 中国水利水电出版社, 2006. LI W. Hydraulic calculation manual[M]. 2nd ed. Beijing: China Water & Power Press, 2006. (in Chinese) |
| [21] | 刘晓燕, 赵军, 石成, 等. 土壤恒温层温度及深度研究[J]. 太阳能学报, 2007, 28(5): 494-498. LIU X Y, ZHAO J, SHI C, et al. Study on soil layer of constant temperature[J]. Acta Energiae Solaris Sinica, 2007, 28(5): 494-498. DOI:10.3321/j.issn:0254-0096.2007.05.007 (in Chinese) |
| [22] | 张建荣, 刘照球. 混凝土对流换热系数的风洞实验研究[J]. 土木工程学报, 2006, 39(9): 39-42, 61. ZHANG J R, LIU Z Q. A study on the convective heat transfer coefficient of concrete in wind tunnel experiment[J]. China Civil Engineering Journal, 2006, 39(9): 39-42, 61. DOI:10.3321/j.issn:1000-131X.2006.09.006 (in Chinese) |
| [23] | LI S Y, ZHANG M Y, PEI W S, et al. Experimental and numerical simulations on heat-water-mechanics interaction mechanism in a freezing soil[J]. Applied Thermal Engineering, 2018, 132: 209-220. DOI:10.1016/j.applthermaleng.2017.12.061 |
