主题037：土壤-结构相互作用

土体本构模型

3.1 弹性模型

3.1.1 线弹性模型

最简单的土体本构模型是线弹性模型，假设应力-应变关系符合广义胡克定律：

$${\sigma }_{ij}=D_{ijkl}{\epsilon }_{kl}$$

对于各向同性材料，弹性矩阵可以用两个独立的弹性常数表示：

用弹性模量E和泊松比ν表示：

$$\left[\begin{array}{c}{\sigma }_x\\ {\sigma }_y\\ {\sigma }_z\\ {\tau }_{xy}\\ {\tau }_{yz}\\ {\tau }_{zx}\end{array}\right]=\frac{E}{(1+\nu )(1-2\nu )}\left[\begin{array}{cccccc}1-\nu & \nu & \nu & 0& 0& 0\\ \nu & 1-\nu & \nu & 0& 0& 0\\ \nu & \nu & 1-\nu & 0& 0& 0\\ 0& 0& 0& \frac{1-2\nu }{2}& 0& 0\\ 0& 0& 0& 0& \frac{1-2\nu }{2}& 0\\ 0& 0& 0& 0& 0& \frac{1-2\nu }{2}\end{array}\right]\left[\begin{array}{c}{\epsilon }_x\\ {\epsilon }_y\\ {\epsilon }_z\\ {\gamma }_{xy}\\ {\gamma }_{yz}\\ {\gamma }_{zx}\end{array}\right]$$

用体积模量K和剪切模量G表示：

$$p=K{\epsilon }_v,s_{ij}=2G{e}_{ij}$$

其中：

• 体积模量：$K=\frac{E}{3(1-2\nu )}$
• 剪切模量：$G=\frac{E}{2(1+\nu )}$
• 体积应变：${\epsilon }_v={\epsilon }_x+{\epsilon }_y+{\epsilon }_z$
• 偏应变：$e_{ij}={\epsilon }_{ij}-\frac{{\epsilon }_v}{3}{\delta }_{ij}$

3.1.2 邓肯-张模型（Duncan-Chang Model）

邓肯-张模型是一种非线性弹性模型，能够反映土体应力-应变关系的非线性特性。

基本假设：

• 土的应力-应变关系符合双曲线方程
• 卸载-再加载模量与初始加载模量不同

双曲线应力-应变关系：

$${\sigma }_1-{\sigma }_3=\frac{{\epsilon }_a}{a+b{\epsilon }_a}$$

其中，${\epsilon }_a$ 为轴向应变，$a$ 和 $b$ 为试验参数。

切线弹性模量：

$$E_t={\left(1-\frac{R_f(1-\sin \varphi )({\sigma }_1-{\sigma }_3)}{2c\cos \varphi +2{\sigma }_3\sin \varphi }\right)}^{2}K{p}_a{\left(\frac{{\sigma }_3}{p_a}\right)}^n$$

其中：

• $R_f$ 为破坏比
• $c,\varphi$ 为土的抗剪强度指标
• $K,n$ 为模量系数
• $p_a$ 为大气压力

切线泊松比：

$${\nu }_t=\frac{G-F\lg ({\sigma }_3/p_a)}{(1-A{)}^2}$$

其中，$A=\frac{({\sigma }_1-{\sigma }_3)R_f}{2c\cos \varphi +2{\sigma }_3\sin \varphi }$，$G,F$ 为试验参数。

3.2 弹塑性模型

3.2.1 Mohr-Coulomb模型

Mohr-Coulomb模型是最经典的土体弹塑性模型，基于Mohr-Coulomb强度准则。

屈服准则：

$$f={\sigma }_1-{\sigma }_3-({\sigma }_1+{\sigma }_3)\sin \varphi -2c\cos \varphi =0$$

或用应力不变量表示：

$$f=\frac{I_1}{3}\sin \varphi +\sqrt{J_2}\left(\cos \theta -\frac{\sin \theta \sin \varphi }{\sqrt{3}}\right)-c\cos \varphi =0$$

其中，$\theta$ 为Lode角。

流动法则：

• 关联流动法则：$g=f$，塑性势面与屈服面相同
• 非关联流动法则：$g\ne f$，通常取剪胀角$\psi <\varphi$

优点：

• 参数物理意义明确，易于通过常规三轴试验确定
• 能够较好地描述土体的剪切破坏

局限性：

• 不能反映土体的硬化/软化特性
• 屈服面在π平面上的形状与实际不符
• 不能考虑应力路径的影响

3.2.2 Drucker-Prager模型

Drucker-Prager模型是对Mohr-Coulomb模型的光滑化处理，用圆锥面代替六棱锥屈服面。

屈服准则：

$$f=\alpha I_1+\sqrt{J_2}-k=0$$

其中，$\alpha$ 和 $k$ 为材料参数，与 $c$ 和 $\varphi$ 的关系为：

• 外接圆锥：$\alpha =\frac{2\sin \varphi }{\sqrt{3}(3-\sin \varphi )}$，$k=\frac{6c\cos \varphi }{\sqrt{3}(3-\sin \varphi )}$
• 内切圆锥：$\alpha =\frac{2\sin \varphi }{\sqrt{3}(3+\sin \varphi )}$，$k=\frac{6c\cos \varphi }{\sqrt{3}(3+\sin \varphi )}$

优点：

• 屈服面光滑，数值计算稳定性好
• 便于与有限元程序结合

局限性：

• 高估平面应变条件下的土体强度
• 不能反映土的剪胀性

3.2.3 Cam-Clay模型

Cam-Clay模型是描述正常固结黏土和弱超固结黏土应力-应变特性的经典模型，基于临界状态土力学理论。

临界状态概念：

临界状态是指土体在持续剪切变形下达到的一种状态，此时：

• 体积不再变化（${\dot{\epsilon }}_v=0$）
• 有效平均应力保持不变（${\dot{p}}^{\prime }=0$）
• 偏应力保持不变（$\dot{q}=0$）
• 剪应变持续增加

屈服面方程（修正Cam-Clay）：

$$f=q^2-M^2{p}^{\prime }(p_c^{\prime }-p^{\prime })=0$$

其中：

• $M$ 为临界状态应力比，$M=6\sin {\varphi }^{\prime }/(3-\sin {\varphi }^{\prime })$
• $p_c^{\prime }$ 为预固结压力，是硬化参数

硬化规律：

$$\frac{d{p}_c^{\prime }}{p_c^{\prime }}=\frac{1+e}{\lambda -\kappa }d{\epsilon }_v^p$$

其中：

• $\lambda$ 为正常固结线斜率
• $\kappa$ 为回弹线斜率
• $e$ 为孔隙比

流动法则：

修正Cam-Clay模型采用关联流动法则：

$$d{\epsilon }_v^p=d\lambda \frac{\partial f}{\partial p^{\prime }}=d\lambda (2p^{\prime }-p_c^{\prime })$$

$$d{\epsilon }_s^p=d\lambda \frac{\partial f}{\partial q}=d\lambda \cdot 2q$$

状态边界面：

Cam-Clay模型定义了正常固结线（NCL）和临界状态线（CSL）：

• 正常固结线：$e=e_N-\lambda \ln p^{\prime }$
• 临界状态线：$e=e_{\Gamma }-\lambda \ln p^{\prime }$
• 回弹线：$e=e_{\kappa }-\kappa \ln p^{\prime }$

3.3 高级本构模型

3.3.1 硬化土模型（Hardening Soil Model）

硬化土模型（HSM）是Plaxis软件中采用的先进本构模型，能够较好地描述砂土和黏土的应力-应变特性。

屈服面：

采用双屈服面形式：

1. 剪切屈服面：
   $$f_s=\frac{q}{p^{\prime }+c\cot \varphi }-\frac{2\sin \varphi }{1-\sin \varphi }\cdot \frac{m}{m+a}=0$$

2. 帽盖屈服面（体积屈服面）：
   $$f_c=\frac{q^2}{M^2}+p^{\prime 2}-p_c^{\prime 2}=0$$

硬化规律：

• 剪切硬化：与塑性剪应变相关
• 压缩硬化：与塑性体积应变相关

刚度参数：

• 三轴排水加载模量 $E_{50}$
• 三轴固结仪加载模量 $E_{ood}$
• 卸载-再加载模量 $E_{ur}$

3.3.2 小应变刚度模型

土体在小应变范围内（<0.001%）表现出很高的刚度，随应变增加刚度迅速衰减。小应变刚度模型（如HSS模型）能够描述这一特性。

刚度衰减曲线：

$$\frac{G}{G_0}=\frac{1}{1+a{\left(\frac{\gamma }{{\gamma }_{0.7}}\right)}^b}$$

其中：

• $G_0$ 为初始剪切模量
• ${\gamma }_{0.7}$ 为剪切模量衰减至0.7$G_0$时的剪应变
• $a,b$ 为拟合参数

初始剪切模量：

$$G_0=A\cdot p_a\cdot {\left(\frac{p^{\prime }}{p_a}\right)}^m$$

其中，$A$ 和 $m$ 为材料参数。

---

---

渗流-应力耦合理论

4.1 Biot固结理论

Biot于1941年建立了饱和多孔介质的 fully coupled 固结理论，考虑了渗流场与应力场的双向耦合。

4.1.1 基本假设

1. 土体是均质、各向同性的饱和多孔介质
2. 土颗粒和孔隙水均不可压缩
3. 渗流符合Darcy定律
4. 变形符合小变形假设
5. 孔隙水渗流是连续的

4.1.2 控制方程

平衡方程（土骨架）：

$$G{\nabla }^2\mathbf{u}+(G+\lambda )\nabla (\nabla \cdot \mathbf{u})-\nabla u_w+\mathbf{f}=0$$

其中：

• $\mathbf{u}$ 为土骨架位移矢量
• $u_w$ 为孔隙水压力
• $\lambda =\frac{E\nu }{(1+\nu )(1-2\nu )}$ 为Lame常数
• $G$ 为剪切模量
• $\mathbf{f}$ 为体积力

连续性方程（孔隙流体）：

$$\nabla \cdot \left(\frac{k}{{\gamma }_w}\nabla u_w\right)=\frac{\partial }{\partial t}(\nabla \cdot \mathbf{u})+\frac{n}{K_w}\frac{\partial u_w}{\partial t}$$

对于不可压缩流体（$K_w\to \infty$）：

$$\nabla \cdot \left(\frac{k}{{\gamma }_w}\nabla u_w\right)=\frac{\partial {\epsilon }_v}{\partial t}$$

其中，${\epsilon }_v=\nabla \cdot \mathbf{u}$ 为体积应变。

4.1.3 有限元离散

采用Galerkin有限元法离散控制方程：

位移场离散：

$$\mathbf{u}={\mathbf{N}}_u\mathbf{d},\epsilon ={\mathbf{B}}_u\mathbf{d}$$

孔压场离散：

$$u_w={\mathbf{N}}_p\mathbf{p}$$

耦合有限元方程：

$$\left[\begin{array}{cc}\mathbf{K}& \mathbf{L}\\ {\mathbf{L}}^T& \mathbf{S}\end{array}\right]\left[\begin{array}{c}\mathbf{d}\\ \mathbf{p}\end{array}\right]+\left[\begin{array}{cc}0& 0\\ 0& \mathbf{H}\end{array}\right]\left[\begin{array}{c}\dot{\mathbf{d}}\\ \dot{\mathbf{p}}\end{array}\right]=\left[\begin{array}{c}{\mathbf{F}}_u\\ {\mathbf{F}}_p\end{array}\right]$$

其中：

• $\mathbf{K}={\int }_V{\mathbf{B}}_u^T\mathbf{D}{\mathbf{B}}_{u}dV$ 为刚度矩阵
• $\mathbf{L}={\int }_V{\mathbf{B}}_u^T\mathbf{m}{\mathbf{N}}_{p}dV$ 为耦合矩阵
• $\mathbf{S}={\int }_V{\mathbf{N}}_p^T\frac{1}{Q}{\mathbf{N}}_{p}dV$ 为压缩矩阵
• $\mathbf{H}={\int }_V\frac{k}{{\gamma }_w}\nabla {\mathbf{N}}_p^T\nabla {\mathbf{N}}_{p}dV$ 为渗透矩阵
• $Q$ 为Biot模量

4.2 非饱和土渗流-应力耦合

4.2.1 非饱和土的基本特性

非饱和土中存在水气两相，其力学性质比饱和土复杂得多：

基质吸力：

$$s=u_a-u_w$$

其中，$u_a$ 为孔隙气压力，$u_w$ 为孔隙水压力。

有效应力公式（Bishop公式）：

$${\sigma }^{\prime }=(\sigma -u_a)+\chi (u_a-u_w)=\sigma -u_a+\chi s$$

其中，$\chi$ 为有效应力参数，与饱和度有关（$0\le \chi \le 1$）。

4.2.2 土-水特征曲线

土-水特征曲线（SWCC）描述了基质吸力与饱和度或含水量的关系：

van Genuchten模型：

$$S_e=\frac{\theta -{\theta }_r}{{\theta }_s-{\theta }_r}={\left[1+(\alpha s{)}^n\right]}^{-m}$$

其中：

• $S_e$ 为有效饱和度
• $\theta$ 为体积含水量
• ${\theta }_r,{\theta }_s$ 为残余和饱和体积含水量
• $\alpha ,n,m$ 为拟合参数（通常 $m=1-1/n$）

渗透系数函数：

$$k(S_e)=k_s{S}_e^{0.5}{\left[1-(1-S_e^{1/m}{)}^m\right]}^2$$

其中，$k_s$ 为饱和渗透系数。

4.2.3 非饱和土控制方程

水相连续性方程：

$$\frac{\partial }{\partial t}(n{S}_r{\rho }_w)+\nabla \cdot ({\rho }_w{\mathbf{v}}_w)=0$$

气相连续性方程：

$$\frac{\partial }{\partial t}[n(1-S_r){\rho }_a]+\nabla \cdot ({\rho }_a{\mathbf{v}}_a)=0$$

耦合平衡方程：

$$\nabla \cdot {\sigma }^{\prime }+\nabla \cdot (u_a\mathbf{I})-\nabla \cdot (S_{r}s\mathbf{I})+\mathbf{f}=0$$

4.3 热-水-力耦合（THM）

在冻土工程、地热开发、核废料处置等工程中，需要考虑温度-渗流-应力的 fully coupled 耦合。

4.3.1 温度场控制方程

$$\rho C_p\frac{\partial T}{\partial t}=\nabla \cdot (\lambda \nabla T)+L_f\frac{\partial {\theta }_i}{\partial t}+Q$$

其中：

• $\rho C_p$ 为体积热容
• $\lambda$ 为导热系数
• $L_f$ 为相变潜热
• ${\theta }_i$ 为含冰量
• $Q$ 为热源项

4.3.2 冻土的本构关系

冻土的力学性质与温度密切相关：

温度对强度的影响：

$${\tau }_f=c(T)+{\sigma }_n\tan \varphi (T)$$

其中，$c(T)$ 和 $\varphi (T)$ 随温度降低而增大。

冻胀模型：

$${\dot{\epsilon }}_v^{frost}={\alpha }_f\frac{\partial T}{\partial t}$$

其中，${\alpha }_f$ 为冻胀系数。

---

数值方法

5.1 有限元法在SSI中的应用

5.1.1 有限元离散

土壤-结构相互作用问题通常采用有限元法求解。基本步骤包括：

1. 域离散：将土体和结构离散为有限单元
2. 形函数选择：
   • 位移场：通常采用二次形函数（如15节点三角形单元）
   • 孔压场：通常采用线性形函数（如6节点三角形单元）
3. 单元刚度矩阵计算：
   • 土体单元：考虑非线性本构关系
   • 结构单元：梁、板、壳单元
4. 整体刚度矩阵组装
5. 边界条件处理
6. 方程组求解

5.1.2 接触面处理

土与结构界面是SSI分析的关键，需要特殊处理：

接触面单元类型：

1. ** Goodman单元**：零厚度接触面单元
2. 界面单元：薄层实体单元
3. 接触算法：罚函数法、拉格朗日乘子法、增广拉格朗日法

界面本构模型：

剪切行为：

$$\tau =\{\begin{array}{ll}{K}_s\Delta u& \text{if }|\tau |<{\tau }_f\\ {\tau }_f\cdot \text{sgn}(\Delta \dot{u})& \text{if }|\tau |={\tau }_f\end{array}$$

其中，$K_s$ 为界面剪切刚度，${\tau }_f$ 为抗剪强度。

抗剪强度准则：

$${\tau }_f=c_i+{\sigma }_n^{\prime }\tan {\varphi }_i$$

其中，$c_i$ 和 ${\varphi }_i$ 为界面黏聚力和摩擦角，通常 ${\varphi }_i=(0.5\sim 0.8){\varphi }_{soil}$。

法向行为：

• 受拉：允许分离（${\sigma }_n=0$）
• 受压：弹性或塑性压缩

5.2 时间积分方法

5.2.1 固结问题的时间离散

固结问题需要求解时间相关的偏微分方程，常用的时间积分方法包括：

隐式向后Euler法：

$${\mathbf{K}}_{n+1}{\mathbf{u}}_{n+1}={\mathbf{F}}_{n+1}+\mathbf{C}\frac{{\mathbf{u}}_n-{\mathbf{u}}_{n+1}}{\Delta t}$$

优点：无条件稳定，可采用较大时间步长
缺点：需要求解非线性方程组

Crank-Nicolson法（梯形法则）：

$${\mathbf{K}}_{n+1/2}\frac{{\mathbf{u}}_n+{\mathbf{u}}_{n+1}}{2}={\mathbf{F}}_{n+1/2}+\mathbf{C}\frac{{\mathbf{u}}_n-{\mathbf{u}}_{n+1}}{\Delta t}$$

优点：二阶精度，数值耗散小

5.2.2 时间步长选择

固结分析中时间步长的选择至关重要：

Courant稳定性条件：

$$\Delta t\le \frac{(\Delta x{)}^2}{2c_v}$$

工程经验：

• 固结初期：采用小时间步长（$\Delta t/t_{50}<0.01$）
• 固结后期：可采用较大时间步长
• 建议采用自适应时间步长策略

5.3 求解策略

5.3.1 耦合求解方法

完全耦合求解（Monolithic）：

同时求解位移和孔压方程组：

$$\left[\begin{array}{cc}\mathbf{K}& \mathbf{L}\\ {\mathbf{L}}^T& \mathbf{S}+\theta \Delta t\mathbf{H}\end{array}\right]\left[\begin{array}{c}\Delta \mathbf{d}\\ \Delta \mathbf{p}\end{array}\right]=\left[\begin{array}{c}\Delta {\mathbf{F}}_u\\ \Delta {\mathbf{F}}_p\end{array}\right]$$

优点：精度高，无条件稳定
缺点：计算量大，系数矩阵非对称

顺序耦合求解（Staggered）：

交替求解位移和孔压：

1. 固定孔压，求解位移：$\mathbf{K}\Delta \mathbf{d}=\Delta {\mathbf{F}}_u-\mathbf{Lp}$
2. 固定位移，求解孔压：$(\mathbf{S}+\theta \Delta t\mathbf{H})\Delta \mathbf{p}=\Delta {\mathbf{F}}_p-{\mathbf{L}}^T\Delta \mathbf{d}$

优点：计算效率高，可利用现有程序
缺点：有条件稳定，可能需要迭代

5.3.2 非线性方程组求解

Newton-Raphson迭代：

$${\mathbf{K}}_T^{(i)}\Delta {\mathbf{u}}^{(i)}={\mathbf{F}}_{ext}-{\mathbf{F}}_{int}^{(i)}$$

$${\mathbf{u}}^{(i+1)}={\mathbf{u}}^{(i)}+\Delta {\mathbf{u}}^{(i)}$$

收敛判据：

• 力收敛：$\frac{\Vert {\mathbf{F}}_{ext}-{\mathbf{F}}_{int}\Vert }{\Vert {\mathbf{F}}_{ext}\Vert }<{\epsilon }_F$
• 位移收敛：$\frac{\Vert \Delta \mathbf{u}\Vert }{\Vert \mathbf{u}\Vert }<{\epsilon }_u$
• 能量收敛：$\frac{|\Delta {\mathbf{u}}^T({\mathbf{F}}_{ext}-{\mathbf{F}}_{int})|}{\Vert {\mathbf{F}}_{ext}\Vert \Vert \mathbf{u}\Vert }<{\epsilon }_E$

---

工程案例分析

案例一：一维固结分析

问题描述

某饱和黏土层厚10m，双面排水，地表瞬时施加均布荷载100kPa。已知：

• 土的渗透系数 $k=1\times 10^{-8}$ m/s
• 压缩系数 $a_v=0.5$ MPa${}^{-1}$
• 初始孔隙比 $e_0=1.0$
• 水的重度 ${\gamma }_w=10$ kN/m³

分析该土层的固结过程，计算：

1. 最终沉降量
2. 固结系数
3. 时间因子为0.5时的沉降量和孔压分布
4. 固结度达到90%所需时间

理论分析

最终沉降量：

$$S_{\infty }=\frac{a_v}{1+e_0}{\sigma }_{z}H=\frac{0.5}{1+1.0}\times 100\times 10=250\text{ mm}$$

固结系数：

$$c_v=\frac{k(1+e_0)}{a_v{\gamma }_w}=\frac{1\times 10^{-8}\times 2}{0.5\times 10^{-3}\times 10}=4\times 10^{-6}{\text{ m}}^2/\text{s}$$

时间因子：

$$T_v=\frac{c_{v}t}{H_{dr}^2}$$

对于双面排水，$H_{dr}=H/2=5$ m。

固结度与时间因子关系（Terzaghi一维固结理论）：

$$U_v=1-\sum _{m=0}^{\infty }\frac{2}{M^2}\exp (-M^2{T}_v)$$

其中，$M=\frac{\pi }{2}(2m+1)$。

近似公式（$U_v<0.6$）：

$$U_v=\sqrt{\frac{4T_v}{\pi }}$$

近似公式（$U_v>0.6$）：

$$U_v=1-\frac{8}{{\pi }^2}\exp \left(-\frac{{\pi }^2{T}_v}{4}\right)$$

固结度90%所需时间：

当 $U_v=0.9$ 时，$T_v\approx 0.848$

$$t_{90}=\frac{T_v{H}_{dr}^2}{c_v}=\frac{0.848\times 25}{4\times 10^{-6}}=5.3\times 10^6\text{ s}\approx 61\text{ days}$$

Python仿真代码

"""
案例一：一维固结分析
饱和黏土层双面排水固结问题
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
import matplotlib
matplotlib.use('Agg')

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

# 土体参数
H = 10.0  # 土层厚度 m
k = 1e-8  # 渗透系数 m/s
a_v = 0.5e-6  # 压缩系数 Pa^-1
e_0 = 1.0  # 初始孔隙比
gamma_w = 10e3  # 水的重度 N/m³
sigma_z = 100e3  # 附加应力 Pa

# 计算参数
H_dr = H / 2  # 排水路径长度（双面排水）
c_v = k * (1 + e_0) / (a_v * gamma_w)  # 固结系数
print(f"固结系数 c_v = {c_v:.2e} m²/s")

# 最终沉降量
S_inf = a_v / (1 + e_0) * sigma_z * H
print(f"最终沉降量 S_∞ = {S_inf*1000:.1f} mm")

# 空间离散
nz = 50
dz = H / nz
z = np.linspace(0, H, nz+1)

# 时间参数
T_v_target = 2.0  # 最大时间因子
t_max = T_v_target * H_dr**2 / c_v
t = np.linspace(0, t_max, 100)

# 解析解：孔压分布
def pore_pressure_analytical(z, t, c_v, H, sigma_z):
    """Terzaghi一维固结解析解"""
    nz = len(z)
    nt = len(t)
    u = np.zeros((nt, nz))
    
    for it, ti in enumerate(t):
        T_v = c_v * ti / (H/2)**2
        for iz, zi in enumerate(z):
            # 双面排水边界条件
            Z = zi / H
            u_sum = 0
            for m in range(20):  # 取前20项
                M = np.pi / 2 * (2*m + 1)
                u_sum += (2 / M) * np.sin(M * (1 - 2*Z)) * np.exp(-M**2 * T_v)
            u[it, iz] = sigma_z * u_sum
    
    return u

# 计算孔压
u = pore_pressure_analytical(z, t, c_v, H, sigma_z)

# 计算沉降量（通过孔压积分）
def consolidation_degree(T_v):
    """计算固结度"""
    if T_v < 0.2:
        return np.sqrt(4 * T_v / np.pi)
    else:
        return 1 - (8 / np.pi**2) * np.exp(-np.pi**2 * T_v / 4)

# 计算各时刻沉降量
S_t = np.zeros(len(t))
for i, ti in enumerate(t):
    T_v = c_v * ti / H_dr**2
    U_v = consolidation_degree(T_v)
    S_t[i] = U_v * S_inf

# 计算固结度90%所需时间
from scipy.optimize import fsolve
def equation(T_v):
    return consolidation_degree(T_v) - 0.9

T_v_90 = fsolve(equation, 1.0)[0]
t_90 = T_v_90 * H_dr**2 / c_v
print(f"固结度90%所需时间 t_90 = {t_90/86400:.1f} days")

# 可视化
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 图1：不同时刻的孔压分布
ax1 = axes[0, 0]
time_indices = [0, 10, 25, 50, 99]
colors = plt.cm.viridis(np.linspace(0, 1, len(time_indices)))
for idx, color in zip(time_indices, colors):
    T_v = c_v * t[idx] / H_dr**2
    ax1.plot(u[idx, :]/1000, z, color=color, linewidth=2, 
             label=f'T_v = {T_v:.3f}')
ax1.set_xlabel('孔压 u (kPa)', fontsize=11)
ax1.set_ylabel('深度 z (m)', fontsize=11)
ax1.set_title('一维固结：孔压分布随时间变化', fontsize=12)
ax1.legend(loc='upper right')
ax1.grid(True, alpha=0.3)
ax1.invert_yaxis()

# 图2：沉降-时间曲线
ax2 = axes[0, 1]
T_v_array = c_v * t / H_dr**2
ax2.plot(T_v_array, S_t * 1000, 'b-', linewidth=2, label='沉降量')
ax2.axhline(y=S_inf*1000, color='r', linestyle='--', label=f'最终沉降 = {S_inf*1000:.1f} mm')
ax2.axvline(x=T_v_90, color='g', linestyle=':', label=f'U=90% at T_v={T_v_90:.3f}')
ax2.set_xlabel('时间因子 T_v', fontsize=11)
ax2.set_ylabel('沉降量 S (mm)', fontsize=11)
ax2.set_title('沉降-时间关系曲线', fontsize=12)
ax2.legend(loc='lower right')
ax2.grid(True, alpha=0.3)

# 图3：固结度-时间因子关系
ax3 = axes[1, 0]
T_v_range = np.linspace(0, 2, 100)
U_v_range = [consolidation_degree(Tv) for Tv in T_v_range]
ax3.plot(T_v_range, U_v_range, 'b-', linewidth=2)
ax3.axhline(y=0.9, color='r', linestyle='--', label='U = 90%')
ax3.axvline(x=T_v_90, color='g', linestyle=':', label=f'T_v = {T_v_90:.3f}')
ax3.set_xlabel('时间因子 T_v', fontsize=11)
ax3.set_ylabel('固结度 U_v', fontsize=11)
ax3.set_title('固结度-时间因子关系', fontsize=12)
ax3.legend(loc='lower right')
ax3.grid(True, alpha=0.3)

# 图4：孔压等时线图
ax4 = axes[1, 1]
T_v_plot = c_v * t / H_dr**2
Z_plot = z / H
T_v_mesh, Z_mesh = np.meshgrid(T_v_plot, Z_plot)
contour = ax4.contourf(T_v_mesh, Z_mesh, u.T/1000, levels=20, cmap='RdYlBu_r')
plt.colorbar(contour, ax=ax4, label='孔压 u (kPa)')
ax4.set_xlabel('时间因子 T_v', fontsize=11)
ax4.set_ylabel('相对深度 z/H', fontsize=11)
ax4.set_title('孔压等时线分布', fontsize=12)
ax4.invert_yaxis()

plt.tight_layout()
plt.savefig('case1_consolidation.png', dpi=150, bbox_inches='tight')
print("已保存 case1_consolidation.png")
plt.close()

# 创建动画：孔压消散过程
fig, ax = plt.subplots(figsize=(10, 6))
line, = ax.plot([], [], 'b-', linewidth=2)
ax.set_xlim(0, sigma_z/1000 * 1.1)
ax.set_xlabel('孔压 u (kPa)', fontsize=11)
ax.set_ylabel('深度 z (m)', fontsize=11)
ax.set_title('一维固结过程动画', fontsize=12)
ax.grid(True, alpha=0.3)
ax.invert_yaxis()
time_text = ax.text(0.02, 0.95, '', transform=ax.transAxes, fontsize=10)

# 选择动画帧
frame_indices = np.linspace(0, len(t)-1, 50, dtype=int)

def init():
    line.set_data([], [])
    time_text.set_text('')
    return line, time_text

def update(frame):
    idx = frame_indices[frame]
    line.set_data(u[idx, :]/1000, z)
    T_v = c_v * t[idx] / H_dr**2
    U_v = consolidation_degree(T_v)
    time_text.set_text(f'T_v = {T_v:.3f}, U = {U_v*100:.1f}%')
    return line, time_text

anim = FuncAnimation(fig, update, init_func=init, frames=len(frame_indices),
                     interval=200, blit=True)
anim.save('case1_consolidation.gif', writer='pillow', fps=5, dpi=100)
print("已保存 case1_consolidation.gif")
plt.close()

print("\n案例一分析完成！")

案例二：条形基础地基沉降分析

问题描述

某条形基础宽度 $B=2$ m，埋深 $D=1$ m，承受均布荷载 $q=200$ kPa。地基为均质黏土，土性参数如下：

• 黏聚力 $c=20$ kPa
• 内摩擦角 $\varphi =20°$
• 弹性模量 $E=10$ MPa
• 泊松比 $\nu =0.35$
• 天然重度 $\gamma =18$ kN/m³

采用有限元法分析：

1. 地基中的应力分布
2. 基础沉降量
3. 地基承载力安全系数
4. 塑性区发展

理论分析

Boussinesq应力解：

对于条形基础，地基中任意点的附加应力可采用Boussinesq解或经验公式计算。

沉降计算（弹性理论法）：

$$S=\frac{qB(1-{\nu }^2)}{E}{I}_s$$

其中，$I_s$ 为影响系数，与基础形状和刚度有关。

对于柔性条形基础，$I_s\approx 2.0$。

极限承载力（Terzaghi公式）：

$$q_u=c{N}_c+\gamma D{N}_q+0.5\gamma B{N}_{\gamma }$$

对于 $\varphi =20°$：

• $N_c=17.7$
• $N_q=7.4$
• $N_{\gamma }=5.0$

$$q_u=20\times 17.7+18\times 1\times 7.4+0.5\times 18\times 2\times 5.0=354+133+90=577\text{ kPa}$$

安全系数：$F_s=q_u/q=577/200=2.89$

Python仿真代码

"""
案例二：条形基础地基沉降分析
采用有限元法分析地基应力、变形和承载力
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse import csr_matrix
from scipy.sparse.linalg import spsolve
import matplotlib
matplotlib.use('Agg')

plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

# 几何参数
B = 2.0  # 基础宽度 m
D = 1.0  # 埋深 m
L_domain = 10 * B  # 计算域宽度
H_domain = 8 * B   # 计算域深度

# 土体参数
c = 20e3  # 黏聚力 Pa
phi = np.radians(20)  # 内摩擦角 rad
E = 10e6  # 弹性模量 Pa
nu = 0.35  # 泊松比
gamma = 18e3  # 重度 N/m³

# 基础荷载
q = 200e3  # 均布荷载 Pa

# 创建网格
nx, ny = 60, 40
dx = L_domain / nx
dy = H_domain / ny

x = np.linspace(0, L_domain, nx+1)
y = np.linspace(0, H_domain, ny+1)
X, Y = np.meshgrid(x, y)

# 节点编号
def node_id(i, j):
    return j * (nx + 1) + i

n_nodes = (nx + 1) * (ny + 1)
n_dof = 2 * n_nodes

# 初始化全局刚度矩阵和载荷向量
K = np.zeros((n_dof, n_dof))
F = np.zeros(n_dof)

# 弹性矩阵（平面应变）
D = E / ((1 + nu) * (1 - 2*nu)) * np.array([
    [1 - nu, nu, 0],
    [nu, 1 - nu, 0],
    [0, 0, (1 - 2*nu)/2]
])

# 组装刚度矩阵
for j in range(ny):
    for i in range(nx):
        # 单元节点编号
        n1 = node_id(i, j)
        n2 = node_id(i+1, j)
        n3 = node_id(i+1, j+1)
        n4 = node_id(i, j+1)
        
        # 单元节点坐标
        x_nodes = np.array([X[j, i], X[j, i+1], X[j+1, i+1], X[j+1, i]])
        y_nodes = np.array([Y[j, i], Y[j, i+1], Y[j+1, i+1], Y[j+1, i]])
        
        # 4节点四边形单元刚度矩阵（简化计算）
        # 使用单点高斯积分
        xi, eta = 0, 0
        
        # 形函数导数
        dN_dxi = np.array([-0.25, 0.25, 0.25, -0.25])
        dN_deta = np.array([-0.25, -0.25, 0.25, 0.25])
        
        # Jacobian矩阵
        J = np.array([
            [np.dot(dN_dxi, x_nodes), np.dot(dN_dxi, y_nodes)],
            [np.dot(dN_deta, x_nodes), np.dot(dN_deta, y_nodes)]
        ])
        detJ = np.linalg.det(J)
        invJ = np.linalg.inv(J)
        
        # 全局坐标下的形函数导数
        dN_dx = invJ[0, 0] * dN_dxi + invJ[0, 1] * dN_deta
        dN_dy = invJ[1, 0] * dN_dxi + invJ[1, 1] * dN_deta
        
        # B矩阵
        B = np.zeros((3, 8))
        for k in range(4):
            B[0, 2*k] = dN_dx[k]
            B[1, 2*k+1] = dN_dy[k]
            B[2, 2*k] = dN_dy[k]
            B[2, 2*k+1] = dN_dx[k]
        
        # 单元刚度矩阵
        Ke = B.T @ D @ B * detJ * 4  # 4是单点积分的权重
        
        # 组装到全局矩阵
        nodes = [n1, n2, n3, n4]
        for a in range(4):
            for b in range(4):
                for p in range(2):
                    for q in range(2):
                        I = 2*nodes[a] + p
                        J = 2*nodes[b] + q
                        K[I, J] += Ke[2*a+p, 2*b+q]

# 施加重力荷载
for j in range(ny+1):
    for i in range(nx+1):
        node = node_id(i, j)
        F[2*node + 1] = -gamma * dy * dx  # y方向重力

# 施加基础荷载
base_start = int((L_domain/2 - B/2) / dx)
base_end = int((L_domain/2 + B/2) / dx)
for i in range(base_start, base_end+1):
    node = node_id(i, 0)
    F[2*node + 1] -= q * dx  # 基础均布荷载

# 边界条件
# 底部固定
for i in range(nx+1):
    node = node_id(i, ny)
    K[2*node, :] = 0
    K[:, 2*node] = 0
    K[2*node, 2*node] = 1
    K[2*node+1, :] = 0
    K[:, 2*node+1] = 0
    K[2*node+1, 2*node+1] = 1
    F[2*node] = 0
    F[2*node+1] = 0

# 左右边界水平约束
for j in range(ny+1):
    for i in [0, nx]:
        node = node_id(i, j)
        K[2*node, :] = 0
        K[:, 2*node] = 0
        K[2*node, 2*node] = 1
        F[2*node] = 0

# 求解方程组
print("求解线性方程组...")
K_sparse = csr_matrix(K)
U = spsolve(K_sparse, F)

# 提取位移
ux = U[0::2].reshape(ny+1, nx+1)
uy = U[1::2].reshape(ny+1, nx+1)

# 计算应力（在单元中心）
sigma_x = np.zeros((ny, nx))
sigma_y = np.zeros((ny, nx))
tau_xy = np.zeros((ny, nx))
sigma_vm = np.zeros((ny, nx))

for j in range(ny):
    for i in range(nx):
        n1 = node_id(i, j)
        n2 = node_id(i+1, j)
        n3 = node_id(i+1, j+1)
        n4 = node_id(i, j+1)
        
        # 单元平均位移
        u_elem = np.array([
            ux[j, i], uy[j, i],
            ux[j, i+1], uy[j, i+1],
            ux[j+1, i+1], uy[j+1, i+1],
            ux[j+1, i], uy[j+1, i]
        ])
        
        # 简化的B矩阵和应力计算
        x_nodes = np.array([X[j, i], X[j, i+1], X[j+1, i+1], X[j+1, i]])
        y_nodes = np.array([Y[j, i], Y[j, i+1], Y[j+1, i+1], Y[j+1, i]])
        
        dN_dxi = np.array([-0.25, 0.25, 0.25, -0.25])
        dN_deta = np.array([-0.25, -0.25, 0.25, 0.25])
        
        J = np.array([
            [np.dot(dN_dxi, x_nodes), np.dot(dN_dxi, y_nodes)],
            [np.dot(dN_deta, x_nodes), np.dot(dN_deta, y_nodes)]
        ])
        invJ = np.linalg.inv(J)
        
        dN_dx = invJ[0, 0] * dN_dxi + invJ[0, 1] * dN_deta
        dN_dy = invJ[1, 0] * dN_dxi + invJ[1, 1] * dN_deta
        
        B = np.zeros((3, 8))
        for k in range(4):
            B[0, 2*k] = dN_dx[k]
            B[1, 2*k+1] = dN_dy[k]
            B[2, 2*k] = dN_dy[k]
            B[2, 2*k+1] = dN_dx[k]
        
        epsilon = B @ u_elem
        sigma = D @ epsilon
        
        sigma_x[j, i] = sigma[0]
        sigma_y[j, i] = sigma[1]
        tau_xy[j, i] = sigma[2]
        sigma_vm[j, i] = np.sqrt(sigma[0]**2 - sigma[0]*sigma[1] + sigma[1]**2 + 3*sigma[2]**2)

# 计算极限承载力（Terzaghi公式）
N_c = (np.exp(np.pi*np.tan(phi)) * np.tan(np.pi/4 + phi/2)**2 - 1) / np.tan(phi)
N_q = np.exp(np.pi*np.tan(phi)) * np.tan(np.pi/4 + phi/2)**2
N_gamma = 2 * (N_q + 1) * np.tan(phi)

q_u = c * N_c + gamma * D * N_q + 0.5 * gamma * B * N_gamma
F_s = q_u / q

print(f"\n地基承载力分析：")
print(f"承载力系数：N_c = {N_c:.2f}, N_q = {N_q:.2f}, N_gamma = {N_gamma:.2f}")
print(f"极限承载力 q_u = {q_u/1000:.1f} kPa")
print(f"安全系数 F_s = {F_s:.2f}")

# 基础中心沉降
 center_idx = nx // 2
S_center = -uy[0, center_idx] * 1000  # mm
print(f"\n基础中心沉降 S = {S_center:.2f} mm")

# 可视化
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 图1：位移云图
ax1 = axes[0, 0]
displacement = np.sqrt(ux**2 + uy**2)
im1 = ax1.contourf(X, Y, displacement*1000, levels=20, cmap='viridis')
plt.colorbar(im1, ax=ax1, label='位移 (mm)')
ax1.set_xlabel('x (m)', fontsize=11)
ax1.set_ylabel('y (m)', fontsize=11)
ax1.set_title('地基位移分布', fontsize=12)
ax1.invert_yaxis()

# 图2：竖向位移剖面
ax2 = axes[0, 1]
ax2.plot(-uy[:, center_idx]*1000, Y[:, center_idx], 'b-', linewidth=2)
ax2.set_xlabel('竖向位移 (mm)', fontsize=11)
ax2.set_ylabel('深度 y (m)', fontsize=11)
ax2.set_title('基础中心线下竖向位移', fontsize=12)
ax2.grid(True, alpha=0.3)
ax2.invert_yaxis()

# 图3：应力云图
ax3 = axes[1, 0]
X_center = (X[:-1, :-1] + X[:-1, 1:] + X[1:, :-1] + X[1:, 1:]) / 4
Y_center = (Y[:-1, :-1] + Y[:-1, 1:] + Y[1:, :-1] + Y[1:, 1:]) / 4
im3 = ax3.contourf(X_center, Y_center, sigma_y/1000, levels=20, cmap='RdYlBu_r')
plt.colorbar(im3, ax=ax3, label='竖向应力 σ_y (kPa)')
ax3.set_xlabel('x (m)', fontsize=11)
ax3.set_ylabel('y (m)', fontsize=11)
ax3.set_title('竖向应力分布', fontsize=12)
ax3.invert_yaxis()

# 图4：Von Mises应力
ax4 = axes[1, 1]
im4 = ax4.contourf(X_center, Y_center, sigma_vm/1000, levels=20, cmap='hot')
plt.colorbar(im4, ax=ax4, label='Von Mises应力 (kPa)')
ax4.set_xlabel('x (m)', fontsize=11)
ax4.set_ylabel('y (m)', fontsize=11)
ax4.set_title('Von Mises应力分布', fontsize=12)
ax4.invert_yaxis()

plt.tight_layout()
plt.savefig('case2_foundation.png', dpi=150, bbox_inches='tight')
print("已保存 case2_foundation.png")
plt.close()

print("\n案例二分析完成！")

案例三：基坑开挖稳定性分析

问题描述

某基坑开挖深度 $H=8$ m，采用悬臂式支护桩，桩长 $L=16$ m，桩径 $d=0.8$ m，桩间距 $s=2$ m。地基土为均质砂土，参数如下：

• 有效内摩擦角 ${\varphi }^{\prime }=30°$
• 有效黏聚力 $c^{\prime }=5$ kPa
• 饱和重度 ${\gamma }_{sat}=20$ kN/m³
• 渗透系数 $k=1\times 10^{-5}$ m/s
• 弹性模量 $E=30$ MPa
• 泊松比 $\nu =0.3$

地下水位位于地表下2m。分析：

1. 基坑开挖后的土压力分布
2. 支护桩的弯矩和变形
3. 坑底抗隆起稳定性
4. 渗流对稳定性的影响

理论分析

Rankine土压力理论：

主动土压力系数：
$$K_a={\tan }^2\left(45°-\frac{{\varphi }^{\prime }}{2}\right)={\tan }^2(30°)=0.333$$

被动土压力系数：
$$K_p={\tan }^2\left(45°+\frac{{\varphi }^{\prime }}{2}\right)={\tan }^2(60°)=3.0$$

主动土压力分布（考虑地下水）：

• 0-2m（地下水位以上）：$p_a=K_a\gamma z$
• 2-8m（地下水位以下）：$p_a=K_a{\gamma }^{\prime }z+{\gamma }_w{z}_w$

其中，有效重度 ${\gamma }^{\prime }={\gamma }_{sat}-{\gamma }_w=10$ kN/m³。

坑底抗隆起稳定性（Prandtl公式）：

$$F_s=\frac{N_{c}c+\gamma D{N}_q+0.5\gamma B{N}_{\gamma }}{\gamma H+q}$$

对于 $\varphi =30°$：$N_c=30.1$, $N_q=18.4$, $N_{\gamma }=22.4$

Python仿真代码

"""
案例三：基坑开挖稳定性分析
悬臂式支护桩的受力变形分析
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from scipy.optimize import brentq
import matplotlib
matplotlib.use('Agg')

plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

# 基坑参数
H_exc = 8.0  # 开挖深度 m
L_pile = 16.0  # 桩长 m
d_pile = 0.8  # 桩径 m
s_pile = 2.0  # 桩间距 m

# 土体参数
phi = np.radians(30)  # 内摩擦角
c = 5e3  # 黏聚力 Pa
gamma_sat = 20e3  # 饱和重度 N/m³
gamma_w = 10e3  # 水重度 N/m³
gamma_d = gamma_sat - gamma_w  # 浮重度
k = 1e-5  # 渗透系数 m/s
E_soil = 30e6  # 土体弹性模量 Pa
nu_soil = 0.3  # 土体泊松比

# 支护桩参数
E_pile = 30e9  # 混凝土弹性模量 Pa
I_pile = np.pi * d_pile**4 / 64  # 桩截面惯性矩 m⁴
EI = E_pile * I_pile  # 抗弯刚度

# 地下水位
z_wt = 2.0  # 地下水位深度 m

# 土压力系数
K_a = np.tan(np.pi/4 - phi/2)**2  # 主动土压力系数
K_p = np.tan(np.pi/4 + phi/2)**2  # 被动土压力系数
print(f"主动土压力系数 K_a = {K_a:.3f}")
print(f"被动土压力系数 K_p = {K_p:.3f}")

# 深度离散
z = np.linspace(0, L_pile, 200)
dz = z[1] - z[0]

# 计算土压力分布
def earth_pressure(z_pos, H_exc):
    """计算深度z处的土压力"""
    p_a = np.zeros_like(z_pos)
    p_p = np.zeros_like(z_pos)
    
    for i, zi in enumerate(z_pos):
        if zi <= z_wt:
            # 地下水位以上
            sigma_v = gamma_sat * zi
            u = 0
        else:
            # 地下水位以下
            sigma_v = gamma_sat * z_wt + gamma_d * (zi - z_wt)
            u = gamma_w * (zi - z_wt)
        
        sigma_v_eff = sigma_v - u
        
        # 主动土压力（开挖侧）
        if zi <= H_exc:
            p_a[i] = 0  # 开挖面以上无土压力
        else:
            p_a[i] = K_a * sigma_v_eff - 2 * c * np.sqrt(K_a) + u
            if p_a[i] < 0:
                p_a[i] = 0
        
        # 被动土压力（坑内侧）
        if zi <= H_exc:
            # 坑内土被挖除，被动土压力从坑底开始
            depth_below_exc = H_exc - zi
            if depth_below_exc > 0:
                sigma_v_p = gamma_sat * min(depth_below_exc, z_wt) + \
                           gamma_d * max(0, depth_below_exc - z_wt)
                u_p = gamma_w * max(0, depth_below_exc - z_wt)
                sigma_v_p_eff = sigma_v_p - u_p
                p_p[i] = K_p * sigma_v_p_eff + 2 * c * np.sqrt(K_p) + u_p
        else:
            # 坑底以下
            depth_below_exc = 0
            sigma_v_p = gamma_sat * min(H_exc, z_wt) + gamma_d * max(0, H_exc - z_wt)
            u_p = gamma_w * max(0, H_exc - z_wt)
            sigma_v_p_eff = sigma_v_p - u_p
            p_p[i] = K_p * sigma_v_p_eff + 2 * c * np.sqrt(K_p) + u_p
    
    return p_a, p_p

# 计算土压力
p_active, p_passive = earth_pressure(z, H_exc)

# 净土压力
p_net = p_active - p_passive

# 桩的挠曲微分方程求解
# EI * d⁴y/dz⁴ = p(z)
# 采用有限差分法

n = len(z)
h = dz

# 构建四阶差分矩阵
A = np.zeros((n, n))
for i in range(2, n-2):
    A[i, i-2] = 1
    A[i, i-1] = -4
    A[i, i] = 6
    A[i, i+1] = -4
    A[i, i+2] = 1

A = A / h**4

# 边界条件
# 桩顶自由：弯矩=0，剪力=0 -> d²y/dz²=0, d³y/dz³=0
# 桩底固定：位移=0，转角=0 -> y=0, dy/dz=0

# 顶部边界（i=0,1）
A[0, 0] = 1  # 简化处理
A[1, 1] = 1

# 底部边界（i=n-2, n-1）
A[n-2, n-2] = 1
A[n-2, n-1] = 0
A[n-1, n-2] = 0
A[n-1, n-1] = 1

# 右端项
b = p_net / EI
b[0] = 0
b[1] = 0
b[n-2] = 0
b[n-1] = 0

# 求解挠度
y_pile = np.linalg.solve(A, b)

# 计算弯矩 M = -EI * d²y/dz²
M_pile = np.zeros(n)
for i in range(1, n-1):
    M_pile[i] = -EI * (y_pile[i-1] - 2*y_pile[i] + y_pile[i+1]) / h**2

# 计算剪力 V = dM/dz
V_pile = np.zeros(n)
for i in range(1, n-1):
    V_pile[i] = (M_pile[i+1] - M_pile[i-1]) / (2*h)

# 坑底抗隆起稳定性分析
# 采用Prandtl公式
N_c = (np.exp(np.pi*np.tan(phi)) * np.tan(np.pi/4 + phi/2)**2 - 1) / np.tan(phi)
N_q = np.exp(np.pi*np.tan(phi)) * np.tan(np.pi/4 + phi/2)**2
N_gamma = 2 * (N_q + 1) * np.tan(phi)

# 坑底隆起阻力
q_bearing = c * N_c + gamma_sat * H_exc * N_q

# 坑底隆起荷载（土压力）
q_load = gamma_sat * H_exc

# 安全系数
F_s_heave = q_bearing / q_load

print(f"\n基坑稳定性分析：")
print(f"坑底抗隆起安全系数 F_s = {F_s_heave:.2f}")
print(f"桩顶水平位移 = {y_pile[0]*1000:.2f} mm")
print(f"最大弯矩 = {np.max(np.abs(M_pile))/1000:.1f} kN·m")

# 渗流分析（简化）
# 渗流力 j = i * gamma_w
# 水力梯度 i = H_w / L_seepage
H_w = H_exc - z_wt  # 水头差
L_seepage = H_exc  # 渗流路径（简化）
i_seepage = H_w / L_seepage
j_seepage = i_seepage * gamma_w

print(f"\n渗流分析：")
print(f"水力梯度 i = {i_seepage:.3f}")
print(f"渗流力 j = {j_seepage/1000:.2f} kN/m³")

# 可视化
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 图1：土压力分布
ax1 = axes[0, 0]
ax1.fill_betweenx(z, p_active/1000, alpha=0.3, color='red', label='主动土压力')
ax1.fill_betweenx(z, -p_passive/1000, alpha=0.3, color='blue', label='被动土压力')
ax1.plot(p_active/1000, z, 'r-', linewidth=2)
ax1.plot(-p_passive/1000, z, 'b-', linewidth=2)
ax1.axhline(y=H_exc, color='k', linestyle='--', label='开挖面')
ax1.axhline(y=z_wt, color='cyan', linestyle=':', label='地下水位')
ax1.set_xlabel('土压力 (kPa)', fontsize=11)
ax1.set_ylabel('深度 z (m)', fontsize=11)
ax1.set_title('土压力分布', fontsize=12)
ax1.legend(loc='upper right')
ax1.grid(True, alpha=0.3)
ax1.invert_yaxis()

# 图2：桩身挠度
ax2 = axes[0, 1]
ax2.plot(y_pile*1000, z, 'g-', linewidth=2)
ax2.fill_betweenx(z, y_pile*1000, alpha=0.3, color='green')
ax2.axhline(y=H_exc, color='k', linestyle='--', label='开挖面')
ax2.set_xlabel('水平位移 (mm)', fontsize=11)
ax2.set_ylabel('深度 z (m)', fontsize=11)
ax2.set_title('支护桩挠度', fontsize=12)
ax2.legend(loc='upper right')
ax2.grid(True, alpha=0.3)
ax2.invert_yaxis()

# 图3：桩身弯矩
ax3 = axes[1, 0]
ax3.plot(M_pile/1000, z, 'purple', linewidth=2)
ax3.fill_betweenx(z, M_pile/1000, alpha=0.3, color='purple')
ax3.axhline(y=H_exc, color='k', linestyle='--', label='开挖面')
ax3.axvline(x=0, color='k', linewidth=0.5)
ax3.set_xlabel('弯矩 (kN·m)', fontsize=11)
ax3.set_ylabel('深度 z (m)', fontsize=11)
ax3.set_title('支护桩弯矩分布', fontsize=12)
ax3.legend(loc='upper right')
ax3.grid(True, alpha=0.3)
ax3.invert_yaxis()

# 图4：基坑剖面示意图
ax4 = axes[1, 1]
# 绘制土体
ax4.fill_between([0, 10], [0, 0], [L_pile, L_pile], alpha=0.3, color='brown', label='土体')
# 绘制开挖区域
ax4.fill_between([3, 7], [0, 0], [H_exc, H_exc], alpha=0.8, color='white', edgecolor='black')
# 绘制支护桩
ax4.plot([5, 5], [0, L_pile], 'k-', linewidth=8, label='支护桩')
# 绘制地下水位
ax4.axhline(y=z_wt, color='cyan', linestyle='--', linewidth=2, label='地下水位')
# 绘制开挖面
ax4.axhline(y=H_exc, color='gray', linestyle='-', linewidth=2)
ax4.set_xlim(0, 10)
ax4.set_ylim(L_pile, 0)
ax4.set_xlabel('水平距离 (m)', fontsize=11)
ax4.set_ylabel('深度 (m)', fontsize=11)
ax4.set_title('基坑剖面示意图', fontsize=12)
ax4.legend(loc='upper right')
ax4.set_aspect('equal')

plt.tight_layout()
plt.savefig('case3_excavation.png', dpi=150, bbox_inches='tight')
print("已保存 case3_excavation.png")
plt.close()

print("\n案例三分析完成！")

案例四：边坡稳定性分析

问题描述

某均质土坡，坡高 $H=10$ m，坡角 $\beta =30°$。土体参数：

• 黏聚力 $c=15$ kPa
• 内摩擦角 $\varphi =25°$
• 重度 $\gamma =19$ kN/m³

分析：

1. 采用瑞典条分法计算安全系数
2. 采用Bishop简化法计算安全系数
3. 确定最危险滑动面
4. 降雨入渗对边坡稳定性的影响

理论分析

瑞典条分法（Fellenius法）：

安全系数定义为抗滑力矩与滑动力矩之比：

$$F_s=\frac{\sum (c_i{l}_i+W_i\cos {\alpha }_i\tan {\varphi }_i)}{\sum W_i\sin {\alpha }_i}$$

其中：

• $W_i$ 为第i土条的重量
• ${\alpha }_i$ 为第i土条底面倾角
• $l_i$ 为第i土条底面长度

Bishop简化法：

$$F_s=\frac{\sum \frac{c_i{b}_i+W_i\tan {\varphi }_i}{m_{\alpha i}}}{\sum W_i\sin {\alpha }_i}$$

其中：

$$m_{\alpha i}=\cos {\alpha }_i\left(1+\frac{\tan {\alpha }_i\tan {\varphi }_i}{F_s}\right)$$

由于 $m_{\alpha i}$ 中包含 $F_s$，需要迭代求解。

Python仿真代码

"""
案例四：边坡稳定性分析
采用条分法计算边坡安全系数
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve
import matplotlib
matplotlib.use('Agg')

plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

# 边坡参数
H_slope = 10.0  # 坡高 m
beta = np.radians(30)  # 坡角

# 土体参数
c = 15e3  # 黏聚力 Pa
phi = np.radians(25)  # 内摩擦角
gamma = 19e3  # 重度 N/m³

# 滑动面参数（圆弧滑动面）
# 圆心坐标 (x0, y0)，半径 R
x0 = 5.0
y0 = 15.0
R = 15.0

# 土条划分
n_slices = 20

# 生成滑动面
alpha_range = np.linspace(np.radians(10), np.radians(80), n_slices+1)
x_slip = x0 - R * np.sin(alpha_range)
y_slip = y0 - R * np.cos(alpha_range)

# 计算每个土条的参数
W = np.zeros(n_slices)  # 土条重量
alpha = np.zeros(n_slices)  # 底面倾角
b = np.zeros(n_slices)  # 土条宽度
l = np.zeros(n_slices)  # 底面长度
u = np.zeros(n_slices)  # 孔隙水压力（初始为0）

for i in range(n_slices):
    # 土条几何参数
    x_left = x_slip[i]
    x_right = x_slip[i+1]
    y_left = y_slip[i]
    y_right = y_slip[i+1]
    
    b[i] = x_right - x_left
    l[i] = np.sqrt((x_right - x_left)**2 + (y_right - y_left)**2)
    alpha[i] = np.arctan((y_left - y_right) / (x_right - x_left))
    
    # 计算土条重量
    # 简化的坡面线
    if x_left < H_slope / np.tan(beta):
        y_ground_left = H_slope - x_left * np.tan(beta)
    else:
        y_ground_left = 0
    
    if x_right < H_slope / np.tan(beta):
        y_ground_right = H_slope - x_right * np.tan(beta)
    else:
        y_ground_right = 0
    
    # 土条平均高度
    h_avg = ((y_ground_left + y_ground_right) / 2) - ((y_left + y_right) / 2)
    if h_avg < 0:
        h_avg = 0
    
    W[i] = gamma * h_avg * b[i]

# 瑞典条分法
def swedish_method(W, alpha, l, c, phi):
    """瑞典条分法计算安全系数"""
    numerator = np.sum(c * l + W * np.cos(alpha) * np.tan(phi))
    denominator = np.sum(W * np.sin(alpha))
    return numerator / denominator

F_swedish = swedish_method(W, alpha, b/np.cos(alpha), c, phi)
print(f"瑞典条分法安全系数: F_s = {F_swedish:.3f}")

# Bishop简化法（需要迭代）
def bishop_method(F_guess, W, alpha, b, c, phi):
    """Bishop简化法计算安全系数"""
    m_alpha = np.cos(alpha) * (1 + np.tan(alpha) * np.tan(phi) / F_guess)
    numerator = np.sum((c * b + W * np.tan(phi)) / m_alpha)
    denominator = np.sum(W * np.sin(alpha))
    return numerator / denominator

# 迭代求解
def bishop_iteration(F_guess):
    return bishop_method(F_guess, W, alpha, b, c, phi)

F_bishop = fsolve(lambda F: bishop_iteration(F) - F, 1.5)[0]
print(f"Bishop简化法安全系数: F_s = {F_bishop:.3f}")

# 考虑降雨入渗的影响
# 假设降雨导致孔隙水压力增加
# 简化模型：孔隙水压力与深度成正比
r_factor = 0.5  # 降雨影响系数（0-1）
u_rain = r_factor * gamma * H_slope * 0.3  # 平均孔隙水压力

# 有效应力法计算
W_eff = W - u_rain * b / np.cos(alpha)
F_swedish_rain = swedish_method(W_eff, alpha, b/np.cos(alpha), c, phi)
F_bishop_rain = fsolve(lambda F: bishop_method(F, W_eff, alpha, b, c, phi) - F, 1.0)[0]

print(f"\n降雨后安全系数:")
print(f"瑞典条分法: F_s = {F_swedish_rain:.3f}")
print(f"Bishop简化法: F_s = {F_bishop_rain:.3f}")

# 搜索最危险滑动面（简化）
# 改变圆心位置和半径，寻找最小安全系数
F_min = float('inf')
best_params = None

for x0_test in np.linspace(3, 8, 10):
    for y0_test in np.linspace(12, 20, 10):
        R_test = np.sqrt(x0_test**2 + (y0_test - H_slope)**2)
        if R_test < 8 or R_test > 25:
            continue
        
        # 重新计算滑动面
        alpha_test = np.linspace(np.radians(10), np.radians(80), n_slices+1)
        x_test = x0_test - R_test * np.sin(alpha_test)
        y_test = y0_test - R_test * np.cos(alpha_test)
        
        # 检查滑动面是否合理
        if np.any(y_test < 0) or np.any(x_test < 0):
            continue
        
        # 计算土条参数（简化）
        W_test = np.zeros(n_slices)
        alpha_slice = np.zeros(n_slices)
        b_test = np.zeros(n_slices)
        
        for i in range(n_slices):
            b_test[i] = x_test[i+1] - x_test[i]
            alpha_slice[i] = np.arctan((y_test[i] - y_test[i+1]) / (x_test[i+1] - x_test[i]))
            # 简化的重量计算
            h_avg = max(0, H_slope - (y_test[i] + y_test[i+1])/2)
            W_test[i] = gamma * h_avg * b_test[i]
        
        F_test = swedish_method(W_test, alpha_slice, b_test/np.cos(alpha_slice), c, phi)
        
        if F_test < F_min and F_test > 0.5:
            F_min = F_test
            best_params = (x0_test, y0_test, R_test)

print(f"\n最危险滑动面参数:")
print(f"圆心: ({best_params[0]:.2f}, {best_params[1]:.2f})")
print(f"半径: {best_params[2]:.2f} m")
print(f"最小安全系数: F_s = {F_min:.3f}")

# 可视化
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# 图1：边坡剖面与滑动面
ax1 = axes[0]
# 绘制坡面
x_ground = np.array([0, H_slope/np.tan(beta), H_slope/np.tan(beta) + 5])
y_ground = np.array([H_slope, 0, 0])
ax1.fill_between(x_ground, y_ground, alpha=0.3, color='brown', label='土坡')
ax1.plot(x_ground, y_ground, 'brown', linewidth=2)

# 绘制滑动面
ax1.plot(x_slip, y_slip, 'r--', linewidth=2, label='滑动面')
ax1.plot(x0, y0, 'r+', markersize=10, markeredgewidth=2, label='圆心')

# 绘制土条（部分）
for i in range(0, n_slices, 3):
    x_rect = [x_slip[i], x_slip[i+1], x_slip[i+1], x_slip[i]]
    y_ground_i = H_slope - (x_slip[i] + x_slip[i+1])/2 * np.tan(beta) if (x_slip[i] + x_slip[i+1])/2 < H_slope/np.tan(beta) else 0
    y_rect = [y_slip[i], y_slip[i+1], y_ground_i, y_ground_i]
    ax1.fill(x_rect, y_rect, alpha=0.2, color='yellow', edgecolor='orange')

ax1.set_xlim(-1, 25)
ax1.set_ylim(-1, 20)
ax1.set_xlabel('水平距离 (m)', fontsize=11)
ax1.set_ylabel('高程 (m)', fontsize=11)
ax1.set_title('边坡稳定性分析 - 条分法', fontsize=12)
ax1.legend(loc='upper right')
ax1.set_aspect('equal')
ax1.grid(True, alpha=0.3)

# 图2：安全系数对比
ax2 = axes[1]
methods = ['瑞典条分法\n(干燥)', 'Bishop法\n(干燥)', '瑞典条分法\n(降雨)', 'Bishop法\n(降雨)']
F_values = [F_swedish, F_bishop, F_swedish_rain, F_bishop_rain]
colors = ['green', 'blue', 'orange', 'red']

bars = ax2.bar(methods, F_values, color=colors, alpha=0.7, edgecolor='black')
ax2.axhline(y=1.0, color='red', linestyle='--', linewidth=2, label='临界值 F_s=1.0')
ax2.axhline(y=1.3, color='orange', linestyle='--', linewidth=2, label='最小安全值 F_s=1.3')
ax2.set_ylabel('安全系数 F_s', fontsize=11)
ax2.set_title('边坡安全系数对比', fontsize=12)
ax2.legend(loc='upper right')
ax2.grid(True, alpha=0.3, axis='y')

# 在柱状图上标注数值
for bar, F_val in zip(bars, F_values):
    height = bar.get_height()
    ax2.text(bar.get_x() + bar.get_width()/2., height,
             f'{F_val:.3f}', ha='center', va='bottom', fontsize=10)

plt.tight_layout()
plt.savefig('case4_slope_stability.png', dpi=150, bbox_inches='tight')
print("已保存 case4_slope_stability.png")
plt.close()

print("\n案例四分析完成！")

案例五：桩-土相互作用分析

问题描述

某钻孔灌注桩，桩径 $d=1.0$ m，桩长 $L=20$ m，承受水平荷载 $H=200$ kN。地基为分层土：

• 0-5m：软黏土，$c_u=30$ kPa，$E_s=5$ MPa
• 5-12m：粉质黏土，$c_u=50$ kPa，$E_s=15$ MPa
• 12-20m：密实砂土，$N=30$，$E_s=50$ MPa

分析：

1. 桩身水平位移分布
2. 桩身弯矩和剪力分布
3. 桩侧土抗力分布
4. p-y曲线分析

理论分析

水平受荷桩的挠曲方程：

$$EI\frac{d^{4}y}{d{z}^4}+p(y,z)=0$$

其中，$p(y,z)$ 为桩侧土抗力，通常表示为：

$$p=k_h\cdot y$$

$k_h$ 为地基水平反力系数。

p-y曲线方法：

p-y曲线描述了桩侧土抗力 $p$ 与桩身位移 $y$ 之间的非线性关系。对于黏土，常用Matlock提出的p-y曲线：

$$\frac{p}{p_u}=0.5{\left(\frac{y}{y_{50}}\right)}^{1/3}$$

其中：

• $p_u$ 为极限土抗力
• $y_{50}=2.5{\epsilon }_{50}d$ 为对应50%极限抗力的位移

Python仿真代码

"""
案例五：桩-土相互作用分析
水平受荷桩的受力变形分析
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import solve_banded
import matplotlib
matplotlib.use('Agg')

plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

# 桩参数
d_pile = 1.0  # 桩径 m
L_pile = 20.0  # 桩长 m
E_concrete = 30e9  # 混凝土弹性模量 Pa
I_pile = np.pi * d_pile**4 / 64  # 截面惯性矩 m⁴
EI = E_concrete * I_pile  # 抗弯刚度

# 水平荷载
H_load = 200e3  # 水平力 N
M_load = 0  # 弯矩（桩顶自由）

# 土层参数（分层）
# 层1：0-5m 软黏土
# 层2：5-12m 粉质黏土
# 层3：12-20m 密实砂土

layer_depths = [0, 5, 12, 20]  # 层界面深度
layer_c_u = [30e3, 50e3, 0]  # 不排水抗剪强度（砂土为0）
layer_E_s = [5e6, 15e6, 50e6]  # 土体弹性模量
layer_N = [0, 0, 30]  # 标贯击数
layer_gamma = [18e3, 19e3, 20e3]  # 重度

# 空间离散
nz = 100
dz = L_pile / nz
z = np.linspace(0, L_pile, nz+1)

# 确定每层土的性质
def get_soil_properties(z_pos):
    """获取深度z处的土体性质"""
    for i in range(len(layer_depths)-1):
        if layer_depths[i] <= z_pos <= layer_depths[i+1]:
            return {
                'c_u': layer_c_u[i],
                'E_s': layer_E_s[i],
                'N': layer_N[i],
                'gamma': layer_gamma[i]
            }
    return {
        'c_u': layer_c_u[-1],
        'E_s': layer_E_s[-1],
        'N': layer_N[-1],
        'gamma': layer_gamma[-1]
    }

# 计算地基反力系数 k_h
# 采用API RP 2A推荐的方法
k_h = np.zeros(nz+1)
for i, zi in enumerate(z):
    soil = get_soil_properties(zi)
    if soil['c_u'] > 0:  # 黏土
        # 黏土的地基反力系数
        k_h[i] = 67 * soil['c_u'] / d_pile
    else:  # 砂土
        # 砂土的地基反力系数
        k_h[i] = 0.5 * soil['gamma'] * zi * np.sqrt(soil['N'])

# 有限差分法求解挠曲方程
# EI * d⁴y/dz⁴ + k_h * y = 0

# 构建系数矩阵（带状矩阵）
# 使用五对角矩阵
n = nz + 1
AB = np.zeros((5, n))

# 内部节点
for i in range(2, n-2):
    AB[0, i+2] = EI / dz**4  # 上二对角
    AB[1, i+1] = -4 * EI / dz**4  # 上一对角
    AB[2, i] = 6 * EI / dz**4 + k_h[i]  # 主对角
    AB[3, i-1] = -4 * EI / dz**4  # 下一对角
    AB[4, i-2] = EI / dz**4  # 下二对角

# 边界条件
# 桩顶 (z=0)：剪力 = H_load，弯矩 = M_load
# 简化为：d³y/dz³ = H_load/EI, d²y/dz² = M_load/EI
AB[2, 0] = 1
AB[2, 1] = 0
AB[1, 2] = 0
AB[2, 2] = 1

# 桩底 (z=L)：固定端（简化）或自由端
# 采用固定端：y = 0, dy/dz = 0
AB[2, n-1] = 1
AB[3, n-2] = 0
AB[4, n-3] = 0
AB[2, n-2] = 1

# 右端项
b = np.zeros(n)
b[0] = H_load * dz**3 / EI  # 剪力边界条件（简化）
b[1] = M_load * dz**2 / EI  # 弯矩边界条件（简化）

# 求解
y_pile = solve_banded((2, 2), AB, b)

# 计算弯矩 M = -EI * d²y/dz²
M_pile = np.zeros(n)
for i in range(1, n-1):
    M_pile[i] = -EI * (y_pile[i-1] - 2*y_pile[i] + y_pile[i+1]) / dz**2

# 计算剪力 V = dM/dz
V_pile = np.zeros(n)
for i in range(1, n-1):
    V_pile[i] = (M_pile[i+1] - M_pile[i-1]) / (2*dz)
V_pile[0] = H_load  # 桩顶剪力

# 计算土抗力 p = k_h * y
p_soil = k_h * y_pile

# 计算最大弯矩和位移
M_max = np.max(np.abs(M_pile))
y_max = np.max(np.abs(y_pile))
M_max_idx = np.argmax(np.abs(M_pile))
z_M_max = z[M_max_idx]

print(f"桩身最大水平位移: {y_max*1000:.2f} mm")
print(f"最大弯矩: {M_max/1000:.1f} kN·m")
print(f"最大弯矩位置: 深度 {z_M_max:.2f} m")

# p-y曲线分析（特定深度）
z_py = 3.0  # 分析深度
soil_py = get_soil_properties(z_py)

# Matlock p-y曲线（黏土）
if soil_py['c_u'] > 0:
    epsilon_50 = 0.02  # 典型值
    y_50 = 2.5 * epsilon_50 * d_pile
    p_u = 9 * soil_py['c_u'] * d_pile  # 极限土抗力（简化）
    
    y_range = np.linspace(0, 10*y_50, 100)
    p_py = np.zeros_like(y_range)
    for i, yi in enumerate(y_range):
        if yi < 8 * y_50:
            p_py[i] = 0.5 * p_u * (yi / y_50)**0.33
        else:
            p_py[i] = p_u

# 可视化
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 图1：桩身水平位移
ax1 = axes[0, 0]
ax1.plot(y_pile*1000, z, 'b-', linewidth=2)
ax1.fill_betweenx(z, y_pile*1000, alpha=0.3, color='blue')
ax1.axvline(x=0, color='k', linewidth=0.5)
# 标记土层界面
for ld in layer_depths[1:-1]:
    ax1.axhline(y=ld, color='gray', linestyle='--', alpha=0.5)
ax1.set_xlabel('水平位移 (mm)', fontsize=11)
ax1.set_ylabel('深度 z (m)', fontsize=11)
ax1.set_title('桩身水平位移分布', fontsize=12)
ax1.invert_yaxis()
ax1.grid(True, alpha=0.3)

# 图2：桩身弯矩
ax2 = axes[0, 1]
ax2.plot(M_pile/1000, z, 'r-', linewidth=2)
ax2.fill_betweenx(z, M_pile/1000, alpha=0.3, color='red')
ax2.axvline(x=0, color='k', linewidth=0.5)
for ld in layer_depths[1:-1]:
    ax2.axhline(y=ld, color='gray', linestyle='--', alpha=0.5)
ax2.set_xlabel('弯矩 (kN·m)', fontsize=11)
ax2.set_ylabel('深度 z (m)', fontsize=11)
ax2.set_title('桩身弯矩分布', fontsize=12)
ax2.invert_yaxis()
ax2.grid(True, alpha=0.3)

# 图3：桩身剪力
ax3 = axes[1, 0]
ax3.plot(V_pile/1000, z, 'g-', linewidth=2)
ax3.fill_betweenx(z, V_pile/1000, alpha=0.3, color='green')
ax3.axvline(x=0, color='k', linewidth=0.5)
for ld in layer_depths[1:-1]:
    ax3.axhline(y=ld, color='gray', linestyle='--', alpha=0.5)
ax3.set_xlabel('剪力 (kN)', fontsize=11)
ax3.set_ylabel('深度 z (m)', fontsize=11)
ax3.set_title('桩身剪力分布', fontsize=12)
ax3.invert_yaxis()
ax3.grid(True, alpha=0.3)

# 图4：p-y曲线
ax4 = axes[1, 1]
if soil_py['c_u'] > 0:
    ax4.plot(y_range*1000, p_py/1000, 'purple', linewidth=2, label='p-y曲线')
    ax4.axhline(y=p_u/1000, color='r', linestyle='--', label=f'极限抗力 p_u={p_u/1000:.1f} kN/m')
    ax4.axvline(x=y_50*1000, color='orange', linestyle=':', label=f'y_50={y_50*1000:.1f} mm')
    ax4.set_xlabel('位移 y (mm)', fontsize=11)
    ax4.set_ylabel('土抗力 p (kN/m)', fontsize=11)
    ax4.set_title(f'p-y曲线 (深度 {z_py}m)', fontsize=12)
    ax4.legend(loc='lower right')
    ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('case5_pile_soil.png', dpi=150, bbox_inches='tight')
print("已保存 case5_pile_soil.png")
plt.close()

print("\n案例五分析完成！")

案例六：地震作用下土-结构动力相互作用

问题描述

某高层建筑基础为筏板基础，尺寸 $30m\times 30m$，埋深 $3m$。地基为砂土，地下水位 $2m$。分析地震作用下：

1. 场地地震反应
2. 基础地震响应
3. 液化可能性评估
4. 地震沉降估算

理论分析

场地地震反应分析：

采用一维波动理论分析场地地震反应：

$$\rho \frac{{\partial }^{2}u}{\partial t^2}=G\frac{{\partial }^{2}u}{\partial z^2}+\eta \frac{{\partial }^{3}u}{\partial z^2\partial t}$$

其中：

• $\rho$ 为土体密度
• $G$ 为剪切模量
• $\eta$ 为黏滞阻尼系数

液化判别（简化Seed方法）：

$$CRR=a{\left(\frac{N_1}{b}\right)}^c$$

其中，$N_1$ 为修正标贯击数，$a,b,c$ 为经验参数。

Python仿真代码

"""
案例六：地震作用下土-结构动力相互作用
场地地震反应与液化分析
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
import matplotlib
matplotlib.use('Agg')

plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

# 场地参数
H_site = 30.0  # 场地深度 m
z_wt = 2.0  # 地下水位 m

# 土体参数（分层）
# 0-10m：松散砂土
# 10-20m：中密砂土
# 20-30m：密实砂土

layer_depths = [0, 10, 20, 30]
layer_N = [10, 20, 35]  # 标贯击数
layer_gamma = [18e3, 19e3, 20e3]  # 重度 N/m³
layer_G = [50e6, 100e6, 200e6]  # 剪切模量 Pa

# 地震参数
M_earthquake = 7.5  # 震级
 PGA = 0.3  # 峰值地面加速度 g
a_max = PGA * 9.81  # m/s²

# 输入地震波（简化的正弦波）
t = np.linspace(0, 20, 2000)  # 20秒
dt = t[1] - t[0]
f_eq = 2.0  # 地震波频率 Hz
omega = 2 * np.pi * f_eq

# 包络函数
envelope = np.ones_like(t)
envelope[t < 2] = t[t < 2] / 2
envelope[t > 15] = np.exp(-(t[t > 15] - 15) / 2)

# 基岩输入加速度
a_rock = a_max * np.sin(omega * t) * envelope

# 一维场地反应分析（简化）
# 使用等效线性化方法

nz = 50
dz = H_site / nz
z = np.linspace(0, H_site, nz+1)

# 获取每层土的性质
def get_layer_property(z_pos, property_array):
    for i in range(len(layer_depths)-1):
        if layer_depths[i] <= z_pos <= layer_depths[i+1]:
            return property_array[i]
    return property_array[-1]

# 计算每层的剪切模量和阻尼比
G_profile = np.array([get_layer_property(zi, layer_G) for zi in z])
N_profile = np.array([get_layer_property(zi, layer_N) for zi in z])
gamma_profile = np.array([get_layer_property(zi, layer_gamma) for zi in z])

# 计算剪切波速
rho_profile = gamma_profile / 9.81
Vs_profile = np.sqrt(G_profile / rho_profile)

# 特征周期
T_site = 4 * H_site / np.mean(Vs_profile)
print(f"场地特征周期 T = {T_site:.3f} s")

# 简化的一维反应分析
# 假设加速度随深度线性衰减（简化模型）
amplification = np.zeros((len(t), nz+1))
for i, zi in enumerate(z):
    # 简化的放大系数
    if zi <= z_wt:
        amp_factor = 1.0 + 0.5 * (1 - zi/z_wt)
    else:
        amp_factor = 1.5 + 0.5 * (zi - z_wt) / (H_site - z_wt)
    amplification[:, i] = a_rock * amp_factor

# 计算每层的剪应变和剪应力
gamma_max = np.zeros(nz+1)
tau_max = np.zeros(nz+1)

for i, zi in enumerate(z):
    a_max_z = np.max(np.abs(amplification[:, i]))
    # 简化的剪应变计算
    gamma_max[i] = a_max_z / (2 * np.pi * f_eq)**2 / dz if i < nz else 0
    tau_max[i] = G_profile[i] * gamma_max[i]

# 液化判别（简化Seed方法）
# CSR = 0.65 * (a_max/g) * (sigma_v / sigma_v') * rd
# CRR 根据N1_60确定

CSR = np.zeros(nz+1)  # 循环应力比
CRR = np.zeros(nz+1)  # 循环抗力比

for i, zi in enumerate(z):
    if zi >= z_wt:  # 地下水位以下才可能液化
        sigma_v = np.sum(gamma_profile[:i+1]) * dz if i > 0 else gamma_profile[0] * dz
        sigma_v_eff = sigma_v - gamma_w * (zi - z_wt) if zi > z_wt else sigma_v
        
        # 深度修正系数 rd
        rd = 1.0 - 0.00765 * zi if zi < 9.15 else 1.174 - 0.0267 * zi
        
        # 循环应力比
        a_max_surface = np.max(np.abs(amplification[:, -1]))
        CSR[i] = 0.65 * (a_max_surface / 9.81) * (sigma_v / sigma_v_eff) * rd
        
        # 修正标贯击数 N1_60
        N1_60 = N_profile[i] * np.sqrt(100 / sigma_v_eff * 1000) if sigma_v_eff > 0 else N_profile[i]
        
        # 循环抗力比（简化公式）
        if N1_60 <= 0:
            CRR[i] = 0
        else:
            CRR[i] = 0.05 + 0.015 * N1_60
    else:
        CSR[i] = 0
        CRR[i] = 1.0  # 地下水位以上不液化

# 安全系数
F_liquefaction = CRR / (CSR + 1e-10)
liquefaction_potential = F_liquefaction < 1.0

# 液化深度
liquefaction_depths = z[liquefaction_potential]
if len(liquefaction_depths) > 0:
    print(f"\n液化可能发生深度范围: {liquefaction_depths[0]:.1f}m - {liquefaction_depths[-1]:.1f}m")
else:
    print("\n无明显液化可能性")

# 地震沉降估算（简化）
# 基于体积剪应变的方法
volumetric_strain = np.zeros(nz+1)
for i, zi in enumerate(z):
    if zi >= z_wt and F_liquefaction[i] < 1.5:  # 接近液化的区域
        # 简化的体积应变估算
        Dr = N_profile[i] / 50  # 相对密实度估算
        gamma_vol = 0.01 * (1 - Dr) * max(0, 1 - F_liquefaction[i])
        volumetric_strain[i] = gamma_vol

settlement_seismic = np.sum(volumetric_strain) * dz * 1000  # mm
print(f"估算地震沉降: {settlement_seismic:.1f} mm")

# 可视化
fig, axes = plt.subplots(2, 3, figsize=(16, 10))

# 图1：输入地震波
ax1 = axes[0, 0]
ax1.plot(t, a_rock/9.81, 'b-', linewidth=1)
ax1.set_xlabel('时间 (s)', fontsize=11)
ax1.set_ylabel('加速度 (g)', fontsize=11)
ax1.set_title('输入地震波', fontsize=12)
ax1.grid(True, alpha=0.3)
ax1.set_xlim(0, 20)

# 图2：加速度反应谱（简化）
ax2 = axes[0, 1]
T_range = np.linspace(0.1, 5, 100)
Sa = np.zeros_like(T_range)
for i, Ti in enumerate(T_range):
    # 简化的反应谱
    if Ti < T_site:
        Sa[i] = PGA * (1 + 2 * (T_site - Ti) / T_site)
    else:
        Sa[i] = PGA * (T_site / Ti)**0.5
ax2.plot(T_range, Sa, 'r-', linewidth=2)
ax2.axvline(x=T_site, color='g', linestyle='--', label=f'场地周期 T={T_site:.2f}s')
ax2.set_xlabel('周期 T (s)', fontsize=11)
ax2.set_ylabel('谱加速度 Sa (g)', fontsize=11)
ax2.set_title('加速度反应谱', fontsize=12)
ax2.legend(loc='upper right')
ax2.grid(True, alpha=0.3)

# 图3：土层剖面与性质
ax3 = axes[0, 2]
ax3_twin = ax3.twinx()
ax3.fill_betweenx(z, N_profile, alpha=0.3, color='brown', label='标贯击数')
ax3_twin.plot(Vs_profile, z, 'b-', linewidth=2, label='剪切波速')
ax3.axhline(y=z_wt, color='cyan', linestyle='--', label='地下水位')
ax3.set_xlabel('标贯击数 N', fontsize=11)
ax3_twin.set_xlabel('剪切波速 Vs (m/s)', fontsize=11, color='blue')
ax3.set_ylabel('深度 (m)', fontsize=11)
ax3.set_title('土层剖面', fontsize=12)
ax3.invert_yaxis()
ax3.legend(loc='upper right')
ax3_twin.legend(loc='lower right')

# 图4：液化判别
ax4 = axes[1, 0]
ax4.plot(CSR, z, 'r-', linewidth=2, label='CSR')
ax4.plot(CRR, z, 'b-', linewidth=2, label='CRR')
ax4.fill_betweenx(z, CSR, CRR, where=(CSR > CRR), alpha=0.3, color='red', label='液化区')
ax4.axhline(y=z_wt, color='cyan', linestyle='--', alpha=0.5)
ax4.set_xlabel('应力比', fontsize=11)
ax4.set_ylabel('深度 (m)', fontsize=11)
ax4.set_title('液化判别 (CSR vs CRR)', fontsize=12)
ax4.legend(loc='upper right')
ax4.invert_yaxis()
ax4.grid(True, alpha=0.3)

# 图5：液化安全系数
ax5 = axes[1, 1]
colors = ['red' if f < 1 else 'green' for f in F_liquefaction]
ax5.barh(z, F_liquefaction, color=colors, alpha=0.7, height=0.5)
ax5.axvline(x=1.0, color='black', linestyle='--', linewidth=2, label='F=1.0')
ax5.axhline(y=z_wt, color='cyan', linestyle='--', alpha=0.5)
ax5.set_xlabel('安全系数 F', fontsize=11)
ax5.set_ylabel('深度 (m)', fontsize=11)
ax5.set_title('液化安全系数分布', fontsize=12)
ax5.legend(loc='upper right')
ax5.invert_yaxis()
ax5.grid(True, alpha=0.3)

# 图6：时程分析（特定深度）
ax6 = axes[1, 2]
z_plot = 5.0  # 显示5m深度的加速度时程
idx_z = int(z_plot / dz)
ax6.plot(t, amplification[:, idx_z]/9.81, 'purple', linewidth=1, label=f'深度 {z_plot}m')
ax6.plot(t, a_rock/9.81, 'gray', linewidth=1, linestyle='--', label='基岩输入')
ax6.set_xlabel('时间 (s)', fontsize=11)
ax6.set_ylabel('加速度 (g)', fontsize=11)
ax6.set_title('加速度时程对比', fontsize=12)
ax6.legend(loc='upper right')
ax6.grid(True, alpha=0.3)
ax6.set_xlim(0, 20)

plt.tight_layout()
plt.savefig('case6_seismic_ssi.png', dpi=150, bbox_inches='tight')
print("已保存 case6_seismic_ssi.png")
plt.close()

print("\n案例六分析完成！")

---

Python仿真实现

7.1 完整仿真代码

以下是整合所有案例的完整Python仿真代码：

"""
主题037：土壤-结构相互作用 - 完整仿真代码
整合所有案例的Python仿真程序
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from scipy.optimize import fsolve
from scipy.sparse import csr_matrix
from scipy.sparse.linalg import spsolve
from scipy.linalg import solve_banded
from scipy.integrate import odeint
import matplotlib
matplotlib.use('Agg')

plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

print("="*60)
print("主题037：土壤-结构相互作用仿真分析")
print("="*60)

# ============================================================================
# 案例一：一维固结分析
# ============================================================================
print("\n" + "="*60)
print("案例一：一维固结分析")
print("="*60)

# 土体参数
H = 10.0  # 土层厚度 m
k = 1e-8  # 渗透系数 m/s
a_v = 0.5e-6  # 压缩系数 Pa^-1
e_0 = 1.0  # 初始孔隙比
gamma_w = 10e3  # 水的重度 N/m³
sigma_z = 100e3  # 附加应力 Pa

# 计算参数
H_dr = H / 2  # 排水路径长度（双面排水）
c_v = k * (1 + e_0) / (a_v * gamma_w)  # 固结系数
print(f"固结系数 c_v = {c_v:.2e} m²/s")

# 最终沉降量
S_inf = a_v / (1 + e_0) * sigma_z * H
print(f"最终沉降量 S_∞ = {S_inf*1000:.1f} mm")

# 空间离散
nz = 50
dz = H / nz
z = np.linspace(0, H, nz+1)

# 时间参数
T_v_target = 2.0  # 最大时间因子
t_max = T_v_target * H_dr**2 / c_v
t = np.linspace(0, t_max, 100)

# 解析解：孔压分布
def pore_pressure_analytical(z, t, c_v, H, sigma_z):
    """Terzaghi一维固结解析解"""
    nz = len(z)
    nt = len(t)
    u = np.zeros((nt, nz))
    
    for it, ti in enumerate(t):
        T_v = c_v * ti / (H/2)**2
        for iz, zi in enumerate(z):
            Z = zi / H
            u_sum = 0
            for m in range(20):
                M = np.pi / 2 * (2*m + 1)
                u_sum += (2 / M) * np.sin(M * (1 - 2*Z)) * np.exp(-M**2 * T_v)
            u[it, iz] = sigma_z * u_sum
    
    return u

# 计算孔压
u = pore_pressure_analytical(z, t, c_v, H, sigma_z)

# 计算固结度
def consolidation_degree(T_v):
    if T_v < 0.2:
        return np.sqrt(4 * T_v / np.pi)
    else:
        return 1 - (8 / np.pi**2) * np.exp(-np.pi**2 * T_v / 4)

# 计算各时刻沉降量
S_t = np.zeros(len(t))
for i, ti in enumerate(t):
    T_v = c_v * ti / H_dr**2
    U_v = consolidation_degree(T_v)
    S_t[i] = U_v * S_inf

# 计算固结度90%所需时间
def equation(T_v):
    return consolidation_degree(T_v) - 0.9

T_v_90 = fsolve(equation, 1.0)[0]
t_90 = T_v_90 * H_dr**2 / c_v
print(f"固结度90%所需时间 t_90 = {t_90/86400:.1f} days")

# 可视化
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

ax1 = axes[0, 0]
time_indices = [0, 10, 25, 50, 99]
colors = plt.cm.viridis(np.linspace(0, 1, len(time_indices)))
for idx, color in zip(time_indices, colors):
    T_v = c_v * t[idx] / H_dr**2
    ax1.plot(u[idx, :]/1000, z, color=color, linewidth=2, label=f'T_v = {T_v:.3f}')
ax1.set_xlabel('孔压 u (kPa)', fontsize=11)
ax1.set_ylabel('深度 z (m)', fontsize=11)
ax1.set_title('一维固结：孔压分布随时间变化', fontsize=12)
ax1.legend(loc='upper right')
ax1.grid(True, alpha=0.3)
ax1.invert_yaxis()

ax2 = axes[0, 1]
T_v_array = c_v * t / H_dr**2
ax2.plot(T_v_array, S_t * 1000, 'b-', linewidth=2, label='沉降量')
ax2.axhline(y=S_inf*1000, color='r', linestyle='--', label=f'最终沉降 = {S_inf*1000:.1f} mm')
ax2.axvline(x=T_v_90, color='g', linestyle=':', label=f'U=90% at T_v={T_v_90:.3f}')
ax2.set_xlabel('时间因子 T_v', fontsize=11)
ax2.set_ylabel('沉降量 S (mm)', fontsize=11)
ax2.set_title('沉降-时间关系曲线', fontsize=12)
ax2.legend(loc='lower right')
ax2.grid(True, alpha=0.3)

ax3 = axes[1, 0]
T_v_range = np.linspace(0, 2, 100)
U_v_range = [consolidation_degree(Tv) for Tv in T_v_range]
ax3.plot(T_v_range, U_v_range, 'b-', linewidth=2)
ax3.axhline(y=0.9, color='r', linestyle='--', label='U = 90%')
ax3.axvline(x=T_v_90, color='g', linestyle=':', label=f'T_v = {T_v_90:.3f}')
ax3.set_xlabel('时间因子 T_v', fontsize=11)
ax3.set_ylabel('固结度 U_v', fontsize=11)
ax3.set_title('固结度-时间因子关系', fontsize=12)
ax3.legend(loc='lower right')
ax3.grid(True, alpha=0.3)

ax4 = axes[1, 1]
T_v_plot = c_v * t / H_dr**2
Z_plot = z / H
T_v_mesh, Z_mesh = np.meshgrid(T_v_plot, Z_plot)
contour = ax4.contourf(T_v_mesh, Z_mesh, u.T/1000, levels=20, cmap='RdYlBu_r')
plt.colorbar(contour, ax=ax4, label='孔压 u (kPa)')
ax4.set_xlabel('时间因子 T_v', fontsize=11)
ax4.set_ylabel('相对深度 z/H', fontsize=11)
ax4.set_title('孔压等时线分布', fontsize=12)
ax4.invert_yaxis()

plt.tight_layout()
plt.savefig('case1_consolidation.png', dpi=150, bbox_inches='tight')
print("已保存 case1_consolidation.png")
plt.close()

# 创建动画
fig, ax = plt.subplots(figsize=(10, 6))
line, = ax.plot([], [], 'b-', linewidth=2)
ax.set_xlim(0, sigma_z/1000 * 1.1)
ax.set_xlabel('孔压 u (kPa)', fontsize=11)
ax.set_ylabel('深度 z (m)', fontsize=11)
ax.set_title('一维固结过程动画', fontsize=12)
ax.grid(True, alpha=0.3)
ax.invert_yaxis()
time_text = ax.text(0.02, 0.95, '', transform=ax.transAxes, fontsize=10)

frame_indices = np.linspace(0, len(t)-1, 50, dtype=int)

def init():
    line.set_data([], [])
    time_text.set_text('')
    return line, time_text

def update(frame):
    idx = frame_indices[frame]
    line.set_data(u[idx, :]/1000, z)
    T_v = c_v * t[idx] / H_dr**2
    U_v = consolidation_degree(T_v)
    time_text.set_text(f'T_v = {T_v:.3f}, U = {U_v*100:.1f}%')
    return line, time_text

anim = FuncAnimation(fig, update, init_func=init, frames=len(frame_indices),
                     interval=200, blit=True)
anim.save('case1_consolidation.gif', writer='pillow', fps=5, dpi=100)
print("已保存 case1_consolidation.gif")
plt.close()

print("案例一分析完成！")

# ============================================================================
# 其他案例的代码将在后续补充...
# ============================================================================

print("\n" + "="*60)
print("所有案例仿真完成！")
print("="*60)

---