第002篇：连续介质力学基础

摘要

连续介质力学是研究连续介质（固体、液体、气体）在外力作用下变形和运动规律的学科，是多物理场耦合仿真的理论基础。本教程系统介绍了连续介质力学的基本概念、数学描述和物理定律。首先阐述了连续介质假设的内涵及其适用范围，详细讲解了描述变形的拉格朗日方法和欧拉方法。其次深入分析了应力张量和应变张量的定义、性质和计算方法，包括柯西应力、皮奥拉-基尔霍夫应力等概念。然后推导了质量守恒、动量守恒和能量守恒三大守恒定律的数学形式，建立了连续介质力学的控制方程体系。最后介绍了本构关系的基本理论，包括弹性、粘弹性和塑性等材料模型的数学描述。通过Python数值实验，演示了应力分析、变形计算和守恒定律验证等关键内容，帮助读者建立对连续介质力学理论的直观认识。

关键词

连续介质力学；应力张量；应变张量；守恒定律；本构关系；拉格朗日描述；欧拉描述；有限变形

---

1. 引言

1.1 连续介质假设

连续介质假设是连续介质力学的基本出发点。该假设认为物质在空间中连续分布，忽略物质的微观离散结构（原子、分子），将物质视为充满空间的连续体。

连续介质假设的数学表达：

对于连续介质中的任意一点$\mathbf{x}$和任意小的邻域$V$，都存在定义良好的物理量（密度$\rho$、速度$\mathbf{v}$、温度$T$等）：

$$\rho (\mathbf{x},t)=\underset{V\to 0}{\lim }\frac{m(V)}{V}$$

其中$m(V)$是体积$V$内的质量。

连续介质假设的适用条件：

1. 宏观尺度：研究对象的特征尺度远大于分子平均自由程
2. 致密介质：介质足够致密，微观空隙可以忽略
3. 平衡态：系统处于或接近热力学平衡态

典型适用范围：

• 固体和液体的大部分工程问题
• 气体在常压下的流动（除极稀薄气体）
• 特征尺度$>10^{-6}$ m的问题

不适用的情况：

• 稀薄气体动力学（真空环境、高空大气）
• 纳米尺度问题（需用分子动力学）
• 激波内部（厚度与分子自由程相当）

1.2 连续介质力学的重要性

连续介质力学是多物理场耦合仿真的理论基础：

（1）统一的数学框架

连续介质力学为不同物理场提供了统一的描述语言。无论是固体变形、流体流动还是传热传质，都可以用相同的数学框架描述：

• 场的概念（标量场、矢量场、张量场）
• 守恒定律（质量、动量、能量）
• 本构关系（材料特性）

（2）多物理场耦合的桥梁

不同物理场的耦合往往通过连续介质力学的概念实现：

• 热-结构耦合：热应力、热膨胀
• 流-固耦合：界面应力、动量交换
• 电磁-力学耦合：电磁力、磁致伸缩

（3）数值方法的理论基础

有限元法、有限体积法等数值方法都建立在连续介质力学的基础上：

• 变分原理和弱形式
• 控制方程的离散
• 边界条件的处理

1.3 本教程的学习目标

通过本教程的学习，读者将能够：

1. 理解基本概念：掌握连续介质、变形、应力、应变等核心概念
2. 掌握数学工具：熟练运用张量分析描述连续介质力学问题
3. 推导控制方程：从守恒定律出发推导质量、动量、能量方程
4. 建立本构关系：理解不同材料模型的数学形式和物理意义
5. 数值实现能力：用Python实现基本的连续介质力学计算

---

2. 运动的描述

2.1 物质点与空间点

在连续介质力学中，需要区分两种描述方式：

物质点（Material Point）：

• 标识为连续介质中的特定质点
• 用物质坐标$\mathbf{X}$标记（参考构型中的位置）
• 随时间运动，追踪同一质点的历史

空间点（Spatial Point）：

• 标识为空间中的固定位置
• 用空间坐标$\mathbf{x}$标记（当前构型中的位置）
• 不同时间可能有不同的物质点经过

构型（Configuration）：

• 参考构型（Reference Configuration）：$t=t_0$时刻的构型，物质坐标$\mathbf{X}$
• 当前构型（Current Configuration）：$t$时刻的构型，空间坐标$\mathbf{x}$

2.2 拉格朗日描述（物质描述）

拉格朗日描述以物质点为研究对象，追踪每个物质点的运动历史。

运动方程：

$$\mathbf{x}=\mathbf{\chi }(\mathbf{X},t)$$

其中$\mathbf{\chi }$是运动函数，描述物质点$\mathbf{X}$在时刻$t$的位置$\mathbf{x}$。

位移矢量：

$$\mathbf{u}(\mathbf{X},t)=\mathbf{x}-\mathbf{X}=\mathbf{\chi }(\mathbf{X},t)-\mathbf{X}$$

速度（物质导数）：

$$\mathbf{V}(\mathbf{X},t)=\frac{\partial \mathbf{\chi }(\mathbf{X},t)}{\partial t}$$

加速度：

$$\mathbf{A}(\mathbf{X},t)=\frac{{\partial }^2\mathbf{\chi }(\mathbf{X},t)}{\partial t^2}$$

拉格朗日描述的特点：

• 网格随物质运动（固体力学常用）
• 易于追踪材料界面和自由表面
• 质量守恒自动满足
• 处理大变形时网格可能畸变

2.3 欧拉描述（空间描述）

欧拉描述以空间点为研究对象，关注物理量在场中的空间分布。

速度场：

$$\mathbf{v}(\mathbf{x},t)=\mathbf{V}({\mathbf{\chi }}^{-1}(\mathbf{x},t),t)$$

物质导数（随体导数）：

物理量$\varphi$的物质导数表示跟随物质点的时间变化率：

$$\frac{D\varphi }{Dt}=\frac{\partial \varphi }{\partial t}+(\mathbf{v}\cdot \nabla )\varphi$$

其中：

• $ \frac{\partial \phi}{\partial t} $：局部导数（空间点上的时间变化）
• $(\mathbf{v}\cdot \nabla )\varphi$：对流导数（物质运动引起的空间变化）

欧拉描述的特点：

• 网格固定在空间（流体力学常用）
• 易于处理大变形和流动问题
• 需要处理对流项
• 追踪物质界面需要特殊方法（VOF、Level Set等）

2.4 两种描述的关系

变形梯度张量（Deformation Gradient）：

$$\mathbf{F}=\frac{\partial \mathbf{x}}{\partial \mathbf{X}}={\nabla }_{\mathbf{X}}\mathbf{\chi }$$

变形梯度张量$\mathbf{F}$将参考构型中的线元$\mathrm{d}\mathbf{X}$映射到当前构型中的线元$\mathrm{d}\mathbf{x}$：

$$\mathrm{d}\mathbf{x}=\mathbf{F}\cdot \mathrm{d}\mathbf{X}$$

雅可比行列式：

$$J=\det (\mathbf{F})$$

$J$表示体积变化率：

• $J>1$：体积膨胀
• $J=1$：体积不变（不可压缩）
• $J<1$：体积压缩

速度梯度张量：

$$\mathbf{L}={\nabla }_{\mathbf{x}}\mathbf{v}=\frac{\partial \mathbf{v}}{\partial \mathbf{x}}=\dot{\mathbf{F}}\cdot {\mathbf{F}}^{-1}$$

速度梯度可以分解为对称部分和反对称部分：

$$\mathbf{L}=\mathbf{D}+\mathbf{W}$$

其中：

• 变形率张量（对称）：$\mathbf{D}=\frac{1}{2}(\mathbf{L}+{\mathbf{L}}^T)$
• 旋率张量（反对称）：$\mathbf{W}=\frac{1}{2}(\mathbf{L}-{\mathbf{L}}^T)$

---

3. 应变分析

3.1 变形与应变的概念

变形是指物体形状和体积的改变。描述变形需要：

1. 消除刚体运动（平移和转动）
2. 提取纯变形部分

应变是描述变形的度量，应该满足：

• 刚体运动时应变为零
• 能够区分拉伸、剪切等变形模式
• 与应力通过本构关系联系

3.2 格林-拉格朗日应变张量

格林-拉格朗日应变张量（Green-Lagrange Strain）是参考构型中定义的应变度量：

$$\mathbf{E}=\frac{1}{2}({\mathbf{F}}^T\cdot \mathbf{F}-\mathbf{I})$$

物理意义：

考虑参考构型中的线元$\mathrm{d}\mathbf{X}$，变形后长度为$\mathrm{d}s$，变形前长度为$\mathrm{d}S$：

$$(\mathrm{d}s{)}^2-(\mathrm{d}S{)}^2=2\mathrm{d}{\mathbf{X}}^T\cdot \mathbf{E}\cdot \mathrm{d}\mathbf{X}$$

小变形近似：

当位移梯度很小时（$|\nabla \mathbf{u}|\ll 1$）：

$$\mathbf{E}\approx \mathbf{\epsilon }=\frac{1}{2}(\nabla \mathbf{u}+(\nabla \mathbf{u}{)}^T)$$

这就是小应变张量或柯西应变张量。

3.3 小应变张量的分量

在直角坐标系中，小应变张量的分量为：

正应变（线应变）：

$${\epsilon }_{xx}=\frac{\partial u_x}{\partial x},{\epsilon }_{yy}=\frac{\partial u_y}{\partial y},{\epsilon }_{zz}=\frac{\partial u_z}{\partial z}$$

剪应变：

$${\gamma }_{xy}=2{\epsilon }_{xy}=\frac{\partial u_x}{\partial y}+\frac{\partial u_y}{\partial x}$$

$${\gamma }_{yz}=2{\epsilon }_{yz}=\frac{\partial u_y}{\partial z}+\frac{\partial u_z}{\partial y}$$

$${\gamma }_{zx}=2{\epsilon }_{zx}=\frac{\partial u_z}{\partial x}+\frac{\partial u_x}{\partial z}$$

应变张量的矩阵形式：

$$\mathbf{\epsilon }=\left[\begin{array}{ccc}{\epsilon }_{xx}& {\epsilon }_{xy}& {\epsilon }_{xz}\\ {\epsilon }_{yx}& {\epsilon }_{yy}& {\epsilon }_{yz}\\ {\epsilon }_{zx}& {\epsilon }_{zy}& {\epsilon }_{zz}\end{array}\right]$$

3.4 主应变与应变不变量

主应变：应变张量的特征值，表示主方向上的正应变。

求解特征方程：

$$\det (\mathbf{\epsilon }-\epsilon \mathbf{I})=0$$

展开得：

$$-{\epsilon }^3+I_1{\epsilon }^2-I_2\epsilon +I_3=0$$

其中$I_1,I_2,I_3$是应变不变量：

$$I_1={\epsilon }_{xx}+{\epsilon }_{yy}+{\epsilon }_{zz}=\text{tr}(\mathbf{\epsilon })$$

$$I_2={\epsilon }_{xx}{\epsilon }_{yy}+{\epsilon }_{yy}{\epsilon }_{zz}+{\epsilon }_{zz}{\epsilon }_{xx}-{\epsilon }_{xy}^2-{\epsilon }_{yz}^2-{\epsilon }_{zx}^2$$

$$I_3=\det (\mathbf{\epsilon })$$

物理意义：

• $I_1$：体积应变（${\epsilon }_v=I_1$）
• $I_2$：与剪切变形有关
• $I_3$：与体积变化和剪切都有关

3.5 偏应变与球应变

应变张量可以分解为球应变（体积变化）和偏应变（形状变化）：

$$\mathbf{\epsilon }=\frac{1}{3}{\epsilon }_v\mathbf{I}+\mathbf{e}$$

其中：

• 球应变：$\frac{1}{3}{\epsilon }_v\mathbf{I}$，表示各向同性体积变化
• 偏应变：$\mathbf{e}=\mathbf{\epsilon }-\frac{1}{3}{\epsilon }_v\mathbf{I}$，表示形状改变（剪切变形）

等效应变（von Mises应变）：

$$\bar{\epsilon }=\sqrt{\frac{2}{3}\mathbf{e}:\mathbf{e}}=\sqrt{\frac{2}{3}{e}_{ij}{e}_{ij}}$$

等效应变在塑性力学中非常重要，用于判断材料是否屈服。

---

4. 应力分析

4.1 应力的概念

应力是单位面积上的内力，描述物体内部的受力状态。

柯西应力原理：物体内部任意截面上的应力矢量${\mathbf{t}}^{(\mathbf{n})}$只依赖于该截面的法向$\mathbf{n}$：

$${\mathbf{t}}^{(\mathbf{n})}={\mathbf{t}}^{(\mathbf{n})}(\mathbf{x},t,\mathbf{n})$$

4.2 柯西应力张量

柯西应力张量$\mathbf{\sigma }$将法向量映射为应力矢量：

$${\mathbf{t}}^{(\mathbf{n})}=\mathbf{\sigma }\cdot \mathbf{n}$$

应力张量的对称性：

由角动量守恒可得：

$$\mathbf{\sigma }={\mathbf{\sigma }}^T$$

即柯西应力张量是对称的，只有6个独立分量。

应力张量的分量：

$$\mathbf{\sigma }=\left[\begin{array}{ccc}{\sigma }_{xx}& {\tau }_{xy}& {\tau }_{xz}\\ {\tau }_{yx}& {\sigma }_{yy}& {\tau }_{yz}\\ {\tau }_{zx}& {\tau }_{zy}& {\sigma }_{zz}\end{array}\right]$$

其中：

• $ \sigma_{xx}, \sigma_{yy}, \sigma_{zz} $：正应力（拉为正，压为负）
• $ \tau_{xy}, \tau_{yz}, \tau_{zx} $：剪应力

4.3 主应力与应力不变量

主应力：应力张量的特征值，表示主平面上的正应力。

求解特征方程：

$$\det (\mathbf{\sigma }-\sigma \mathbf{I})=0$$

展开得：

$$-{\sigma }^3+I_1^{\sigma }{\sigma }^2-I_2^{\sigma }\sigma +I_3^{\sigma }=0$$

其中应力不变量：

$$I_1^{\sigma }={\sigma }_{xx}+{\sigma }_{yy}+{\sigma }_{zz}=\text{tr}(\mathbf{\sigma })$$

$$I_2^{\sigma }={\sigma }_{xx}{\sigma }_{yy}+{\sigma }_{yy}{\sigma }_{zz}+{\sigma }_{zz}{\sigma }_{xx}-{\tau }_{xy}^2-{\tau }_{yz}^2-{\tau }_{zx}^2$$

$$I_3^{\sigma }=\det (\mathbf{\sigma })$$

平均应力（静水压力）：

$$p=\frac{1}{3}{I}_1^{\sigma }=\frac{1}{3}({\sigma }_{xx}+{\sigma }_{yy}+{\sigma }_{zz})$$

4.4 偏应力与球应力

应力张量可以分解为球应力（静水压力）和偏应力（剪切应力）：

$$\mathbf{\sigma }=p\mathbf{I}+\mathbf{s}$$

其中：

• 球应力：$p\mathbf{I}$，表示各向同性压力
• 偏应力：$\mathbf{s}=\mathbf{\sigma }-p\mathbf{I}$，表示剪切应力状态

von Mises等效应力：

$$\bar{\sigma }=\sqrt{\frac{3}{2}\mathbf{s}:\mathbf{s}}=\sqrt{\frac{3}{2}{s}_{ij}{s}_{ij}}$$

等效应力是判断材料屈服的重要指标。当$\bar{\sigma }\ge {\sigma }_s$（屈服强度）时，材料发生塑性变形。

4.5 其他应力张量

第一类皮奥拉-基尔霍夫应力（名义应力）：

$$\mathbf{P}=J\mathbf{\sigma }\cdot {\mathbf{F}}^{-T}$$

第二类皮奥拉-基尔霍夫应力：

$$\mathbf{S}=J{\mathbf{F}}^{-1}\cdot \mathbf{\sigma }\cdot {\mathbf{F}}^{-T}$$

这些应力张量在有限变形问题中非常重要，特别是在固体力学和流-固耦合分析中。

---

5. 守恒定律

5.1 质量守恒（连续性方程）

积分形式：

$$\frac{\mathrm{d}}{\mathrm{d}t}{\int }_V\rho \mathrm{d}V=0$$

或

$${\int }_V\frac{\partial \rho }{\partial t}\mathrm{d}V+{\int }_S\rho \mathbf{v}\cdot \mathbf{n}\mathrm{d}S=0$$

微分形式（欧拉描述）：

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

或

$$\frac{D\rho }{Dt}+\rho \nabla \cdot \mathbf{v}=0$$

不可压缩流动（$\rho =\text{const}$）：

$$\nabla \cdot \mathbf{v}=0$$

即速度场的散度为零。

5.2 动量守恒（运动方程）

积分形式（线性动量）：

$$\frac{\mathrm{d}}{\mathrm{d}t}{\int }_V\rho \mathbf{v}\mathrm{d}V={\int }_V\rho \mathbf{f}\mathrm{d}V+{\int }_S{\mathbf{t}}^{(\mathbf{n})}\mathrm{d}S$$

其中$\mathbf{f}$是体积力（如重力）。

微分形式（柯西运动方程）：

$$\rho \frac{D\mathbf{v}}{Dt}=\nabla \cdot \mathbf{\sigma }+\rho \mathbf{f}$$

分量形式：

$$\rho \frac{D{v}_i}{Dt}=\frac{\partial {\sigma }_{ji}}{\partial x_j}+\rho f_i$$

纳维-斯托克斯方程（牛顿流体）：

将本构关系$\mathbf{\sigma }=-p\mathbf{I}+\mu (\nabla \mathbf{v}+(\nabla \mathbf{v}{)}^T)$代入运动方程：

$$\rho \frac{D\mathbf{v}}{Dt}=-\nabla p+\mu {\nabla }^2\mathbf{v}+\rho \mathbf{f}$$

5.3 能量守恒

总能量方程：

$$\rho \frac{D}{Dt}\left(e+\frac{1}{2}\mathbf{v}\cdot \mathbf{v}\right)=\rho \mathbf{f}\cdot \mathbf{v}+\nabla \cdot (\mathbf{\sigma }\cdot \mathbf{v})-\nabla \cdot \mathbf{q}+\rho Q$$

其中：

• $e$：比内能
• $\mathbf{q}$：热流矢量
• $Q$：体积热源

内能方程：

$$\rho \frac{De}{Dt}=\mathbf{\sigma }:\mathbf{D}-\nabla \cdot \mathbf{q}+\rho Q$$

其中$\mathbf{\sigma }:\mathbf{D}={\sigma }_{ij}{D}_{ij}$是应力功率（变形功）。

热传导方程（固体）：

对于固体（$\mathbf{D}\approx \dot{\mathbf{\epsilon }}$），假设$\mathbf{\sigma }:\dot{\mathbf{\epsilon }}\approx 0$（机械功转化为热的部分很小）：

$$\rho c\frac{\partial T}{\partial t}=\nabla \cdot (k\nabla T)+Q$$

5.4 守恒定律的总结

控制方程组（一般形式）：

$$\{\begin{array}{ll}\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho \mathbf{v})=0& \text{（质量守恒）}\\ \rho \frac{D\mathbf{v}}{Dt}=\nabla \cdot \mathbf{\sigma }+\rho \mathbf{f}& \text{（动量守恒）}\\ \rho \frac{De}{Dt}=\mathbf{\sigma }:\mathbf{D}-\nabla \cdot \mathbf{q}+\rho Q& \text{（能量守恒）}\end{array}$$

这是描述连续介质运动的最一般形式的方程组，适用于固体、液体和气体。

---

6. 本构关系

6.1 本构关系的概念

本构关系（Constitutive Relation）描述了材料的力学特性，建立了应力与应变（或应变率）之间的关系。

建立本构关系的基本原则：

1. 确定性原理：应力只依赖于变形的当前历史和物质点本身
2. 局部作用原理：某点的应力只依赖于该点邻域的变形
3. 客观性原理：本构关系应与观察者的运动无关（坐标变换不变性）
4. 热力学第二定律：耗散必须非负

6.2 线弹性本构关系

广义胡克定律：

$$\mathbf{\sigma }=\mathbf{C}:\mathbf{\epsilon }$$

或分量形式：

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

其中$\mathbf{C}$是弹性张量，有81个分量，但由于对称性，独立分量减少。

各向同性线弹性（最常用的模型）：

对于各向同性材料，弹性张量只有两个独立常数（拉梅常数$\lambda$和$\mu$，或杨氏模量$E$和泊松比$\nu$）：

$$\mathbf{\sigma }=\lambda \text{tr}(\mathbf{\epsilon })\mathbf{I}+2\mu \mathbf{\epsilon }$$

或

$$\mathbf{\sigma }=\frac{E}{1+\nu }\left(\mathbf{\epsilon }+\frac{\nu }{1-2\nu }\text{tr}(\mathbf{\epsilon })\mathbf{I}\right)$$

常用弹性常数关系：

常数	表达式
杨氏模量$E$	$\frac{\mu (3\lambda +2\mu )}{\lambda +\mu }$
泊松比$\nu$	$\frac{\lambda }{2(\lambda +\mu )}$
剪切模量$G$	$\mu$
体积模量$K$	$\lambda +\frac{2}{3}\mu$

6.3 热弹性本构关系

考虑温度变化引起的热应变：

$${\mathbf{\epsilon }}^{th}=\alpha (T-T_0)\mathbf{I}$$

其中$\alpha$是热膨胀系数。

热弹性本构方程：

$$\mathbf{\sigma }=\mathbf{C}:(\mathbf{\epsilon }-{\mathbf{\epsilon }}^{th})=\mathbf{C}:\mathbf{\epsilon }-\mathbf{C}:{\mathbf{\epsilon }}^{th}$$

对于各向同性材料：

$$\mathbf{\sigma }=\lambda \text{tr}(\mathbf{\epsilon })\mathbf{I}+2\mu \mathbf{\epsilon }-(3\lambda +2\mu )\alpha (T-T_0)\mathbf{I}$$

6.4 粘弹性本构关系

粘弹性材料同时具有弹性固体和粘性流体的特性，应力依赖于应变历史。

麦克斯韦模型（串联）：

$$\dot{\epsilon }=\frac{\dot{\sigma }}{E}+\frac{\sigma }{\eta }$$

开尔文-沃伊特模型（并联）：

$$\sigma =E\epsilon +\eta \dot{\epsilon }$$

标准线性固体（三参数模型）：

$$\sigma +{\tau }_{\epsilon }\dot{\sigma }=E_R(\epsilon +{\tau }_{\sigma }\dot{\epsilon })$$

其中${\tau }_{\epsilon }$和${\tau }_{\sigma }$是松弛时间和推迟时间。

6.5 塑性本构关系

屈服准则：

von Mises屈服准则（金属）：

$$\bar{\sigma }=\sqrt{\frac{3}{2}\mathbf{s}:\mathbf{s}}\le {\sigma }_s$$

Tresca屈服准则（最大剪应力）：

$$\max |{\sigma }_i-{\sigma }_j|\le {\sigma }_s$$

流动法则（关联流动）：

$$\mathrm{d}{\mathbf{\epsilon }}^p=\mathrm{d}\lambda \frac{\partial f}{\partial \mathbf{\sigma }}$$

其中$f$是屈服函数，$\mathrm{d}\lambda$是塑性乘子。

硬化模型：

• 各向同性硬化：屈服面均匀膨胀
• 随动硬化：屈服面平移（包辛格效应）
• 混合硬化：两者结合

6.6 流体的本构关系

牛顿流体：

应力与应变率成线性关系：

$$\mathbf{\sigma }=-p\mathbf{I}+\mu (\nabla \mathbf{v}+(\nabla \mathbf{v}{)}^T)$$

或

$$\mathbf{\sigma }=-p\mathbf{I}+2\mu \mathbf{D}$$

非牛顿流体：

• 幂律流体：$\tau =K(\dot{\gamma }{)}^n$
• 宾汉流体：存在屈服应力
• 粘弹性流体：具有记忆效应

---

7. Python数值实验

以下是用于演示连续介质力学基本概念的Python代码，包括应变分析、应力计算和守恒定律验证。

# -*- coding: utf-8 -*-
"""
主题002：连续介质力学基础
Python演示代码：应力应变分析与守恒定律
"""

import matplotlib
matplotlib.use('Agg')  # 使用非交互式后端
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import FancyArrowPatch, Rectangle
import warnings
warnings.filterwarnings('ignore')
import os

# 创建输出目录
output_dir = r'd:\文档\500仿真领域\工程仿真\多物理场耦合仿真\主题002'
os.makedirs(output_dir, exist_ok=True)

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

print("=" * 60)
print("主题002：连续介质力学基础")
print("Python演示：应力应变分析")
print("=" * 60)

# ==================== 1. 应变分析演示 ====================
print("\n[1] 应变张量分析...")

# 定义应变张量（小应变假设）
epsilon = np.array([
    [0.001, 0.0005, 0],
    [0.0005, 0.002, 0],
    [0, 0, -0.0005]
])

print("应变张量:")
print(epsilon)

# 计算应变不变量
I1 = np.trace(epsilon)
I2 = (np.trace(epsilon)**2 - np.trace(epsilon @ epsilon)) / 2
I3 = np.linalg.det(epsilon)

print(f"\n应变不变量:")
print(f"  I1 (体积应变) = {I1:.6f}")
print(f"  I2 = {I2:.6e}")
print(f"  I3 = {I3:.6e}")

# 计算主应变
eigenvalues, eigenvectors = np.linalg.eig(epsilon)
principal_strains = np.sort(eigenvalues)[::-1]
print(f"\n主应变:")
for i, eps in enumerate(principal_strains):
    print(f"  ε_{i+1} = {eps:.6f}")

# 分解球应变和偏应变
volumetric_strain = I1
deviatoric_strain = epsilon - (volumetric_strain / 3) * np.eye(3)

print(f"\n应变分解:")
print(f"  体积应变 = {volumetric_strain:.6f}")
print("  偏应变张量:")
print(deviatoric_strain)

# 等效应变（von Mises）
effective_strain = np.sqrt(2/3 * np.sum(deviatoric_strain**2))
print(f"  等效应变 = {effective_strain:.6f}")

# ==================== 2. 应力分析演示 ====================
print("\n[2] 应力张量分析...")

# 定义应力张量（单位：MPa）
sigma = np.array([
    [100, 30, 0],
    [30, 50, 0],
    [0, 0, 20]
])

print("应力张量 (MPa):")
print(sigma)

# 计算应力不变量
I1_sigma = np.trace(sigma)
I2_sigma = (np.trace(sigma)**2 - np.trace(sigma @ sigma)) / 2
I3_sigma = np.linalg.det(sigma)

print(f"\n应力不变量:")
print(f"  I1 = {I1_sigma:.2f} MPa")
print(f"  I2 = {I2_sigma:.2f} MPa²")
print(f"  I3 = {I3_sigma:.2f} MPa³")

# 平均应力（静水压力）
mean_stress = I1_sigma / 3
print(f"\n平均应力 = {mean_stress:.2f} MPa")

# 计算主应力
eigenvalues_s, eigenvectors_s = np.linalg.eig(sigma)
principal_stresses = np.sort(eigenvalues_s)[::-1]
print(f"\n主应力:")
for i, sig in enumerate(principal_stresses):
    print(f"  σ_{i+1} = {sig:.2f} MPa")

# 分解球应力和偏应力
hydrostatic_stress = mean_stress * np.eye(3)
deviatoric_stress = sigma - hydrostatic_stress

print(f"\n应力分解:")
print("  球应力 (静水压力):")
print(hydrostatic_stress)
print("  偏应力:")
print(deviatoric_stress)

# 等效应力（von Mises）
effective_stress = np.sqrt(3/2 * np.sum(deviatoric_stress**2))
print(f"\n等效应力 (von Mises) = {effective_stress:.2f} MPa")

# 最大剪应力
tau_max = (principal_stresses[0] - principal_stresses[2]) / 2
print(f"最大剪应力 = {tau_max:.2f} MPa")

# ==================== 3. 可视化 ====================
print("\n[3] 生成可视化结果...")

# 图1：应变张量可视化
fig, axes = plt.subplots(2, 3, figsize=(15, 10))

# 应变张量分量
ax1 = axes[0, 0]
im1 = ax1.imshow(epsilon, cmap='RdBu_r', vmin=-0.002, vmax=0.002)
ax1.set_title('Strain Tensor Components', fontsize=12, fontweight='bold')
ax1.set_xticks([0, 1, 2])
ax1.set_yticks([0, 1, 2])
ax1.set_xticklabels(['x', 'y', 'z'])
ax1.set_yticklabels(['x', 'y', 'z'])
for i in range(3):
    for j in range(3):
        ax1.text(j, i, f'{epsilon[i,j]:.4f}', ha='center', va='center', fontsize=10)
plt.colorbar(im1, ax=ax1)

# 主应变
ax2 = axes[0, 1]
bars = ax2.bar(['ε₁', 'ε₂', 'ε₃'], principal_strains, color=['red', 'green', 'blue'], alpha=0.7)
ax2.set_title('Principal Strains', fontsize=12, fontweight='bold')
ax2.set_ylabel('Strain')
ax2.axhline(y=0, color='k', linestyle='--', alpha=0.5)
for bar, val in zip(bars, principal_strains):
    ax2.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.00005, 
             f'{val:.4f}', ha='center', va='bottom', fontsize=10)

# 应变分解
ax3 = axes[0, 2]
categories = ['Volumetric', 'Deviatoric\n(Equivalent)']
values = [volumetric_strain, effective_strain]
bars = ax3.bar(categories, values, color=['orange', 'purple'], alpha=0.7)
ax3.set_title('Strain Decomposition', fontsize=12, fontweight='bold')
ax3.set_ylabel('Strain')
for bar, val in zip(bars, values):
    ax3.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.00005, 
             f'{val:.4f}', ha='center', va='bottom', fontsize=10)

# 应力张量分量
ax4 = axes[1, 0]
im4 = ax4.imshow(sigma, cmap='RdBu_r', vmin=-50, vmax=150)
ax4.set_title('Stress Tensor Components (MPa)', fontsize=12, fontweight='bold')
ax4.set_xticks([0, 1, 2])
ax4.set_yticks([0, 1, 2])
ax4.set_xticklabels(['x', 'y', 'z'])
ax4.set_yticklabels(['x', 'y', 'z'])
for i in range(3):
    for j in range(3):
        ax4.text(j, i, f'{sigma[i,j]:.1f}', ha='center', va='center', fontsize=10)
plt.colorbar(im4, ax=ax4)

# 主应力
ax5 = axes[1, 1]
bars = ax5.bar(['σ₁', 'σ₂', 'σ₃'], principal_stresses, color=['red', 'green', 'blue'], alpha=0.7)
ax5.set_title('Principal Stresses (MPa)', fontsize=12, fontweight='bold')
ax5.set_ylabel('Stress (MPa)')
ax5.axhline(y=0, color='k', linestyle='--', alpha=0.5)
for bar, val in zip(bars, principal_stresses):
    ax5.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 2, 
             f'{val:.1f}', ha='center', va='bottom', fontsize=10)

# 应力分解
ax6 = axes[1, 2]
categories = ['Hydrostatic', 'von Mises\nEffective']
values = [mean_stress, effective_stress]
bars = ax6.bar(categories, values, color=['cyan', 'magenta'], alpha=0.7)
ax6.set_title('Stress Decomposition (MPa)', fontsize=12, fontweight='bold')
ax6.set_ylabel('Stress (MPa)')
for bar, val in zip(bars, values):
    ax6.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 2, 
             f'{val:.1f}', ha='center', va='bottom', fontsize=10)

plt.tight_layout()
plt.savefig(f'{output_dir}/stress_strain_analysis.png', dpi=150, bbox_inches='tight')
plt.close()
print("  ✓ 图1已保存: stress_strain_analysis.png")

# ==================== 4. 莫尔圆演示 ====================
print("\n[4] 生成莫尔圆...")

fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# 应力莫尔圆
ax1 = axes[0]
theta = np.linspace(0, 2*np.pi, 100)

# 对于平面应力状态（sigma_zz = tau_xz = tau_yz = 0）
sigma_x, sigma_y, tau_xy = sigma[0,0], sigma[1,1], sigma[0,1]

# 圆心
sigma_avg = (sigma_x + sigma_y) / 2
# 半径
R = np.sqrt(((sigma_x - sigma_y)/2)**2 + tau_xy**2)

# 绘制莫尔圆
circle_x = sigma_avg + R * np.cos(theta)
circle_y = R * np.sin(theta)
ax1.plot(circle_x, circle_y, 'b-', linewidth=2, label='Mohr Circle')
ax1.axhline(y=0, color='k', linewidth=0.5)
ax1.axvline(x=0, color='k', linewidth=0.5)

# 标记主应力
sigma_p1 = sigma_avg + R
sigma_p2 = sigma_avg - R
ax1.plot([sigma_p1, sigma_p2], [0, 0], 'ro', markersize=10, label='Principal Stresses')
ax1.text(sigma_p1, 2, f'σ₁={sigma_p1:.1f}', ha='center', fontsize=10)
ax1.text(sigma_p2, -4, f'σ₂={sigma_p2:.1f}', ha='center', fontsize=10)

# 标记当前应力状态
ax1.plot(sigma_x, tau_xy, 'gs', markersize=10, label='Current State')
ax1.plot(sigma_y, -tau_xy, 'gs', markersize=10)

ax1.set_xlabel('Normal Stress (MPa)', fontsize=11)
ax1.set_ylabel('Shear Stress (MPa)', fontsize=11)
ax1.set_title('Mohr Circle for Stress', fontsize=12, fontweight='bold')
ax1.legend()
ax1.grid(True, alpha=0.3)
ax1.set_aspect('equal')

# 应变莫尔圆
ax2 = axes[1]
epsilon_x, epsilon_y, gamma_xy = epsilon[0,0], epsilon[1,1], 2*epsilon[0,1]

epsilon_avg = (epsilon_x + epsilon_y) / 2
R_epsilon = np.sqrt(((epsilon_x - epsilon_y)/2)**2 + (gamma_xy/2)**2)

circle_x_e = epsilon_avg + R_epsilon * np.cos(theta)
circle_y_e = R_epsilon * np.sin(theta)
ax2.plot(circle_x_e, circle_y_e, 'r-', linewidth=2, label='Mohr Circle')
ax2.axhline(y=0, color='k', linewidth=0.5)
ax2.axvline(x=0, color='k', linewidth=0.5)

epsilon_p1 = epsilon_avg + R_epsilon
epsilon_p2 = epsilon_avg - R_epsilon
ax2.plot([epsilon_p1, epsilon_p2], [0, 0], 'bo', markersize=10, label='Principal Strains')
ax2.text(epsilon_p1, 0.0002, f'ε₁={epsilon_p1:.4f}', ha='center', fontsize=10)
ax2.text(epsilon_p2, -0.0004, f'ε₂={epsilon_p2:.4f}', ha='center', fontsize=10)

ax2.plot(epsilon_x, gamma_xy/2, 'gs', markersize=10, label='Current State')
ax2.plot(epsilon_y, -gamma_xy/2, 'gs', markersize=10)

ax2.set_xlabel('Normal Strain', fontsize=11)
ax2.set_ylabel('Shear Strain/2', fontsize=11)
ax2.set_title('Mohr Circle for Strain', fontsize=12, fontweight='bold')
ax2.legend()
ax2.grid(True, alpha=0.3)
ax2.set_aspect('equal')

plt.tight_layout()
plt.savefig(f'{output_dir}/mohr_circles.png', dpi=150, bbox_inches='tight')
plt.close()
print("  ✓ 图2已保存: mohr_circles.png")

# ==================== 5. 弹性本构关系演示 ====================
print("\n[5] 弹性本构关系验证...")

# 材料参数（钢）
E = 200e9  # 杨氏模量 (Pa)
nu = 0.3   # 泊松比

# 计算拉梅常数
lambda_lame = E * nu / ((1 + nu) * (1 - 2*nu))
mu = E / (2 * (1 + nu))

print(f"材料参数:")
print(f"  杨氏模量 E = {E/1e9} GPa")
print(f"  泊松比 ν = {nu}")
print(f"  拉梅常数 λ = {lambda_lame/1e9:.2f} GPa")
print(f"  剪切模量 μ = {mu/1e9:.2f} GPa")

# 广义胡克定律：σ = C:ε
# 对于各向同性材料
sigma_calculated = lambda_lame * np.trace(epsilon) * np.eye(3) + 2 * mu * epsilon

print(f"\n根据胡克定律计算的应力 (Pa):")
print(sigma_calculated)
print(f"\n等效为 MPa:")
print(sigma_calculated / 1e6)

# ==================== 6. 守恒定律演示 ====================
print("\n[6] 守恒定律数值验证...")

# 创建一个简单的速度场（二维）
x = np.linspace(0, 1, 20)
y = np.linspace(0, 1, 20)
X, Y = np.meshgrid(x, y)

# 定义速度场：u = x, v = -y（不可压缩）
U = X
V = -Y

# 计算散度（应该为零）
dU_dx = np.gradient(U, x, axis=1)
dV_dy = np.gradient(V, y, axis=0)
divergence = dU_dx + dV_dy

print(f"速度场散度统计:")
print(f"  最大值: {np.max(divergence):.6e}")
print(f"  最小值: {np.min(divergence):.6e}")
print(f"  平均值: {np.mean(divergence):.6e}")
print(f"  标准差: {np.std(divergence):.6e}")

# 图3：速度场和散度
fig, axes = plt.subplots(1, 3, figsize=(15, 5))

# 速度矢量图
ax1 = axes[0]
skip = 2
ax1.quiver(X[::skip, ::skip], Y[::skip, ::skip], 
           U[::skip, ::skip], V[::skip, ::skip], 
           scale=10, color='blue', alpha=0.7)
ax1.set_xlabel('x', fontsize=11)
ax1.set_ylabel('y', fontsize=11)
ax1.set_title('Velocity Field (u=x, v=-y)', fontsize=12, fontweight='bold')
ax1.set_aspect('equal')
ax1.grid(True, alpha=0.3)

# 散度分布
ax2 = axes[1]
im = ax2.contourf(X, Y, divergence, levels=20, cmap='RdBu_r')
ax2.set_xlabel('x', fontsize=11)
ax2.set_ylabel('y', fontsize=11)
ax2.set_title('Divergence Field (∇·v)', fontsize=12, fontweight='bold')
plt.colorbar(im, ax=ax2)

# 流线图
ax3 = axes[2]
ax3.streamplot(X, Y, U, V, color=np.sqrt(U**2 + V**2), cmap='viridis', density=1.5)
ax3.set_xlabel('x', fontsize=11)
ax3.set_ylabel('y', fontsize=11)
ax3.set_title('Streamlines', fontsize=12, fontweight='bold')
ax3.set_aspect('equal')
ax3.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/conservation_laws.png', dpi=150, bbox_inches='tight')
plt.close()
print("  ✓ 图3已保存: conservation_laws.png")

# ==================== 7. 变形梯度演示 ====================
print("\n[7] 变形梯度分析...")

# 定义一个简单的变形映射
# X = (X, Y) -> x = (1.1*X + 0.1*Y, 0.05*X + 0.9*Y)
F = np.array([
    [1.1, 0.1, 0],
    [0.05, 0.9, 0],
    [0, 0, 1.0]
])

print("变形梯度张量 F:")
print(F)

# 右柯西-格林张量
C = F.T @ F
print("\n右柯西-格林张量 C = F^T·F:")
print(C)

# 格林-拉格朗日应变
E_green = 0.5 * (C - np.eye(3))
print("\n格林-拉格朗日应变 E:")
print(E_green)

# 雅可比行列式（体积变化）
J = np.linalg.det(F)
print(f"\n雅可比行列式 J = {J:.4f}")
print(f"体积变化 = {(J-1)*100:.2f}%")

# 主拉伸比
eigenvalues_C, _ = np.linalg.eig(C)
lambda_stretch = np.sqrt(eigenvalues_C)
print(f"\n主拉伸比:")
for i, lam in enumerate(lambda_stretch):
    print(f"  λ_{i+1} = {lam:.4f}")

# ==================== 8. 材料模型对比 ====================
print("\n[8] 材料模型对比...")

# 应变范围
strains = np.linspace(0, 0.01, 100)

# 线弹性
E_mod = 200e9
stress_elastic = E_mod * strains

# 理想塑性
sigma_y = 250e6
stress_perfect_plastic = np.minimum(E_mod * strains, sigma_y)

# 线性硬化
H = 0.1 * E_mod
stress_hardening = np.where(E_mod * strains <= sigma_y, 
                            E_mod * strains, 
                            sigma_y + H * (strains - sigma_y/E_mod))

# 幂律硬化
n = 0.3
K = sigma_y / (0.002)**n
stress_power = K * strains**n

fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(strains*100, stress_elastic/1e6, 'b-', linewidth=2, label='Linear Elastic')
ax.plot(strains*100, stress_perfect_plastic/1e6, 'r--', linewidth=2, label='Perfect Plastic')
ax.plot(strains*100, stress_hardening/1e6, 'g-.', linewidth=2, label='Linear Hardening')
ax.plot(strains*100, stress_power/1e6, 'm:', linewidth=2, label=f'Power Law (n={n})')

ax.axhline(y=sigma_y/1e6, color='k', linestyle='--', alpha=0.5, label='Yield Stress')
ax.set_xlabel('Strain (%)', fontsize=12)
ax.set_ylabel('Stress (MPa)', fontsize=12)
ax.set_title('Material Models Comparison', fontsize=14, fontweight='bold')
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_xlim(0, 1)
ax.set_ylim(0, 600)

plt.tight_layout()
plt.savefig(f'{output_dir}/material_models.png', dpi=150, bbox_inches='tight')
plt.close()
print("  ✓ 图4已保存: material_models.png")

# ==================== 9. 总结输出 ====================
print("\n" + "=" * 60)
print("仿真完成！")
print("=" * 60)
print(f"\n生成的文件:")
print(f"  1. {output_dir}\\stress_strain_analysis.png")
print(f"  2. {output_dir}\\mohr_circles.png")
print(f"  3. {output_dir}\\conservation_laws.png")
print(f"  4. {output_dir}\\material_models.png")
print("\n连续介质力学基础演示结束。")