做碲锌镉CdZnTeCZT这类化合物半导体晶体生长模拟的人早晚会面对同一个问题凝固过程根本不是一个纯传热问题也不只是一个流场问题更不是单纯画一条固液界面移动线就能交代的几何变化——它是变形几何、固体传热和流场三个过程搅在一起的整体行为。市面上能摸到的晶体生长仿真方案很多但真正肯用COMSOL Multiphysics把这三样东西同时放进一个模型里跑的人都会经历一段怎么搭都别扭、怎么调都发散的时期。这篇文章不谈教科书推导只讲我在CZT结晶过程模拟里实际做过的事情模型怎么搭、三场耦合怎么设、网格怎么救以及那些让COMSOL报错但帮助文档里根本不会写清楚的细节。给正在做晶体生长仿真、特别是准备用COMSOL动网格的朋友一份可以直接抄作业的参考。1. 为什么碲锌镉结晶模拟必须同时盯着传热、流动和几何变形1.1 碲锌镉晶体的实际应用与生长工艺背景碲锌镉CdZnTe常简写CZT是碲化镉基体中固溶一定比例锌形成的三元化合物半导体典型组分为Cd₀.₉Zn₀.₁Te。它之所以在半导体材料里地位特殊是因为平均原子序数高、密度接近6 g/cm³对X射线和伽马射线有很强的阻止能力是目前少有的能在室温下工作的核辐射探测器材料。同时它还是HgCdTe红外探测器的主流衬底探测器级CZT晶体的市场价值一直很高但生长难度也常年排在前列。生长CZT晶体的主流工艺是垂直布里奇曼法VBM和移动加热器法THM。VBM的操作逻辑很直观把多晶原料密封在高纯石英安瓿里整只安瓿以1~3 mm/h的速度在具有特定温度分布的炉膛中移动熔体从冷端开始定向凝固。这个速度慢到什么程度一根150 mm长的晶锭往往要连续生长两三天。工艺窗口又窄炉温梯度、安瓿移动速度、熔体对流强度只要有一点偏差轻则界面形状不对重则长出多晶或出现组分偏析。靠纯实验去试错成本和周期都扛不住所以数值模拟成了工艺前移的主要手段。而这类模拟最麻烦之处在于凝固界面本身在移动移动的界面改变固体域和液体域的比例进而改变热阻分布热阻分布又反过来决定界面位置的推进速度。整个问题天然是强耦合的想偷懒只算其中一个场结果基本没有参考价值。1.2 凝固过程中三个物理场的耦合逻辑先把物理链条理清楚。热量从高温炉区经安瓿壁传入熔体在固液界面处熔体凝固并释放潜热这部分热量通过已凝固的晶体向冷端导出。与此同时熔体内部因为轴向和径向都存在温度差密度分布不均匀浮力驱动出自然对流如果熔体有自由表面表面张力随温度变化还会叠加一层热毛细对流也就是Marangoni对流。流场一旦起来熔体里的热量输运就不再只是导热对流项会把温度场搅得面目全非。温度场一变固液界面的位置和形状也跟着变。界面所在位置决定了高导热固体域和低导热液体域各占多大空间这个几何比例又反过来改变整个计算域的传热路径。更麻烦的是界面移动本身有自己的时间尺度整个生长过程看似几天里缓慢推进但实际上每一步都接近稳态、又永远不完全稳态这种准静态推进对数值方法非常不友好。这三个场的耦合强度并不对称。传热和界面几何之间是强耦合界面位置由热平衡决定热阻分布又由界面位置决定流场与传热之间在格拉晓夫数Gr比较大时也极强流场与几何之间的耦合相对弱一些但在动网格框架下熔体域收缩时流场边界形态改变对涡结构和界面稳定性影响很直接。理解了这个主次关系后面设计求解策略才有依据。1.3 为什么选COMSOL Multiphysics而非自编程序早期做晶体生长模拟的人大多自编有限体积或有限差分程序。一套二维轴对称代码要从N-S方程写起再叠加热方程、固液界面追踪、潜热释放没有大半年很难跑出像样的结果而且每改一种工艺条件就要改一次边界条件维护成本极高。COMSOL Multiphysics的优势在于把这些物理场封装成了标准接口固体传热、层流、变形几何/移动网格都可以分别添加再用多物理场节点把它们绑定在一起初版模型一周内就能跑通。对晶体生长模拟来说COMSOL最有价值的其实是变形几何Deformed Geometry和移动网格Moving Mesh这一族工具。前者处理欧拉框架下计算域形状随时间变化的问题允许材料流过变形网格后者处理ALE框架下边界移动追踪问题网格跟着材料边界走。CZT凝固过程里这两种思路都有用武之地后面我会专门说怎么选。当然COMSOL也不是没有脾气。全耦合求解器面对强非线性时经常不收敛动网格畸变会毫无征兆地报错而且很多报错信息根本不说人话。接下来这几章重点就是怎么避开和解决这些问题。2. 模型搭建从材料参数到几何与边界条件2.1 几何与坐标系选择一上来就做全尺寸3D是新手最容易踩的坑。CZT安瓿直径几十毫米、高度两百毫米界面附近又需要微米级分辨率才能分辨温度梯度和流动结构3D模型随便一划就是几十万甚至上百万单元再加上非稳态动网格普通工作站根本转不动。我强烈建议从2D轴对称开始垂直布里奇曼的安瓿绕轴旋转几何和外场条件基本是轴对称的用2D轴对称建模就能抓住主部物理网格量省一个数量级。几何上我习惯建三个域熔体域初始为液相、晶体域初始为固相、石英安瓿壁。以直径D80 mm、总长200 mm的安瓿为例熔体初始高度取100 mm晶体初始高度取100 mm。在COMSOL里画成矩形分区即可。注意2D轴对称模型的宽度是半径而不是直径旋转对称轴放在左边。安瓿壁要单独画成一个域后面给材料属性和外壁换热条件都方便。这里有一个关键决定初始的熔体/晶体分界面到底怎么处理。如果后续用等效热容法两个域其实使用同一套温度相关材料属性界面位置完全由温度场决定不是几何上的真实边界。如果走显式移动网格追踪路线则两个域在几何建模上就是不同的域界面作为真实边界存在固液属性分别指定。我的经验是初版模型用等效热容法快速摸规律确认流场和温度场都稳定后再升到显式移动网格做精细化复核。2.2 CdZnTe材料热物性参数的选取与处理材料参数是CZT模拟最刁钻的一关。不同文献给出的热物性数据很分散尤其是液态热导率、动力粘度和表面张力温度系数这三个差别能达到两三倍。我通常先按主流文献值建一套基线参数再对最不确定的参数做灵敏度扫描绝不把一个孤立的数值当场真理。参数固相晶体液相熔体简要说明密度kg/m³约5900约5680液相比固相约低4%浮力项要配合参考密度使用比热容J/(kg·K)约160约200文献分散建议做灵敏度热导率W/(m·K)约1.0~1.5约1.0~2.0液态有效热导率受对流影响要区分表观值和本征值动力粘度Pa·s—约1.5×10⁻³~3×10⁻³温度相关常用指数律近似热膨胀系数1/K约5×10⁻⁶约2×10⁻⁴体积膨胀系数浮力项关键参数液相线温度/熔点℃约1092约1092Zn含量升高则液相线略升熔化焓kJ/kg约200约200等效热容法必须用表面张力温度系数N/(m·K)—约-2×10⁻⁴自由液面模型才需要别小看固液密度那4%的差别。在布里奇曼生长中熔体受热膨胀而晶体致密密度差既是自然对流的驱动力也是后期热应力的来源。做流动模拟时我在非等温流多物理场节点里勾选Boussinesq近似参考密度取ρ₀5680 kg/m³浮力项写作ρ₀·β·(T-T₀)·gβ取2×10⁻⁴ K⁻¹。这个量级下半径方向几十K的温差就足以驱动出显著的自然对流流场不能忽视。比热和热导率的温度依赖如果有实测数据强烈建议做成插值函数而不是常数。CZT固态热导率随温度变化很明显从高温降到室温能差好几倍。不过如果只关注凝固阶段的界面附近温度场也就是1050~1150℃这个窄窗口常数假设勉强能接受前提是初版模型。石英安瓿的参数也要认真给热导率约2 W/(m·K)比热约700 J/(kg·K)密度约2200 kg/m³。很多人把注意力全放在CZT上却忘了安瓿壁本身的热阻对界面温度场的影响。厚度几个毫米的石英壁在热流路径上造成的影响经常比想象中大得多。2.3 边界条件与初始条件设置炉温分布是整个模型最重要的边界条件。实测的Bridgman炉轴向温度分布是一条典型S型曲线高温区约1150℃低温区约1000℃中间绝热区温度梯度最大。在COMSOL里我习惯用定义→解析函数→分段定义轴向温度曲线T_outer(z)再把它作为安瓿外壁的温度边界条件。如果工艺是炉体移动而不是安瓿下降这个温度曲线在空间上要随时间平移写成T_outer(z - v_growth·t)即可。安瓿外壁不建议一上来就用表面对表面辐射。辐射项是T⁴非线性会直接把初次运行的稳态和瞬态求解器拖入泥潭。先用第一类给定温度或带对流系数的第三类边界把模型跑通之后再开辐射做精细化验证成功率会高很多。熔体接触安瓿壁的内壁设置为无滑移边界。如果保留自由上表面上表面设置为滑移边界同时施加表面张力温度系数产生的Marangoni剪切应力。这个应力不能忽略CZT熔体粘度虽然不算低但自由液面附近的Marangoni流动对局部温度场影响很明显。初始条件设置同样有讲究。我通常先单独求一次纯传热稳态也就是只开导热不开流场和动网格得到温度初值T_0流场初值设0压力初值设0。这比直接让一个非均匀温度场在耦合求解器里硬算要稳定得多。把整个模型初始温度都设为炉温曲线对应值也比均匀设成1120℃再让边界慢慢拉平衡要靠谱。2.4 潜热处理的注意点潜热是能量方程里最凶的源项。最常见的做法是等效热容法把相变潜热折算进一个等效比热容。核心公式就是cp_eq cp Lf * d(fL)/dT fL 0.5 * (1 - tanh((T - Tm) / dT))其中fL是液相分数Tm是熔点dT是人为设定的相变过渡半宽度。用光滑的tanh函数代替阶跃是为了避免两个网格节点之间发生温度跳变否则求解器会把时间步压缩到无法接受。dT的取值直接影响数值表现。ΔT取2~4 K是一个比较稳的起步区间。ΔT越小越接近真实相变面但数值越硬ΔT太大又会人为扩大两相区厚度导致算出来的界面位置模糊。CZT这种需要精细获取界面形状的场景我一般从ΔT3 K开始数值稳定后再逐步收紧。等效热容法还有一个副作用相变区间内等效比热会出现一个高尖峰让能量方程变成强非线性跟流场耦合时尤其敏感。所以正确的操作顺序是先固定网格、固定流场单独把温度场和相变收敛跑稳再开流场最后才引入动网格。分阶段逼近比一次性全耦合省太多时间。3. 三场耦合的实现细节COMSOL接口设置与求解策略3.1 物理场接口搭建具体操作上在添加物理场里选择以下接口传热→固体传热Heat Transfer in Solidsht作用域选熔体域、晶体域和石英壁流体流动→单相流→层流Laminar Flowspf作用域只选熔体域数学→变形网格→变形几何Deformed Geometrydg或移动网格Moving Meshale作用域选所有发生变形的域。如果采用等效热容网格变形的组合物理本质上还是材料域在变形固体传热和流场都定义在变形后的坐标里。这里有一个很容易漏掉的细节用移动网格时层流接口的物理模型记得勾选与网格速度相关的选项。COMSOL在移动网格下默认会用ALE形式的N-S方程自动带出网格速度修正项但如果你是从固定网格模型复制修改过来的一定要去物理模型面板确认一遍否则移动边界会在动量方程里产生虚假源项结果怎么调都怪。多物理场节点按这个思路搭一是非等温流动Nonisothermal Flow把层流与流体传热耦合它会在流场里加入浮力项、在能量方程里加入对流项二是流体-固体边界传热让熔体与晶体在界面上满足温度连续和热流连续。用等效热容法后这个条件基本由材料属性隐式完成但显式建一个节点方便后续查看界面热流。材料节点里固液两套属性分别指定到两个域。在等效热容法框架下更常规的做法是只在界面可能扫过的区域启用含潜热尖峰的等效比热晶体远端和熔体远端保持单一相这样能避免不必要的非线性。3.2 界面移动的驱动方式Stefan条件与变形几何这部分是核心中的核心。凝固界面移动的本构关系是Stefan条件物理含义一句话界面法向单位面积上释放的潜热等于固相侧导走的热流与液相侧补进来的热流之差。公式写出来是ρ_s · Lf · v_n k_s · ∂T_s/∂n - k_l · ∂T_l/∂n界面法向速度v_n不是人为指定的而是由界面两侧的温度梯度差决定。这是整个变形几何模型里最重要的生长法律。等效热容法的好处就在于不需要显式写出这个界面速度表达式因为能量方程已经把潜热消耗吞进去了温度场会自动让液相分数等值线移动。后处理时用fL0.5等值线提取界面位置再对时间求导就能得到生长速度。这是集中参数式的界面追踪实现难度低适合先把整条物理链路打通。如果必须显式追踪界面则在移动网格接口里把界面边界上的法向速度写成上面那个公式的形式。这个速度不是常数而是温度梯度的函数需要在边界变量或边界偏微分方程里承载。COMSOL支持用边界偏微分方程来解决这类问题但调试难度明显上了一个台阶符号写反、法向方向取错都会直接发散成垃圾结果。所以我不建议一上来就硬刚显式Stefan条件先拿等效热容法把物理趋势摸对再用移动网格做复核更省钱。这两种方法如果拿到的界面形状差在3%以内说明模型不仅稳定而且正确。差很多的话大概率不是数值方法的问题而是材料参数或边界条件在两个模型里没有保持一致。3.3 求解顺序与时间步控制求解器配置是COMSOL动网格仿真里最大的坑之一。默认的全耦合求解器在三个强非线性物理场同时启动时十次有九次不收敛。我的操作顺序是第一步关闭流场和变形几何单独解稳态传热拿到初始温度场。 第二步打开层流固定网格做瞬态非等温流让流动达到准稳态。 第三步打开变形几何或移动网格保留上一步结果作为初值。 第四步改用分离式Segregated求解器把压力-速度放一组温度场单独一组网格位移单独一组。每个子问题都是相对成熟的小框架整体反而比强耦合更快更稳。时间步进用BDF向后差分公式同时设置最大时间步。凝固速率按1 μm/s估算如果界面附近最小网格尺寸是0.1 mm界面穿过一个单元需要100秒。时间步要是超过几十秒界面状态可能一下跳过好几个单元网格变形会突然剧烈起来。我习惯给BDF设初始步长不超过0.1 s最大步长不超过30~60 s让它在这个框架内做自适应。容差建议收紧到1e-4而不是默认的1e-3。默认容差对瞬态多物理场往往偏松某些时刻已经偏离物理规律但求解器还在闷头跑。容差收紧后界面位置的锯齿明显减少虽然会慢几倍但翻车率低一个数量级。另外一定要设停止条件保护比如界面移动超过设定值或出现非物理温度就提前终止避免烧掉一整天算出一堆垃圾。4. 网格畸变、数值振荡与收敛性踩坑记录4.1 变形几何中网格畸变的典型表现与成因动网格模型报错最吓人的一句话是网格反移或倒置单元有时候直接报negative Jacobian。看到这些别急着改网格先回去看几何为什么变成那样。CZT布里奇曼生长里网格畸变的主要导火索有三个界面推进导致熔体域轴向变薄温度梯度大、界面处网格被剧烈拉伸流体连续性方程与动网格速度耦合产生压力尖峰。界面向熔体方向推进几十毫米后熔体域最初100 mm的高度最后可能只剩几个毫米网格单元一边被压缩一边被拉长长宽比飙升。COMSOL变形几何内置了Laplace和Winslow等网格平滑算法但能力有限网格相对厚度低于约3~5层时界面附近的压力离散就会开始病态。我的习惯是在整个模型里给变形几何设置一个初始网格位移为零的初始条件然后在所有外边界用滑移位移而不是固定位移让网格能沿安瓿壁滑移来容纳变形。很多人一出问题就把边界全设成固定结果界面动一下整片网格像拉一块硬板似的全被拽弯倒置单元必然出现。4.2 一次不收敛问题的完整定位过程分享一次真实的排查过程可以当作模板用。那轮参数扫描里基线工况能跑通但把冷端冷却速率从0.5 K/min改成1.5 K/min后t≈3000 s时求解中止报错信息是最后步长小于最小步长瞬态求解失败。我的第一步不是调网格而是关掉流场只保留传热变形几何。结果这个组合顺利跑过了3000 s说明问题不在动网格本身。第二步关掉变形几何保留传热流场也能跑说明问题不在流热耦合。问题就这样被锁定在流场动网格的联合效应上。第三步打开流场和动网格把时间步固定在一个极小值跑几步看中间输出。果然在熔体域厚度压缩到大约还剩8 mm时界面附近压力出现高频振荡紧接着最低层的网格单元开始反转。原因清楚了熔体域变薄、界面又持续释放潜热流场在越来越窄的通道里被挤压力-速度耦合在动网格框架下没有足够阻尼振荡超过阈值后网格直接塌掉。解法分三路第一在熔体域顶部与气相接触的边界设置一个压力约束给封闭空间一个泄压通道第二在界面边界上对网格位移施加额外的切向滑移自由度让网格能顺着壁面滑而不被压扁第三检查层流接口的稳定化设置COMSOL默认的流线扩散没有被自定义设置偷关掉就保留。改完后同一个扫描点顺利跑通后面更极端的参数点只是变慢没有再崩。4.3 稳定收敛的几招储备整理几条我在CZT模拟里沉淀下来的保命经验。设置网格重剖分触发条件。COMSOL变形几何可以在网格质量低于阈值时自动重新划分网格这个功能很诱人但要慎用任何重剖分都会带来历史场变量的插值误差Crank-Nicolson类格式对插值尤其敏感。实在要用就设一个最大单元扭曲度阈值让重剖分尽量少触发。给界面法向速度加限制。把界面推进速度限制在一个物理合理范围比如0~5 μm/s防止初期微小数值振荡被放大成不合理的界面位移。这个限制在COMSOL变量定义里用min和max函数就能实现成本极低但效果明显。用参数扫描接力逼近目标工况。把目标工艺参数拆成多个子工况从简单到复杂逐个计算每一步以上一步的解作为初值。这比直接怼一个极端工况稳定得多本质上是用前一个解当了热启动。定期检查能量守恒。每跑完一段用COMSOL的报告功能统计界面释放潜热总量和通过边界累计导热量两者偏差超过5%就优先怀疑潜热过渡区太宽或网格分辨率不足。能量账不平后面所有结果都不能信。5. 结果后处理与物理可信度检验5.1 界面形态演化与生长速率判据等效热容法的直接输出是温度场和液相分数场。界面位置可以通过派生值来定义在每条轴向采样线上找fL0.5的z坐标。把界面各径向位置的z坐标与中心位置的z坐标做差就是界面凸度。经验值上界面中心略微凸向熔体1~3 mm是比较理想的工艺状态这个形状有利于把杂质推向晶体外周而不是滞留在中心对探测器级CZT尤其重要。生长速率后处理也很直接记录中心点界面z坐标随时间的变化取斜率就是局部生长速度。理想情况下经过一段起步过渡后生长速度应接近安瓿移动速度比如1 mm/h。如果差太多说明热平衡没有建立或者冷却速率设置太快界面在追赶温度场。这种追赶在实验中表现为界面过冷或局部回熔模拟的价值就在这里——不用动炉子就能提前摸清工艺参数区间。把界面推进速率、熔体最大流速和界面凸度三者放在同一张时间曲线上观察几乎一眼就能看出某组参数是否和谐。如果流速随时间剧烈波动、界面推进速度跟着抖那多半是对流进入了非稳态区间需要进一步检查格拉晓夫数。5.2 流场与温度场的交叉验证解出流场后先别急着截图先做无量纲判断。格拉晓夫数Gr g·β·ΔT·H³/ν²以熔体高度H50 mm、温差ΔT30 K、体积膨胀系数2×10⁻⁴、运动粘度5.3×10⁻⁷ m²/s估算Gr大概在10⁶量级这已经是非常确定的自然对流不是可以忽略的次级效应。如果你模拟出来的流速小到几乎看不见第一时间怀疑Boussinesq项的符号或参考密度设置错了。温度场和流场的交叉检查我从三个维度做。一是界面附近的等温线是否被对流拉变形纯导热情况下等温线应该平滑地穿过界面对流则会带来明显扭曲。二是对称轴附近是否出现规则环流环流的方向应该与浮力方向一致热区上升、冷区下降。三是通过界面的总热流是否等于潜热消耗两侧导热的收支差热流账不平基本就是潜热处理的问题。我还喜欢看一个细节熔体顶部或自由液面附近是否存在第二个反向涡。垂直布里奇曼里只出现一个主环流是常见情况但如果工艺参数比较极端可能出现上下层双涡结构这种结构对溶质分布影响很大。做参数扫描时可以专门输出涡心位置随参数的变化这比单看界面形貌更早暴露潜在不稳定。5.3 参数扫描与工艺优化视角等模型跑通COMSOL的参数化扫描就变成了效率工具。用辅助扫描把炉温梯度、安瓿下降速度、熔体初始高度作为参数列表一次性提交多个工况。每个工况都把上一个工况的解作为初值相当于给每个极端参数点搭了一座桥收敛率大幅提升。扫描结果的判读要回到物理。提高炉温梯度会让界面更凸向熔体但伴随而来的是更大的热应力降低下降速度会让界面更接近稳态、流场更平静但生长周期被拉长。用派生值表把每个参数点的界面凸度、最大流速、轴向温度梯度、界面过冷度整理成表基本就能画出工艺窗口的轮廓。有一点必须提醒模拟得到的趋势比具体数值更可靠。比如梯度增加10 K/cm界面凸度增大约0.5 mm这样的趋势对工艺调整有直接指导意义但梯度30 K/cm时界面凸度是1.52 mm这种绝对值只能当参考因为材料参数和边界条件总有不精确的地方。我在对外输出报告时永远把结果写成一个区间而不是一个点这是做仿真的人最基本的专业素养。6. 继续做下去的几个方向与我的实际体会6.1 从二维到三维、从热流到场2D轴对称模型能回答的问题已经很多但它天然过滤掉了周向非对称扰动。真实安瓿在炉膛里不是完美对中炉体也在缓慢旋转这些因素会在晶体里留下螺旋状杂质条纹。想研究这些就必须上3D模型并考虑安瓿旋转无论是旋转坐标系还是滑移网格复杂度都是2D的几倍建议在2D物理彻底摸清后再碰。更值得优先做的是把第三个场加进来——溶质场。CdZnTe熔体中的Zn和Cd在凝固时会发生再分配造成组分偏析。用稀物质传递接口耦合浓度场配合界面处的分凝系数可以定量预测沿晶锭方向的电阻率均匀性。这个结果对探测器材料来说比纯温度分布更接近用户真正关心的指标。6.2 我的体会回头做CZT结晶模拟最有价值的时刻不是模型第一次跑通而是拿模拟得到的温度曲线和实测热电偶数据对比、发现两者相差在2℃以内那次。从那以后我确定了一件事材料参数校准和边界条件确认的优先级高于任何高级物理算法。很多人一上来就追求相场、追求3D、追求辐射换热结果连基线传热都没校出来这是最划不来的。我的最后一条建议是先用最朴素的方法把能量守恒做对再慢慢往模型里加细节。变形几何、固体传热、流场三者的耦合本质上是把晶体生长最核心的物理搬进电脑里。基础骨架稳了后面加什么细节都只是锦上添花。
做碲锌镉CdZnTeCZT这类化合物半导体晶体生长模拟的人早晚会面对同一个问题凝固过程根本不是一个纯传热问题也不只是一个流场问题更不是单纯画一条固液界面移动线就能交代的几何变化——它是变形几何、固体传热和流场三个过程搅在一起的整体行为。市面上能摸到的晶体生长仿真方案很多但真正肯用COMSOL Multiphysics把这三样东西同时放进一个模型里跑的人都会经历一段怎么搭都别扭、怎么调都发散的时期。这篇文章不谈教科书推导只讲我在CZT结晶过程模拟里实际做过的事情模型怎么搭、三场耦合怎么设、网格怎么救以及那些让COMSOL报错但帮助文档里根本不会写清楚的细节。给正在做晶体生长仿真、特别是准备用COMSOL动网格的朋友一份可以直接抄作业的参考。1. 为什么碲锌镉结晶模拟必须同时盯着传热、流动和几何变形1.1 碲锌镉晶体的实际应用与生长工艺背景碲锌镉CdZnTe常简写CZT是碲化镉基体中固溶一定比例锌形成的三元化合物半导体典型组分为Cd₀.₉Zn₀.₁Te。它之所以在半导体材料里地位特殊是因为平均原子序数高、密度接近6 g/cm³对X射线和伽马射线有很强的阻止能力是目前少有的能在室温下工作的核辐射探测器材料。同时它还是HgCdTe红外探测器的主流衬底探测器级CZT晶体的市场价值一直很高但生长难度也常年排在前列。生长CZT晶体的主流工艺是垂直布里奇曼法VBM和移动加热器法THM。VBM的操作逻辑很直观把多晶原料密封在高纯石英安瓿里整只安瓿以1~3 mm/h的速度在具有特定温度分布的炉膛中移动熔体从冷端开始定向凝固。这个速度慢到什么程度一根150 mm长的晶锭往往要连续生长两三天。工艺窗口又窄炉温梯度、安瓿移动速度、熔体对流强度只要有一点偏差轻则界面形状不对重则长出多晶或出现组分偏析。靠纯实验去试错成本和周期都扛不住所以数值模拟成了工艺前移的主要手段。而这类模拟最麻烦之处在于凝固界面本身在移动移动的界面改变固体域和液体域的比例进而改变热阻分布热阻分布又反过来决定界面位置的推进速度。整个问题天然是强耦合的想偷懒只算其中一个场结果基本没有参考价值。1.2 凝固过程中三个物理场的耦合逻辑先把物理链条理清楚。热量从高温炉区经安瓿壁传入熔体在固液界面处熔体凝固并释放潜热这部分热量通过已凝固的晶体向冷端导出。与此同时熔体内部因为轴向和径向都存在温度差密度分布不均匀浮力驱动出自然对流如果熔体有自由表面表面张力随温度变化还会叠加一层热毛细对流也就是Marangoni对流。流场一旦起来熔体里的热量输运就不再只是导热对流项会把温度场搅得面目全非。温度场一变固液界面的位置和形状也跟着变。界面所在位置决定了高导热固体域和低导热液体域各占多大空间这个几何比例又反过来改变整个计算域的传热路径。更麻烦的是界面移动本身有自己的时间尺度整个生长过程看似几天里缓慢推进但实际上每一步都接近稳态、又永远不完全稳态这种准静态推进对数值方法非常不友好。这三个场的耦合强度并不对称。传热和界面几何之间是强耦合界面位置由热平衡决定热阻分布又由界面位置决定流场与传热之间在格拉晓夫数Gr比较大时也极强流场与几何之间的耦合相对弱一些但在动网格框架下熔体域收缩时流场边界形态改变对涡结构和界面稳定性影响很直接。理解了这个主次关系后面设计求解策略才有依据。1.3 为什么选COMSOL Multiphysics而非自编程序早期做晶体生长模拟的人大多自编有限体积或有限差分程序。一套二维轴对称代码要从N-S方程写起再叠加热方程、固液界面追踪、潜热释放没有大半年很难跑出像样的结果而且每改一种工艺条件就要改一次边界条件维护成本极高。COMSOL Multiphysics的优势在于把这些物理场封装成了标准接口固体传热、层流、变形几何/移动网格都可以分别添加再用多物理场节点把它们绑定在一起初版模型一周内就能跑通。对晶体生长模拟来说COMSOL最有价值的其实是变形几何Deformed Geometry和移动网格Moving Mesh这一族工具。前者处理欧拉框架下计算域形状随时间变化的问题允许材料流过变形网格后者处理ALE框架下边界移动追踪问题网格跟着材料边界走。CZT凝固过程里这两种思路都有用武之地后面我会专门说怎么选。当然COMSOL也不是没有脾气。全耦合求解器面对强非线性时经常不收敛动网格畸变会毫无征兆地报错而且很多报错信息根本不说人话。接下来这几章重点就是怎么避开和解决这些问题。2. 模型搭建从材料参数到几何与边界条件2.1 几何与坐标系选择一上来就做全尺寸3D是新手最容易踩的坑。CZT安瓿直径几十毫米、高度两百毫米界面附近又需要微米级分辨率才能分辨温度梯度和流动结构3D模型随便一划就是几十万甚至上百万单元再加上非稳态动网格普通工作站根本转不动。我强烈建议从2D轴对称开始垂直布里奇曼的安瓿绕轴旋转几何和外场条件基本是轴对称的用2D轴对称建模就能抓住主部物理网格量省一个数量级。几何上我习惯建三个域熔体域初始为液相、晶体域初始为固相、石英安瓿壁。以直径D80 mm、总长200 mm的安瓿为例熔体初始高度取100 mm晶体初始高度取100 mm。在COMSOL里画成矩形分区即可。注意2D轴对称模型的宽度是半径而不是直径旋转对称轴放在左边。安瓿壁要单独画成一个域后面给材料属性和外壁换热条件都方便。这里有一个关键决定初始的熔体/晶体分界面到底怎么处理。如果后续用等效热容法两个域其实使用同一套温度相关材料属性界面位置完全由温度场决定不是几何上的真实边界。如果走显式移动网格追踪路线则两个域在几何建模上就是不同的域界面作为真实边界存在固液属性分别指定。我的经验是初版模型用等效热容法快速摸规律确认流场和温度场都稳定后再升到显式移动网格做精细化复核。2.2 CdZnTe材料热物性参数的选取与处理材料参数是CZT模拟最刁钻的一关。不同文献给出的热物性数据很分散尤其是液态热导率、动力粘度和表面张力温度系数这三个差别能达到两三倍。我通常先按主流文献值建一套基线参数再对最不确定的参数做灵敏度扫描绝不把一个孤立的数值当场真理。参数固相晶体液相熔体简要说明密度kg/m³约5900约5680液相比固相约低4%浮力项要配合参考密度使用比热容J/(kg·K)约160约200文献分散建议做灵敏度热导率W/(m·K)约1.0~1.5约1.0~2.0液态有效热导率受对流影响要区分表观值和本征值动力粘度Pa·s—约1.5×10⁻³~3×10⁻³温度相关常用指数律近似热膨胀系数1/K约5×10⁻⁶约2×10⁻⁴体积膨胀系数浮力项关键参数液相线温度/熔点℃约1092约1092Zn含量升高则液相线略升熔化焓kJ/kg约200约200等效热容法必须用表面张力温度系数N/(m·K)—约-2×10⁻⁴自由液面模型才需要别小看固液密度那4%的差别。在布里奇曼生长中熔体受热膨胀而晶体致密密度差既是自然对流的驱动力也是后期热应力的来源。做流动模拟时我在非等温流多物理场节点里勾选Boussinesq近似参考密度取ρ₀5680 kg/m³浮力项写作ρ₀·β·(T-T₀)·gβ取2×10⁻⁴ K⁻¹。这个量级下半径方向几十K的温差就足以驱动出显著的自然对流流场不能忽视。比热和热导率的温度依赖如果有实测数据强烈建议做成插值函数而不是常数。CZT固态热导率随温度变化很明显从高温降到室温能差好几倍。不过如果只关注凝固阶段的界面附近温度场也就是1050~1150℃这个窄窗口常数假设勉强能接受前提是初版模型。石英安瓿的参数也要认真给热导率约2 W/(m·K)比热约700 J/(kg·K)密度约2200 kg/m³。很多人把注意力全放在CZT上却忘了安瓿壁本身的热阻对界面温度场的影响。厚度几个毫米的石英壁在热流路径上造成的影响经常比想象中大得多。2.3 边界条件与初始条件设置炉温分布是整个模型最重要的边界条件。实测的Bridgman炉轴向温度分布是一条典型S型曲线高温区约1150℃低温区约1000℃中间绝热区温度梯度最大。在COMSOL里我习惯用定义→解析函数→分段定义轴向温度曲线T_outer(z)再把它作为安瓿外壁的温度边界条件。如果工艺是炉体移动而不是安瓿下降这个温度曲线在空间上要随时间平移写成T_outer(z - v_growth·t)即可。安瓿外壁不建议一上来就用表面对表面辐射。辐射项是T⁴非线性会直接把初次运行的稳态和瞬态求解器拖入泥潭。先用第一类给定温度或带对流系数的第三类边界把模型跑通之后再开辐射做精细化验证成功率会高很多。熔体接触安瓿壁的内壁设置为无滑移边界。如果保留自由上表面上表面设置为滑移边界同时施加表面张力温度系数产生的Marangoni剪切应力。这个应力不能忽略CZT熔体粘度虽然不算低但自由液面附近的Marangoni流动对局部温度场影响很明显。初始条件设置同样有讲究。我通常先单独求一次纯传热稳态也就是只开导热不开流场和动网格得到温度初值T_0流场初值设0压力初值设0。这比直接让一个非均匀温度场在耦合求解器里硬算要稳定得多。把整个模型初始温度都设为炉温曲线对应值也比均匀设成1120℃再让边界慢慢拉平衡要靠谱。2.4 潜热处理的注意点潜热是能量方程里最凶的源项。最常见的做法是等效热容法把相变潜热折算进一个等效比热容。核心公式就是cp_eq cp Lf * d(fL)/dT fL 0.5 * (1 - tanh((T - Tm) / dT))其中fL是液相分数Tm是熔点dT是人为设定的相变过渡半宽度。用光滑的tanh函数代替阶跃是为了避免两个网格节点之间发生温度跳变否则求解器会把时间步压缩到无法接受。dT的取值直接影响数值表现。ΔT取2~4 K是一个比较稳的起步区间。ΔT越小越接近真实相变面但数值越硬ΔT太大又会人为扩大两相区厚度导致算出来的界面位置模糊。CZT这种需要精细获取界面形状的场景我一般从ΔT3 K开始数值稳定后再逐步收紧。等效热容法还有一个副作用相变区间内等效比热会出现一个高尖峰让能量方程变成强非线性跟流场耦合时尤其敏感。所以正确的操作顺序是先固定网格、固定流场单独把温度场和相变收敛跑稳再开流场最后才引入动网格。分阶段逼近比一次性全耦合省太多时间。3. 三场耦合的实现细节COMSOL接口设置与求解策略3.1 物理场接口搭建具体操作上在添加物理场里选择以下接口传热→固体传热Heat Transfer in Solidsht作用域选熔体域、晶体域和石英壁流体流动→单相流→层流Laminar Flowspf作用域只选熔体域数学→变形网格→变形几何Deformed Geometrydg或移动网格Moving Meshale作用域选所有发生变形的域。如果采用等效热容网格变形的组合物理本质上还是材料域在变形固体传热和流场都定义在变形后的坐标里。这里有一个很容易漏掉的细节用移动网格时层流接口的物理模型记得勾选与网格速度相关的选项。COMSOL在移动网格下默认会用ALE形式的N-S方程自动带出网格速度修正项但如果你是从固定网格模型复制修改过来的一定要去物理模型面板确认一遍否则移动边界会在动量方程里产生虚假源项结果怎么调都怪。多物理场节点按这个思路搭一是非等温流动Nonisothermal Flow把层流与流体传热耦合它会在流场里加入浮力项、在能量方程里加入对流项二是流体-固体边界传热让熔体与晶体在界面上满足温度连续和热流连续。用等效热容法后这个条件基本由材料属性隐式完成但显式建一个节点方便后续查看界面热流。材料节点里固液两套属性分别指定到两个域。在等效热容法框架下更常规的做法是只在界面可能扫过的区域启用含潜热尖峰的等效比热晶体远端和熔体远端保持单一相这样能避免不必要的非线性。3.2 界面移动的驱动方式Stefan条件与变形几何这部分是核心中的核心。凝固界面移动的本构关系是Stefan条件物理含义一句话界面法向单位面积上释放的潜热等于固相侧导走的热流与液相侧补进来的热流之差。公式写出来是ρ_s · Lf · v_n k_s · ∂T_s/∂n - k_l · ∂T_l/∂n界面法向速度v_n不是人为指定的而是由界面两侧的温度梯度差决定。这是整个变形几何模型里最重要的生长法律。等效热容法的好处就在于不需要显式写出这个界面速度表达式因为能量方程已经把潜热消耗吞进去了温度场会自动让液相分数等值线移动。后处理时用fL0.5等值线提取界面位置再对时间求导就能得到生长速度。这是集中参数式的界面追踪实现难度低适合先把整条物理链路打通。如果必须显式追踪界面则在移动网格接口里把界面边界上的法向速度写成上面那个公式的形式。这个速度不是常数而是温度梯度的函数需要在边界变量或边界偏微分方程里承载。COMSOL支持用边界偏微分方程来解决这类问题但调试难度明显上了一个台阶符号写反、法向方向取错都会直接发散成垃圾结果。所以我不建议一上来就硬刚显式Stefan条件先拿等效热容法把物理趋势摸对再用移动网格做复核更省钱。这两种方法如果拿到的界面形状差在3%以内说明模型不仅稳定而且正确。差很多的话大概率不是数值方法的问题而是材料参数或边界条件在两个模型里没有保持一致。3.3 求解顺序与时间步控制求解器配置是COMSOL动网格仿真里最大的坑之一。默认的全耦合求解器在三个强非线性物理场同时启动时十次有九次不收敛。我的操作顺序是第一步关闭流场和变形几何单独解稳态传热拿到初始温度场。 第二步打开层流固定网格做瞬态非等温流让流动达到准稳态。 第三步打开变形几何或移动网格保留上一步结果作为初值。 第四步改用分离式Segregated求解器把压力-速度放一组温度场单独一组网格位移单独一组。每个子问题都是相对成熟的小框架整体反而比强耦合更快更稳。时间步进用BDF向后差分公式同时设置最大时间步。凝固速率按1 μm/s估算如果界面附近最小网格尺寸是0.1 mm界面穿过一个单元需要100秒。时间步要是超过几十秒界面状态可能一下跳过好几个单元网格变形会突然剧烈起来。我习惯给BDF设初始步长不超过0.1 s最大步长不超过30~60 s让它在这个框架内做自适应。容差建议收紧到1e-4而不是默认的1e-3。默认容差对瞬态多物理场往往偏松某些时刻已经偏离物理规律但求解器还在闷头跑。容差收紧后界面位置的锯齿明显减少虽然会慢几倍但翻车率低一个数量级。另外一定要设停止条件保护比如界面移动超过设定值或出现非物理温度就提前终止避免烧掉一整天算出一堆垃圾。4. 网格畸变、数值振荡与收敛性踩坑记录4.1 变形几何中网格畸变的典型表现与成因动网格模型报错最吓人的一句话是网格反移或倒置单元有时候直接报negative Jacobian。看到这些别急着改网格先回去看几何为什么变成那样。CZT布里奇曼生长里网格畸变的主要导火索有三个界面推进导致熔体域轴向变薄温度梯度大、界面处网格被剧烈拉伸流体连续性方程与动网格速度耦合产生压力尖峰。界面向熔体方向推进几十毫米后熔体域最初100 mm的高度最后可能只剩几个毫米网格单元一边被压缩一边被拉长长宽比飙升。COMSOL变形几何内置了Laplace和Winslow等网格平滑算法但能力有限网格相对厚度低于约3~5层时界面附近的压力离散就会开始病态。我的习惯是在整个模型里给变形几何设置一个初始网格位移为零的初始条件然后在所有外边界用滑移位移而不是固定位移让网格能沿安瓿壁滑移来容纳变形。很多人一出问题就把边界全设成固定结果界面动一下整片网格像拉一块硬板似的全被拽弯倒置单元必然出现。4.2 一次不收敛问题的完整定位过程分享一次真实的排查过程可以当作模板用。那轮参数扫描里基线工况能跑通但把冷端冷却速率从0.5 K/min改成1.5 K/min后t≈3000 s时求解中止报错信息是最后步长小于最小步长瞬态求解失败。我的第一步不是调网格而是关掉流场只保留传热变形几何。结果这个组合顺利跑过了3000 s说明问题不在动网格本身。第二步关掉变形几何保留传热流场也能跑说明问题不在流热耦合。问题就这样被锁定在流场动网格的联合效应上。第三步打开流场和动网格把时间步固定在一个极小值跑几步看中间输出。果然在熔体域厚度压缩到大约还剩8 mm时界面附近压力出现高频振荡紧接着最低层的网格单元开始反转。原因清楚了熔体域变薄、界面又持续释放潜热流场在越来越窄的通道里被挤压力-速度耦合在动网格框架下没有足够阻尼振荡超过阈值后网格直接塌掉。解法分三路第一在熔体域顶部与气相接触的边界设置一个压力约束给封闭空间一个泄压通道第二在界面边界上对网格位移施加额外的切向滑移自由度让网格能顺着壁面滑而不被压扁第三检查层流接口的稳定化设置COMSOL默认的流线扩散没有被自定义设置偷关掉就保留。改完后同一个扫描点顺利跑通后面更极端的参数点只是变慢没有再崩。4.3 稳定收敛的几招储备整理几条我在CZT模拟里沉淀下来的保命经验。设置网格重剖分触发条件。COMSOL变形几何可以在网格质量低于阈值时自动重新划分网格这个功能很诱人但要慎用任何重剖分都会带来历史场变量的插值误差Crank-Nicolson类格式对插值尤其敏感。实在要用就设一个最大单元扭曲度阈值让重剖分尽量少触发。给界面法向速度加限制。把界面推进速度限制在一个物理合理范围比如0~5 μm/s防止初期微小数值振荡被放大成不合理的界面位移。这个限制在COMSOL变量定义里用min和max函数就能实现成本极低但效果明显。用参数扫描接力逼近目标工况。把目标工艺参数拆成多个子工况从简单到复杂逐个计算每一步以上一步的解作为初值。这比直接怼一个极端工况稳定得多本质上是用前一个解当了热启动。定期检查能量守恒。每跑完一段用COMSOL的报告功能统计界面释放潜热总量和通过边界累计导热量两者偏差超过5%就优先怀疑潜热过渡区太宽或网格分辨率不足。能量账不平后面所有结果都不能信。5. 结果后处理与物理可信度检验5.1 界面形态演化与生长速率判据等效热容法的直接输出是温度场和液相分数场。界面位置可以通过派生值来定义在每条轴向采样线上找fL0.5的z坐标。把界面各径向位置的z坐标与中心位置的z坐标做差就是界面凸度。经验值上界面中心略微凸向熔体1~3 mm是比较理想的工艺状态这个形状有利于把杂质推向晶体外周而不是滞留在中心对探测器级CZT尤其重要。生长速率后处理也很直接记录中心点界面z坐标随时间的变化取斜率就是局部生长速度。理想情况下经过一段起步过渡后生长速度应接近安瓿移动速度比如1 mm/h。如果差太多说明热平衡没有建立或者冷却速率设置太快界面在追赶温度场。这种追赶在实验中表现为界面过冷或局部回熔模拟的价值就在这里——不用动炉子就能提前摸清工艺参数区间。把界面推进速率、熔体最大流速和界面凸度三者放在同一张时间曲线上观察几乎一眼就能看出某组参数是否和谐。如果流速随时间剧烈波动、界面推进速度跟着抖那多半是对流进入了非稳态区间需要进一步检查格拉晓夫数。5.2 流场与温度场的交叉验证解出流场后先别急着截图先做无量纲判断。格拉晓夫数Gr g·β·ΔT·H³/ν²以熔体高度H50 mm、温差ΔT30 K、体积膨胀系数2×10⁻⁴、运动粘度5.3×10⁻⁷ m²/s估算Gr大概在10⁶量级这已经是非常确定的自然对流不是可以忽略的次级效应。如果你模拟出来的流速小到几乎看不见第一时间怀疑Boussinesq项的符号或参考密度设置错了。温度场和流场的交叉检查我从三个维度做。一是界面附近的等温线是否被对流拉变形纯导热情况下等温线应该平滑地穿过界面对流则会带来明显扭曲。二是对称轴附近是否出现规则环流环流的方向应该与浮力方向一致热区上升、冷区下降。三是通过界面的总热流是否等于潜热消耗两侧导热的收支差热流账不平基本就是潜热处理的问题。我还喜欢看一个细节熔体顶部或自由液面附近是否存在第二个反向涡。垂直布里奇曼里只出现一个主环流是常见情况但如果工艺参数比较极端可能出现上下层双涡结构这种结构对溶质分布影响很大。做参数扫描时可以专门输出涡心位置随参数的变化这比单看界面形貌更早暴露潜在不稳定。5.3 参数扫描与工艺优化视角等模型跑通COMSOL的参数化扫描就变成了效率工具。用辅助扫描把炉温梯度、安瓿下降速度、熔体初始高度作为参数列表一次性提交多个工况。每个工况都把上一个工况的解作为初值相当于给每个极端参数点搭了一座桥收敛率大幅提升。扫描结果的判读要回到物理。提高炉温梯度会让界面更凸向熔体但伴随而来的是更大的热应力降低下降速度会让界面更接近稳态、流场更平静但生长周期被拉长。用派生值表把每个参数点的界面凸度、最大流速、轴向温度梯度、界面过冷度整理成表基本就能画出工艺窗口的轮廓。有一点必须提醒模拟得到的趋势比具体数值更可靠。比如梯度增加10 K/cm界面凸度增大约0.5 mm这样的趋势对工艺调整有直接指导意义但梯度30 K/cm时界面凸度是1.52 mm这种绝对值只能当参考因为材料参数和边界条件总有不精确的地方。我在对外输出报告时永远把结果写成一个区间而不是一个点这是做仿真的人最基本的专业素养。6. 继续做下去的几个方向与我的实际体会6.1 从二维到三维、从热流到场2D轴对称模型能回答的问题已经很多但它天然过滤掉了周向非对称扰动。真实安瓿在炉膛里不是完美对中炉体也在缓慢旋转这些因素会在晶体里留下螺旋状杂质条纹。想研究这些就必须上3D模型并考虑安瓿旋转无论是旋转坐标系还是滑移网格复杂度都是2D的几倍建议在2D物理彻底摸清后再碰。更值得优先做的是把第三个场加进来——溶质场。CdZnTe熔体中的Zn和Cd在凝固时会发生再分配造成组分偏析。用稀物质传递接口耦合浓度场配合界面处的分凝系数可以定量预测沿晶锭方向的电阻率均匀性。这个结果对探测器材料来说比纯温度分布更接近用户真正关心的指标。6.2 我的体会回头做CZT结晶模拟最有价值的时刻不是模型第一次跑通而是拿模拟得到的温度曲线和实测热电偶数据对比、发现两者相差在2℃以内那次。从那以后我确定了一件事材料参数校准和边界条件确认的优先级高于任何高级物理算法。很多人一上来就追求相场、追求3D、追求辐射换热结果连基线传热都没校出来这是最划不来的。我的最后一条建议是先用最朴素的方法把能量守恒做对再慢慢往模型里加细节。变形几何、固体传热、流场三者的耦合本质上是把晶体生长最核心的物理搬进电脑里。基础骨架稳了后面加什么细节都只是锦上添花。
