即：
上述本构方程中的弹性模量E即是利用我们推导出的由 土和冰组成的饱和冻土弹性模量：E (1 i ) Es i Ei 。 冰的体积含量 i 取0.05，砂土的弹性模量 Es 为64MPa，冰 的 弹性模量随温度变化曲线由调研得到的数据拟合得到：
Ei 3.70659 0.63004 T 0.02842 T 2
剪应力 分布：
3.7 小结
从复合材料的细观机理出发，利用土力学 原理，先后建立了土、冰耦合的冻土本构模型， 土、冰、水耦合的冻土本构模型，并对温度影 响下的冻土本构关系进行了初步探讨，建立了 基本分析模型。另外，用飞箭软件对一些前人 所作的实际问题进行了数值模拟，为编制适用 于我们提出的本构模型的计算程序做了前期准 备工作。
1－Eshelby等效夹杂方法 2－微分法 3－Mori-Tanaka方法 4－Reuss等应力假设
1－Eshelby等效夹杂方法 2－微分法 3－Mori-Tanaka方法 4－Reuss等应力假设 5－Voigt等应变假设
图1 不同细观方法得到冻结砂 土弹性模量随温度变化图
图2 不同细观方法得到冻结砂 土弹性模量随温度变化图
冻结50 小时后：
冻结100 小时后：
冻结200 小时后：
冻结500 小时后：
冻结1000 小时后：
冻结2000 小时后：
冻结5000 小时后：
前述渠道静力平衡方程为： 几何方程为：
x L 0 0
L 0 Lu
0 y 0 0 0 z y x 0 0 z y z 0 x
3.2.5 混合率方法
Voigt等应变假设：
K* Ci K i
i 1
1
N
G Ci Gi
* i 1
1
N
Reuss等应力假设:
Ci K i 1 K i
* R N

Ci G i 1 Gi
* R N
Leabharlann 3.2.6 几种弹性模量的比较
其中：
D 物理方程（本构方程）为：
E x 1 x y z E y z x y 1 E z z x y 1 xy G xy G yz yz zx G zx
右边界温度初始条件原始数据：10.09
初始温度场云图：
上、下边界温度边界条件：
F1 (t ) 11.2 sin( t / 3120 .0) (h) 0.014( t 336.0 ) F ( t ) 10 . 9 e ( h) 2
冻结10小时后的渠道温度分布图：
冻结20 小时后：
(1 ) i ln(T / T0 ) n pw (1 ) i
3.6 数值模拟
3.6.1 飞箭有限元程序自动生成系统 FEPG
飞箭有限元程序自动生成系统 FEPG (Finite Element Program Generator)是一个开放、灵活的有限元程序开 发平台，它可以根据用户输入的有限元方法所需的各种 表达式和公式，自动产生解决该问题所需的全部有限元 计算的 FORTRAN 源程序。
3.3 与实际冻土的比较
不同的冰体积含量 i ，由式 E (1 i ) Es i Ei 计算得 到砂土冻结后的等效弹性模量随温度变化曲线：
冻结砂土的弹性模量随温度变化曲线
图a的应变速率为 1.1×10-4 (s-1 )， 图b的应变速率为
1.1×10-3 (s-1 )。实验曲线的围压为： 1－ 3 =8Mpa， 2－ 3 =12MPa， 3－ 3 =22MPa
四、下一阶段的工作计划
进一步完善和修正饱和冻土的本构模型，如可以 考虑粘弹性、塑性、粘塑性等不同的本构模型。 从细观的角度研究非饱和冻土的力学特性和本构 模型。对土体冻融过程中未冻水和冰、气之间的 相变机理进行研究。
用我们推导出的饱和冻土、非饱和冻土的本构模 型，用有限元法对一些工程实际问题如路基、隧 道、地基等进行数值模拟，以验证我们模型的合 理性和可靠性。
G G0
C1 G1 G0 1 (1 C1 )(G1 G0 ) (G0 G * )
3.2.4 Mori-Tanaka方法
C1 0 1 0 0 0 1 C1 1 0
G G0 C1 G0 G1 G0 G0 1 C1 G1 G0
2．参数的确定
C 3076 .7
C 2279.0
5.53
K1 1.8E 0.7
6.58
K2 0.0036 exp(0.551 (T 0.30))
K3 3.6E 3
T p 0.30
L 334.7
Tb 0.75
0.0664 exp( 0.0551 T ) T 0.75 Wu 0.3058 0.596(T 0.30) 0.75 T 0.30 0.3058 T 0.30
' x ' y ' z
3.5 温度对饱和冻土本构关系的影响
冻土体内一点处的水压、冰压以及温度之间的关 系用Clapeyron方程描述：
T L ln i T0
水压、冰压与该点的孔隙正压力 n 的关系为：
p
pi
n p (1 ) pi
水压力 pw 为：
L s K Ta
T Tp 1 ( ) Tb T T p 2 T Tb
K 3 K K 2 exp(A T ) K 1
T Tp Tb T T p T Tb
冻结砂土的应力－应变曲线及其比较（ i =0.05, D =0.5）
冻结砂土的应力－应变曲线及其比较（ i =0.05,
D
=0.8）
冻结砂土的应力－应变曲线及其比较（ i =0.1,
D
=0.5）
冻结砂土的应力－应变曲线及其比较（ i =0.1,
D
=0.8）
冻结砂土的应力－应变曲线及其比较（ i =0.2,
冰的弹性模量 随温度变化的 拟合曲线：
模型的位移边界条件为：下边界固定，上边界自由，左 右边界只有y方向的位移。 应力边界条件为：左右边界剪应力均为0。模型上边界受 相当于渠道衬切板重量和与衬切板垂直的冻结力的均布压力 P=1400N，方向垂直于边界线。计算截面模型跟前面热传导 问题的模型一样，截面的有限元网格图如下：
3．数值计算结果
计算截面的有限元网格图
下边界温度线性拟合
上边界温度线性拟合
左边界温度原 始实测数据:
左边界温度的分 段线性拟合函数:
0 y 0.25 0.1 5.6 y 1.3 0.25 y 0.3784375 T ( y ) 4 y 2.333 0.3784375 y 0.619017857 7.2 y 4.6 0.619017857 y 0.86111 0.86111 y 1.0 1.6
Wu f (T )
初始条件和边界条件为：
T ( x, y,0) ( x, y)
T ( x, L1 , t ) F2 (t )
T ( L2 , y, t ) 0 x
T ( x,0, t ) F1 (t )
T (0, y, t ) 0 x
上述各式构成渠道基土二维相变温度场的定解问题
3.2 土、冰混合本构关系的细观分析
3.2.1 Eshelby等效夹杂理论
C1 0 ( K1 K 0 ) K K0 0 ( 1 0 ) C1 G0 (G1 G0 ) G G0 G0 (G1 G0 )
3.2.3 微分法
C1 1 0 0 1 (1 C1 )( 1 0 ) ( 0 * )
1. 功能特点：FEPG 可自动产生计算单元刚度矩 阵、质量矩阵、阻尼矩阵和单元荷载向量的单 元子程序。算法程序可按照用户给出的代数方 程组的矩阵表达式及右端项表达式，自动产生 算法主程序。
2. 适用领域：FEPG系统可以提供静力分析、热 传导分析、电磁场分析、流场分析、结构力学 分析、渗流场分析、围岩稳定性分析、高边坡 稳定性分析、蜗壳结构受力分析、岔管结构分 析、水力学计算等。 3. 前、后处理系统：前后处理器FEPG.GID通过 体、面的布尔加、减、交等操作得到模型，可 快速将几何模型自动离散成各种单元，并可根 据用户需要对网格进行局部加密。后处理支持 的结果显示方式有：带状云图显示、等直线显 示、切片显示、矢量显示、变形等。
取渠道的初始温度分布为冰的弹性模量关系式中的T 值，计算得到渠道上边界受均布外载作用时内部位移场 和应立场分布。如果温度取不同冻结时刻的温度值，则 应力场可与温度场相互耦合，得到不同冻结时刻渠道内 部的温度场和位移场、应立场分布。
x方向位 移分布
y方向位 移分布
x方向应 力分布
y方向应 力分布
D
=0.5）
冻结砂土的应力－应变曲线及其比较（ i =0.2,
D
=0.8）
3.4 饱和冻土的本构关系
借用土力学中饱和土体的 Terzaghi 有效 应力原理，我们对饱和冻土体的 Terzaghi 有 效应力原理表述为：饱和冻土体中任一点的 总应力为该点有效力与孔隙水压力之和，即：
x p y p p z
3.6.2 数值模拟实例
用飞箭软件对饱水渠道基土冻结过程中的热传导问题进 行数值模拟，得到二维渠道基土冻结过程中的温度分布。 假定渠道无限长，截面为平面问题，认为土质均匀、水 热蒸发热量及其它势场忽略不计。渠道模型下边界宽3.0m； 渠底距下边界1.0m，宽0.5m；渠堤距下边界2.0m,宽1.0m。 计算截面示意图如下：