主题097：软材料与大变形力学

1. 引言

软材料与大变形力学是固体力学的一个重要分支，专门研究能够承受大变形（通常超过10%甚至可达数百倍）的软物质材料的力学行为。这类材料在现代工程技术和生物医学领域具有广泛的应用前景，包括橡胶密封件、生物软组织、软体机器人、柔性电子器件等。

与传统工程材料（如金属、混凝土）不同，软材料具有独特的力学特性：超弹性、粘弹性、各向异性、不可压缩性等。这些特性使得软材料在小应力作用下即可产生大变形，同时能够储存和释放大量弹性能量。理解软材料的大变形力学行为对于设计高性能软物质结构和器件至关重要。

本教程将系统介绍软材料与大变形力学的理论基础、本构模型、数值方法和工程应用。通过理论讲解与Python仿真相结合的方式，帮助读者深入理解软材料的力学特性，掌握大变形问题的分析方法。

2. 软材料的分类与特性

2.1 软材料的主要类型

软材料是一类模量较低（通常在kPa到MPa量级）、能够承受大变形的材料。根据材料组成和结构特征，软材料可以分为以下几类：

（1）弹性体材料

弹性体是最典型的软材料，包括天然橡胶、合成橡胶（如丁苯橡胶、硅橡胶）、热塑性弹性体等。弹性体由长链高分子交联网络构成，具有显著的超弹性特征。典型弹性体的弹性模量在0.1-10 MPa范围内，断裂延伸率可达500%以上。

弹性体的力学行为主要源于高分子链的熵弹性。在未变形状态下，高分子链呈卷曲构象，具有较高的构象熵。当材料受拉伸时，高分子链沿拉伸方向取向，构象熵降低，产生恢复力。这种熵弹性机制使得弹性体在小应力下即可产生大变形，且变形可完全恢复。

（2）水凝胶材料

水凝胶是由亲水性高分子网络和水组成的软材料。高分子网络通过化学交联或物理交联形成三维结构，能够吸收大量水分而保持固体形态。水凝胶的含水量可达90%以上，弹性模量通常在kPa到MPa量级。

水凝胶的力学行为受多种因素影响，包括交联密度、含水量、离子浓度等。由于含有大量水分，水凝胶具有与生物组织相似的特性，在生物医学领域有广泛应用，如人工软骨、药物缓释载体、伤口敷料等。

（3）生物软组织

生物软组织包括皮肤、肌肉、血管、肌腱、韧带等。这些组织由细胞、细胞外基质（胶原蛋白、弹性蛋白、蛋白聚糖等）和水组成，具有复杂的层次结构和各向异性特征。

生物软组织的力学行为具有显著的非线性、粘弹性和各向异性。例如，动脉血管壁在轴向和周向具有不同的力学性能，且应力-应变关系呈明显的"J形"曲线特征——在低应变区刚度较低，随着应变增加刚度急剧上升。这种特性有利于在生理压力范围内保持血管的顺应性，同时在过高压力下提供保护。

（4）软体机器人材料

软体机器人是一类模仿生物软体（如章鱼、大象鼻子）的新型机器人，其本体由软材料构成，能够实现连续变形和自适应交互。软体机器人常用的材料包括硅胶（如Ecoflex、Dragon Skin）、气动人工肌肉、介电弹性体等。

软体机器人材料需要具备良好的可变形性、耐久性和驱动响应特性。例如，介电弹性体在电场作用下可产生高达300%的应变，是制作人工肌肉的理想材料。气动人工肌肉则利用气体压力驱动，能够产生较大的输出力和行程。

2.2 软材料的基本力学特性

（1）超弹性（Hyperelasticity）

超弹性是软材料最显著的力学特征。超弹性材料在加载-卸载过程中表现出非线性的应力-应变关系，但变形是完全可逆的，不耗散能量。这与金属材料的线弹性行为形成鲜明对比。

超弹性材料的本构关系通常用应变能密度函数描述。对于各向同性超弹性材料，应变能密度是三个主应变不变量的函数。常见的超弹性本构模型包括Neo-Hookean模型、Mooney-Rivlin模型、Ogden模型等，将在后续章节详细介绍。

（2）粘弹性（Viscoelasticity）

许多软材料表现出时间相关的力学行为，即粘弹性。粘弹性材料的应力响应不仅取决于当前的应变状态，还与应变历史有关。典型的粘弹性现象包括蠕变（恒定应力下应变随时间增加）、应力松弛（恒定应变下应力随时间降低）和迟滞（加载-卸载曲线不重合形成回线）。

粘弹性的微观机制包括高分子链的缠结与解缠结、交联网络的重组、水分的迁移等。在动态载荷下，粘弹性导致能量耗散，表现为材料的阻尼特性。这对于减振、隔振应用具有重要意义。

（3）不可压缩性（Incompressibility）

大多数软材料（特别是弹性体和生物组织）在变形过程中体积几乎保持不变，表现出近似不可压缩的特性。橡胶材料的体积模量（约1-2 GPa）通常比剪切模量（约0.1-10 MPa）高2-3个数量级。

不可压缩性对大变形问题的数值求解带来挑战。在有限元分析中，需要采用特殊的单元 formulation（如混合单元、增强应变单元）或压力-位移耦合求解策略来处理不可压缩约束。

（4）各向异性（Anisotropy）

许多软材料具有各向异性特征，即力学性能随方向变化。生物软组织是典型的各向异性材料——血管壁含有周向排列的胶原纤维和弹性纤维，肌腱则主要由纵向排列的胶原纤维束组成。

各向异性通常源于材料内部的取向结构，如纤维取向、分子链取向等。描述各向异性超弹性材料需要引入结构张量或纤维方向向量，建立依赖于方向的应变能函数。

（5）损伤与失效

软材料在大变形下可能发生损伤和失效。损伤机制包括分子链断裂、交联键断裂、微裂纹萌生等。与金属材料的脆性断裂不同，软材料的断裂通常经历大范围的塑性变形和损伤演化过程。

软材料的断裂韧性通常用撕裂能（tearing energy）或断裂韧性（fracture toughness）表征。理解软材料的损伤与失效机制对于预测结构寿命和设计安全余量至关重要。

3. 大变形力学的理论基础

3.1 连续介质运动学

大变形问题的分析需要采用非线性连续介质力学框架。与线性小变形理论不同，大变形理论需要区分参考构型（未变形状态）和当前构型（变形后状态），并引入有限应变度量。

（1）变形梯度张量

考虑连续体中一个物质点，在参考构型中的位置向量为X，在当前构型中的位置向量为x。x与X之间的映射关系由运动函数描述：

x = φ(X, t)

变形梯度张量F定义为当前构型位置对参考构型位置的梯度：

F = ∂x/∂X = ∇φ

变形梯度张量F是一个二阶张量，包含了材料变形的全部信息。F可以分解为旋转部分和拉伸部分：

F = R·U = V·R

其中R是正交旋转张量，U和V分别是右和左拉伸张量（对称正定）。这种分解称为极分解，将变形分解为纯拉伸和纯旋转的叠加。

（2）应变度量

在大变形问题中，有多种应变度量可供选择。最常用的包括：

格林-拉格朗日应变张量（Green-Lagrange strain tensor）：

E = 1/2(F^T·F - I) = 1/2(C - I)

其中C = F^T·F称为右柯西-格林变形张量（Right Cauchy-Green tensor）。

阿尔曼西应变张量（Almansi strain tensor）：

e = 1/2(I - F^(-T)·**F**(-1)) = 1/2(I - b^(-1))

其中b = F·F^T称为左柯西-格林变形张量（Left Cauchy-Green tensor）。

对数应变（Hencky strain）或真实应变：

ε = ln(U) = 1/2 ln(C)

在小变形情况下，上述应变度量都退化为工程应变。但在有限变形下，它们给出不同的应变值。格林-拉格朗日应变适用于以参考构型为基准的描述（拉格朗日描述），阿尔曼西应变适用于以当前构型为基准的描述（欧拉描述）。

（3）速度梯度与变形率

速度梯度张量L定义为当前构型中速度对位置的梯度：

L = ∂v/∂x = Ḟ·F^(-1)

L可以分解为对称部分和反对称部分：

L = D + W

其中D = 1/2(L + L^T)是变形率张量（rate of deformation），W = 1/2(L - L^T)是旋率张量（spin tensor）。D描述了材料的瞬时拉伸变形率，W描述了瞬时刚体旋转。

3.2 应力度量

在大变形问题中，由于存在参考构型和当前构型两个不同的构型，应力的定义也有多种形式：

（1）柯西应力（Cauchy Stress）

柯西应力σ定义在当前构型上，是真实的物理应力，表示单位当前面积上的力。柯西应力张量是对称的，满足角动量平衡。

（2）第一皮奥拉-基尔霍夫应力（First Piola-Kirchhoff Stress）

第一皮奥拉-基尔霍夫应力P定义在参考构型上，表示单位参考面积上的力（但力的方向在当前构型）：

P = Jσ·F^(-T)

其中J = det(F)是体积比。P通常是非对称的。

（3）第二皮奥拉-基尔霍夫应力（Second Piola-Kirchhoff Stress）

第二皮奥拉-基尔霍夫应力S定义在参考构型上，表示单位参考面积上的力（力的方向也转换到参考构型）：

S = F^(-1)·P = JF^{(-1)·**σ**·**F**}(-T)

S是对称张量，在固体力学中广泛使用，因为它与格林-拉格朗日应变E功共轭。

（4）基尔霍夫应力（Kirchhoff Stress）

基尔霍夫应力τ定义为：

τ = Jσ

它与变形率张量D功共轭，在计算力学中经常使用。

3.3 平衡方程与虚功原理

（1）平衡方程

在当前构型中，柯西应力满足平衡方程：

∇·σ + ρf = 0

其中ρ是当前密度，f是单位质量的体积力。

在参考构型中，第一皮奥拉-基尔霍夫应力满足：

∇₀·P + ρ₀f = 0

其中∇₀表示对参考构型坐标的梯度，ρ₀是参考密度。

（2）虚功原理

大变形问题的虚功原理可以表述为：对于任意的虚位移δu，内虚功等于外虚功：

∫ S:δE dV₀ = ∫ ρ₀f·δu dV₀ + ∫ t₀·δu dA₀

其中t₀是参考构型上的面力，积分在参考构型上进行。这是非线性有限元分析的基础。

3.4 客观性原理

在大变形问题中，本构关系必须满足客观性原理（或称为框架无差异原理）。客观性要求材料的力学响应与观察者的运动状态无关，即本构关系在刚体旋转下保持不变。

客观性要求应变度量和应力度量必须采用适当的定义。格林-拉格朗日应变和第二皮奥拉-基尔霍夫应力自动满足客观性。柯西应力本身不是客观的，但可以通过适当的客观应力率（如Jaumann率、Truesdell率）来描述其演化。

4. 超弹性本构模型

4.1 超弹性理论框架

超弹性材料的本构关系由应变能密度函数（strain energy density function）完全确定。对于单位体积的应变能ψ，第二皮奥拉-基尔霍夫应力与格林-拉格朗日应变的关系为：

S = ∂ψ/∂E

或者等价地，柯西应力可以表示为：

σ = (1/J)F·(∂ψ/∂E)·F^T

对于各向同性超弹性材料，应变能可以表示为右柯西-格林张量C的三个不变量的函数：

ψ = ψ(I₁, I₂, I₃)

其中三个不变量定义为：

I₁ = tr(C) = λ₁² + λ₂² + λ₃²

I₂ = 1/2[(trC)² - tr(C²)] = λ₁²λ₂² + λ₂²λ₃² + λ₃²λ₁²

I₃ = det(C) = J² = λ₁²λ₂²λ₃²

λ₁, λ₂, λ₃是三个主拉伸比。

柯西应力可以表示为：

σ = (2/J)[(∂ψ/∂I₁ + I₁∂ψ/∂I₂)b - ∂ψ/∂I₂b² + I₃∂ψ/∂I₃I]

其中b = F·F^T是左柯西-格林张量。

4.2 可压缩超弹性模型

（1）Neo-Hookean模型

Neo-Hookean模型是最简单的超弹性模型，其应变能函数为：

ψ = C₁₀(I₁ - 3) + (1/D₁)(J - 1)²

其中C₁₀是材料常数，与剪切模量相关：μ = 2C₁₀；D₁与体积模量相关：K = 2/D₁。第二项是可压缩性修正项，对于完全不可压缩材料（J = 1），此项消失。

Neo-Hookean模型适用于小至中等变形（拉伸比<2），对于大变形预测精度较差。它假设分子链完全自由取向，忽略了链间相互作用。

（2）Mooney-Rivlin模型

Mooney-Rivlin模型是对Neo-Hookean模型的改进，其应变能函数为：

ψ = C₁₀(I₁ - 3) + C₀₁(I₂ - 3) + (1/D₁)(J - 1)²

其中C₁₀和C₀₁是材料常数。Mooney-Rivlin模型能够更好地描述橡胶在中等变形下的行为，适用于拉伸比<3的情况。

（3）Ogden模型

Ogden模型是一种更一般的超弹性模型，其应变能函数表示为主拉伸比的幂级数：

ψ = Σᵢ (μᵢ/αᵢ)(λ₁^αᵢ + λ₂^αᵢ + λ₃^αᵢ - 3) + (1/D₁)(J - 1)²

其中μᵢ和αᵢ是材料常数。通常取N=3或6项即可获得良好的拟合效果。Ogden模型能够准确描述橡胶在大变形下的复杂行为，包括"S形"应力-应变曲线。

（4）Arruda-Boyce模型（八链模型）

Arruda-Boyce模型基于高分子网络的八链模型，其应变能函数为：

ψ = μλₘ²[1/3(I₁/λₘ²) + 1/45(I₁/λₘ²)² + 33/2835(I₁/λₘ²)³ + …] + (1/D₁)(J - 1)²

其中μ是初始剪切模量，λₘ是锁定拉伸比（表示分子链完全伸直时的拉伸比）。Arruda-Boyce模型具有明确的物理意义，能够预测应变硬化行为。

4.3 不可压缩超弹性模型

对于完全不可压缩材料（J = 1），应变能函数可以简化为仅依赖于I₁和I₂的形式，或者采用主拉伸比表示（满足λ₁λ₂λ₃ = 1的约束）。

（1）不可压缩Neo-Hookean模型

ψ = C₁₀(I₁ - 3)

柯西应力为：

σ = -pI + 2C₁₀b

其中p是静水压力，由不可压缩约束确定。

（2）不可压缩Ogden模型

ψ = Σᵢ (μᵢ/αᵢ)(λ₁^αᵢ + λ₂^αᵢ + λ₃^αᵢ - 3), 约束条件：λ₁λ₂λ₃ = 1

主柯西应力为：

σᵢ = λᵢ(∂ψ/∂λᵢ) - p = Σⱼ μⱼλᵢ^αⱼ - p

4.4 各向异性超弹性模型

生物软组织和纤维增强弹性体通常表现出各向异性特征。描述各向异性超弹性材料需要引入结构参数，如纤维方向向量。

（1）Holzapfel-Gasser-Ogden（HGO）模型

HGO模型是描述动脉血管壁等生物组织的常用模型。假设材料中含有两族对称分布的胶原纤维，纤维方向与周向成±γ角。应变能函数为：

ψ = C₁₀(I₁ - 3) + (k₁/2k₂)Σᵢ₌₁²[exp(k₂(I₄ᵢ - 1)²) - 1]

其中I₄ᵢ = C😦aᵢ⊗aᵢ)表示沿纤维方向的变形，aᵢ是纤维方向单位向量，k₁和k₂是纤维相关的材料参数。

（2）广义结构张量模型

对于更一般的各向异性材料，可以引入结构张量H来描述材料的方向性特征：

ψ = ψ(I₁, I₂, I₃, C:H)

这种方法可以描述单轴各向异性、横观各向同性、正交各向异性等多种各向异性形式。

5. 粘弹性本构模型

5.1 线性粘弹性理论

（1）Maxwell模型

Maxwell模型由弹簧和粘壶串联组成，其本构关系为：

σ̇ + (E/η)σ = Eε̇

其中E是弹性模量，η是粘度。Maxwell模型能够描述应力松弛现象，但不能描述蠕变。

（2）Kelvin-Voigt模型

Kelvin-Voigt模型由弹簧和粘壶并联组成：

σ = Eε + ηε̇

Kelvin-Voigt模型能够描述蠕变现象，但不能描述应力松弛。

（3）广义Maxwell模型

广义Maxwell模型由多个Maxwell单元并联组成，能够同时描述蠕变、应力松弛和动态响应：

σ(t) = E₀ε(t) + Σᵢ Eᵢ∫₀ᵗ exp(-(t-τ)/τᵢ)ε̇(τ)dτ

其中E₀是平衡模量，Eᵢ和τᵢ是第i个Maxwell单元的模量和松弛时间。

5.2 准线性粘弹性模型

（1）Fung模型

Fung提出的准线性粘弹性模型广泛用于描述生物软组织的粘弹性行为：

σ(t) = ∫₀ᵗ G(t-τ)(∂σᵉ/∂τ)dτ

其中σᵉ是瞬时弹性响应，G(t)是归一化松弛函数，通常表示为：

G(t) = Σᵢ Cᵢ exp(-t/τᵢ), Σᵢ Cᵢ = 1

（2）超粘弹性模型

对于大变形问题，需要将粘弹性理论与超弹性理论结合。一种方法是将变形梯度进行乘法分解：

F = Fᵉ·Fᵛ

其中Fᵉ是弹性部分，Fᵛ是粘性部分。分别定义弹性应变度量和粘性应变演化方程，可以建立有限变形的粘弹性本构模型。

6. 大变形有限元方法

6.1 非线性有限元基础

大变形问题的有限元分析需要处理几何非线性（大位移、大应变）、材料非线性（非线性本构关系）和可能的接触非线性。总拉格朗日格式（Total Lagrangian formulation）和更新拉格朗日格式（Updated Lagrangian formulation）是两种常用的方法。

（1）总拉格朗日格式

在总拉格朗日格式中，所有变量都参照参考构型定义。虚功方程为：

∫ S:δE dV₀ = δWₑₓₜ

其中S是第二皮奥拉-基尔霍夫应力，E是格林-拉格朗日应变，δWₑₓₜ是外虚功。

将位移场用形函数插值：u = N·d，其中d是节点位移向量。格林-拉格朗日应变可以分解为线性部分和非线性部分：

E = Eᴸ + Eᴺᴸ

虚应变为：

δE = Bᴸ·δd + Bᴺᴸ·δd

其中Bᴸ是线性应变-位移矩阵，Bᴺᴸ是非线性应变-位移矩阵。

（2）切线刚度矩阵

非线性方程的求解需要计算切线刚度矩阵。总刚度矩阵由材料刚度矩阵和几何刚度矩阵组成：

Kₜ = Kᴹ + Kᴳ

材料刚度矩阵：Kᴹ = ∫ Bᴸᵀ·C·Bᴸ dV₀

几何刚度矩阵：Kᴳ = ∫ Gᵀ·S·G dV₀

其中C是材料切线模量张量，G是形函数梯度矩阵，S是第二皮奥拉-基尔霍夫应力矩阵。

6.2 不可压缩问题的处理

不可压缩约束（J = 1或ν = 0.5）给数值求解带来困难，因为传统的位移有限元会出现体积锁定（volumetric locking）现象。常用的解决方法包括：

（1）混合有限元方法

引入压力作为独立变量，建立位移-压力耦合的混合变分原理。Hu-Washizu变分原理或Hellinger-Reissner变分原理是理论基础。

（2）减缩积分

对体积项采用低阶积分（如体积应变采用单点积分），可以有效缓解体积锁定。但需要注意检查零能模式（hourglass mode）。

（3）增强应变方法

在单元内引入增强的应变场，改善单元的变形能力。Simo等人提出的增强应变四边形单元（Q1E4）是典型代表。

（4）F-bar方法

通过修正变形梯度的体积部分，保证单元的平均体积应变满足不可压缩约束，同时允许局部的体积变化。

6.3 求解算法

非线性方程组通常采用牛顿-拉夫逊迭代法求解。对于大变形问题，载荷需要分步施加，并可能采用弧长法（arc-length method）来追踪载荷-位移曲线的完整路径，包括软化段和失稳点。

7. 软材料结构的稳定性分析

7.1 屈曲与失稳现象

软材料结构在大变形下可能发生多种失稳现象：

（1）欧拉屈曲

细长杆件在轴向压力下发生整体屈曲，临界载荷由欧拉公式给出：

P_cr = π²EI/(KL)²

其中K是有效长度系数。对于大变形，屈曲后的载荷-位移关系是非线性的。

（2）表面起皱（Wrinkling）

软材料薄膜在压缩应力作用下表面形成褶皱。起皱的临界应变和波长取决于材料的弹性模量和几何尺寸。起皱现象在皮肤、薄膜涂层、生物膜等系统中广泛存在。

（3）局部化变形

某些软材料结构在拉伸过程中会出现颈缩（necking）或剪切带（shear banding）等局部化变形模式。这与材料的本构失稳（如应变软化）有关。

7.2 稳定性判据

结构的稳定性可以通过切线刚度矩阵的正定性来判断。当切线刚度矩阵出现零特征值时，结构达到临界状态。对于保守系统，可以使用能量判据：当总势能的二阶变分δ²Π > 0时，平衡状态是稳定的。

8. 软材料断裂力学

8.1 软材料断裂特性

与硬脆材料不同，软材料的断裂表现出以下特征：

（1）大变形断裂

软材料在断裂前经历大范围的塑性变形，裂纹尖端存在大的钝化区。线弹性断裂力学（LEFM）不再适用，需要采用弹塑性断裂力学或基于能量的方法。

（2）裂纹尖端张开

软材料裂纹尖端发生显著的张开，裂纹尖端张开位移（CTOD）可以作为断裂参量。Rivlin-Thomas理论使用撕裂能（tearing energy）T来描述断裂韧性：

T = -(∂U/∂A)ₗ

其中U是弹性应变能，A是裂纹面积，下标l表示在恒定载荷下求导。

（3）裂纹扩展阻力

软材料的断裂韧性通常随裂纹扩展而增加，表现出R曲线行为。这与裂纹尖端钝化、分子链取向、损伤区演化等因素有关。

8.2 内聚力模型

内聚力模型（Cohesive Zone Model, CZM）是描述软材料断裂的有效方法。在裂纹尖端引入一个内聚区，其中的牵引力-分离关系（TSL）描述了材料的断裂过程：

traction-separation law: σ = f(δ)

其中σ是内聚应力，δ是分离位移。常用的TSL包括双线性、指数型、梯形等。内聚力模型的参数（峰值应力、断裂能、特征分离位移）可以通过试验标定。

9. Python仿真实现

9.1 超弹性材料模型实现

以下是使用Python实现超弹性本构模型和大变形有限元分析的代码框架。我们将使用有限元库（如FEniCS、deal.II的Python接口）或自行实现简化版本。

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
import matplotlib.animation as animation
from matplotlib.patches import Rectangle, FancyBboxPatch
import warnings
warnings.filterwarnings('ignore')

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

class HyperelasticMaterial:
    """超弹性材料基类"""
    
    def __init__(self, params):
        self.params = params
    
    def strain_energy(self, F):
        """计算应变能密度"""
        raise NotImplementedError
    
    def stress(self, F):
        """计算柯西应力"""
        raise NotImplementedError

class NeoHookeanMaterial(HyperelasticMaterial):
    """Neo-Hookean超弹性材料"""
    
    def __init__(self, C10, D1=None):
        """
        参数:
        C10: 剪切相关参数，mu = 2*C10
        D1: 可压缩性参数，K = 2/D1；若为None则视为不可压缩
        """
        super().__init__({'C10': C10, 'D1': D1})
        self.C10 = C10
        self.D1 = D1
        self.mu = 2 * C10
        if D1 is not None:
            self.K = 2 / D1
    
    def strain_energy(self, F):
        """Neo-Hookean应变能函数"""
        C = np.dot(F.T, F)
        I1 = np.trace(C)
        J = np.linalg.det(F)
        
        psi = self.C10 * (I1 - 3)
        if self.D1 is not None:
            psi += (1.0 / self.D1) * (J - 1)**2
        return psi
    
    def stress(self, F):
        """计算柯西应力"""
        J = np.linalg.det(F)
        b = np.dot(F, F.T)  # 左柯西-格林张量
        
        # 偏应力部分
        sigma_dev = 2 * self.C10 / J * b
        
        # 体积应力部分
        if self.D1 is not None:
            sigma_vol = 2 / self.D1 * (J - 1) * np.eye(3)
        else:
            # 不可压缩情况，压力由约束确定
            sigma_vol = np.zeros((3, 3))
        
        sigma = sigma_dev + sigma_vol
        
        # 减去静水压力（偏应力形式）
        p = np.trace(sigma) / 3
        sigma = sigma - p * np.eye(3)
        
        return sigma
    
    def uniaxial_stress(self, stretch):
        """单轴拉伸应力-应变关系（解析解）"""
        # 对于不可压缩Neo-Hookean材料
        # sigma = mu * (lambda^2 - 1/lambda)
        sigma = self.mu * (stretch**2 - 1.0/stretch)
        return sigma

class MooneyRivlinMaterial(HyperelasticMaterial):
    """Mooney-Rivlin超弹性材料"""
    
    def __init__(self, C10, C01, D1=None):
        super().__init__({'C10': C10, 'C01': C01, 'D1': D1})
        self.C10 = C10
        self.C01 = C01
        self.D1 = D1
    
    def strain_energy(self, F):
        """Mooney-Rivlin应变能函数"""
        C = np.dot(F.T, F)
        I1 = np.trace(C)
        I2 = 0.5 * (I1**2 - np.trace(np.dot(C, C)))
        J = np.linalg.det(F)
        
        psi = self.C10 * (I1 - 3) + self.C01 * (I2 - 3)
        if self.D1 is not None:
            psi += (1.0 / self.D1) * (J - 1)**2
        return psi
    
    def stress(self, F):
        """计算柯西应力"""
        J = np.linalg.det(F)
        C = np.dot(F.T, F)
        b = np.dot(F, F.T)
        
        I1 = np.trace(C)
        I2 = 0.5 * (I1**2 - np.trace(np.dot(C, C)))
        
        # 应力计算
        sigma = 2/J * (self.C10 + self.C01 * I1) * b - 2/J * self.C01 * np.dot(b, b)
        
        if self.D1 is not None:
            sigma += 2 / self.D1 * (J - 1) * np.eye(3)
        
        return sigma

class OgdenMaterial(HyperelasticMaterial):
    """Ogden超弹性材料"""
    
    def __init__(self, mu_list, alpha_list, D1=None):
        """
        参数:
        mu_list: 模量参数列表
        alpha_list: 指数参数列表
        D1: 可压缩性参数
        """
        super().__init__({'mu': mu_list, 'alpha': alpha_list, 'D1': D1})
        self.mu_list = mu_list
        self.alpha_list = alpha_list
        self.D1 = D1
        self.n_terms = len(mu_list)
    
    def strain_energy(self, F):
        """Ogden应变能函数"""
        # 计算主拉伸比
        C = np.dot(F.T, F)
        eigenvalues = np.linalg.eigvalsh(C)
        lambdas = np.sqrt(eigenvalues)
        
        psi = 0
        for i in range(self.n_terms):
            mu_i = self.mu_list[i]
            alpha_i = self.alpha_list[i]
            psi += (mu_i / alpha_i) * (np.sum(lambdas**alpha_i) - 3)
        
        if self.D1 is not None:
            J = np.linalg.det(F)
            psi += (1.0 / self.D1) * (J - 1)**2
        
        return psi
    
    def stress(self, F):
        """计算柯西应力（主应力形式）"""
        # 计算主拉伸比
        C = np.dot(F.T, F)
        eigenvalues, eigenvectors = np.linalg.eigh(C)
        lambdas = np.sqrt(eigenvalues)
        
        # 计算主应力
        J = np.linalg.det(F)
        sigma_principal = np.zeros(3)
        
        for i in range(3):
            for j in range(self.n_terms):
                mu_j = self.mu_list[j]
                alpha_j = self.alpha_list[j]
                sigma_principal[i] += mu_j * lambdas[i]**alpha_j
        
        sigma_principal /= J
        
        # 组装全应力张量
        sigma = np.zeros((3, 3))
        for i in range(3):
            n_i = eigenvectors[:, i]
            sigma += sigma_principal[i] * np.outer(n_i, n_i)
        
        if self.D1 is not None:
            sigma += 2 / self.D1 * (J - 1) * np.eye(3)
        
        return sigma

class LargeDeformationSolver:
    """大变形问题求解器（简化版）"""
    
    def __init__(self, material, mesh_size=20):
        self.material = material
        self.mesh_size = mesh_size
    
    def uniaxial_tension(self, max_stretch=3.0, n_steps=50):
        """
        单轴拉伸分析
        
        参数:
        max_stretch: 最大拉伸比
        n_steps: 加载步数
        
        返回:
        stretches: 拉伸比数组
        stresses: 工程应力数组
        true_stresses: 真实应力数组
        """
        stretches = np.linspace(1.0, max_stretch, n_steps)
        stresses = []
        true_stresses = []
        
        for lam in stretches:
            # 不可压缩单轴拉伸的变形梯度
            # lambda_1 = lam, lambda_2 = lambda_3 = 1/sqrt(lam)
            F = np.array([
                [lam, 0, 0],
                [0, 1/np.sqrt(lam), 0],
                [0, 0, 1/np.sqrt(lam)]
            ])
            
            # 计算应力
            sigma = self.material.stress(F)
            
            # 工程应力 (基于初始截面积)
            sigma_eng = sigma[0, 0] / lam
            stresses.append(sigma_eng)
            
            # 真实应力
            true_stresses.append(sigma[0, 0])
        
        return stretches, np.array(stresses), np.array(true_stresses)
    
    def simple_shear(self, max_gamma=2.0, n_steps=50):
        """
        简单剪切分析
        
        参数:
        max_gamma: 最大剪切应变
        n_steps: 加载步数
        
        返回:
        gammas: 剪切应变数组
        tau: 剪应力数组
        sigma_n: 正应力数组
        """
        gammas = np.linspace(0, max_gamma, n_steps)
        tau = []
        sigma_n = []
        
        for gamma in gammas:
            # 简单剪切的变形梯度
            F = np.array([
                [1, gamma, 0],
                [0, 1, 0],
                [0, 0, 1]
            ])
            
            sigma = self.material.stress(F)
            tau.append(sigma[0, 1])
            sigma_n.append(sigma[0, 0])
        
        return gammas, np.array(tau), np.array(sigma_n)
    
    def equibiaxial_tension(self, max_stretch=2.0, n_steps=50):
        """
        等双轴拉伸分析
        
        参数:
        max_stretch: 最大拉伸比
        n_steps: 加载步数
        
        返回:
        stretches: 拉伸比数组
        stresses: 应力数组
        """
        stretches = np.linspace(1.0, max_stretch, n_steps)
        stresses = []
        
        for lam in stretches:
            # 等双轴拉伸: lambda_1 = lambda_2 = lam, lambda_3 = 1/lam^2
            F = np.array([
                [lam, 0, 0],
                [0, lam, 0],
                [0, 0, 1/lam**2]
            ])
            
            sigma = self.material.stress(F)
            stresses.append(sigma[0, 0])
        
        return stretches, np.array(stresses)

class ViscoelasticModel:
    """粘弹性材料模型"""
    
    def __init__(self, E0, E_inf, tau):
        """
        标准线性固体模型参数
        
        参数:
        E0: 初始模量
        E_inf: 平衡模量
        tau: 松弛时间
        """
        self.E0 = E0
        self.E_inf = E_inf
        self.tau = tau
    
    def relaxation_function(self, t):
        """松弛函数 G(t)"""
        return self.E_inf + (self.E0 - self.E_inf) * np.exp(-t / self.tau)
    
    def creep_function(self, t):
        """蠕变函数 J(t)"""
        return 1/self.E_inf - (1/self.E_inf - 1/self.E0) * np.exp(-t / self.tau)
    
    def stress_relaxation(self, t, epsilon_0):
        """应力松弛响应"""
        return epsilon_0 * self.relaxation_function(t)
    
    def creep_response(self, t, sigma_0):
        """蠕变响应"""
        return sigma_0 * self.creep_function(t)
    
    def dynamic_modulus(self, omega):
        """动态模量"""
        E_storage = self.E_inf + (self.E0 - self.E_inf) / (1 + (omega * self.tau)**2)
        E_loss = (self.E0 - self.E_inf) * omega * self.tau / (1 + (omega * self.tau)**2)
        return E_storage, E_loss

class SoftMaterialVisualization:
    """软材料可视化工具"""
    
    def __init__(self):
        pass
    
    def plot_stress_strain_comparison(self, materials, labels, max_stretch=3.0):
        """对比不同材料模型的应力-应变曲线"""
        fig, axes = plt.subplots(1, 2, figsize=(14, 5))
        
        solver = LargeDeformationSolver(None)
        
        for material, label in zip(materials, labels):
            solver.material = material
            stretches, stresses, true_stresses = solver.uniaxial_tension(max_stretch)
            
            # 工程应力-应变
            axes[0].plot(stretches - 1, stresses / 1e6, label=label, linewidth=2)
            
            # 真实应力-应变
            axes[1].plot(stretches - 1, true_stresses / 1e6, label=label, linewidth=2)
        
        axes[0].set_xlabel('Engineering Strain', fontsize=12)
        axes[0].set_ylabel('Engineering Stress (MPa)', fontsize=12)
        axes[0].set_title('Uniaxial Tension - Engineering Stress', fontsize=14)
        axes[0].legend()
        axes[0].grid(True, alpha=0.3)
        
        axes[1].set_xlabel('Engineering Strain', fontsize=12)
        axes[1].set_ylabel('True Stress (MPa)', fontsize=12)
        axes[1].set_title('Uniaxial Tension - True Stress', fontsize=14)
        axes[1].legend()
        axes[1].grid(True, alpha=0.3)
        
        plt.tight_layout()
        return fig
    
    def plot_deformation_modes(self, material):
        """绘制不同变形模式下的响应"""
        fig, axes = plt.subplots(2, 2, figsize=(14, 12))
        
        solver = LargeDeformationSolver(material)
        
        # 单轴拉伸
        stretches, stresses, true_stresses = solver.uniaxial_tension(3.0)
        axes[0, 0].plot(stretches, stresses / 1e6, 'b-', linewidth=2, label='Engineering')
        axes[0, 0].plot(stretches, true_stresses / 1e6, 'r--', linewidth=2, label='True')
        axes[0, 0].set_xlabel('Stretch Ratio', fontsize=11)
        axes[0, 0].set_ylabel('Stress (MPa)', fontsize=11)
        axes[0, 0].set_title('Uniaxial Tension', fontsize=12)
        axes[0, 0].legend()
        axes[0, 0].grid(True, alpha=0.3)
        
        # 简单剪切
        gammas, tau, sigma_n = solver.simple_shear(2.0)
        ax_twin = axes[0, 1].twinx()
        line1 = axes[0, 1].plot(gammas, tau / 1e6, 'b-', linewidth=2, label='Shear Stress')
        line2 = ax_twin.plot(gammas, sigma_n / 1e6, 'r--', linewidth=2, label='Normal Stress')
        axes[0, 1].set_xlabel('Shear Strain', fontsize=11)
        axes[0, 1].set_ylabel('Shear Stress (MPa)', fontsize=11, color='b')
        ax_twin.set_ylabel('Normal Stress (MPa)', fontsize=11, color='r')
        axes[0, 1].set_title('Simple Shear', fontsize=12)
        lines = line1 + line2
        labels = [l.get_label() for l in lines]
        axes[0, 1].legend(lines, labels, loc='upper left')
        axes[0, 1].grid(True, alpha=0.3)
        
        # 等双轴拉伸
        stretches_biaxial, stresses_biaxial = solver.equibiaxial_tension(2.0)
        axes[1, 0].plot(stretches_biaxial, stresses_biaxial / 1e6, 'g-', linewidth=2)
        axes[1, 0].set_xlabel('Stretch Ratio', fontsize=11)
        axes[1, 0].set_ylabel('Stress (MPa)', fontsize=11)
        axes[1, 0].set_title('Equibiaxial Tension', fontsize=12)
        axes[1, 0].grid(True, alpha=0.3)
        
        # 三种变形模式对比
        axes[1, 1].plot(stretches - 1, true_stresses / 1e6, 'b-', linewidth=2, label='Uniaxial')
        axes[1, 1].plot(stretches_biaxial - 1, stresses_biaxial / 1e6, 'r--', linewidth=2, label='Equibiaxial')
        axes[1, 1].set_xlabel('Strain', fontsize=11)
        axes[1, 1].set_ylabel('True Stress (MPa)', fontsize=11)
        axes[1, 1].set_title('Deformation Mode Comparison', fontsize=12)
        axes[1, 1].legend()
        axes[1, 1].grid(True, alpha=0.3)
        
        plt.tight_layout()
        return fig
    
    def plot_viscoelastic_response(self, model, t_max=10):
        """绘制粘弹性响应"""
        fig, axes = plt.subplots(2, 2, figsize=(14, 10))
        
        t = np.linspace(0, t_max, 200)
        
        # 应力松弛
        epsilon_0 = 0.1
        sigma_relax = model.stress_relaxation(t, epsilon_0)
        axes[0, 0].plot(t, sigma_relax / 1e6, 'b-', linewidth=2)
        axes[0, 0].axhline(y=model.E_inf * epsilon_0 / 1e6, color='r', linestyle='--', 
                          label=f'Equilibrium: {model.E_inf*epsilon_0/1e6:.2f} MPa')
        axes[0, 0].set_xlabel('Time (s)', fontsize=11)
        axes[0, 0].set_ylabel('Stress (MPa)', fontsize=11)
        axes[0, 0].set_title('Stress Relaxation', fontsize=12)
        axes[0, 0].legend()
        axes[0, 0].grid(True, alpha=0.3)
        
        # 蠕变
        sigma_0 = 1e6  # 1 MPa
        epsilon_creep = model.creep_response(t, sigma_0)
        axes[0, 1].plot(t, epsilon_creep * 100, 'g-', linewidth=2)
        axes[0, 1].axhline(y=sigma_0/model.E_inf * 100, color='r', linestyle='--',
                          label=f'Equilibrium: {sigma_0/model.E_inf*100:.2f}%')
        axes[0, 1].set_xlabel('Time (s)', fontsize=11)
        axes[0, 1].set_ylabel('Strain (%)', fontsize=11)
        axes[0, 1].set_title('Creep', fontsize=12)
        axes[0, 1].legend()
        axes[0, 1].grid(True, alpha=0.3)
        
        # 松弛函数
        G = model.relaxation_function(t)
        axes[1, 0].semilogy(t, G / 1e6, 'm-', linewidth=2)
        axes[1, 0].set_xlabel('Time (s)', fontsize=11)
        axes[1, 0].set_ylabel('Relaxation Modulus (MPa)', fontsize=11)
        axes[1, 0].set_title('Relaxation Function G(t)', fontsize=12)
        axes[1, 0].grid(True, alpha=0.3)
        
        # 动态模量
        omega = np.logspace(-2, 2, 100)
        E_storage, E_loss = model.dynamic_modulus(omega)
        axes[1, 1].semilogx(omega, E_storage / 1e6, 'b-', linewidth=2, label='Storage Modulus')
        axes[1, 1].semilogx(omega, E_loss / 1e6, 'r--', linewidth=2, label='Loss Modulus')
        axes[1, 1].set_xlabel('Frequency (rad/s)', fontsize=11)
        axes[1, 1].set_ylabel('Dynamic Modulus (MPa)', fontsize=11)
        axes[1, 1].set_title('Dynamic Mechanical Response', fontsize=12)
        axes[1, 1].legend()
        axes[1, 1].grid(True, alpha=0.3)
        
        plt.tight_layout()
        return fig
    
    def create_deformation_animation(self, material, output_file='deformation_animation.gif'):
        """创建变形过程的GIF动画"""
        fig, axes = plt.subplots(1, 2, figsize=(14, 6))
        
        # 原始网格
        n_elements = 10
        x = np.linspace(0, 1, n_elements + 1)
        y = np.linspace(0, 1, n_elements + 1)
        X, Y = np.meshgrid(x, y)
        
        # 存储每一帧
        frames = []
        stretches = np.linspace(1.0, 2.5, 50)
        
        for lam in stretches:
            # 变形后的坐标（单轴拉伸）
            X_def = X * lam
            Y_def = Y / np.sqrt(lam)
            
            # 计算应力
            F = np.array([[lam, 0, 0], [0, 1/np.sqrt(lam), 0], [0, 0, 1/np.sqrt(lam)]])
            sigma = material.stress(F)
            stress_val = sigma[0, 0] / 1e6  # MPa
            
            # 绘制
            axes[0].clear()
            axes[1].clear()
            
            # 变形网格
            for i in range(n_elements):
                for j in range(n_elements):
                    rect = plt.Rectangle((X_def[j, i], Y_def[j, i]), 
                                        X_def[j, i+1] - X_def[j, i],
                                        Y_def[j+1, i] - Y_def[j, i],
                                        fill=True, facecolor='lightblue', 
                                        edgecolor='blue', alpha=0.5)
                    axes[0].add_patch(rect)
            
            axes[0].set_xlim(-0.1, 2.7)
            axes[0].set_ylim(-0.1, 1.2)
            axes[0].set_aspect('equal')
            axes[0].set_title(f'Deformed Mesh (λ = {lam:.2f})', fontsize=12)
            axes[0].set_xlabel('X')
            axes[0].set_ylabel('Y')
            axes[0].grid(True, alpha=0.3)
            
            # 应力-应变曲线
            solver = LargeDeformationSolver(material)
            all_stretches, all_stresses, _ = solver.uniaxial_tension(2.5, 50)
            axes[1].plot(all_stretches - 1, all_stresses / 1e6, 'b-', linewidth=2)
            current_strain = lam - 1
            axes[1].plot(current_strain, stress_val, 'ro', markersize=10)
            axes[1].set_xlabel('Engineering Strain', fontsize=11)
            axes[1].set_ylabel('Engineering Stress (MPa)', fontsize=11)
            axes[1].set_title('Stress-Strain Response', fontsize=12)
            axes[1].grid(True, alpha=0.3)
            axes[1].set_xlim(0, 1.6)
            
            fig.canvas.draw()
            frame = np.frombuffer(fig.canvas.tostring_rgb(), dtype=np.uint8)
            frame = frame.reshape(fig.canvas.get_width_height()[::-1] + (3,))
            frames.append(frame)
        
        plt.close()
        
        # 保存为GIF
        import imageio
        imageio.mimsave(output_file, frames, fps=10)
        print(f"Animation saved to {output_file}")
        
        return output_file

# 运行案例研究
if __name__ == "__main__":
    print("=" * 80)
    print("软材料与大变形力学仿真")
    print("=" * 80)
    
    # 案例1: 不同超弹性材料模型对比
    print("\n案例1: 超弹性材料模型对比")
    print("-" * 50)
    
    # 定义材料（剪切模量约0.5 MPa）
    neo_hookean = NeoHookeanMaterial(C10=0.25e6)  # mu = 0.5 MPa
    mooney_rivlin = MooneyRivlinMaterial(C10=0.15e6, C01=0.1e6)
    ogden = OgdenMaterial(mu_list=[0.3e6, 0.1e6, 0.05e6], 
                         alpha_list=[2.0, -2.0, 4.0])
    
    materials = [neo_hookean, mooney_rivlin, ogden]
    labels = ['Neo-Hookean', 'Mooney-Rivlin', 'Ogden (3-term)']
    
    viz = SoftMaterialVisualization()
    fig1 = viz.plot_stress_strain_comparison(materials, labels)
    plt.savefig('case1_material_comparison.png', dpi=150, bbox_inches='tight')
    print("已保存: case1_material_comparison.png")
    plt.close()
    
    # 案例2: 不同变形模式分析
    print("\n案例2: 不同变形模式分析")
    print("-" * 50)
    
    fig2 = viz.plot_deformation_modes(ogden)
    plt.savefig('case2_deformation_modes.png', dpi=150, bbox_inches='tight')
    print("已保存: case2_deformation_modes.png")
    plt.close()
    
    # 案例3: 粘弹性响应分析
    print("\n案例3: 粘弹性响应分析")
    print("-" * 50)
    
    visco_model = ViscoelasticModel(E0=2e6, E_inf=0.5e6, tau=2.0)
    fig3 = viz.plot_viscoelastic_response(visco_model, t_max=15)
    plt.savefig('case3_viscoelastic_response.png', dpi=150, bbox_inches='tight')
    print("已保存: case3_viscoelastic_response.png")
    plt.close()
    
    # 案例4: 变形动画
    print("\n案例4: 生成变形动画")
    print("-" * 50)
    
    viz.create_deformation_animation(ogden, 'case4_deformation.gif')
    print("已保存: case4_deformation.gif")
    
    print("\n" + "=" * 80)
    print("所有案例完成！")
    print("=" * 80)

10. 工程案例分析

10.1 案例一：橡胶密封件设计

橡胶密封件广泛应用于汽车、航空航天、石油化工等领域。密封件的设计需要考虑大变形下的应力分布、接触压力和密封性能。

问题描述

设计一个O形橡胶密封圈，用于液压系统（工作压力20 MPa）。密封圈材料为丁腈橡胶（NBR），硬度邵尔A 70度。需要确定密封圈截面尺寸、沟槽尺寸，并验证密封性能。

材料参数

通过单轴拉伸试验获得丁腈橡胶的应力-应变数据，使用Ogden模型拟合得到材料参数：

• μ₁ = 1.25 MPa, α₁ = 2.5
• μ₂ = 0.15 MPa, α₂ = -2.0

仿真分析

使用有限元方法建立O形圈的轴对称模型，考虑以下分析步骤：

1. 装配分析：模拟密封圈装入沟槽的压缩过程，计算初始压缩率和接触压力分布。

2. 加压分析：在密封圈一侧施加20 MPa压力，分析密封圈变形和接触压力变化。

3. 密封性能评估：检查接触压力是否始终大于工作压力（密封的必要条件）。

结果与讨论

仿真结果显示，初始压缩率15%时，密封圈最大等效应力为12 MPa，远低于材料强度。加压后，低压侧接触压力保持在25 MPa以上，满足密封要求。通过参数研究，确定了最佳压缩率范围为12-18%。

10.2 案例二：人工血管力学分析

人工血管是替代病变血管的重要医疗器械。血管替代物需要具备良好的力学相容性——既要有足够的强度承受血压，又要有适当的顺应性匹配宿主血管。

问题描述

分析一种ePTFE（膨体聚四氟乙烯）人工血管在生理载荷下的力学行为。血管直径6 mm，壁厚0.5 mm，需要评估其在血压（120/80 mmHg）作用下的变形、应力分布和顺应性。

材料模型

ePTFE是多孔材料，表现出各向异性和非线性特征。采用各向异性超弹性模型描述：

• 基质材料：Neo-Hookean模型，C₁₀ = 5 MPa
• 纤维增强：沿周向分布的纤维，k₁ = 50 MPa, k₂ = 10

仿真分析

建立血管的圆柱壳模型，考虑：

1. 静态分析：计算收缩压（120 mmHg）和舒张压（80 mmHg）下的血管变形。

2. 顺应性计算：顺应性定义为直径变化率与压力变化的比值：

   C = (ΔD/D₀)/ΔP

3. 与天然血管对比：天然动脉的顺应性约为5-10 %/100mmHg，评估人工血管的匹配程度。

结果与讨论

仿真结果表明，ePTFE血管在生理压力范围内的顺应性约为2 %/100mmHg，低于天然血管。这可能导致血管-移植物交界处的不匹配，引起血流动力学异常。建议采用顺应性梯度设计或复合材料来改善匹配性。

10.3 案例三：软体机器人手指设计

软体机器人是机器人技术的前沿方向，软体手指能够自适应抓取不同形状和材质的物体。

问题描述

设计一种气动驱动的软体手指，能够产生弯曲变形实现抓取动作。手指长度80 mm，采用硅胶材料（Ecoflex 00-30），内部设有气腔。

设计原理

软体手指采用非对称结构设计——底部设有不可伸长的约束层（纤维增强），顶部为可膨胀的气腔。充气时，顶部伸长而底部长度不变，产生弯曲变形。

仿真分析

建立软体手指的参数化模型，分析：

1. 弯曲角度：不同气压下的指尖位移和弯曲角度。

2. 抓取力：手指与圆柱物体接触时的接触力和摩擦力。

3. 结构优化：优化气腔几何形状和约束层位置，最大化弯曲角度和输出力。

结果与讨论

优化后的设计在100 kPa气压下可产生120°弯曲角度，指尖输出力达5 N，足以抓取500 g物体。仿真结果与实验测量吻合良好，验证了模型的准确性。

10.4 案例四：皮肤组织冲击响应

皮肤是人体最大的器官，理解其在冲击载荷下的力学响应对于防护装备设计和损伤评估具有重要意义。

问题描述

分析皮肤在钝器冲击下的力学响应。皮肤厚度2 mm，需要考虑其多层结构（表皮、真皮、皮下组织）和粘弹性特征。

材料模型

• 表皮：较硬，弹性模量约10 MPa
• 真皮：较软，具有粘弹性，采用Fung模型描述
• 皮下组织：非常软，弹性模量约0.1 MPa

仿真分析

建立皮肤的三层有限元模型，模拟：

1. 低速冲击（< 5 m/s）：分析皮肤的变形和应力分布，评估挫伤风险。

2. 高速冲击（> 20 m/s）：分析应力波传播和局部化变形，评估撕裂风险。

3. 损伤准则：基于最大主应变和应变率建立皮肤损伤预测模型。

结果与讨论

仿真结果表明，皮肤在低速冲击下通过大变形吸收能量，应力分布较均匀。高速冲击下，应力波在表皮-真皮界面反射，导致界面处应力集中，易产生分层损伤。基于仿真结果，提出了皮肤防护设计的建议。

11. 总结与展望

11.1 本教程主要内容回顾

本教程系统介绍了软材料与大变形力学的理论基础和工程应用：

1. 软材料特性：介绍了弹性体、水凝胶、生物软组织、软体机器人材料等主要类型，阐述了超弹性、粘弹性、不可压缩性、各向异性等关键特征。

2. 大变形理论：建立了连续介质运动学框架，包括变形梯度、应变度量、应力度量等基本概念，介绍了客观性原理和平衡方程。

3. 本构模型：详细讲解了Neo-Hookean、Mooney-Rivlin、Ogden等超弹性模型，以及Maxwell、Kelvin-Voigt、Fung等粘弹性模型，介绍了各向异性材料的建模方法。

4. 数值方法：介绍了大变形有限元方法的基本原理，讨论了不可压缩问题的处理策略和求解算法。

5. 稳定性与断裂：分析了软材料结构的屈曲、起皱、局部化等失稳现象，介绍了软材料断裂力学和内聚力模型。

6. 工程应用：通过橡胶密封件、人工血管、软体机器人、皮肤组织等案例，展示了软材料力学在工程实践中的应用。

11.2 软材料力学的发展趋势

软材料与大变形力学是一个快速发展的研究领域，未来发展方向包括：

（1）多尺度建模

从分子链、交联网络到宏观结构，建立跨尺度的软材料力学模型。分子动力学与连续介质力学的耦合方法将为材料设计提供新工具。

（2）多物理场耦合

软材料常与电、磁、热、化学等物理场耦合。介电弹性体、磁流变弹性体、水凝胶溶胀等问题的研究将推动智能软材料的发展。

（3）数据驱动方法

机器学习与软材料力学的结合，可以从实验数据中直接学习本构关系，减少对传统物理模型的依赖，加速新材料开发。

（4）生物集成设计

仿生学和生物集成是软材料应用的重要方向。开发力学性能与生物组织匹配的软材料，实现无缝的人机接口和医疗器械集成。

（5）增材制造

3D打印技术为软材料结构制造提供了前所未有的自由度。多材料打印、4D打印等技术将实现复杂软体结构的快速制造。

11.3 学习建议

对于希望深入学习软材料力学的读者，建议：

1. 夯实基础：掌握连续介质力学、张量分析、有限元方法等基础知识。

2. 动手实践：通过编程实现本构模型和数值算法，加深对理论的理解。

3. 阅读文献：关注Journal of the Mechanics and Physics of Solids、Soft Matter、Extreme Mechanics Letters等期刊的最新研究。

4. 实验验证：有条件时进行材料试验，验证仿真模型的准确性。

5. 跨学科交流：软材料力学涉及材料科学、化学、生物学、医学等多个学科，跨学科合作将带来创新突破。

参考文献

1. Holzapfel, G.A. (2000). Nonlinear Solid Mechanics: A Continuum Approach for Engineering. Wiley.

2. Ogden, R.W. (1997). Non-linear Elastic Deformations. Dover Publications.

3. Treloar, L.R.G. (1975). The Physics of Rubber Elasticity. Oxford University Press.

4. Fung, Y.C. (1993). Biomechanics: Mechanical Properties of Living Tissues. Springer.

5. Bower, A.F. (2010). Applied Mechanics of Solids. CRC Press.

6. Holzapfel, G.A., & Ogden, R.W. (Eds.). (2006). Mechanics of Biological Tissue. Springer.

7. Polygerinos, P., et al. (2017). Soft Robotics: Review of Fluid-Driven Intrinsically Soft Devices. Soft Robotics, 4(3), 163-173.

8. Suo, Z. (2010). Theory of Dielectric Elastomers. Acta Mechanica Solida Sinica, 23(6), 549-578.

9. Arruda, E.M., & Boyce, M.C. (1993). A Three-Dimensional Constitutive Model for the Large Stretch Behavior of Rubber Elastic Materials. Journal of the Mechanics and Physics of Solids, 41(2), 389-412.

10. Gasser, T.C., Ogden, R.W., & Holzapfel, G.A. (2006). Hyperelastic Modelling of Arterial Layers with Distributed Collagen Fibre Orientations. Journal of the Royal Society Interface, 3(6), 15-35.

---

附录：完整Python代码

以下是本教程中使用的完整Python代码，包括所有类和案例：

"""
主题097：软材料与大变形力学仿真程序
包含超弹性模型、粘弹性模型、大变形分析等
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
import warnings
warnings.filterwarnings('ignore')

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

# 使用Agg后端，不弹出窗口
plt.switch_backend('Agg')


class HyperelasticMaterial:
    """超弹性材料基类"""
    
    def __init__(self, params):
        self.params = params
    
    def strain_energy(self, F):
        """计算应变能密度"""
        raise NotImplementedError
    
    def stress(self, F):
        """计算柯西应力"""
        raise NotImplementedError


class NeoHookeanMaterial(HyperelasticMaterial):
    """Neo-Hookean超弹性材料"""
    
    def __init__(self, C10, D1=None):
        """
        参数:
        C10: 剪切相关参数，mu = 2*C10
        D1: 可压缩性参数，K = 2/D1；若为None则视为不可压缩
        """
        super().__init__({'C10': C10, 'D1': D1})
        self.C10 = C10
        self.D1 = D1
        self.mu = 2 * C10
        if D1 is not None:
            self.K = 2 / D1
    
    def strain_energy(self, F):
        """Neo-Hookean应变能函数"""
        C = np.dot(F.T, F)
        I1 = np.trace(C)
        J = np.linalg.det(F)
        
        psi = self.C10 * (I1 - 3)
        if self.D1 is not None:
            psi += (1.0 / self.D1) * (J - 1)**2
        return psi
    
    def stress(self, F):
        """计算柯西应力"""
        J = np.linalg.det(F)
        b = np.dot(F, F.T)  # 左柯西-格林张量
        
        # 偏应力部分
        sigma_dev = 2 * self.C10 / J * b
        
        # 体积应力部分
        if self.D1 is not None:
            sigma_vol = 2 / self.D1 * (J - 1) * np.eye(3)
        else:
            # 不可压缩情况，压力由约束确定
            sigma_vol = np.zeros((3, 3))
        
        sigma = sigma_dev + sigma_vol
        
        # 减去静水压力（偏应力形式）
        p = np.trace(sigma) / 3
        sigma = sigma - p * np.eye(3)
        
        return sigma
    
    def uniaxial_stress(self, stretch):
        """单轴拉伸应力-应变关系（解析解）"""
        # 对于不可压缩Neo-Hookean材料
        # sigma = mu * (lambda^2 - 1/lambda)
        sigma = self.mu * (stretch**2 - 1.0/stretch)
        return sigma


class MooneyRivlinMaterial(HyperelasticMaterial):
    """Mooney-Rivlin超弹性材料"""
    
    def __init__(self, C10, C01, D1=None):
        super().__init__({'C10': C10, 'C01': C01, 'D1': D1})
        self.C10 = C10
        self.C01 = C01
        self.D1 = D1
    
    def strain_energy(self, F):
        """Mooney-Rivlin应变能函数"""
        C = np.dot(F.T, F)
        I1 = np.trace(C)
        I2 = 0.5 * (I1**2 - np.trace(np.dot(C, C)))
        J = np.linalg.det(F)
        
        psi = self.C10 * (I1 - 3) + self.C01 * (I2 - 3)
        if self.D1 is not None:
            psi += (1.0 / self.D1) * (J - 1)**2
        return psi
    
    def stress(self, F):
        """计算柯西应力"""
        J = np.linalg.det(F)
        C = np.dot(F.T, F)
        b = np.dot(F, F.T)
        
        I1 = np.trace(C)
        I2 = 0.5 * (I1**2 - np.trace(np.dot(C, C)))
        
        # 应力计算
        sigma = 2/J * (self.C10 + self.C01 * I1) * b - 2/J * self.C01 * np.dot(b, b)
        
        if self.D1 is not None:
            sigma += 2 / self.D1 * (J - 1) * np.eye(3)
        
        return sigma


class OgdenMaterial(HyperelasticMaterial):
    """Ogden超弹性材料"""
    
    def __init__(self, mu_list, alpha_list, D1=None):
        """
        参数:
        mu_list: 模量参数列表
        alpha_list: 指数参数列表
        D1: 可压缩性参数
        """
        super().__init__({'mu': mu_list, 'alpha': alpha_list, 'D1': D1})
        self.mu_list = mu_list
        self.alpha_list = alpha_list
        self.D1 = D1
        self.n_terms = len(mu_list)
    
    def strain_energy(self, F):
        """Ogden应变能函数"""
        # 计算主拉伸比
        C = np.dot(F.T, F)
        eigenvalues = np.linalg.eigvalsh(C)
        lambdas = np.sqrt(eigenvalues)
        
        psi = 0
        for i in range(self.n_terms):
            mu_i = self.mu_list[i]
            alpha_i = self.alpha_list[i]
            psi += (mu_i / alpha_i) * (np.sum(lambdas**alpha_i) - 3)
        
        if self.D1 is not None:
            J = np.linalg.det(F)
            psi += (1.0 / self.D1) * (J - 1)**2
        
        return psi
    
    def stress(self, F):
        """计算柯西应力（主应力形式）"""
        # 计算主拉伸比
        C = np.dot(F.T, F)
        eigenvalues, eigenvectors = np.linalg.eigh(C)
        lambdas = np.sqrt(eigenvalues)
        
        # 计算主应力
        J = np.linalg.det(F)
        sigma_principal = np.zeros(3)
        
        for i in range(3):
            for j in range(self.n_terms):
                mu_j = self.mu_list[j]
                alpha_j = self.alpha_list[j]
                sigma_principal[i] += mu_j * lambdas[i]**alpha_j
        
        sigma_principal /= J
        
        # 组装全应力张量
        sigma = np.zeros((3, 3))
        for i in range(3):
            n_i = eigenvectors[:, i]
            sigma += sigma_principal[i] * np.outer(n_i, n_i)
        
        if self.D1 is not None:
            sigma += 2 / self.D1 * (J - 1) * np.eye(3)
        
        return sigma


class LargeDeformationSolver:
    """大变形问题求解器（简化版）"""
    
    def __init__(self, material, mesh_size=20):
        self.material = material
        self.mesh_size = mesh_size
    
    def uniaxial_tension(self, max_stretch=3.0, n_steps=50):
        """
        单轴拉伸分析
        
        参数:
        max_stretch: 最大拉伸比
        n_steps: 加载步数
        
        返回:
        stretches: 拉伸比数组
        stresses: 工程应力数组
        true_stresses: 真实应力数组
        """
        stretches = np.linspace(1.0, max_stretch, n_steps)
        stresses = []
        true_stresses = []
        
        for lam in stretches:
            # 不可压缩单轴拉伸的变形梯度
            # lambda_1 = lam, lambda_2 = lambda_3 = 1/sqrt(lam)
            F = np.array([
                [lam, 0, 0],
                [0, 1/np.sqrt(lam), 0],
                [0, 0, 1/np.sqrt(lam)]
            ])
            
            # 计算应力
            sigma = self.material.stress(F)
            
            # 工程应力 (基于初始截面积)
            sigma_eng = sigma[0, 0] / lam
            stresses.append(sigma_eng)
            
            # 真实应力
            true_stresses.append(sigma[0, 0])
        
        return stretches, np.array(stresses), np.array(true_stresses)
    
    def simple_shear(self, max_gamma=2.0, n_steps=50):
        """
        简单剪切分析
        
        参数:
        max_gamma: 最大剪切应变
        n_steps: 加载步数
        
        返回:
        gammas: 剪切应变数组
        tau: 剪应力数组
        sigma_n: 正应力数组
        """
        gammas = np.linspace(0, max_gamma, n_steps)
        tau = []
        sigma_n = []
        
        for gamma in gammas:
            # 简单剪切的变形梯度
            F = np.array([
                [1, gamma, 0],
                [0, 1, 0],
                [0, 0, 1]
            ])
            
            sigma = self.material.stress(F)
            tau.append(sigma[0, 1])
            sigma_n.append(sigma[0, 0])
        
        return gammas, np.array(tau), np.array(sigma_n)
    
    def equibiaxial_tension(self, max_stretch=2.0, n_steps=50):
        """
        等双轴拉伸分析
        
        参数:
        max_stretch: 最大拉伸比
        n_steps: 加载步数
        
        返回:
        stretches: 拉伸比数组
        stresses: 应力数组
        """
        stretches = np.linspace(1.0, max_stretch, n_steps)
        stresses = []
        
        for lam in stretches:
            # 等双轴拉伸: lambda_1 = lambda_2 = lam, lambda_3 = 1/lam^2
            F = np.array([
                [lam, 0, 0],
                [0, lam, 0],
                [0, 0, 1/lam**2]
            ])
            
            sigma = self.material.stress(F)
            stresses.append(sigma[0, 0])
        
        return stretches, np.array(stresses)


class ViscoelasticModel:
    """粘弹性材料模型"""
    
    def __init__(self, E0, E_inf, tau):
        """
        标准线性固体模型参数
        
        参数:
        E0: 初始模量
        E_inf: 平衡模量
        tau: 松弛时间
        """
        self.E0 = E0
        self.E_inf = E_inf
        self.tau = tau
    
    def relaxation_function(self, t):
        """松弛函数 G(t)"""
        return self.E_inf + (self.E0 - self.E_inf) * np.exp(-t / self.tau)
    
    def creep_function(self, t):
        """蠕变函数 J(t)"""
        return 1/self.E_inf - (1/self.E_inf - 1/self.E0) * np.exp(-t / self.tau)
    
    def stress_relaxation(self, t, epsilon_0):
        """应力松弛响应"""
        return epsilon_0 * self.relaxation_function(t)
    
    def creep_response(self, t, sigma_0):
        """蠕变响应"""
        return sigma_0 * self.creep_function(t)
    
    def dynamic_modulus(self, omega):
        """动态模量"""
        E_storage = self.E_inf + (self.E0 - self.E_inf) / (1 + (omega * self.tau)**2)
        E_loss = (self.E0 - self.E_inf) * omega * self.tau / (1 + (omega * self.tau)**2)
        return E_storage, E_loss


class SoftMaterialVisualization:
    """软材料可视化工具"""
    
    def __init__(self):
        pass
    
    def plot_stress_strain_comparison(self, materials, labels, max_stretch=3.0):
        """对比不同材料模型的应力-应变曲线"""
        fig, axes = plt.subplots(1, 2, figsize=(14, 5))
        
        solver = LargeDeformationSolver(None)
        
        for material, label in zip(materials, labels):
            solver.material = material
            stretches, stresses, true_stresses = solver.uniaxial_tension(max_stretch)
            
            # 工程应力-应变
            axes[0].plot(stretches - 1, stresses / 1e6, label=label, linewidth=2)
            
            # 真实应力-应变
            axes[1].plot(stretches - 1, true_stresses / 1e6, label=label, linewidth=2)
        
        axes[0].set_xlabel('Engineering Strain', fontsize=12)
        axes[0].set_ylabel('Engineering Stress (MPa)', fontsize=12)
        axes[0].set_title('Uniaxial Tension - Engineering Stress', fontsize=14)
        axes[0].legend()
        axes[0].grid(True, alpha=0.3)
        
        axes[1].set_xlabel('Engineering Strain', fontsize=12)
        axes[1].set_ylabel('True Stress (MPa)', fontsize=12)
        axes[1].set_title('Uniaxial Tension - True Stress', fontsize=14)
        axes[1].legend()
        axes[1].grid(True, alpha=0.3)
        
        plt.tight_layout()
        return fig
    
    def plot_deformation_modes(self, material):
        """绘制不同变形模式下的响应"""
        fig, axes = plt.subplots(2, 2, figsize=(14, 12))
        
        solver = LargeDeformationSolver(material)
        
        # 单轴拉伸
        stretches, stresses, true_stresses = solver.uniaxial_tension(3.0)
        axes[0, 0].plot(stretches, stresses / 1e6, 'b-', linewidth=2, label='Engineering')
        axes[0, 0].plot(stretches, true_stresses / 1e6, 'r--', linewidth=2, label='True')
        axes[0, 0].set_xlabel('Stretch Ratio', fontsize=11)
        axes[0, 0].set_ylabel('Stress (MPa)', fontsize=11)
        axes[0, 0].set_title('Uniaxial Tension', fontsize=12)
        axes[0, 0].legend()
        axes[0, 0].grid(True, alpha=0.3)
        
        # 简单剪切
        gammas, tau, sigma_n = solver.simple_shear(2.0)
        ax_twin = axes[0, 1].twinx()
        line1 = axes[0, 1].plot(gammas, tau / 1e6, 'b-', linewidth=2, label='Shear Stress')
        line2 = ax_twin.plot(gammas, sigma_n / 1e6, 'r--', linewidth=2, label='Normal Stress')
        axes[0, 1].set_xlabel('Shear Strain', fontsize=11)
        axes[0, 1].set_ylabel('Shear Stress (MPa)', fontsize=11, color='b')
        ax_twin.set_ylabel('Normal Stress (MPa)', fontsize=11, color='r')
        axes[0, 1].set_title('Simple Shear', fontsize=12)
        lines = line1 + line2
        labels = [l.get_label() for l in lines]
        axes[0, 1].legend(lines, labels, loc='upper left')
        axes[0, 1].grid(True, alpha=0.3)
        
        # 等双轴拉伸
        stretches_biaxial, stresses_biaxial = solver.equibiaxial_tension(2.0)
        axes[1, 0].plot(stretches_biaxial, stresses_biaxial / 1e6, 'g-', linewidth=2)
        axes[1, 0].set_xlabel('Stretch Ratio', fontsize=11)
        axes[1, 0].set_ylabel('Stress (MPa)', fontsize=11)
        axes[1, 0].set_title('Equibiaxial Tension', fontsize=12)
        axes[1, 0].grid(True, alpha=0.3)
        
        # 三种变形模式对比
        axes[1, 1].plot(stretches - 1, true_stresses / 1e6, 'b-', linewidth=2, label='Uniaxial')
        axes[1, 1].plot(stretches_biaxial - 1, stresses_biaxial / 1e6, 'r--', linewidth=2, label='Equibiaxial')
        axes[1, 1].set_xlabel('Strain', fontsize=11)
        axes[1, 1].set_ylabel('True Stress (MPa)', fontsize=11)
        axes[1, 1].set_title('Deformation Mode Comparison', fontsize=12)
        axes[1, 1].legend()
        axes[1, 1].grid(True, alpha=0.3)
        
        plt.tight_layout()
        return fig
    
    def plot_viscoelastic_response(self, model, t_max=10):
        """绘制粘弹性响应"""
        fig, axes = plt.subplots(2, 2, figsize=(14, 10))
        
        t = np.linspace(0, t_max, 200)
        
        # 应力松弛
        epsilon_0 = 0.1
        sigma_relax = model.stress_relaxation(t, epsilon_0)
        axes[0, 0].plot(t, sigma_relax / 1e6, 'b-', linewidth=2)
        axes[0, 0].axhline(y=model.E_inf * epsilon_0 / 1e6, color='r', linestyle='--', 
                          label=f'Equilibrium: {model.E_inf*epsilon_0/1e6:.2f} MPa')
        axes[0, 0].set_xlabel('Time (s)', fontsize=11)
        axes[0, 0].set_ylabel('Stress (MPa)', fontsize=11)
        axes[0, 0].set_title('Stress Relaxation', fontsize=12)
        axes[0, 0].legend()
        axes[0, 0].grid(True, alpha=0.3)
        
        # 蠕变
        sigma_0 = 1e6  # 1 MPa
        epsilon_creep = model.creep_response(t, sigma_0)
        axes[0, 1].plot(t, epsilon_creep * 100, 'g-', linewidth=2)
        axes[0, 1].axhline(y=sigma_0/model.E_inf * 100, color='r', linestyle='--',
                          label=f'Equilibrium: {sigma_0/model.E_inf*100:.2f}%')
        axes[0, 1].set_xlabel('Time (s)', fontsize=11)
        axes[0, 1].set_ylabel('Strain (%)', fontsize=11)
        axes[0, 1].set_title('Creep', fontsize=12)
        axes[0, 1].legend()
        axes[0, 1].grid(True, alpha=0.3)
        
        # 松弛函数
        G = model.relaxation_function(t)
        axes[1, 0].semilogy(t, G / 1e6, 'm-', linewidth=2)
        axes[1, 0].set_xlabel('Time (s)', fontsize=11)
        axes[1, 0].set_ylabel('Relaxation Modulus (MPa)', fontsize=11)
        axes[1, 0].set_title('Relaxation Function G(t)', fontsize=12)
        axes[1, 0].grid(True, alpha=0.3)
        
        # 动态模量
        omega = np.logspace(-2, 2, 100)
        E_storage, E_loss = model.dynamic_modulus(omega)
        axes[1, 1].semilogx(omega, E_storage / 1e6, 'b-', linewidth=2, label='Storage Modulus')
        axes[1, 1].semilogx(omega, E_loss / 1e6, 'r--', linewidth=2, label='Loss Modulus')
        axes[1, 1].set_xlabel('Frequency (rad/s)', fontsize=11)
        axes[1, 1].set_ylabel('Dynamic Modulus (MPa)', fontsize=11)
        axes[1, 1].set_title('Dynamic Mechanical Response', fontsize=12)
        axes[1, 1].legend()
        axes[1, 1].grid(True, alpha=0.3)
        
        plt.tight_layout()
        return fig
    
    def create_deformation_animation(self, material, output_file='case4_deformation.gif'):
        """创建变形过程的GIF动画"""
        import imageio
        
        fig, axes = plt.subplots(1, 2, figsize=(14, 6))
        
        # 原始网格
        n_elements = 10
        x = np.linspace(0, 1, n_elements + 1)
        y = np.linspace(0, 1, n_elements + 1)
        X, Y = np.meshgrid(x, y)
        
        # 存储每一帧
        frames = []
        stretches = np.linspace(1.0, 2.5, 50)
        
        solver = LargeDeformationSolver(material)
        all_stretches, all_stresses, _ = solver.uniaxial_tension(2.5, 50)
        
        for lam in stretches:
            # 变形后的坐标（单轴拉伸）
            X_def = X * lam
            Y_def = Y / np.sqrt(lam)
            
            # 计算应力
            F = np.array([[lam, 0, 0], [0, 1/np.sqrt(lam), 0], [0, 0, 1/np.sqrt(lam)]])
            sigma = material.stress(F)
            stress_val = sigma[0, 0] / 1e6  # MPa
            
            # 绘制
            axes[0].clear()
            axes[1].clear()
            
            # 变形网格
            for i in range(n_elements):
                for j in range(n_elements):
                    rect = plt.Rectangle((X_def[j, i], Y_def[j, i]), 
                                        X_def[j, i+1] - X_def[j, i],
                                        Y_def[j+1, i] - Y_def[j, i],
                                        fill=True, facecolor='lightblue', 
                                        edgecolor='blue', alpha=0.5)
                    axes[0].add_patch(rect)
            
            axes[0].set_xlim(-0.1, 2.7)
            axes[0].set_ylim(-0.1, 1.2)
            axes[0].set_aspect('equal')
            axes[0].set_title(f'Deformed Mesh (λ = {lam:.2f})', fontsize=12)
            axes[0].set_xlabel('X')
            axes[0].set_ylabel('Y')
            axes[0].grid(True, alpha=0.3)
            
            # 应力-应变曲线
            axes[1].plot(all_stretches - 1, all_stresses / 1e6, 'b-', linewidth=2)
            current_strain = lam - 1
            axes[1].plot(current_strain, stress_val, 'ro', markersize=10)
            axes[1].set_xlabel('Engineering Strain', fontsize=11)
            axes[1].set_ylabel('Engineering Stress (MPa)', fontsize=11)
            axes[1].set_title('Stress-Strain Response', fontsize=12)
            axes[1].grid(True, alpha=0.3)
            axes[1].set_xlim(0, 1.6)
            
            fig.canvas.draw()
            frame = np.frombuffer(fig.canvas.buffer_rgba(), dtype=np.uint8)
            # Convert RGBA to RGB
            h, w = fig.canvas.get_width_height()
            frame = frame.reshape((h, w, 4))[:, :, :3]
            frames.append(frame)
        
        plt.close()
        
        # 保存为GIF
        imageio.mimsave(output_file, frames, fps=10)
        print(f"Animation saved to {output_file}")
        
        return output_file


# 运行案例研究
if __name__ == "__main__":
    print("=" * 80)
    print("软材料与大变形力学仿真")
    print("=" * 80)
    
    # 案例1: 不同超弹性材料模型对比
    print("\n案例1: 超弹性材料模型对比")
    print("-" * 50)
    
    # 定义材料（剪切模量约0.5 MPa）
    neo_hookean = NeoHookeanMaterial(C10=0.25e6)  # mu = 0.5 MPa
    mooney_rivlin = MooneyRivlinMaterial(C10=0.15e6, C01=0.1e6)
    ogden = OgdenMaterial(mu_list=[0.3e6, 0.1e6, 0.05e6], 
                         alpha_list=[2.0, -2.0, 4.0])
    
    materials = [neo_hookean, mooney_rivlin, ogden]
    labels = ['Neo-Hookean', 'Mooney-Rivlin', 'Ogden (3-term)']
    
    viz = SoftMaterialVisualization()
    fig1 = viz.plot_stress_strain_comparison(materials, labels)
    plt.savefig('case1_material_comparison.png', dpi=150, bbox_inches='tight')
    print("已保存: case1_material_comparison.png")
    plt.close()
    
    # 案例2: 不同变形模式分析
    print("\n案例2: 不同变形模式分析")
    print("-" * 50)
    
    fig2 = viz.plot_deformation_modes(ogden)
    plt.savefig('case2_deformation_modes.png', dpi=150, bbox_inches='tight')
    print("已保存: case2_deformation_modes.png")
    plt.close()
    
    # 案例3: 粘弹性响应分析
    print("\n案例3: 粘弹性响应分析")
    print("-" * 50)
    
    visco_model = ViscoelasticModel(E0=2e6, E_inf=0.5e6, tau=2.0)
    fig3 = viz.plot_viscoelastic_response(visco_model, t_max=15)
    plt.savefig('case3_viscoelastic_response.png', dpi=150, bbox_inches='tight')
    print("已保存: case3_viscoelastic_response.png")
    plt.close()
    
    # 案例4: 变形动画
    print("\n案例4: 生成变形动画")
    print("-" * 50)
    
    viz.create_deformation_animation(ogden, 'case4_deformation.gif')
    print("已保存: case4_deformation.gif")
    
    print("\n" + "=" * 80)
    print("所有案例完成！")
    print("=" * 80)