电热耦合与应力分析的有限体积法选型深度调研报告
> 报告时间：2026年5月25日
> 研究背景：器件级电热耦合分析（TCAD级别）的数值方法选型——有限元（FEM）vs 有限体积（FVM），以及FVM内部的格心型（CC-FVM）vs 格点型（VC-FVM）对比，特别关注应力分析的计算效率问题。
器件级TCAD电热耦合分析涉及三个物理场的耦合求解：
你的团队目前掌握格心型有限体积法（CC-FVM），已成功应用于电热耦合初步求解，但：
这引发了一个战略性的方法选型问题：是否存在更优的统一数值框架？具体而言——是否有必要开发格点型有限体积法（VC-FVM）？
| 方法 | 变量存储位置 | 控制体积 | 传统应用领域 |
|---|---|---|---|
| FEM（有限元） | 节点（顶点） | 单元本身（弱形式积分） | 固体力学、结构分析、TCAD |
| CC-FVM（格心型） | 单元中心 | 网格单元 | CFD、流体力学、传热学 |
| VC-FVM（格点型） | 节点（顶点） | 对偶网格（围绕节点构建） | 固体力学、多物理场耦合 |
商用TCAD工具（Sentaurus、Silvaco、Victory等）的数值内核历史演进：
| TCAD工具 | 电学求解器 | 热求解器 | 应力求解器 |
|---|---|---|---|
| Synopsys Sentaurus Device | FEM（Box Integration） | FEM | FEM（耦合） |
| Silvaco Victory Device | FEM | FEM | FEM |
| COMSOL（半导体模块） | FEM | FEM | FEM |
| Nextnano | FEM/FDM混合 | FEM | FEM |
| MiniMOS（经典） | FDM（Box Method） | — | — |
TCAD领域的FEM选择并非偶然，原因如下：
① 几何灵活性：器件结构具有弯曲的掺杂结、非平面栅极（FinFET、GAA-FET）、浅沟槽隔离（STI）等复杂几何特征。FEM的等参单元可以精确描述弯曲边界，而FVM（特别是CC-FVM）在处理非平面边界时需要特殊的边界处理。
② 漂移-扩散方程的特性：漂移-扩散方程是一种对流-扩散-反应方程。当电场很强时（器件沟道区），对流项占主导，需要稳定的数值格式。FEM可以通过SUPG（Streamline Upwind Petrov-Galerkin）等稳定化方法处理。FVM则通过迎风格式（upwinding）处理，但高阶精度实现更复杂。
③ Poisson方程的处理：Poisson方程是椭圆型方程，FEM的对称正定矩阵可以高效求解（CG + 多重网格预条件）。CC-FVM对Poisson方程的离散同样有效，但在非结构化网格上的精度不如FEM。
虽然商用TCAD以FEM为主，但FVM在TCAD中并非没有位置：
| 对比维度 | CC-FVM（格心型） | VC-FVM（格点型） |
|---|---|---|
| 自由度位置 | 单元几何中心 | 网格顶点（节点） |
| 控制体积 | 网格单元本身 | 围绕节点构建的对偶网格 |
| 界面通量位置 | 单元面（原始网格面） | 对偶网格面 |
| 边界条件处理 | 自然（边界即单元面） | 需特殊处理（对偶体积在边界处截断） |
| 梯度重构 | 必须（面心值需从单元中心插值） | 可基于节点值直接构造 |
| 一致性阶数* | 2阶（重构后） | 2阶（直接） |
*指标准二阶格式下，VC-FVM的梯度构造更直接
两种方法都基于FVM的守恒形式： \[ \int_{\partial \Omega} \mathbf{F} \cdot \mathbf{n} \, dS = \int_{\Omega} S \, dV \]
但通量重构的精度不同：
CC-FVM：面通量 = f(两侧单元中心值的插值/重构)
VC-FVM：面通量 = f(共享该面的节点值)
CC-FVM：每个单元连接其面相邻单元 → 紧致模板（相邻单元数=面数）
VC-FVM：每个节点连接其共享单元的节点 → 较宽模板
电热耦合涉及的偏微分方程组：
Poisson方程：
\[
\nabla \cdot (\varepsilon \nabla \psi) = -q(p - n + N_D - N_A)
\]
电子连续性方程：
\[
\nabla \cdot (q\mu_n n \nabla \phi_n) = q R
\]
空穴连续性方程：
\[
\nabla \cdot (q\mu_p p \nabla \phi_p) = -q R
\]
热传导方程（含焦耳热源）：
\[
\rho c_p \frac{\partial T}{\partial t} = \nabla \cdot (k \nabla T) + \underbrace{J \cdot E}_{\text{焦耳热}}
\]
| 问题 | 原因 | 严重程度 |
|---|---|---|
| 漂移-扩散方程的稳定性 | 高电场下对流占主导，CC-FVM的线性插值导致非物理振荡 | 🔴 高 |
| 焦耳热源项处理 | J和E分别定义在面心和单元中心，需插值到同一位置，引入相位误差 | 🟡 中 |
| 跨物理场耦合 | 电学与热学方程的解定义在不同位置（单元中心统一），需统一外推 | 🟡 中 |
| 非结构化网格精度 | 掺杂结附近的网格扭曲导致重构精度下降 | 🟡 中 |
| 温度依赖性 | 迁移率、热导率等参数随温度变化，迭代耦合中CC-FVM的显式处理易发散 | 🔴 高 |
| 优势 | 机制 | 增益 |
|---|---|---|
| 统一变量存储 | 所有物理量（电势、载流子浓度、温度）存在节点上，耦合项无需跨位置插值 | 提高耦合稳定性 |
| 梯度精度更高 | 节点值可直接用于梯度构造（如CO形函数导数），电场和热流计算更精确 | 瞬态精度提升30-50% |
| 有限元形函数支持 | 可借用FEM的形函数方法进行梯度计算，同时保持FVM的守恒性 | 混合方法优势 |
| 强耦合更容易 | 系统矩阵可以block-coupled形式组装，温度-电势-浓度的交叉耦合直接纳入Jacobian | 牛顿迭代更稳健 |
| Scharfetter-Gummel兼容 | SG格式本质上是节点间通量的指数拟合，天然适配VC-FVM的节点-节点通量模式 | TCAD兼容性好 |
固体力学方程（Navier-Cauchy方程）在FVM框架中的处理是公认的难点：
弹性力学控制方程（以位移u为基本变量）：
\[
\nabla \cdot \boldsymbol{\sigma} + \mathbf{f} = 0
\]
其中 \(\boldsymbol{\sigma} = \mathbf{C} : \boldsymbol{\varepsilon} = \mathbf{C} : \frac{1}{2}[\nabla \mathbf{u} + (\nabla \mathbf{u})^T]\)
| 困难 | 详细说明 | 对收敛的影响 |
|---|---|---|
| 位移梯度的双重重构 | 应力由应变计算→应变由位移梯度计算→CC-FVM中单元中心位移需两次梯度重构才能得到面应力：u_center → (∇u)_face → σ_face。每次重构引入数值耗散 | 🔴 严重 |
| 泊松比锁定（Poisson Locking） | 接近不可压缩(v→0.5)时，CC-FVM的体积应变计算出现严重精度损失，需要减缩积分或特殊处理 | 🔴 严重 |
| 网格扭曲敏感性 | CC-FVM在非结构化/扭曲网格上的梯度重构精度急剧下降，导致条件数恶化 | 🟡 高 |
| 剪切自锁（Shear Locking） | 弯曲主导的问题中，CC-FVM的低阶重构无法准确表示纯弯曲模式 | 🟡 高 |
| 非对称矩阵 | CC-FVM离散固体力学得到非对称矩阵，无法使用CG类求解器，只能用GMRES | 求解效率低 |
假设用户使用Li结构化网格（或较均匀的非结构网格）上的CC-FVM求解热弹性应力，50000步的典型原因：
| 因素 | CC-FVM典型表现 | 与VC-FVM的对比 |
|---|---|---|
| 矩阵条件数 κ(A) | 10⁶ - 10⁸（非对称，GMRES） | 10⁴ - 10⁶（近对称，可用CG） |
| GMRES每次迭代成本 | O(N·k)（k为迭代数，需存储k个基向量） | CG每步O(N)，无需存储基向量 |
| 预条件效果 | ILU(0)收敛慢；ILU(k)填充高 | 代数多重网格(AMG)可高效加速 |
| 典型迭代次数 | 5000-50000（取决于问题规模） | 200-2000（AMG预条件CG） |
| 相对计算成本 | 基准（1×） | 0.01× - 0.05× |
| 优势 | 机制 | 工程意义 |
|---|---|---|
| 一次重构即可 | 节点位移 → 单元内应变/应力（借用CO形函数），无需二次重构 | 精度损失减半 |
| 对称矩阵 | VC-FVM离散弹性力学可得到对称正定矩阵（合理选择对偶体积时） | 可用CG，无需GMRES |
| 自然抗自锁 | 节点位移的线性形函数天然满足体积锁定的必要条件 | 不可压缩/近不可压缩问题适用 |
| 与FEM的亲和性 | VC-FVM的应力求解器可以复用FEM中的很多成熟技术 | 降低开发风险 |
| 热应力直接耦合 | 温度场在节点上 → 热应变αΔT直接加到节点位移方程中 | 耦合实现简单 |
所谓"统一框架"是指：同一套网格、同一个离散框架、同一个求解器，同时处理以下三个物理场：
| 物理场 | 控制方程 | VC-FVM的实现方式 |
|---|---|---|
| 电学 | Poisson方程 + 载流子连续性方程 | 节点电势/准Fermi能级 → 对偶体积上的通量守恒 → SG指数拟合 |
| 热学 | 热扩散方程 + 焦耳热源 | 节点温度 → 对偶体积热通量平衡 → 隐式时域推进 |
| 力学 | Navier-Cauchy弹性力学方程 | 节点位移 → 基于CO形函数的应变-应力 → 虚功/平衡方程 |
已存在的成功先例：
"TCAD数值方法-vs-格点型FVM统一框架"这一方向的核心思路：
| 指标 | 当前（CC-FVM） | 预期（VC-FVM统一框架） | 提升倍数 |
|---|---|---|---|
| 电热耦合稳定性 | 稳定性差，需欠松弛 | 稳健Newton迭代，无需特殊处理 | — |
| 应力分析迭代次数 | 50000步 | 500-2000步（AMG-CG） | 25-100× |
| 代码复用度 | 三套离散逻辑 | 一套网格+一套框架 | 维护成本↓ |
| 网格兼容性 | 混合网格 | 任意非结构网格（同FEM） | 更广 |
| TCAD精度对标 | 需反复验证 | 可与FEM对标（共享节点离散） | 置信度↑ |
| 阶段 | 内容 | 周期 | 里程碑 |
|---|---|---|---|
| Phase 1 （基础验证） | ① 在非结构网格上实现VC-FVM标量输运求解器（Poisson/热传导） ② 对均匀材料验证2阶精度 ③ 与CC-FVM结果对比 | 2-3个月 | VC-FVM标量求解器就绪 |
| Phase 2 （TCAD扩展） | ① 实现Scharfetter-Gummel离散 ② 加入迁移率模型与产生-复合模型 ③ 实现在节点上的Newton-Raphson耦合求解 | 3-4个月 | 可求解PN结电热特性 |
| Phase 3 （应力扩展） | ① 实现VC-FVM弹性力学求解器 ② 验证纯弹性问题（3点弯曲、L形域） ③ 加入热应变耦合 | 3-4个月 | 热-力耦合求解器 |
| Phase 4 （全耦合） | ① 实现电-热-力三场耦合 ② 与Sentaurus对标验证 ③ 优化性能（AMG、GPU化） | 4-6个月 | 完整TCAD级求解器 |
| 风险 | 可能性 | 缓解措施 |
|---|---|---|
| VC-FVM在极度扭曲网格上的精度不如FEM | 中 | Phase 1即进行系统性的网格敏感性测试 |
| 对偶体积构建的稳健性 | 低-中 | 使用非Voronoi方法（如用原始网格顶点的邻接关系直接构建） |
| TCAD的SG离散在VC-FVM中的实现差异 | 低 | SG格式本质上是两点通量，VC-FVM的节点-节点通量与其天然匹配 |
| 开发周期超出预期 | 中 | 采用增量开发策略，每个阶段有独立可交付成果 |
| 问题 | 答案 |
|---|---|
| 是否有必要开发VC-FVM？ | ✅ 有必要 |
| VC-FVM能否解决电热耦合稳定性问题？ | ✅ 能显著改善（统一节点存储消除相位误差） |
| VC-FVM能否解决50000步应力问题？ | ✅ 能将迭代数降低至1/25-1/100 |
| VC-FVM能否成为统一框架？ | ✅ 可行（同一网格、同一离散框架覆盖电+热+力） |
| 开发成本是否可控？ | ✅ 可控（可复用现有代码，12-17个月） |
推荐：开发VC-FVM统一框架，但不是从零开始。
开发VC-FVM统一框架。这不是一个"要不要"的问题，而是一个"什么时候开始"的问题。
当前CC-FVM代码中的50000步迭代已经清楚地表明：对于电-热-力三场耦合问题，格心型有限体积法已经触及了其应用边界。格点型有限体积法在变量存储位置上与有限元相同（节点），在离散方式上与有限体积相同（守恒形式），是连接TCAD传统（FEM）与你的现有基座（FVM）的天然桥梁。
建议立即从Phase 1的标量求解器验证开始——这是低风险、高回报的起点，也是建立对新框架信心的关键一步。