坝基渗流

在本节中，分析了大坝下面的地下水流动。大坝的几何形状如图17所示，取材自[34]第1章练习30。所有长度单位为英尺（0.3048米）。大坝上游水深20英尺，下游水深5英尺。大坝下面的土壤是各向异性的。上游的渗透系数为 $k_1= 4 k_2 = 10^{-2}$ cms，下游为 $25 k_3 = 100 k_4 = 10^{-2}$ cms。我们主要关心水力梯度，即 $\nabla h$ ，因为这是判断是否会发生管涌的指标。管涌是指土壤被地下水冲走（通常发生在下游），是一种不稳定状态。作为经验法则，当水力梯度约为1时会发生管涌。

根据第6.9.14节，我们知道稳态地下水流动的方程与热传导方程相同。总水头的等效量是温度，流速的等效量是热流。有限元分析采用国际单位制，因此英尺转换为米。此外，假设上游和下游远处有垂直隔水墙（实际上，上游距大坝中点30米，下游距大坝中点30米）。

边界条件如下：

大坝、左上和右下垂直边界以及底部水平边界是不可渗透的。这意味着这些边界上的法向水速为零，或者说热流为零。
取总水头定义中z坐标的参考点位于大坝底部（参见方程582中总水头的定义），并假设大气压为零，则上游总水头为28英尺，下游为13英尺。热等效中这对应于温度边界条件。

输入卡汇总如图18所示。完整的输入卡是示例问题的一部分。该问题是真正的二维问题，因此在CalculiX GraphiX中使用qu8单元生成网格。为了在大坝附近获得更高的分辨率，使用了偏斜网格（网格如图19所示）。

**图17:** 大坝几何形状
$\begin{figure}\epsfig{file=dam.eps,width=11cm}\end{figure}$

**图18:** 大坝问题的输入卡
$\begin{figure}\epsfig{file=daminp.eps,width=11cm}\end{figure}$

**图19:** 总水头
$\begin{figure}\epsfig{file=dam1.eps,width=10cm}\end{figure}$

**图20:** y方向的渗流速度
$\begin{figure}\epsfig{file=dam2.eps,width=10cm}\end{figure}$

在输入卡的开头定义了节点和单元的拓扑结构。CalculiX GraphiX中的qu8单元类型默认通过send命令转换为CalculiX CrunchiX中的S8（壳）单元。然而，这里平面单元更为合适。由于计算是热分析而非机械分析，采用平面应变（CPE8）还是平面应力（CPS8）单元实际上并无区别。使用*ELSET关键字定义两组不同土壤的单元集。应用恒定总水头的节点由*NSET卡定义。土壤的渗透系数对应于热分析中的热导系数。注意，渗透系数定义为正交各向异性，使用*CONDUCTIVITY,TYPE=ORO卡。此卡下面的值是x、y和z方向的渗透系数（国际单位制：m/s）。实际上z方向的值无关紧要，因为该方向没有梯度预期。*SOLID SECTION卡用于将材料分配给相应的土壤区域。*INITIAL CONDITIONS卡实际上不是必需的，因为计算是稳态的，但CalculiX CrunchiX在热传导计算中正式需要它。

在step中选择*HEAT TRANSFER, STEADY STATE计算，无需额外的时间步信息。这意味着将采用step长度（1）和初始增量大小（1）的默认值。使用*BOUNDARY卡定义上游和下游的总水头（11是温度自由度）。最后，使用*NODE PRINT、*NODE FILE和*EL FILE卡定义输出：NT是温度，或者说总水头（图19），HFL是热流，或者说地下水渗流速度（图20中的y分量）。

由于上游渗透系数高，总水头梯度较小。下游情况相反。流速在下游尤为重要。在那里，它达到约 $2.25 \times 10^{-4}$ m/s（图20中的红点），对应的水力梯度约为0.56，因为下游y方向的渗透系数为 $4 \times 10^{-4}$ m/s。这小于1，因此不会发生管涌。注意，速度自然在大坝旁边最高。

本示例展示了如何利用CalculiX GraphiX的热传导功能来求解渗流问题。这同样适用于任何其他受拉普拉斯类型方程控制的现象。