主题036：生物传热与流体

1. 引言

1.1 生物传热与流体的重要性

生物传热与流体是研究生物体内热量传递和流体运动的交叉学科，在医学、生物学和工程学中具有重要应用价值。人体是一个复杂的热力学系统，维持恒定的体温（约37°C）是生命活动的基础。理解生物传热机制对于疾病诊断、治疗计划制定和医疗器械设计至关重要。

主要应用领域：

• 热疗与冷冻治疗：肿瘤热消融、冷冻手术的温度场预测
• 医学影像：红外热成像、磁共振温度监测
• 医疗器械：血液透析、人工心脏、药物输送系统
• 人体舒适度：服装热湿舒适性、建筑环境设计
• 运动医学：体温调节、热应激分析
• 组织工程：生物反应器设计、组织培养温度控制

1.2 生物系统的复杂性

生物传热与流体问题具有以下特点：

多尺度特性：

• 微观：细胞尺度（μm）的代谢产热
• 介观：血管网络（mm）的对流换热
• 宏观：器官/全身（cm-m）的温度分布

多物理场耦合：

• 传热-流动耦合：血液流动与组织换热
• 传热-代谢耦合：组织产热与温度调节
• 传热-相变耦合：冻结/融化过程中的潜热

生物变异性：

• 个体差异：年龄、性别、体型、健康状况
• 组织异质性：不同器官的热物性差异
• 动态变化：代谢率、血流量的时变特性

1.3 本教程内容安排

本教程将系统介绍生物传热与流体的理论基础、数学模型和数值仿真方法，并通过六个典型工程案例展示Python在生物医学校拟中的应用：

1. 皮肤烧伤温度场分析
2. 肿瘤热疗计划优化
3. 心血管血流动力学仿真
4. 低温冷冻保存过程模拟
5. 全身热应激响应分析
6. 药物经皮渗透动力学

---

2. 生物组织热特性

2.1 生物组织的热物性

2.1.1 导热系数

生物组织的导热系数受含水量、脂肪含量和结构影响：

组织类型	导热系数 (W/m·K)	密度 (kg/m³)	比热容 (J/kg·K)
血液	0.50-0.60	1050-1060	3600-3800
肌肉	0.40-0.50	1030-1050	3400-3800
脂肪	0.20-0.25	850-950	2200-2500
皮肤	0.30-0.40	1000-1100	3000-3600
骨骼	0.30-0.70	1700-2000	1300-1700
肝脏	0.50-0.55	1050-1070	3500-3600
脑组织	0.50-0.60	1030-1040	3600-3700

导热系数随温度变化：

$$k(T)=k_{37}[1+{\alpha }_k(T-37)]$$

其中${\alpha }_k$是温度系数，约为-0.001 ~ -0.003 /°C（生物组织导热系数随温度升高略有下降）。

2.1.2 比热容与热扩散率

生物组织的比热容主要取决于含水量：

$$c_p=w_{water}{c}_{p,water}+(1-w_{water})c_{p,dry}$$

其中$w_{water}$是含水率，$c_{p,water}$ = 4186 J/kg·K，$c_{p,dry}$ ≈ 1800 J/kg·K。

热扩散率：

$$\alpha =\frac{k}{\rho c_p}$$

生物组织的热扩散率通常在$0.1\times 10^{-6}$ ~ $0.2\times 10^{-6}$ m²/s范围内。

2.1.3 代谢产热

基础代谢率（BMR）与体表面积的关系（Du Bois公式）：

$$BMR=71.2\times A_{body}^{0.425}\times h^{0.725}$$

其中$A_{body}$是体重（kg），$h$是身高（m）。

组织代谢产热率：

组织	代谢率 (W/kg)
肝脏	10-15
脑	8-12
心脏	8-10
肾脏	15-20
肌肉（静息）	0.5-1.0
肌肉（运动）	5-15
脂肪	0.2-0.5

代谢产热的温度依赖性（Q₁₀效应）：

$$q_m(T)=q_{m,37}\times Q_{10}^{(T-37)/10}$$

其中$Q_{10}$ ≈ 2-3，表示温度每升高10°C代谢率增加的倍数。

2.2 血液的热特性

2.2.1 血液组成与热物性

血液由血浆（55%）和血细胞（45%）组成：

血浆：

• 含水率：90-92%
• 蛋白质：7-8%
• 密度：1025-1030 kg/m³
• 比热容：3800-3900 J/kg·K

全血：

• 密度：1050-1060 kg/m³
• 比热容：3600-3800 J/kg·K
• 导热系数：0.50-0.60 W/m·K

血液粘度与温度关系：

$$\mu (T)={\mu }_{37}\times \exp [E_{visc}/R\times (1/T-1/310)]$$

其中$E_{visc}$ ≈ 15-20 kJ/mol是流动活化能。

2.2.2 血液灌注率

血液灌注率${\omega }_b$（单位：mL/min/100g或1/s）表示单位体积组织中的血流量：

组织	静息灌注率	最大灌注率
皮肤	0.5-2.0	10-20
肌肉	2-5	50-100
脑	50-60	60-70
肝脏	25-30	30-40
肾脏	200-400	400-500

灌注率与温度关系：

$${\omega }_b(T)={\omega }_{b,0}[1+{\alpha }_{vaso}(T-T_{set})]$$

其中${\alpha }_{vaso}$是血管舒张系数，$T_{set}$是体温设定点（约37°C）。

2.3 温度对生物组织的影响

2.3.1 热损伤机制

蛋白质变性动力学（Arrhenius模型）：

$$\Omega ={\int }_0^{t}A\exp \left(-\frac{E_a}{RT}\right)dt$$

其中$\Omega$是损伤积分，$A$是频率因子，$E_a$是活化能。

热损伤阈值：

损伤程度	温度	时间
可逆损伤	42-45°C	数小时
轻度烧伤	45-50°C	数分钟
中度烧伤	50-55°C	数十秒
重度烧伤	>55°C	<10秒
即时坏死	>60°C	瞬间

2.3.2 低温效应

相变温度：

• 纯水：0°C
• 生理盐水：-0.5 ~ -1.0°C
• 细胞外液：-2 ~ -5°C
• 细胞内液：-5 ~ -15°C

冷冻损伤机制：

• 机械损伤：冰晶形成导致细胞膜破裂
• 溶质损伤：细胞外冰晶导致细胞内溶质浓度升高
• 渗透损伤：水分跨膜运输引起的细胞皱缩/肿胀

---

3. 生物传热模型

3.1 Pennes生物传热方程

3.1.1 经典Pennes方程

Pennes方程（1948年）是描述生物组织传热最广泛使用的模型：

$$\rho c_p\frac{\partial T}{\partial t}=\nabla \cdot (k\nabla T)+{\omega }_b{\rho }_b{c}_{p,b}(T_a-T)+q_m$$

其中：

• $\rho ,c_p,k$：组织密度、比热容、导热系数
• ${\omega }_b$：血液灌注率（1/s）
• ${\rho }_b,c_{p,b}$：血液密度、比热容
• $T_a$：动脉血温度（通常取核心温度37°C）
• $q_m$：代谢产热率（W/m³）

3.1.2 方程各项物理意义

热存储项：$\rho c_p\frac{\partial T}{\partial t}$

• 表示组织温度随时间的变化率

热传导项：$\nabla \cdot (k\nabla T)$

• 描述热量在组织中的扩散

血液对流项：${\omega }_b{\rho }_b{c}_{p,b}(T_a-T)$

• 表示血液灌注引起的对流换热
• 当$T<T_a$时，该项为正（加热组织）
• 当$T>T_a$时，该项为负（冷却组织）

代谢产热项：$q_m$

• 组织代谢产生的热量

3.1.3 Pennes方程的局限性

1. 假设血液在毛细血管瞬间达到组织温度
2. 忽略血管走向和尺寸效应
3. 假设动脉血温度恒定
4. 不考虑血管网络结构

3.2 改进的生物传热模型

3.2.1 Chen-Holmes模型

考虑血液灌注的非均匀性：

$$\rho c_p\frac{\partial T}{\partial t}=\nabla \cdot (k_{eff}\nabla T)+{\omega }_b^{\ast }{\rho }_b{c}_{p,b}(T_a^{\ast }-T)+q_m$$

其中$k_{eff}$是有效导热系数，${\omega }_b^{\ast }$是有效灌注率，$T_a^{\ast }$是有效动脉血温度。

有效导热系数：

$$k_{eff}=k+k_{perf}$$

其中$k_{perf}$是灌注引起的增强导热。

3.2.2 Weinbaum-Jiji模型

考虑血管对传热的影响：

$$\rho c_p\frac{\partial T}{\partial t}=\nabla \cdot (k\nabla T)-\nabla \cdot q_p+q_m$$

其中$q_p$是血管对热流的贡献：

$$q_p={\int }_0^{\infty }n(r){\rho }_b{c}_{p,b}v(r)[T_a(r)-T]\pi r^{2}dr$$

$n(r)$是血管数密度分布，$v(r)$是血流速度。

3.2.3 多孔介质模型

将组织视为多孔介质，血液在孔隙中流动：

$$(1-\epsilon ){\rho }_t{c}_{p,t}\frac{\partial T_t}{\partial t}=(1-\epsilon )\nabla \cdot (k_t\nabla T_t)-h_{sf}{a}_{sf}(T_t-T_b)+(1-\epsilon )q_m$$

$$\epsilon {\rho }_b{c}_{p,b}\frac{\partial T_b}{\partial t}=\epsilon {\rho }_b{c}_{p,b}\mathbf{u}\cdot \nabla T_b+h_{sf}{a}_{sf}(T_t-T_b)$$

其中$\epsilon$是血液体积分数，$h_{sf}$是表面换热系数，$a_{sf}$是比表面积。

3.3 相变传热模型

3.3.1 相变问题描述

生物组织冻结/融化涉及相变潜热：

$$\rho L_f\frac{\partial f_l}{\partial t}=\nabla \cdot (k\nabla T)+{\omega }_b{\rho }_b{c}_{p,b}(T_a-T)+q_m$$

其中$L_f$是相变潜热（约334 kJ/kg），$f_l$是液相分数。

3.3.2 焓方法

定义焓：

$$H={\int }_0^T\rho c_{p}dT+\rho L_f(1-f_l)$$

焓形式的传热方程：

$$\frac{\partial H}{\partial t}=\nabla \cdot (k\nabla T)+{\omega }_b{\rho }_b{c}_{p,b}(T_a-T)+q_m$$

3.3.3 相变温度范围

生物组织相变发生在一定温度范围内：

$$f_l=\{\begin{array}{ll}1& T>T_{m,2}\\ \frac{T_{m,2}-T}{T_{m,2}-T_{m,1}}& T_{m,1}\le T\le T_{m,2}\\ 0& T<T_{m,1}\end{array}$$

其中$T_{m,1}$ ≈ -8°C，$T_{m,2}$ ≈ -1°C。

---

4. 血流动力学基础

4.1 血液流变学

4.1.1 非牛顿流体特性

血液是典型的非牛顿流体，表现出剪切稀化特性：

Casson模型：

$$\sqrt{\tau }=\sqrt{{\tau }_y}+\sqrt{{\mu }_c\dot{\gamma }}$$

其中${\tau }_y$是屈服应力，${\mu }_c$是Casson粘度。

幂律模型：

$$\tau =K{\dot{\gamma }}^n$$

其中$K$是稠度系数，$n$是幂律指数（血液$n$ ≈ 0.7-0.8）。

4.1.2 Fahraeus-Lindqvist效应

在微小血管（直径<1 mm）中，血液的表观粘度随管径减小而降低：

$${\mu }_{app}={\mu }_{\infty }{\left(\frac{D}{D+c}\right)}^2$$

其中${\mu }_{\infty }$是宏观粘度，$c$是特征长度（约10-20 μm）。

4.1.3 红细胞聚集

低剪切速率下，红细胞形成缗钱状聚集体：

$$\tau ={\tau }_y+{\mu }_p\dot{\gamma }+\frac{A}{\sqrt{\dot{\gamma }}}$$

其中$A$是聚集系数。

4.2 血管系统血流

4.2.1 动脉血流

大动脉中的脉动血流可用Womersley数表征：

$$Wo=D\sqrt{\frac{\omega \rho }{\mu }}$$

其中$\omega$是脉动频率（心率×2π/60）。

• $Wo<1$：准稳态流动
• $Wo>10$：充分发展的脉动流

4.2.2 微循环血流

毛细血管中的血流特性：

• 红细胞单行通过
• 血浆间隙层效应
• 血管舒缩调节

Starling方程描述跨血管壁的液体交换：

$$J_v=L_{p}A[(P_c-P_i)-\sigma ({\pi }_c-{\pi }_i)]$$

其中$L_p$是水力传导率，$\sigma$是反射系数，$P$和$\pi$分别是静水压和胶体渗透压。

4.2.3 血管网络模型

分形模型：

血管网络遵循Murray定律：

$$r_0^3=r_1^3+r_2^3$$

其中$r_0$是母血管半径，$r_1$和$r_2$是子血管半径。

阻力网络模型：

总血管阻力：

$$R_{total}=\sum _i{R}_i=\sum _i\frac{8\mu L_i}{\pi r_i^4}$$

4.3 心血管系统建模

4.3.1 集总参数模型

Windkessel模型：

$$Q_{in}-Q_{out}=C\frac{dP}{dt}+\frac{P-P_{venous}}{R}$$

其中$C$是血管顺应性，$R$是阻力。

4.3.2 一维血流模型

质量守恒：

$$\frac{\partial A}{\partial t}+\frac{\partial Q}{\partial x}=0$$

动量守恒：

$$\frac{\partial Q}{\partial t}+\frac{\partial }{\partial x}\left(\frac{Q^2}{A}\right)+\frac{A}{\rho }\frac{\partial P}{\partial x}=-f_{wall}$$

其中$A$是血管截面积，$Q$是流量，$f_{wall}$是壁面摩擦力。

4.3.3 三维血流模拟

Navier-Stokes方程：

$$\rho \left(\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\cdot \nabla \mathbf{u}\right)=-\nabla p+\nabla \cdot \mathbf{\tau }+\mathbf{f}$$

$$abla\cdot \mathbf{u}=0$$

对于非牛顿流体：

$$\mathbf{\tau }=\mu (\dot{\gamma })(\nabla \mathbf{u}+\nabla {\mathbf{u}}^T)$$

---

5. 工程应用案例

5.1 案例一：皮肤烧伤温度场分析

5.1.1 问题描述

皮肤是人体最大的器官，也是热交换的主要界面。分析皮肤在热暴露下的温度响应对于烧伤评估和防护设计具有重要意义。

皮肤结构：

• 表皮：厚度0.1-0.2 mm，无血管
• 真皮：厚度1-4 mm，含血管网
• 皮下组织：厚度5-20 mm，含脂肪层

热暴露条件：

• 热流体温度：60-100°C
• 暴露时间：1-10秒
• 对流换热系数：100-1000 W/m²·K

5.1.2 数学模型

一维多层皮肤模型：

$${\rho }_i{c}_{p,i}\frac{\partial T_i}{\partial t}=\frac{\partial }{\partial x}\left(k_i\frac{\partial T_i}{\partial x}\right)+{\omega }_{b,i}{\rho }_b{c}_{p,b}(T_a-T_i)+q_{m,i}$$

边界条件：

• 表面：$-k\frac{\partial T}{\partial x}=h_{ext}(T_{ext}-T_s)+\epsilon \sigma (T_{rad}^4-T_s^4)$
• 深层：$T=T_{core}$（核心温度37°C）

初始条件：$T(x,0)=T_{initial}$

5.1.3 热损伤评估

Arrhenius损伤积分：

$$\Omega (x,t)={\int }_0^{t}A\exp \left(-\frac{E_a}{RT(x,\tau )}\right)d\tau$$

损伤判据：

• $\Omega <0.53$：无损伤
• $0.53\le \Omega <1$：可逆损伤
• $\Omega \ge 1$：不可逆损伤（二级烧伤）
• $\Omega \ge 10^4$：三级烧伤

5.1.4 仿真结果

• 表皮基底层温度达到44°C所需时间：约3秒（60°C水）
• 二级烧伤临界时间：约5秒
• 热损伤深度与暴露时间呈正相关

工程应用：

• 防护服热防护性能评估
• 烧伤深度预测
• 急救时间窗口确定

5.2 案例二：肿瘤热疗计划优化

5.2.1 问题描述

热疗是利用高温（42-45°C）杀死癌细胞的治疗方法。射频消融（RFA）和微波消融（MWA）是常用的局部热疗技术。

系统参数：

• 肿瘤尺寸：直径3 cm球形
• 目标温度：50-60°C（消融区）
• 治疗时间：10-15分钟
• 电极：单针或多针射频电极

5.2.2 射频消融模型

电磁场分布（简化电阻模型）：

$$q^{\prime \prime \prime }=\sigma E^2=\frac{\sigma V^2}{4{\pi }^2{r}^4}$$

其中$\sigma$是组织电导率，$V$是施加电压。

生物传热方程：

$$\rho c_p\frac{\partial T}{\partial t}=\nabla \cdot (k\nabla T)+{\omega }_b{\rho }_b{c}_{p,b}(T_a-T)+q^{\prime \prime \prime }+q_m$$

组织电导率温度依赖性：

$$\sigma (T)={\sigma }_{37}[1+{\alpha }_{\sigma }(T-37)]$$

其中${\alpha }_{\sigma }$ ≈ 0.02 /°C。

5.2.3 热损伤预测

消融区边界（等温线法）：

• 50°C等温线：即时坏死边界
• 54°C等温线：可靠消融边界

热损伤积分法：

• 累积等效时间（CEM）：$CE{M}_{43}={\int }_0^t{R}^{(43-T)}dt$
• $R=0.5$（T<43°C），$R=0.25$（T≥43°C）
• CEM₄₃ > 240分钟：细胞死亡

5.2.4 仿真结果

• 消融区直径：约2.5-3.0 cm
• 消融区形状：椭球形（受血流影响）
• 大血管附近存在"热沉效应"

治疗优化：

• 多针电极布局设计
• 功率调制策略
• 冷却电极防止炭化

5.3 案例三：心血管血流动力学仿真

5.3.1 问题描述

心血管疾病的诊断和治疗需要理解血流动力学特性。颈动脉狭窄是脑卒中的重要危险因素。

几何参数：

• 颈总动脉直径：6-8 mm
• 分叉角度：20-30°
• 狭窄程度：50-75%
• 狭窄长度：5-10 mm

5.3.2 血流模型

脉动入口条件：

$$Q(t)=Q_{mean}\left[1+0.5\sin \left(\frac{2\pi t}{T_{cardiac}}\right)-0.2\sin \left(\frac{4\pi t}{T_{cardiac}}\right)\right]$$

非牛顿粘度模型（Carreau-Yasuda）：

$$\frac{\mu -{\mu }_{\infty }}{{\mu }_0-{\mu }_{\infty }}=[1+(\lambda \dot{\gamma }{)}^a{]}^{(n-1)/a}$$

其中${\mu }_0$ = 0.056 Pa·s，${\mu }_{\infty }$ = 0.0035 Pa·s，$\lambda$ = 3.313 s，$a$ = 1.25，$n$ = 0.356。

5.3.3 壁面剪切应力分析

壁面剪切应力（WSS）：

$${\tau }_w=\mu \frac{\partial u}{\partial n}{|}_{wall}$$

病理相关指标：

• 低WSS（<1 Pa）：促进动脉粥样硬化
• 高WSS（>10 Pa）：内皮损伤风险
• WSS振荡：不稳定血流标志

5.3.4 仿真结果

• 狭窄处最大流速：增加3-5倍
• 狭窄下游出现流动分离和涡流
• 狭窄肩部高WSS区域
• 狭窄远心端低WSS和振荡区域

临床意义：

• 评估卒中风险
• 支架设计优化
• 手术方案规划

5.4 案例四：低温冷冻保存过程模拟

5.4.1 问题描述

细胞和组织冷冻保存是再生医学和辅助生殖技术的关键。优化冷冻/复温过程对于保持细胞活力至关重要。

系统参数：

• 样品：细胞悬浮液（1 mL）
• 冷冻保护剂：10% DMSO
• 冷冻速率：1-100°C/min
• 目标温度：-196°C（液氮）

5.4.2 两参数模型

Mazur两参数模型预测细胞存活率：

$$\frac{dV}{dt}=-L_{p}A\left[\Delta P-\Delta \pi -\frac{RT}{V_w}\ln \left(\frac{V-V_b}{V_0-V_b}\right)\right]$$

其中$V$是细胞体积，$L_p$是水力传导率，$V_b$是不可压缩体积。

最佳冷冻速率：

$$B_{opt}=\frac{L_{p}A\Delta H_f}{k_{cell}{V}_0}\sqrt{\frac{D_{CPA}}{L_p}}$$

其中$\Delta H_f$是相变潜热，$D_{CPA}$是冷冻保护剂扩散系数。

5.4.3 传热分析

考虑相变的传热方程：

$$\rho c_{eff}\frac{\partial T}{\partial t}=\nabla \cdot (k\nabla T)$$

有效比热容：

$$c_{eff}=c_p+L_f\frac{\partial f_l}{\partial T}$$

5.4.4 仿真结果

• 最佳冷冻速率：约5-10°C/min（细胞）
• 快速复温（>100°C/min）减少重结晶损伤
• 玻璃化冷冻可实现无冰保存

应用：

• 干细胞库建设
• 卵子/精子冷冻
• 组织工程产品保存

5.5 案例五：全身热应激响应分析

5.5.1 问题描述

人体在高温环境下的热应激反应涉及复杂的生理调节机制。预测体温变化对于职业安全和运动医学具有重要意义。

环境条件：

• 干球温度：30-45°C
• 相对湿度：40-80%
• 代谢率：静息（100 W）到剧烈运动（800 W）
• 服装热阻：0.5-1.5 clo

5.5.2 人体热调节模型

核心温度变化：

$$\rho c_p{V}_{core}\frac{d{T}_{core}}{dt}=q_{met}-q_{res}-q_{skin}-q_{resp}$$

皮肤温度：

$$\rho c_p{V}_{skin}\frac{d{T}_{skin}}{dt}=q_{blood}-q_{evap}-q_{conv}-q_{rad}$$

出汗调节：

$${\dot{m}}_{sw}={\dot{m}}_{sw,max}\times \frac{T_{core}-T_{set,core}}{T_{set,core}+\Delta T-T_{set,core}}\times \frac{T_{skin}-T_{set,skin}}{T_{set,skin}+\Delta T-T_{set,skin}}$$

5.5.3 热应激指标

湿球黑球温度（WBGT）：

$$WBGT=0.7T_{wet}+0.2T_{globe}+0.1T_{dry}$$

预测热应激指数（PHS）：

基于热平衡方程预测可暴露时间。

5.5.4 仿真结果

• 静息状态下，环境温度40°C时核心温度上升速率：约0.5°C/h
• 剧烈运动时，热平衡可能被打破
• 高湿度显著降低蒸发散热效率

安全建议：

• 工作-休息周期规划
• 补水策略
• 个体防护装备选择

5.6 案例六：药物经皮渗透动力学

5.6.1 问题描述

经皮给药是一种非侵入性药物递送方式。理解药物在皮肤中的传输机制对于贴剂设计至关重要。

皮肤结构模型：

• 角质层：主要屏障，厚度10-20 μm
• 活性表皮：厚度50-100 μm
• 真皮：厚度1-2 mm

5.6.2 扩散模型

Fick第二定律：

$$\frac{\partial C}{\partial t}=D\frac{{\partial }^{2}C}{\partial x^2}$$

角质层有效扩散系数：

$$D_{sc}=D_0\exp \left(-\frac{E_a}{RT}\right)\times {\epsilon }_{lipid}$$

其中${\epsilon }_{lipid}$是脂质通道体积分数。

5.6.3 增强技术

化学增强剂：

• 表面活性剂
• 脂肪酸
• 氮酮（Azone）

物理增强方法：

• 离子导入：$J_{ion}=\frac{zFD}{RT}C\frac{d\varphi }{dx}$
• 超声导入
• 微针阵列

5.6.4 仿真结果

• 稳态通量：$J_{ss}=\frac{DK{C}_0}{h}$
• 滞后时间：$t_{lag}=\frac{h^2}{6D}$
• 增强剂可提高渗透率10-100倍

贴剂设计：

• 储库型vs骨架型
• 控释膜设计
• 粘附性能优化

---

6. Python仿真实现

以下是六个案例的Python仿真代码实现：

# -*- coding: utf-8 -*-
"""
主题036：生物传热与流体 - 案例仿真代码
"""
import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint, solve_ivp
from scipy.sparse import diags
from scipy.sparse.linalg import spsolve
from scipy.special import erf, erfc
import warnings
warnings.filterwarnings('ignore')
import os

output_dir = r'd:\文档\500仿真领域\工程仿真\多物理场耦合仿真\主题036'
os.makedirs(output_dir, exist_ok=True)

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

# 物理常数
sigma_sb = 5.67e-8  # Stefan-Boltzmann常数 (W/m²·K⁴)
R_gas = 8.314  # 气体常数 (J/mol·K)

print("=" * 70)
print("主题036：生物传热与流体 - 案例仿真")
print("=" * 70)

# ==================== 案例一：皮肤烧伤温度场分析 ====================
print("\n" + "=" * 70)
print("案例一：皮肤烧伤温度场分析")
print("=" * 70)

# 皮肤分层参数
layers = {
    'epidermis': {'thickness': 0.15e-3, 'k': 0.21, 'rho': 1200, 'cp': 3600, 'omega': 0, 'qm': 0},
    'dermis': {'thickness': 2.0e-3, 'k': 0.45, 'rho': 1200, 'cp': 3600, 'omega': 0.001, 'qm': 368},
    'subcutaneous': {'thickness': 5.0e-3, 'k': 0.19, 'rho': 850, 'cp': 2500, 'omega': 0.0005, 'qm': 368}
}

# 血液参数
rho_b = 1060  # kg/m³
cp_b = 3800  # J/kg·K
T_a = 37.0  # 动脉血温度 °C

# 热暴露条件
T_ext = 80.0  # 外部热流体温度 °C
h_ext = 500.0  # 对流换热系数 W/m²·K

# 数值参数
t_total = 10.0  # 总时间 s
dt = 0.01  # 时间步长 s

# 建立网格
layer_names = ['epidermis', 'dermis', 'subcutaneous']
x_faces = [0]
for name in layer_names:
    x_faces.append(x_faces[-1] + layers[name]['thickness'])

# 在每个层内创建均匀网格
nx_per_layer = 20
x_grid = []
for i, name in enumerate(layer_names):
    x_start = x_faces[i]
    x_end = x_faces[i+1]
    x_layer = np.linspace(x_start, x_end, nx_per_layer + 1)
    if i > 0:
        x_layer = x_layer[1:]  # 避免重复节点
    x_grid.extend(x_layer)

x_grid = np.array(x_grid)
nx = len(x_grid)
dx = np.diff(x_grid)

# 获取每个节点的物性
def get_properties(x):
    """获取位置x处的物性参数"""
    if x <= layers['epidermis']['thickness']:
        layer = 'epidermis'
    elif x <= layers['epidermis']['thickness'] + layers['dermis']['thickness']:
        layer = 'dermis'
    else:
        layer = 'subcutaneous'
    return layers[layer]

# 初始化温度
T_initial = 37.0
T = np.ones(nx) * T_initial

# 时间步进
nt = int(t_total / dt)
t_history = [0]
T_history = [T.copy()]

# 热损伤参数
A_damage = 3.1e98  # 1/s
Ea_damage = 6.28e5  # J/mol
Omega = np.zeros(nx)  # 损伤积分

print(f"\n皮肤结构:")
print(f"  表皮厚度: {layers['epidermis']['thickness']*1e6:.0f} μm")
print(f"  真皮厚度: {layers['dermis']['thickness']*1e3:.1f} mm")
print(f"  皮下组织厚度: {layers['subcutaneous']['thickness']*1e3:.1f} mm")
print(f"\n热暴露条件:")
print(f"  热流体温度: {T_ext:.1f}°C")
print(f"  对流换热系数: {h_ext:.0f} W/m²·K")
print(f"  暴露时间: {t_total:.1f} s")

# 简化的一维有限差分求解
for n in range(nt):
    T_new = T.copy()
    
    # 内部节点
    for i in range(1, nx-1):
        props = get_properties(x_grid[i])
        k = props['k']
        rho = props['rho']
        cp = props['cp']
        omega = props['omega']
        qm = props['qm']
        
        dx_local = (dx[i-1] + dx[i]) / 2
        
        # 热传导项
        d2T_dx2 = (T[i+1] - 2*T[i] + T[i-1]) / dx_local**2
        q_cond = k * d2T_dx2
        
        # 血液对流项
        q_blood = omega * rho_b * cp_b * (T_a - T[i])
        
        # 代谢产热
        q_met = qm
        
        # 温度更新
        dT_dt = (q_cond + q_blood + q_met) / (rho * cp)
        T_new[i] = T[i] + dt * dT_dt
        
        # 热损伤积分
        if T[i] > 43 + 273.15:
            Omega[i] += A_damage * np.exp(-Ea_damage / (R_gas * (T[i] + 273.15))) * dt
    
    # 边界条件
    # 表面（对流边界）
    props_surf = get_properties(x_grid[0])
    q_conv = h_ext * (T_ext - T[0])
    q_cond_surf = props_surf['k'] * (T[1] - T[0]) / dx[0]
    T_new[0] = T[0] + dt * (q_conv + q_cond_surf) / (props_surf['rho'] * props_surf['cp'] * dx[0]/2)
    
    # 深层（Dirichlet边界）
    T_new[-1] = T_initial
    
    T = T_new.copy()
    
    if (n+1) % 100 == 0:
        t_history.append((n+1) * dt)
        T_history.append(T.copy())

T_history = np.array(T_history)

# 分析结果
T_surface_max = np.max(T_history[:, 0])
T_base = T_history[-1, np.argmin(np.abs(x_grid - layers['epidermis']['thickness'] - layers['dermis']['thickness']))]

print(f"\n仿真结果:")
print(f"  表面最高温度: {T_surface_max:.1f}°C")
print(f"  基底温度: {T_base:.1f}°C")
print(f"  最大热损伤积分: {np.max(Omega):.2e}")

# 判断烧伤程度
if np.max(Omega) < 0.53:
    burn_degree = "无损伤"
elif np.max(Omega) < 1:
    burn_degree = "一级烧伤（可逆）"
elif np.max(Omega) < 1e4:
    burn_degree = "二级烧伤"
else:
    burn_degree = "三级烧伤"
print(f"  烧伤程度: {burn_degree}")

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

ax1 = axes[0, 0]
for i, t in enumerate([0, 2, 5, 10]):
    idx = min(int(t / dt / 100), len(t_history)-1)
    ax1.plot(x_grid * 1e3, T_history[idx], linewidth=2, label=f't = {t_history[idx]:.1f} s')
ax1.axvline(x=layers['epidermis']['thickness']*1e3, color='gray', linestyle='--', alpha=0.5, label='Epidermis/Dermis')
ax1.axvline(x=(layers['epidermis']['thickness']+layers['dermis']['thickness'])*1e3, color='gray', linestyle=':', alpha=0.5, label='Dermis/Subcutaneous')
ax1.set_xlabel('Depth (mm)', fontsize=11)
ax1.set_ylabel('Temperature (°C)', fontsize=11)
ax1.set_title('Temperature Distribution in Skin', fontsize=12, fontweight='bold')
ax1.legend()
ax1.grid(True, alpha=0.3)

ax2 = axes[0, 1]
ax2.semilogy(x_grid * 1e3, Omega + 1e-10, 'r-', linewidth=2)
ax2.axhline(y=0.53, color='orange', linestyle='--', alpha=0.7, label='1st degree')
ax2.axhline(y=1, color='red', linestyle='--', alpha=0.7, label='2nd degree')
ax2.axhline(y=1e4, color='darkred', linestyle='--', alpha=0.7, label='3rd degree')
ax2.set_xlabel('Depth (mm)', fontsize=11)
ax2.set_ylabel('Damage Integral Ω', fontsize=11)
ax2.set_title('Thermal Damage Distribution', fontsize=12, fontweight='bold')
ax2.legend()
ax2.grid(True, alpha=0.3)

ax3 = axes[1, 0]
time_points = np.array(t_history)
surface_temp = T_history[:, 0]
base_temp = T_history[:, np.argmin(np.abs(x_grid - layers['epidermis']['thickness']))]
ax3.plot(time_points, surface_temp, 'b-', linewidth=2, label='Surface')
ax3.plot(time_points, base_temp, 'r-', linewidth=2, label='Dermis base')
ax3.axhline(y=44, color='gray', linestyle='--', alpha=0.5, label='Pain threshold')
ax3.set_xlabel('Time (s)', fontsize=11)
ax3.set_ylabel('Temperature (°C)', fontsize=11)
ax3.set_title('Temperature Transient', fontsize=12, fontweight='bold')
ax3.legend()
ax3.grid(True, alpha=0.3)

ax4 = axes[1, 1]
# 不同外部温度下的烧伤深度
T_ext_range = [60, 70, 80, 90, 100]
burn_depths = []
for T_ext_test in T_ext_range:
    # 简化的烧伤深度估算
    burn_depth = 0.5e-3 * ((T_ext_test - 37) / 43)**2 * t_total / 5
    burn_depths.append(min(burn_depth * 1e3, 7))  # 限制在7mm内
ax4.plot(T_ext_range, burn_depths, 'mo-', linewidth=2, markersize=8)
ax4.set_xlabel('External Temperature (°C)', fontsize=11)
ax4.set_ylabel('Estimated Burn Depth (mm)', fontsize=11)
ax4.set_title('Burn Depth vs Temperature', fontsize=12, fontweight='bold')
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/case1_skin_burn.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 案例一完成")

# ==================== 案例二：肿瘤热疗计划优化 ====================
print("\n" + "=" * 70)
print("案例二：肿瘤热疗计划优化")
print("=" * 70)

# 肿瘤和周围组织参数
r_tumor = 0.015  # 肿瘤半径 m
r_domain = 0.05  # 计算域半径 m

# 组织物性
k_tissue = 0.5  # W/m·K
rho_tissue = 1050  # kg/m³
cp_tissue = 3600  # J/kg·K
omega_b_normal = 0.001  # 1/s
omega_b_tumor = 0.002  # 肿瘤血流更丰富

# 射频消融参数
P_rf = 50  # RF功率 W
r_electrode = 0.001  # 电极半径 m

# 建立径向网格
nr = 100
r_grid = np.linspace(r_electrode, r_domain, nr)
dr = r_grid[1] - r_grid[0]

# 初始条件
T_initial = 37.0
T = np.ones(nr) * T_initial

# 治疗时间
t_treatment = 600  # 10分钟
dt = 0.5  # 时间步长
nt = int(t_treatment / dt)

print(f"\n肿瘤参数:")
print(f"  肿瘤半径: {r_tumor*1e3:.1f} mm")
print(f"  治疗时间: {t_treatment/60:.0f} min")
print(f"  RF功率: {P_rf:.0f} W")

# 热源分布（简化模型）
def heat_source(r, P, r_elec):
    """RF热源分布"""
    sigma = 0.3  # S/m 组织电导率
    V = np.sqrt(P / (2 * np.pi * sigma * r_elec))
    q = sigma * V**2 / (4 * np.pi**2 * r**4) * np.exp(-(r - r_elec) / 0.01)
    return q * 1e-6  # 转换为 W/m³

# 时间步进
for n in range(nt):
    T_new = T.copy()
    
    for i in range(1, nr-1):
        # 判断位置
        if r_grid[i] < r_tumor:
            omega = omega_b_tumor
        else:
            omega = omega_b_normal
        
        # 径向热传导
        dT_dr = (T[i+1] - T[i-1]) / (2 * dr)
        d2T_dr2 = (T[i+1] - 2*T[i] + T[i-1]) / dr**2
        q_cond = k_tissue * (d2T_dr2 + (2/r_grid[i]) * dT_dr)
        
        # 血液对流
        q_blood = omega * rho_b * cp_b * (T_a - T[i])
        
        # RF热源
        q_rf = heat_source(r_grid[i], P_rf, r_electrode)
        
        # 温度更新
        dT_dt = (q_cond + q_blood + q_rf) / (rho_tissue * cp_tissue)
        T_new[i] = T[i] + dt * dT_dt
    
    # 边界条件
    T_new[0] = T[1]  # 电极表面绝热
    T_new[-1] = T_initial  # 远场恒温
    
    T = T_new.copy()

# 分析消融区
ablation_50 = r_grid[T >= 50]  # 50°C等温线
ablation_54 = r_grid[T >= 54]  # 54°C等温线

if len(ablation_50) > 0:
    r_ablation_50 = np.max(ablation_50)
else:
    r_ablation_50 = 0

if len(ablation_54) > 0:
    r_ablation_54 = np.max(ablation_54)
else:
    r_ablation_54 = 0

print(f"\n消融结果:")
print(f"  最高温度: {np.max(T):.1f}°C")
print(f"  50°C消融区直径: {r_ablation_50*2*1e3:.1f} mm")
print(f"  54°C消融区直径: {r_ablation_54*2*1e3:.1f} mm")

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

ax1 = axes[0, 0]
ax1.plot(r_grid * 1e3, T, 'b-', linewidth=2)
ax1.axvline(x=r_tumor*1e3, color='red', linestyle='--', alpha=0.7, label='Tumor boundary')
ax1.axhline(y=50, color='orange', linestyle='--', alpha=0.7, label='50°C (ablation)')
ax1.axhline(y=54, color='red', linestyle='--', alpha=0.7, label='54°C (coagulation)')
ax1.set_xlabel('Radius (mm)', fontsize=11)
ax1.set_ylabel('Temperature (°C)', fontsize=11)
ax1.set_title('Temperature Distribution (RFA)', fontsize=12, fontweight='bold')
ax1.legend()
ax1.grid(True, alpha=0.3)

ax2 = axes[0, 1]
# 不同功率下的消融区
P_range = [20, 30, 40, 50, 60, 80]
r_ablation_range = []
for P_test in P_range:
    T_test = np.ones(nr) * T_initial
    for n in range(nt):
        T_new = T_test.copy()
        for i in range(1, nr-1):
            omega = omega_b_tumor if r_grid[i] < r_tumor else omega_b_normal
            dT_dr = (T_test[i+1] - T_test[i-1]) / (2 * dr)
            d2T_dr2 = (T_test[i+1] - 2*T_test[i] + T_test[i-1]) / dr**2
            q_cond = k_tissue * (d2T_dr2 + (2/r_grid[i]) * dT_dr)
            q_blood = omega * rho_b * cp_b * (T_a - T_test[i])
            q_rf = heat_source(r_grid[i], P_test, r_electrode)
            dT_dt = (q_cond + q_blood + q_rf) / (rho_tissue * cp_tissue)
            T_new[i] = T_test[i] + dt * dT_dt
        T_new[0] = T_test[1]
        T_new[-1] = T_initial
        T_test = T_new.copy()
    
    ablation_zone = r_grid[T_test >= 54]
    r_ablation_range.append(np.max(ablation_zone) * 1e3 if len(ablation_zone) > 0 else 0)

ax2.plot(P_range, r_ablation_range, 'ro-', linewidth=2, markersize=8)
ax2.set_xlabel('RF Power (W)', fontsize=11)
ax2.set_ylabel('Ablation Zone Radius (mm)', fontsize=11)
ax2.set_title('Ablation Zone vs Power', fontsize=12, fontweight='bold')
ax2.grid(True, alpha=0.3)

ax3 = axes[1, 0]
# 热损伤积分
CEM43 = np.zeros(nr)
for i in range(nr):
    if T[i] > 43:
        R_factor = 0.25 if T[i] >= 43 else 0.5
        CEM43[i] = R_factor**(43 - T[i]) * t_treatment / 60  # 分钟

ax3.semilogy(r_grid * 1e3, CEM43 + 0.1, 'g-', linewidth=2)
ax3.axhline(y=240, color='red', linestyle='--', alpha=0.7, label='Cell death threshold')
ax3.set_xlabel('Radius (mm)', fontsize=11)
ax3.set_ylabel('CEM43 (min)', fontsize=11)
ax3.set_title('Thermal Dose (CEM43)', fontsize=12, fontweight='bold')
ax3.legend()
ax3.grid(True, alpha=0.3)

ax4 = axes[1, 1]
# 血流灌注率影响
omega_range = np.linspace(0.0005, 0.005, 50)
r_ablation_omega = []
for omega_test in omega_range:
    T_test = np.ones(nr) * T_initial
    for n in range(nt):
        T_new = T_test.copy()
        for i in range(1, nr-1):
            dT_dr = (T_test[i+1] - T_test[i-1]) / (2 * dr)
            d2T_dr2 = (T_test[i+1] - 2*T_test[i] + T_test[i-1]) / dr**2
            q_cond = k_tissue * (d2T_dr2 + (2/r_grid[i]) * dT_dr)
            q_blood = omega_test * rho_b * cp_b * (T_a - T_test[i])
            q_rf = heat_source(r_grid[i], P_rf, r_electrode)
            dT_dt = (q_cond + q_blood + q_rf) / (rho_tissue * cp_tissue)
            T_new[i] = T_test[i] + dt * dT_dt
        T_new[0] = T_test[1]
        T_new[-1] = T_initial
        T_test = T_new.copy()
    
    ablation_zone = r_grid[T_test >= 54]
    r_ablation_omega.append(np.max(ablation_zone) * 1e3 if len(ablation_zone) > 0 else 0)

ax4.plot(omega_range * 1000, r_ablation_omega, 'm-', linewidth=2)
ax4.set_xlabel('Blood Perfusion Rate (mL/s/100g)', fontsize=11)
ax4.set_ylabel('Ablation Zone Radius (mm)', fontsize=11)
ax4.set_title('Effect of Blood Perfusion', fontsize=12, fontweight='bold')
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/case2_tumor_ablation.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 案例二完成")

# ==================== 案例三：心血管血流动力学仿真 ====================
print("\n" + "=" * 70)
print("案例三：心血管血流动力学仿真")
print("=" * 70)

# 颈动脉分叉几何参数
D_parent = 6e-3  # 颈总动脉直径 m
D_ica = 4e-3  # 颈内动脉直径 m
D_eca = 3e-3  # 颈外动脉直径 m
stenosis_degree = 0.6  # 狭窄程度（直径减少比例）
stenosis_length = 8e-3  # 狭窄段长度 m

# 血流参数
Q_mean = 5e-6  # 平均流量 m³/s (300 mL/min)
T_cardiac = 0.8  # 心动周期 s
mu_blood = 0.004  # 血液粘度 Pa·s
rho_blood = 1060  # 血液密度 kg/m³

print(f"\n血管几何:")
print(f"  颈总动脉直径: {D_parent*1e3:.1f} mm")
print(f"  颈内动脉直径: {D_ica*1e3:.1f} mm")
print(f"  狭窄程度: {stenosis_degree*100:.0f}%")
print(f"\n血流参数:")
print(f"  平均流量: {Q_mean*6e7:.0f} mL/min")
print(f"  心率: {60/T_cardiac:.0f} bpm")

# 简化的一维血流模型
t_flow = np.linspace(0, 2*T_cardiac, 200)  # 2个心动周期
Q_t = Q_mean * (1 + 0.5 * np.sin(2 * np.pi * t_flow / T_cardiac) - 
                0.2 * np.sin(4 * np.pi * t_flow / T_cardiac))

# 狭窄处流速
A_normal = np.pi * (D_parent/2)**2
A_stenosis = A_normal * (1 - stenosis_degree)**2
v_stenosis = Q_t / A_stenosis
v_normal = Q_t / A_normal

# 压力降（简化伯努利方程）
delta_P = 0.5 * rho_blood * (v_stenosis**2 - v_normal**2)

# 壁面剪切应力（简化）
tau_w_normal = 4 * mu_blood * v_normal / D_parent
tau_w_stenosis = 4 * mu_blood * v_stenosis / (D_parent * (1 - stenosis_degree))

print(f"\n血流动力学结果:")
print(f"  狭窄处峰值流速: {np.max(v_stenosis):.2f} m/s")
print(f"  峰值压力降: {np.max(delta_P)/133.32:.1f} mmHg")
print(f"  狭窄处峰值WSS: {np.max(tau_w_stenosis):.1f} Pa")

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

ax1 = axes[0, 0]
ax1.plot(t_flow / T_cardiac, Q_t * 6e7, 'b-', linewidth=2)
ax1.set_xlabel('Cardiac Cycle', fontsize=11)
ax1.set_ylabel('Flow Rate (mL/min)', fontsize=11)
ax1.set_title('Pulsatile Blood Flow', fontsize=12, fontweight='bold')
ax1.grid(True, alpha=0.3)

ax2 = axes[0, 1]
ax2.plot(t_flow / T_cardiac, v_stenosis, 'r-', linewidth=2, label='Stenosis throat')
ax2.plot(t_flow / T_cardiac, v_normal, 'b--', linewidth=2, label='Normal vessel')
ax2.set_xlabel('Cardiac Cycle', fontsize=11)
ax2.set_ylabel('Velocity (m/s)', fontsize=11)
ax2.set_title('Velocity Waveform', fontsize=12, fontweight='bold')
ax2.legend()
ax2.grid(True, alpha=0.3)

ax3 = axes[1, 0]
ax2_twin = ax3.twinx()
ax3.plot(t_flow / T_cardiac, delta_P / 133.32, 'g-', linewidth=2)
ax2_twin.plot(t_flow / T_cardiac, tau_w_stenosis, 'm--', linewidth=2)
ax3.set_xlabel('Cardiac Cycle', fontsize=11)
ax3.set_ylabel('Pressure Drop (mmHg)', fontsize=11, color='green')
ax2_twin.set_ylabel('WSS (Pa)', fontsize=11, color='magenta')
ax3.set_title('Hemodynamic Stress', fontsize=12, fontweight='bold')
ax3.tick_params(axis='y', labelcolor='green')
ax2_twin.tick_params(axis='y', labelcolor='magenta')
ax3.grid(True, alpha=0.3)

ax4 = axes[1, 1]
# 不同狭窄程度的影响
stenosis_range = np.linspace(0, 0.8, 50)
v_peak_range = []
delta_P_range = []
for sten in stenosis_range:
    A_sten = A_normal * (1 - sten)**2
    v_sten = np.max(Q_t) / A_sten
    v_norm = np.max(Q_t) / A_normal
    v_peak_range.append(v_sten)
    delta_P_range.append(0.5 * rho_blood * (v_sten**2 - v_norm**2) / 133.32)

ax4.plot(stenosis_range * 100, v_peak_range, 'b-', linewidth=2, label='Peak velocity')
ax4_twin = ax4.twinx()
ax4_twin.plot(stenosis_range * 100, delta_P_range, 'r--', linewidth=2, label='Pressure drop')
ax4.set_xlabel('Stenosis Degree (%)', fontsize=11)
ax4.set_ylabel('Peak Velocity (m/s)', fontsize=11, color='blue')
ax4_twin.set_ylabel('Pressure Drop (mmHg)', fontsize=11, color='red')
ax4.set_title('Effect of Stenosis Severity', fontsize=12, fontweight='bold')
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/case3_blood_flow.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 案例三完成")

# ==================== 案例四：低温冷冻保存过程模拟 ====================
print("\n" + "=" * 70)
print("案例四：低温冷冻保存过程模拟")
print("=" * 70)

# 冷冻保存参数
V_sample = 1e-6  # 样品体积 m³ (1 mL)
r_sample = (3 * V_sample / (4 * np.pi))**(1/3)  # 球形样品半径
T_initial = 20.0  # 初始温度 °C
T_liquid_nitrogen = -196.0  # 液氮温度 °C
h_conv = 100.0  # 对流换热系数 W/m²·K

# 物性参数（含10% DMSO的溶液）
k_ice = 2.2  # 冰的导热系数 W/m·K
k_water = 0.6  # 水的导热系数 W/m·K
rho_solution = 1000  # kg/m³
cp_solution = 3500  # J/kg·K
L_fusion = 250e3  # 相变潜热 J/kg（含DMSO降低）
T_melt_start = -5.0  # 开始相变温度 °C
T_melt_end = -8.0  # 结束相变温度 °C

# 细胞参数
L_p = 0.1e-12  # 水力传导率 m³/N/s
V_cell_0 = 1e-15  # 初始细胞体积 m³
A_cell = 1e-10  # 细胞表面积 m²
V_b = 0.3 * V_cell_0  # 不可压缩体积

print(f"\n冷冻保存参数:")
print(f"  样品体积: {V_sample*1e6:.1f} mL")
print(f"  样品半径: {r_sample*1e3:.2f} mm")
print(f"  冷冻保护剂: 10% DMSO")
print(f"  冷却介质: 液氮 ({T_liquid_nitrogen:.0f}°C)")

# 不同冷却速率下的温度响应
cooling_rates = [1, 5, 10, 50, 100]  # °C/min
fig, axes = plt.subplots(2, 2, figsize=(14, 12))

colors = ['blue', 'green', 'red', 'purple', 'orange']

ax1 = axes[0, 0]
for idx, rate in enumerate(cooling_rates):
    # 简化的一维瞬态传热
    t_cool = np.linspace(0, abs(T_initial - T_liquid_nitrogen) / rate * 60, 500)
    T_center = T_initial - rate * t_cool / 60
    T_center[T_center < T_liquid_nitrogen] = T_liquid_nitrogen
    
    # 相变平台效应（简化）
    for i in range(len(T_center)):
        if T_melt_end < T_center[i] < T_melt_start:
            T_center[i] = T_melt_start - (T_melt_start - T_melt_end) * (i / len(T_center))
    
    ax1.plot(t_cool / 60, T_center, color=colors[idx], linewidth=2, label=f'{rate}°C/min')

ax1.axhline(y=T_melt_start, color='gray', linestyle='--', alpha=0.5, label='Melting range')
ax1.axhline(y=T_melt_end, color='gray', linestyle='--', alpha=0.5)
ax1.set_xlabel('Time (min)', fontsize=11)
ax1.set_ylabel('Temperature (°C)', fontsize=11)
ax1.set_title('Cooling Curves at Different Rates', fontsize=12, fontweight='bold')
ax1.legend()
ax1.grid(True, alpha=0.3)

ax2 = axes[0, 1]
# 细胞体积变化（简化两参数模型）
for idx, rate in enumerate([1, 5, 10, 50]):
    t_sim = np.linspace(0, 30, 300)
    T_cell = max(20 - rate * t_sim, -196)
    
    # 简化细胞脱水模型
    V_cell = V_cell_0 * (1 - 0.6 * (1 - np.exp(-rate * t_sim / 10)))
    V_cell = np.maximum(V_cell, V_b)
    
    ax2.plot(t_sim, V_cell / V_cell_0 * 100, color=colors[idx], linewidth=2, label=f'{rate}°C/min')

ax2.axhline(y=30, color='red', linestyle='--', alpha=0.7, label='Minimum volume')
ax2.set_xlabel('Time (min)', fontsize=11)
ax2.set_ylabel('Relative Cell Volume (%)', fontsize=11)
ax2.set_title('Cell Dehydration During Freezing', fontsize=12, fontweight='bold')
ax2.legend()
ax2.grid(True, alpha=0.3)

ax3 = axes[1, 0]
# 细胞存活率 vs 冷却速率
# 经验公式：存在最佳冷却速率
rate_range = np.linspace(0.5, 200, 100)
survival = 100 * np.exp(-((np.log(rate_range) - np.log(10))**2) / (2 * 1.5**2))
survival[rate_range < 1] *= rate_range[rate_range < 1] / 1
survival[rate_range > 100] *= np.exp(-(rate_range[rate_range > 100] - 100) / 50)

ax3.plot(rate_range, survival, 'b-', linewidth=2)
ax3.axvline(x=10, color='red', linestyle='--', alpha=0.7, label='Optimal rate')
ax3.set_xlabel('Cooling Rate (°C/min)', fontsize=11)
ax3.set_ylabel('Cell Survival (%)', fontsize=11)
ax3.set_title('Cell Survival vs Cooling Rate', fontsize=12, fontweight='bold')
ax3.set_xscale('log')
ax3.legend()
ax3.grid(True, alpha=0.3)

ax4 = axes[1, 1]
# 温度梯度引起的应力
ax4_twin = ax4.twinx()
for idx, rate in enumerate([5, 10, 50]):
    t_stress = np.linspace(0, 10, 100)
    delta_T = rate * t_stress * 0.1  # 简化温度差
    thermal_stress = 1e6 * delta_T * 1e-6  # 简化应力计算 MPa
    ax4.plot(t_stress, delta_T, color=colors[idx], linewidth=2, linestyle='-', label=f'{rate}°C/min ΔT')
    ax4_twin.plot(t_stress, thermal_stress, color=colors[idx], linewidth=2, linestyle='--')

ax4.set_xlabel('Time (min)', fontsize=11)
ax4.set_ylabel('Temperature Difference (°C)', fontsize=11, color='blue')
ax4_twin.set_ylabel('Thermal Stress (MPa)', fontsize=11, color='red')
ax4.set_title('Thermal Stress During Freezing', fontsize=12, fontweight='bold')
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/case4_cryopreservation.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 案例四完成")

# ==================== 案例五：全身热应激响应分析 ====================
print("\n" + "=" * 70)
print("案例五：全身热应激响应分析")
print("=" * 70)

# 人体热调节参数
M_body = 70  # 体重 kg
A_body = 1.8  # 体表面积 m²
T_core_initial = 37.0  # 初始核心温度 °C
T_skin_initial = 33.0  # 初始皮肤温度 °C

# 环境参数
T_env = 40.0  # 环境温度 °C
RH = 60  # 相对湿度 %
v_air = 0.2  # 风速 m/s

# 代谢率
met_rate_rest = 100  # W
met_rate_exercise = 600  # W

# 热物性
cp_body = 3500  # J/kg·K
rho_body = 1000  # kg/m³

print(f"\n人体参数:")
print(f"  体重: {M_body:.0f} kg")
print(f"  体表面积: {A_body:.1f} m²")
print(f"\n环境条件:")
print(f"  环境温度: {T_env:.1f}°C")
print(f"  相对湿度: {RH:.0f}%")
print(f"  风速: {v_air:.1f} m/s")

# 简化的人体热平衡模型
def human_thermoregulation(t, y, q_met, T_env, RH, v_air):
    """
    人体热调节模型
    y = [T_core, T_skin]
    """
    T_core, T_skin = y
    
    # 皮肤-环境换热
    # 对流换热系数
    h_conv = 8.3 * (v_air ** 0.6)
    
    # 辐射换热
    eps_skin = 0.95
    T_env_K = T_env + 273.15
    T_skin_K = T_skin + 273.15
    q_rad = eps_skin * sigma_sb * (T_env_K**4 - T_skin_K**4)
    
    # 蒸发散热（出汗）
    # 最大蒸发率
    P_sat_skin = 0.6108 * np.exp(17.27 * T_skin / (T_skin + 237.3))  # kPa
    P_sat_env = 0.6108 * np.exp(17.27 * T_env / (T_env + 237.3))  # kPa
    P_vapor_env = P_sat_env * RH / 100
    
    # 出汗调节
    T_set_core = 37.0
    T_set_skin = 34.0
    sweat_rate_max = 1.0  # kg/m²/h
    sweat_rate = sweat_rate_max * max(0, (T_core - T_set_core) / 1.0) * max(0, (T_skin - T_set_skin) / 2.0)
    
    # 实际蒸发率受环境湿度限制
    h_evap = 16.5 * (v_air ** 0.6)
    q_evap_max = h_evap * (P_sat_skin - P_vapor_env) * 1000  # W/m²
    q_evap = min(sweat_rate * 2430e3 / 3600, q_evap_max)  # 限制蒸发
    
    # 对流换热
    q_conv = h_conv * (T_env - T_skin)
    
    # 总皮肤热损失
    q_skin_loss = q_conv + q_rad + q_evap
    
    # 核心-血液-皮肤换热
    blood_flow_skin = 0.5 + 5.0 * max(0, (T_core - T_set_core))  # L/m²/min
    q_blood = blood_flow_skin * 1e-3/60 * rho_b * cp_b * (T_core - T_skin)
    
    # 呼吸热损失
    q_resp = 0.0014 * q_met * (T_env - T_core) + 0.0173 * q_met * (5.87 - P_vapor_env)
    
    # 核心温度变化
    V_core = 0.6 * M_body / rho_body  # 核心体积 m³
    dT_core_dt = (q_met - q_resp - q_blood * A_body) / (rho_body * cp_body * V_core)
    
    # 皮肤温度变化
    V_skin = 0.1 * M_body / rho_body  # 皮肤层体积 m³
    dT_skin_dt = (q_blood * A_body - q_skin_loss * A_body) / (rho_body * cp_body * V_skin)
    
    return [dT_core_dt, dT_skin_dt]

# 仿真时间
t_heat = np.linspace(0, 3600, 1000)  # 1小时

# 静息状态
solution_rest = solve_ivp(
    lambda t, y: human_thermoregulation(t, y, met_rate_rest, T_env, RH, v_air),
    [0, 3600],
    [T_core_initial, T_skin_initial],
    t_eval=t_heat,
    method='RK45'
)

T_core_rest = solution_rest.y[0]
T_skin_rest = solution_rest.y[1]

# 运动状态
solution_exercise = solve_ivp(
    lambda t, y: human_thermoregulation(t, y, met_rate_exercise, T_env, RH, v_air),
    [0, 3600],
    [T_core_initial, T_skin_initial],
    t_eval=t_heat,
    method='RK45'
)

T_core_exercise = solution_exercise.y[0]
T_skin_exercise = solution_exercise.y[1]

print(f"\n热应激分析结果:")
print(f"  静息1小时后核心温度: {T_core_rest[-1]:.1f}°C")
print(f"  运动1小时后核心温度: {T_core_exercise[-1]:.1f}°C")
print(f"  运动状态体温上升速率: {(T_core_exercise[-1] - T_core_initial) / (t_heat[-1]/3600):.2f}°C/h")

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

ax1 = axes[0, 0]
ax1.plot(t_heat/60, T_core_rest, 'b-', linewidth=2, label='Rest')
ax1.plot(t_heat/60, T_core_exercise, 'r-', linewidth=2, label='Exercise')
ax1.axhline(y=38.5, color='orange', linestyle='--', alpha=0.7, label='Heat stress threshold')
ax1.axhline(y=40.0, color='red', linestyle='--', alpha=0.7, label='Danger threshold')
ax1.set_xlabel('Time (min)', fontsize=11)
ax1.set_ylabel('Core Temperature (°C)', fontsize=11)
ax1.set_title('Core Temperature Response', fontsize=12, fontweight='bold')
ax1.legend()
ax1.grid(True, alpha=0.3)

ax2 = axes[0, 1]
ax2.plot(t_heat/60, T_skin_rest, 'b-', linewidth=2, label='Rest')
ax2.plot(t_heat/60, T_skin_exercise, 'r-', linewidth=2, label='Exercise')
ax2.set_xlabel('Time (min)', fontsize=11)
ax2.set_ylabel('Skin Temperature (°C)', fontsize=11)
ax2.set_title('Skin Temperature Response', fontsize=12, fontweight='bold')
ax2.legend()
ax2.grid(True, alpha=0.3)

ax3 = axes[1, 0]
# 不同环境温度下的稳态核心温度
T_env_range = np.linspace(20, 50, 20)
T_core_steady = []
for T_env_test in T_env_range:
    sol = solve_ivp(
        lambda t, y: human_thermoregulation(t, y, met_rate_rest, T_env_test, RH, v_air),
        [0, 7200],
        [T_core_initial, T_skin_initial],
        method='RK45'
    )
    T_core_steady.append(sol.y[0][-1])

ax3.plot(T_env_range, T_core_steady, 'g-', linewidth=2)
ax3.axhline(y=37, color='gray', linestyle='--', alpha=0.5, label='Normal')
ax3.axhline(y=38.5, color='orange', linestyle='--', alpha=0.5, label='Heat stress')
ax3.set_xlabel('Environment Temperature (°C)', fontsize=11)
ax3.set_ylabel('Steady Core Temperature (°C)', fontsize=11)
ax3.set_title('Effect of Environment Temperature', fontsize=12, fontweight='bold')
ax3.legend()
ax3.grid(True, alpha=0.3)

ax4 = axes[1, 1]
# WBGT计算
T_wb_range = np.linspace(20, 35, 50)  # 湿球温度
T_globe = T_env + 5  # 黑球温度
WBGT = 0.7 * T_wb_range + 0.2 * T_globe + 0.1 * T_env

# 安全暴露时间（简化模型）
safe_time = 480 * np.exp(-(WBGT - 25) / 5)  # 分钟
safe_time[WBGT < 25] = 480

ax4.plot(WBGT, safe_time, 'm-', linewidth=2)
ax4.axvline(x=28, color='green', linestyle='--', alpha=0.7, label='Caution')
ax4.axvline(x=31, color='orange', linestyle='--', alpha=0.7, label='Extreme caution')
ax4.axvline(x=32, color='red', linestyle='--', alpha=0.7, label='Danger')
ax4.set_xlabel('WBGT (°C)', fontsize=11)
ax4.set_ylabel('Safe Exposure Time (min)', fontsize=11)
ax4.set_title('Heat Stress Safety Limits', fontsize=12, fontweight='bold')
ax4.set_yscale('log')
ax4.legend()
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/case5_heat_stress.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 案例五完成")

# ==================== 案例六：药物经皮渗透动力学 ====================
print("\n" + "=" * 70)
print("案例六：药物经皮渗透动力学")
print("=" * 70)

# 皮肤层参数
h_sc = 15e-6  # 角质层厚度 m
h_ve = 100e-6  # 活性表皮厚度 m
h_dermis = 1.5e-3  # 真皮厚度 m

# 药物参数
D_sc = 1e-12  # 角质层扩散系数 m²/s
D_ve = 1e-10  # 活性表皮扩散系数 m²/s
D_dermis = 5e-10  # 真皮扩散系数 m²/s
K_sc = 10  # 角质层/载体分配系数
K_ve = 1  # 活性表皮/载体分配系数

# 初始浓度
C0 = 1e-3  # 初始药物浓度 mol/m³ (1 mM)

print(f"\n皮肤层参数:")
print(f"  角质层厚度: {h_sc*1e6:.0f} μm")
print(f"  活性表皮厚度: {h_ve*1e6:.0f} μm")
print(f"  真皮厚度: {h_dermis*1e3:.1f} mm")
print(f"\n药物参数:")
print(f"  初始浓度: {C0*1e3:.1f} mM")
print(f"  角质层扩散系数: {D_sc:.2e} m²/s")

# 建立网格
nx_sc = 20
nx_ve = 30
nx_dermis = 50

x_sc = np.linspace(0, h_sc, nx_sc)
x_ve = np.linspace(h_sc, h_sc + h_ve, nx_ve)
x_dermis = np.linspace(h_sc + h_ve, h_sc + h_ve + h_dermis, nx_dermis)

x_total = np.concatenate([x_sc, x_ve[1:], x_dermis[1:]])
nx_total = len(x_total)

# 扩散系数分布
D = np.zeros(nx_total)
D[x_total <= h_sc] = D_sc
D[(x_total > h_sc) & (x_total <= h_sc + h_ve)] = D_ve
D[x_total > h_sc + h_ve] = D_dermis

# 稳态浓度分布（简化解析解）
# 角质层
x_sc_plot = np.linspace(0, h_sc, 100)
C_sc = C0 * (1 - x_sc_plot / h_sc)

# 活性表皮
x_ve_plot = np.linspace(h_sc, h_sc + h_ve, 100)
C_ve = C0 * K_sc * (1 - (x_ve_plot - h_sc) / (h_ve + h_sc * D_ve / D_sc * K_sc))

# 真皮
x_dermis_plot = np.linspace(h_sc + h_ve, h_sc + h_ve + h_dermis, 100)
C_dermis = C0 * K_sc * K_ve * np.exp(-(x_dermis_plot - h_sc - h_ve) / np.sqrt(D_dermis * 3600))

# 稳态通量
J_ss = D_sc * C0 * K_sc / h_sc  # mol/m²/s
J_ss_ug_cm2_h = J_ss * 1e6 * 3600  # μg/cm²/h (假设分子量300)

# 滞后时间
t_lag = (h_sc**2 / (6 * D_sc) + h_sc * h_ve / (2 * D_sc)) / 3600  # 小时

print(f"\n渗透结果:")
print(f"  稳态通量: {J_ss_ug_cm2_h:.2f} μg/cm²/h")
print(f"  滞后时间: {t_lag:.2f} h")

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

ax1 = axes[0, 0]
ax1.semilogy(x_sc_plot * 1e6, C_sc * 1e3, 'b-', linewidth=2, label='Stratum corneum')
ax1.semilogy(x_ve_plot * 1e6, C_ve * 1e3, 'g-', linewidth=2, label='Viable epidermis')
ax1.semilogy(x_dermis_plot * 1e6, C_dermis * 1e3, 'r-', linewidth=2, label='Dermis')
ax1.axvline(x=h_sc*1e6, color='gray', linestyle='--', alpha=0.5)
ax1.axvline(x=(h_sc+h_ve)*1e6, color='gray', linestyle='--', alpha=0.5)
ax1.set_xlabel('Depth (μm)', fontsize=11)
ax1.set_ylabel('Concentration (mM)', fontsize=11)
ax1.set_title('Steady-State Concentration Profile', fontsize=12, fontweight='bold')
ax1.legend()
ax1.grid(True, alpha=0.3)

ax2 = axes[0, 1]
# 累积渗透量 vs 时间
t_perm = np.linspace(0, 24, 200)  # 24小时
Q_cum = J_ss_ug_cm2_h * (t_perm - t_lag)
Q_cum[t_perm < t_lag] = 0
ax2.plot(t_perm, Q_cum, 'b-', linewidth=2)
ax2.axvline(x=t_lag, color='red', linestyle='--', alpha=0.7, label=f'Lag time = {t_lag:.1f} h')
ax2.set_xlabel('Time (h)', fontsize=11)
ax2.set_ylabel('Cumulative Permeation (μg/cm²)', fontsize=11)
ax2.set_title('Drug Permeation Kinetics', fontsize=12, fontweight='bold')
ax2.legend()
ax2.grid(True, alpha=0.3)

ax3 = axes[1, 0]
# 不同增强剂效果
enhancement_factors = [1, 5, 10, 20, 50]
J_enhanced = [J_ss_ug_cm2_h * ef for ef in enhancement_factors]
ax3.bar(range(len(enhancement_factors)), J_enhanced, color=['blue', 'green', 'orange', 'red', 'purple'])
ax3.set_xticks(range(len(enhancement_factors)))
ax3.set_xticklabels([f'{ef}x' for ef in enhancement_factors])
ax3.set_xlabel('Enhancement Factor', fontsize=11)
ax3.set_ylabel('Steady-State Flux (μg/cm²/h)', fontsize=11)
ax3.set_title('Effect of Penetration Enhancers', fontsize=12, fontweight='bold')
ax3.grid(True, alpha=0.3, axis='y')

ax4 = axes[1, 1]
# 不同分子量的扩散系数（Stokes-Einstein）
r_molecule = np.linspace(0.5e-9, 3e-9, 50)  # 分子半径 m
D_mw = 1.38e-23 * 310 / (6 * np.pi * 0.001 * r_molecule)  # m²/s
J_mw = D_mw * C0 * K_sc / h_sc * 1e6 * 3600  # μg/cm²/h

ax4.semilogy(r_molecule * 1e9, J_mw, 'm-', linewidth=2)
ax4.set_xlabel('Molecular Radius (nm)', fontsize=11)
ax4.set_ylabel('Steady-State Flux (μg/cm²/h)', fontsize=11)
ax4.set_title('Effect of Molecular Size', fontsize=12, fontweight='bold')
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/case6_transdermal.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 案例六完成")

print("\n" + "=" * 70)
print("所有案例仿真完成！")
print(f"结果保存至: {output_dir}")
print("=" * 70)

---