t3code：托卡马克等离子体一维径向输运模拟实战指南

做了这么多年聚变等离子体模拟，电脑里装过的输运代码一只手数不过来，t3code是其中资历最老、却最不该被忽略的一个。如果你盯着“t3code”这个名字发愁，多半是在查托卡马克等离子体的一维径向输运演化怎么做，或者看到某篇老文献里引用它作为基准工具，想拿它复现一个密度剖面或温度剖面，却不知道从哪下手。

T3不是那种开箱即用、文档齐全的现代框架，它更像一台做工扎实的老式机床——没有花哨界面，但物理建模思路非常通透，跑一个典型场景只需要几秒到几分钟，特别适合做参数扫描、概念设计、以及快速验证“如果我把加热功率抬高到多少，台基会不会塌”这类问题。这篇文章我按自己的使用经验，把t3code从物理模型、编译运行、算例设置到常见坑点完整梳理了一遍，给准备入门或被迫接手老代码的你省点时间。

1. t3code到底解决什么问题：一维径向输运的取舍

1.1 为什么聚变模拟偏要选一维

托卡马克等离子体真实状态是三维的，有环向、极向和径向三个空间方向，再加上时间，是个四维系统。但你真想研究“密度分布怎么随时间演化”“温度台基怎么形成”，直接上三维湍流代码是杀鸡用牛刀——成本高、周期长，且大量细节对你要回答的问题并没有帮助。

t3code做的是把问题简化成一条径向线。它假设等离子体在磁面上是均匀的，所有物理量只随小半径r变化，再用通量面平均的方式，把三维效应塞进一维系数里。这样做的核心逻辑是：输运过程主要由径向梯度驱动，环向和极向的快速输运会把不均匀抹平，所以描述长时间尺度演化时，“径向一维+时间演化”是一个非常好的近似。

这个取舍的收益非常直接：运行速度快、参数空间容易扫、稳定性和可复现性都远好于大代码。用T3做一个基准算例，从给定初始剖面到接近稳态，通常就是几分钟的事情。它不会给你湍流的细节，但能给你整个等离子体的“骨架演化”，这在做设计参数扫描和趋势判断时足够用了。

1.2 时间尺度分离：欧姆定律是准静态的

T3背后有个关键物理思想：时间尺度分离。等离子体里有极快的磁流体阿尔芬过程，有微秒量级的湍流涨落，也有毫秒到秒量级的能量约束时间，还有更慢的电流扩散和杂质积累。

如果我们把全部过程放在一起算，数值步长必须追到最快过程，计算成本会爆炸。所以t3code只追踪慢过程：等离子体密度、电子温度、离子温度、电流密度剖面、杂质浓度和极向磁场的演化走“输运方程组”；而快过程如磁流体平衡、粒子源和加热沉积，则用准静态平衡和经验模型即时更新。

实际实现中，这意味着每个时间步内，T3会先调用一个平衡计算，基于当前的电流剖面和压强剖面，当场解出磁面几何，更新通量坐标，然后才推进输运方程。这是一个典型的“算子分裂”思路：快物理作为背景，慢物理作为演化主体，两步交替前进。如果你读T3源码，会发现主循环里就是“求平衡-推进输运-更新源项”这个结构在不断重复，逻辑清楚到即使Fortran老代码也不难看懂。

2. 核心物理模型与代码结构速览

2.1 输运方程组：粒子、热、电流一条龙

T3求解的核心方程组可以分成三块，都源自标准通量面平均输运理论。

密度演化方程写为粒子平衡：

∂n/∂t = -(1/V')(∂/∂r)(V' Γ) + S_n

其中V'是体积对径向坐标的导数，Γ是径向粒子通量，S_n是粒子源项，包括中性束注入、壁解吸和再循环等。Γ由扩散系数D和对流速度V convective构成：Γ = -D ∂n/∂r + V_conv n。

电子和离子的能量方程则是热平衡：

(3/2)∂(nT_e)/∂t = -(1/V')(∂/∂r)(V' q_e) + P_heat_e - P_eq - P_rad (3/2)∂(nT_i)/∂t = -(1/V')(∂/∂r)(V' q_i) + P_heat_i + P_eq

q是径向热流，P_heat是外加加热功率密度，P_eq是电子离子碰撞能量交换项，P_rad是辐射损失。真正精细的地方在于P_rad通常包含轫致辐射、线辐射和杂质辐射模型，其中杂质含量一旦涨起来，辐射会把边缘温度压得很低。

电流密度演化则是一个磁扩散方程：

∂B_θ/∂t = (η/μ₀)(1/r)(∂/∂r)(r ∂B_θ/∂r) + 驱动项

这里η是电阻率，由Spitzer公式给出，驱动项包括欧姆电流、束驱动电流和自举电流。自举电流的模型直接决定了剖面是否“空心”，也影响最终是否进入稳态。

这三套方程用径向网格离散，然后时间推进。T3默认推进方式是显式的，优点是简单、不需要大型矩阵求解；缺点是时间步长会被稳定性条件卡住。所以跑长时间演化时，注意时间步长设到满足库朗条件，具体经验我在后文会讲。

2.2 径向网格与边界条件：近轴点最容易出问题

T3的径向网格一般是等间距或略微加密的近轴分布。网格数不用太多，一二百个点已经能描述剖面形状，多了反而会拖慢速度。反正它是一维代码，单步计算成本本来就低。

边界条件分两个端点。中心r=0处，所有通量归零，所有物理量的径向梯度为零——这是对称性条件，不能乱改。边界r=a处，也就是等离子体边缘，温度密度要被设定为有限的边界值。在实践中，T3的边界条件通常由刮削层模型或直接给定边界温度和密度表格来提供。

踩坑提示：多数运行发散问题都出在边界条件不自洽。你设的边界温度低于辐射和输运决定的自然边界值，边缘立刻出现陡梯度，导致热流数值振荡。如果你看到T_e剖面在最后一个网格点上跳来跳去，十有八九是边界值不匹配，先调边界比调时间步长更管用。

2.3 输运系数：新经典基础加反常提升

输运系数是整个模拟的“灵魂”，T3里通常拆成两部分写。

新经典输运部分基于Hinton-Hazeltine模型，处理碰撞导致的不连续轨道效应，包括香蕉区、平台区和箍缩区的输运。这一项在高温低碰撞率区域较弱，但它是自洽的，由等离子体参数直接推导出来，不需要额外调节。

反常输运部分则是“调参的重头戏”。T3默认提供几种参数化模型，比如定常D、定常χ，以及更复杂的缩放模型。实际操作中绝大多数人都会在D和χ上乘一个系数，去匹配实验约束时间。我的做法一般是先跑一组基准：给定密度和加热，调χ_e和χ_i让温度剖面匹配实验点，然后再用这个标定过的输运系数去扫描其他参数。

需要警惕的是参数过拟合。T3这种参数化输运代码，总是能通过调系数拟合出一两组实验剖面，但到了预测新场景时，系数不一定可迁移。所以做预测前，至少要用两三个不同放电的基准点做交叉验证，不要只靠一个放电定系数。

3. 从拿到t3code到跑出第一张剖面图的完整流程

3.1 获取源码与编译环境准备

t3code是老牌Fortran代码，很多研究所内部都有本地拷贝，如果你能访问项目仓库或代码分发渠道，直接拿到源码包最省事。拿到后先看压缩包里的文件结构：一般包含一个主程序文件、若干物理模块子程序和一个示例输入文件，还可能附带一个用于绘图的后处理脚本。

编译前确认你装了Fortran编译器。这里我推荐gfortran，免费、跨平台、对老代码兼容性好。Intel的ifort也行，有些老代码针对ifort做过向量化优化，但如果没有licence就没什么好纠结的。在Linux环境下，一条命令就能完成编译：

gfortran -O2 -o t3code t3main.f90 t3transport.f90 t3balance.f90 ...

如果你的源码包里有Makefile，就直接：

make

编译报错一般是两类：一类是Fortran语法版本差异，老代码用的固定格式（.f后缀）和自由格式（.f90后缀）混在一起，gfortran默认可能按后缀区分格式，注意把固定格式文件加上 -ffixed-form；另一类是缺少某些工具库或数学函数库引用，用 -lm 或添加对应链接路径即可解决。

一个务实的建议：先用小网格、最少物理模块编译通过，再逐步打开完整物理选项。一次全开，报起错来就要排查很久。

3.2 输入文件结构：namelist风格要耐心

t3code的输入控制通常走Namellist（Fortran的命名列表）方式，你会在示例输入文件里看到类似下面的结构：

&name_list nr = 128 rmax = 0.85 dt = 0.01 ntstep = 10000 ne_bc = 1.0e19 te_bc = 100.0 ti_bc = 100.0 ipedestal = 1 chi_e = 1.2 chi_i = 0.8 d_coef = 0.5 /

每个参数都有明确含义。nr是径向网格数，rmax是等离子体半径（米），dt是时间步长（秒），ntstep是总时间步数，ne_bc、te_bc、ti_bc是边界电子密度、电子温度、离子温度，ipedestal开关决定是否启动台基高度模型，chi_e、chi_i、d_coef分别控制电子热扩散系数、离子热扩散系数和粒子扩散系数的基值。

第一次跑，请务必先保留所有参数默认值，只改输出文件名。T3内置的基准算例大概率能直接跑通，给你一个已知正确的参照。先建立“能跑”的基线，再做物理修改，这个顺序能帮你避免“改了十个参数出了问题却不知道是哪一个改坏的”的尴尬局面。

3.3 正式运行与标准输出解读

编译和输入文件准备完毕后，运行非常简单：

./t3code < t3code_input.dat

终端会实时打印当前时间步、等离子体电流、约束时间、中心温度等关键指标。如果你看到这些量稳定增长并趋于饱和，恭喜，这个算例走上了正轨。T3的每个时间步都会把剖面数据写入输出文件，一般是纯文本格式，多列分别为径向位置、电子密度、电子温度、离子温度、电流密度等，方便你用任何工具处理。

我一般会把终端输出重定向到log文件，方便逐段回看：

./t3code < t3code_input.dat > run.log 2>&1

如果中途跑飞（数值发散），run.log里通常能看到时间步推进失败或剖面出现NaN的痕迹，这是第一手排查资料，丢掉很可惜。

3.4 把输出文件画成剖面演化图

T3的历史输出文件记录的是“每一小间隔时间步的全剖面快照”。文件如果很大，画图前要先抽样，否则gnuplot或matplotlib都能被卡死。

用Python清洗和绘图可以参考如下脚本思路：

import numpy as np import matplotlib.pyplot as plt data = np.loadtxt("t3code_profiles.dat") # 假设第1列是时间轴，第2列是r，第3列是Te tvals = np.unique(data[:, 0]) fig, ax = plt.subplots(figsize=(8, 5)) for t in tvals[:: max(1, len(tvals)//10)]: # 抽样10个时刻 mask = data[:, 0] == t ax.plot(data[mask, 1], data[mask, 2], label=f"t={t:.2f}s") ax.set_xlabel("r (m)") ax.set_ylabel("Te (eV)") ax.legend() plt.tight_layout() plt.savefig("te_evolution.png", dpi=150)

如果你拿到的输出格式是“每个时间步一个文件”的形式，那就更简单了，循环读取文件绘制即可。画图时注意横轴是“小半径”还是“归一化半径”，T3默认常见输出是物理半径，画图前记得看单位，避免出现一根看起来像是被压缩过的奇怪曲线。

3.5 参数扫描实例：加热功率对温度剖面和约束时间的影响

用t3code做参数扫描，是它最擅长的事情。现在演示一个典型任务：扫描总加热功率从5MW到25MW，看中心电子温度和能量约束时间怎么变化。

做法很简单：写一个小循环，每次改输入文件里的加热功率参数，调t3code，抓取输出里的中心温度和约束时间。我习惯用Python脚本控制这个过程：

import subprocess import numpy as np powers = np.linspace(5.0, 25.0, 9) results = [] for p in powers: with open("t3code_input.dat", "r") as f: text = f.read() text = text.replace("p_heat = 10.0", f"p_heat = {p:.2f}") with open("t3code_input.dat", "w") as f: f.write(text) subprocess.run(["./t3code < t3code_input.dat > scan_run.log 2>&1"], shell=True) # 从输出文件读取中心温度和约束时间 data = np.loadtxt("t3code_profiles_last.dat") results.append((p, data[-1, 2], data[-1, 0]))

跑完之后把结果画出来，你会看到中心温度随功率上升、能量约束时间却往下掉——这正是典型L模约束定标率的体现。如果你给代码开了台基模型，可能还会看到功率超过某个阈值后台基高度跳变、约束时间突然变化的“相变”行为。这种参数空间里的拐点，用大代码扫描是奢侈的，用t3code可以轻松探索。

4. 实战中经常翻车的几个点：问题定位与修正

4.1 初始剖面的选择：冷启动还是热启动

T3的初始剖面直接决定前几百步的稳定性。最常见的翻车是初始剖面与世界平衡差别太大，第一大步的梯度就超出了显式时间步的稳定性极限，直接振荡。

解决方式有两条路。第一，使用更平滑的初始剖面，比如从抛物线形状开头，不要上来就塞一个很尖锐的台基状剖面。第二，减小初始时间步长，让代码自己慢慢适应，通常跑过前几百步后就可以适当增大步长。T3允许你在输入文件里做时间步长的分段控制，可以设定前500步用dt=0.001，之后切换到dt=0.01，这样的“预热”能有效避免冷启动发散。

如果你的研究对象是台基演化，那就必须从高约束模开始做热启动。给一个事先跑好的稳态剖面作为初值更稳妥。T3支持读入上次运行的最终剖面作为新运行的初始条件，这个功能强烈建议用起来，别每回都从同一个粗糙初值算起。

4.2 显式时间步长、网格数与源项增量的三角约束

T3采用显式推进，时间步长受限于网格分辨率与最大输运系数的乘积。公式上大概满足：

dt ≤ Δr² / (2 max(D, χ))

如果你把网格数从128加到256，时间步长理论上要缩小4倍才能维持稳定。实操中我遇到过有人只加密网格、忘改时间步长，然后五分钟前还在正常跑的算例突然就爆了。

排查这种问题先看run.log：如果发散点固定出现在同一个时间步附近，并且剖面在某个径向位置开始锯齿状振荡，基本可以判定是稳定性条件被破坏。对策是缩小dt，或者用输运系数更保守的参数值。源项增量也得注意，加热功率密度在网格边界处跳变过大，会让对应网格点的温度更新过猛，局部产生虚假热流。更平滑的沉积剖面能缓解这个问题。

4.3 台基高度与边界温度：两个容易混淆的调参旋钮

很多T3的中文用户上来就调“台基参数”，但台基模型效果不明显，就开始怀疑代码。这里有个容易混淆的点：T3这类代码里的“台基”模型，一般说的是控制台基高度或宽度的经验关系式，它在漂移动理方程里不是自觉出现的，而是被“钉”在某个值附近。

真正的边界温度则直接影响台基底部的位置。边界温度给得太低，台基难形成；给得太高，台基会异常高。我的经验是把边界温度作为一个“锚点”，配合台基标度关系一起调。先固定边界温度不变量，扫描台基模型参数，找到匹配实验高度的数值；再回头微调边界温度，让台基底面位置也match上。两步分开做，不要一起动，否则你永远不知道是谁在起主导作用。

4.4 杂质含量与辐射损失：隐藏的剖面“刺客”

如果你跑参数扫描时发现温度剖面整体偏低，特别是边缘温度上不去，而密度剖面看起来正常，先别急着调热输运系数——看一眼杂质浓度参数。T3支持输入有效电荷数Z_eff剖面或固定杂质浓度，辐射功率密度跟Z_eff近似成正比，Z_eff从1涨到3，辐射损失能吃掉一大块加热功率，边缘温度应声而落。

这是个非常容易“误诊”的问题。我见过有人为了抬高边缘温度，费劲半天把chi_e调了五倍，结果发现只是辐射模型开关忘了关。建议每次换工况时，都先确认杂质相关的参数和辐射模型开关是否符合你的物理预期，再谈其他调整。

5. 快速排查速查表

查问题时我习惯只改一个变量，跑短时间步对比，再改下一个。T3跑得快，试错成本低，不要一上来就大改参数组合，那是把简单问题复杂化。

6. 从t3code往后走：扩展思路与个人体会

t3code这类代码的价值在框架，而不在某个参数的精确值。跑通它、看懂它，等于把一个完整的一维输运思路装进脑子里。后续如果要做更复杂的分析，完全可以以T3的输运方程和推进框架为蓝本，加入自己的模型，比如新的反常输运系数表达式、更精细的边界模型、或者耦合一个简单的刮削层模型。我在实际项目里就用T3做过沙箱测试，把一种新的台基高度标度式写进去，对比基准算例，确认趋势合理后，才决定要不要移植到大代码里跑高成本模拟。

还有个很容易被忽略的小技巧：T3速度快，非常适合作不确定性量化。你可以在输入参数上加随机扰动，跑几十次ensemble，看输出剖面的散布范围。这种事在TRANSP这种级别的大代码上几乎不敢想，但在t3code上只是多花几分钟的事。

如果你是被课题或项目“强制”需要弄懂这个代码的，建议你从基准算例开始，逐行过一遍主循环，搞清楚“平衡—输运—源项”三个模块之间的调用顺序。这个顺序搞懂了，T3就已经是你手里的工具，而不是一团黑箱。后续想往里面加物理或者换数值格式，都有明确的下手点。