一维辐射传递方程的谱方法求解

Ａｂｓｔｒａｃｔ ： Ｕｓｅｉｎｇ ｔｈｅ ｄｉｆｆｅｒｅｎｃｅ ｔｏ ｇｏ ｏｎ ｎｕｍｅｒｉｃａｌ ｃａｌｃｕｌａｔｉｏｎ ｆｏｒ ｏｎｅ－ｄｉｍｅｎｓｉｏｎａｌ ｈｅａｔ ｃｏｎｄｕｃｔｉｏｎ ｅｑｕａｔｉｏｎ ，ａｎｄ ｇｉｖｅ ａｎ ｉｌｌｕｓｔｒａｔｉｏｎ ｏｆ ｔｈｅ ｓｏｌｕｔｉｏｎ ． Ｋｅｙ ｗｏｒｄｓ ： ｄｉｆｆｅｒｅｎｃｅ ｍｅｔｈｏｄ ；ｏｎｅ－ｄｉｍｅｎｓｉｏｎａｌ ｈｅａｔ ｃｏｎｄｕｃｔｉｏｎ ｅｑｕａｔｉｏｎ ；ｂｏｕｎｄａｒｙ ｃｏｎｄｉｔｉｏｎ
ｉ ，ｊ
Δｔ
（５）
以（４） 、（５） 代入（１） 式得
ｕ － ｕ ｉ ，ｊ＋ １
ｉ ，ｊ
Δｔ
＝
ａ２
ｕｉ －１ ，ｊ
－ ２ ｕｉ ，ｊ Δｘ２
＋
ｕｉ＋ １ ，ｊ
＋
ｆ （ ｉ ，ｊ）
（６）
解得
ｕｉ ，ｊ＋ １ ＝ ｃ （ ｕｉ －１ ，ｊ ＋ ｕｉ＋ １ ，ｊ ） ＋ （１ － ２ ｃ） ｕｉ ，ｊ ＋ Δｔｆ （ ｉ ，ｊ ）
式处理即可 ．
２） 如果是波动方程 ，那么初始条件中含有“初始速度” ，即 ｕｔ （ｘ ，ｔ） ｜ ｔ ＝ ０ ＝ 宝（ｘ） ，利用向后差分法得
抄 ｕ（ ｉ ，ｌ） 抄ｔ
＝
ｕｉ ，２ － ｕｉ ，０ ２ Δｔ
＝
宝（ ｘ ）
（３５）
则
ｕｉ ，０ ＝ ｕｉ ，２ － ２ Δｔ宝（ｘ ）
（３６）
利用（３６） 式便得以求解 ．
用中心差分近似代替对空间的偏微分即2u抄用向前差分近似代替对时间的偏微分即uij1uij2uxxfxt0xl0t123抄2xui1j2uijui1jx24抄u抄tt5以45代入1式得uij1uijta2ui1j2uijui1jx2fij6解得uij1cui1jui1j12cuijtfij7其中ctax228根据式7如果已知j不同i坐标每一个格点的温度值并且由11类边界条件可知两边界i1及in上的温度值那么就可以求出j1坐标上每一个格点上的温度值

1 DRESOR 法求解含有边界的各 向异性散射问题
方 向 辐射 信 号时 体 现 了 该方 法 的 优 越性。 DRESO R法是一种求解辐射传递方程的新方法 , 它是基于蒙特卡洛方法提出的 ,继承了蒙特卡洛 方法 适 应 性强 、 能处 理 复 杂 几何 形 状 的辐 射 传 递 问题的特点 ,同时 ,它还能够给出高方向分辨率的 辐射强度 [5 ]。
表 1 无量纲反射热流
T ab. 1 Dimensio nless reflective hea t flux
辐射元法
蒙特卡洛法 有限体积法
0. 355 53 0. 514 59
0. 719 48
0. 357 64 0. 515 33 0. 610 27 0. 674 42 0. 719 22 0. 754 02
Iη( r ,s ) =
1 π
{πX( rw ) I bη( rw
)+
∫s Rd
(
′
rw
,
r
w
,s )
[πX( r′w
)
Ibη(
′
rw
)
]dA′+
W
∫Rsd (r″, rw , s) V
[4πUZ( 1 -
kη) I bη(r″)
]dV″} ×
∫ ∫ exp -
s
_
ηds″
+
0
s 0
1 4π{
4πUη(
具 体做法是: 将通过式 ( 5) 计算得到的辐射 强度 IR ( r , s) 代入到式 ( 3) 和 ( 4) 的右边计算得到 IRw ( rw , s) 和 SR ( r , s) ,再将 IRw ( rw , s) 和 SR ( r , s) 代 入到方程 ( 2) 的右边计算得到验证辐射强度 IV ( r , s)。这样将计算辐射强度 IR ( r, s) 和验证辐 射强度 IV ( r, s) 作比较检验辐射传递方程求解结 果的准确性。

一维热传导方程的差分法一维热传导方程是指在一维空间中，描述材料内部温度分布随时间的变化过程的方程式。

可以表示为：$$\frac{\partial u(x, t)}{\partial t} = \alpha \frac{\partial^2 u(x,t)}{\partial x^2}$$其中，$u(x, t)$ 是空间坐标为 $x$，时间为 $t$ 时的温度，$\alpha$ 是热扩散系数。

在边界条件确定的情况下，可以得到一维热传导方程的解。

然而，在实际应用中，解析解并不总是容易或可行的，因此需要使用数值方法进行近似求解。

常用的数值方法有有限差分法、有限元法和谱方法等，其中有限差分法是最为简单、易于实现的方法之一。

有限差分法的基本思想是将连续的空间区间离散化为若干个节点，将时间轴离散化为若干个时间步。

在空间和时间轴两个方向上，分别对热传导方程进行差分，得到离散的差分方程组，从而可以求得数值解。

在一维热传导方程的差分过程中，我们首先需要将空间区间 $(0, L)$ 划分为 $N$ 个等间距的节点，每个节点间距为 $h = \frac{L}{N}$。

我们使用 $u_i^n$ 表示节点$i$ 在时间步 $n$ 时的温度，其中 $i = 0, 1, ..., N$，$n = 0, 1, ..., M$。

接下来，我们对一维热传导方程进行中心差分，得到：$$\frac{u_i^{n+1} - u_i^n}{\Delta t} = \alpha \frac{u_{i+1}^n - 2u_i^n + u_{i-1}^n}{h^2}$$其中，$\Delta t$ 是时间步长。

可以将上式改写为：其中，$r = \frac{\alpha \Delta t}{h^2}$。

由于在一维空间中，有两个边界，因此需要对边界进行特殊处理。

常见的边界处理方式有三种：1. 固定边界：将边界上的温度固定为某一值；2. 自然边界：假设边界处热通量为零，从而根据傅里叶定律可以求得边界处温度的梯度值，从而推算出边界的温度值；3. 第二类边界：将边界节点的温度根据边界条件与内部的节点做差分，从而计算出边界节点的温度值。

摘 要 ：针 对 三 维 长方 形 炉 内具 有 吸 收 一 射 介 质 的 辐 射 换 热 ，基 于 Ｃ ｅｙ ｈｖ配 置 点 谱 方 法 和 Ｓ ｈ ｒ 解 开 发 了 发 ｈ ｂ ｓｅ ｅｕ 分 直 接求 解 辐 射 离 散 坐标 方 程 的求 解 器 。针 对 离 散 后 所 得 到 的 三 维 矩 阵 方 程 ，分 别 用 两 种 方 法 进 行 求 解 ，一 种 是 用 张量 积 将 三 维 转变 成二 维 然 后 直 接 用 Ｓｈ ｒ 解 求 解 ；另 一 种 是 自行 开 发 三 维 Ｓ ｈ ｒ 解 直接 求 解 。数 值 实 验 ｅｕ 分 ｅｕ 分 表 明 ，在 相 同 的输 入 参 数 下 ，新 求 解 器具 有很 好 的 精 度 ，尤 其 相 比 于 标 准 离 散 坐 标 法 ，新 求 解 器 能 节 省 大 量 计
最近作者将chebyshev配置点谱方法成功应用于求解一维具有各向异性散射的辐射问题25一维梯度折射率辐射问题26半透明梯度折射率内辐射与导热耦合问题27一维瞬态各向异性散射梯度折射率介质内辐射与导热耦合问题28以及同心球内具有参与介质时辐射与导热耦合问题24作为chebyshev配置点谱方法计算辐射换2428的多维扩展本文开发了直接求解三维长方体炉内具有吸收发射介质情况下的离散坐标辐射方程组
算 时间 。特 别是 基 于三 维 Ｓ ｈ ｒ 解 的 直 接 求 解 器 ，在 相 同 的 输 人 参 数 下 ， 计 算 时 间 只 有 标 准 离 散 坐 标 法 的 ｃｕ 分
１Ｏ ～ １ 。
关 键 词 ：辐 射换 热 ；离 散 坐 标 法 ；Ｃ ｅｙ ｈｖ配置 点 谱 方 法 ；Ｓｈ ｒ 解 ｈ ｂ ｓｅ ｃｕ 分
Ａｂ ｔ ａ ｔ Ｂａ ｅ ｎ ｔ ｓｒ ｃ ： ｓ ｄ ｏ ｈｅＣｈｅ ｓ ｅ ｏ ｌ ｃ ｔｏ ｐｅ ｔ ａ ｅ ｈ ｎｄ Ｓｃ ｂｙ ｈ ｖ ｃ ｌｏ ａ ｉ ｎ ｓ ｃｒ ｌｍ ｔ ｏｄ ａ ｈｕｒｄ ｃ ｍｐｏ ｉｉ ｎ，ｔ ｉ ｅ ｔｓ ｌ ｅ ｓ ｅｏ ｓ ｔｏ ＷＯ ｄ ｒ ｃ ｏ ｖ ｒ ｆ ｒ ａ ａ ｉ ｅ ｉｃ ｅ ｅ ｏ ｒ ｄｉｔｖ ｄ ｓ ｒ ｔ ｏｒ ｉ ａ ｅ ｅ ａ ｉ ｓ ｏ ａ ｈｒ ｅ ｄ ｍｅ ｉ ａ ｒ ｃ ａ ｇｕ ａ ｆ ｒ ｃ ｗｉｈ ｂｓ ｒ ｎ — ｄ ｎ ｔ ｓ ｑｕ ｔｏｎ ｆ ｒ ｔ ｅ — ｉ ｎｓｏｎ ｌ ｅ ｔ ｎ ｌ ｒ ｕ ｎａ ｅ ｔ ａ ｏ ｂｉ ｇ ｅ ｉｔｎｇ ｍ ｔ ｉ ｍｅ ｉ ｍ ａ ｅ ｅ ｅ ｏ ｅ ． ｒ ｈｅ ｈｒ ｅ ｄｉ ｅ ｉ ｎａ ｍ ａ ｒｘ ｑｕ ｔｏｎ ｒ ｍ ｄ ｓ ｒ ｔｚ ｔｏ ｏ ｔ ｅ ｄｕ ｒ ｄ ｖ ｌ ｐ ｄ Ｆｏ ｔ ｔ ｅ — ｍ ｎｓｏ ｌ ｔ ｉ ｅ ａ ｉ ｆ ｏ ｉｃ ｅ ｉａ ｉｎ ｆ ｈ ｒ ｄ ａ ｉ ｅ ｔ ａ ｆ ｒ ｅ ｕａ ｉ ｎ，ｏ ｅ ｓ ｖ ｒｉ ａ ｅ ｎ ｔ — ｍｅ ｓｏ ｌＳ ｈｕ ｅ ｏ ｐ ｉｉ ｎ ａ ｔ ｒｔ ｒ ｎｓｏｒ ａ ｉｔｖ ｒ ｎｓ ｅ ｑ ｔｏ ｎ ｏｌ ｅ ｓ ｂ ｓ ｄ ｏ ｗｏ ｄｉ ｎ ｉ ｎａ ｃ ｒｄ ｃ ｍ ｏｓｔｏ ｆｅ ｈｅｔ ａ ｆ ｍ ｏ ｈｒ ｅｄｉ ｅ ｉ ｎａ ａ ｒｘ ｅ ｕ ｔｏｎ ｔ ｗｏ ｄｉ ｎ ｉ ｎａ ｎｅ ｕｓｎｇ ｔ ｎｓ ｒ ｐ ｏ ｕｃ ，ａ ｈ ｈｅ ｓ ｂ ｓ ｄ ｆｔ ｅ — ｍ ｎｓｏ ｌｍ ｔ ｉ ｑ ａ ｉ ｏ ｔ ～ ｍｅ ｓｏ ｌｏ ｉ ｅ ｏ ｒ ｄ ｔ ｎｄ ｔ ｅ ｏｔ ｒｉ ａ ｅ

十三、程序流程说明(FDTD程序)
提取波形注意事项： 1)第一次运行该程序，为了提取入射波，需要将计算空间全部设置 为自由空间，并得到数据文件，如inc.dat 2)第二次运行程序，此时提取的波形为入射波和反射波叠加，并得 数据文件，inc_ref.dat 3)两次数据相减，得到反射波波形，并生成数据文件ref.dat 4)透射波的提取是第二次运行程序时提取。
一维简易时域有限差分法(FDTD) 计算介质板反射和透射系数
汤炜
一、一维Maxwell方程
假设电场极化方向为x，磁场极化方向为y，波传播方向为z。所有 物理量在x、y平面不变，而仅仅为z的函数，一维Maxwell方程为
Ex H y
z
t
H y Ex
z
t
二、时间，空间离散
•核心思想：将计算区域的空间和时间进行划分。 空间：例如：三维空间划分为立方块，二维空间划分为正方柱，一 维空间划分为平面板。划分的区域非常小，以至于可以认为场量在 该区域是不变的。 时间：将电磁波的与目标的作用时间划分为很多时间小段，可以认 为场量在该时间段内是不变的。
0.5
t
s
Exn
k
1
Exn
k
五、差分与Maxwell方程的结合
H
n.5 y
k
0.5
H
n.5 y
k
0.5
t
s
Exn
k
1
Exn
k
同理根据另一Maxwell方程得：
E n1 x
k
Exn
k
t
s
H
n.5 y
k
0.5
H
n.5 y
k
0.5
两个迭代方程中 右边：后时刻场量 左边：前时刻场量 即如果能够得到前一时刻的电场和磁场，根据方程即可得到后一时 刻的场量。

z →−m
lim (z + m)Γ(z ) lim
(z + m)(z + m − 1) · · · (z + 1)Γ(z + 1) z →−m (z + m − 1) · · · (z + 1)z Γ(z + m + 1) = lim z →−m (z + m − 1) · · · (z + 1)z 1 = (−1)m m!
ξ →+ ∞
因而在这种情况下, 原 Schrödinger 方程无有界解. 现讨论 v 是非负整数的情况, 即 v = m, m 为非负整数. 在式 (23), (24)的基础上, 若 m 为偶数, 则式 y1 变为未定式, 而 y2 仍为原无界函数. 这时由递推式 (19), am+2 = 0, y1 为 m 次多项式；同样, 若 m 为奇数, 则 y1 为无界函数, y2 为 m 次多项式. 现研究 Γ 函数在 0 和负整数附近的性质. 由于 Res[Γ(z ), −m] = =
ξ2
(10) (11)
y ′′ − 2ξy ′ + (ν − 1)y = 0. 式 (11)即为 Hermite 微分方程.
2
2 Hermite 方程的幂级数求解
下面用幂级数法对 Hermite 方程进行求解. 重写 Hermite 方程如下, 在其中令 λ = ν − 1： y ′′ − 2ξy ′ + λy = 0. 设 y (ξ ) 为如下幂级数形式： y (ξ ) = 那么有 y
l
(39)
变换求和哑元, 令 k = l − i, 则 y1 =
l ∑ (−1)l−i 22l−2i l! i=0
(2l − 2i)!i!

其中 τ = u l ( z ) dz ，即 dτ ( z ) = ul ( z )dz
∂
∫
z
如果单片叶子的单次散射反照率是一个常数，那么辐射传输方程可变换为另一种形式。
Q
1
π
1
Γ ( Ω' → Ω ) =
1 2π
2π
∫ g l (Ω l ) | Ω l ⋅ Ω' | f (Ω' → Ω, Ω l )dΩ l
−µ
dL( Z , Ω) + σ e ( Z , Ω) L( Z , Ω) = ∫ σ s ( Z , Ω ' → Ω)L( Z , Ω ' )dΩ ' dz φπ
此处 L 代表光亮度，其中
σ e 称为消光系数，它代表光路介质对光子的吸收与散射致使
57
光亮度在传播方向上减弱，
σ s 称为散射削弱系数（包含了相位函数） ，它描述了经多次散射
f s = K ( k , µ ' ) F ( n, µ ' )δ ( µ − µ ' )
其中 K ( k , µ ) = exp −
'
2 kt gθ ' π
K 为描述叶子表面粗糙程度而引入的修正系数（0<K<1） ，其中 k 称为叶毛系数，取值 范围为 0.1～0.3。
1 sin 2 (θ '−θ s ) t g (θ '−θ s ) F ( n, µ ' ) = 2 + 2 sin (θ '+θ s ) t g (θ '+θ s )
− +
↓
−
−
F + 与F − ，这样微分——积分辐射传输方程便可简化为一组线性微分方程。

C p
DT Dt
div(gradT
qr ) qo
T
DP Dt
左端为瞬态能量的储存，称非稳态项
右端第一项为导热、辐射项
右端第二项为源项
右端第三项为粘性耗散项
2021/5/4 右端第四项为压缩膨胀式压力作的功10
辐射传递方程的积分形式 RTE积分法
di
dk
i
(k
,
)
I
(k
,
)
用expkλ乘全式，并从0到kλ(S)积分得
a
dv
, T, Peb , T i d
2021/5/4
18
辐射平衡
qr 0 或
0
a dv
, T, Peb , Td
0 a
dv
, T, Pi d
2021/5/4
19
普朗克平均吸收系数
aP T, P
0
a
,
T,
Peb
,
T
d
0
eb
,
T
d
0
a
,
T,
Peb
,
T
d
T 4
aP是一个用黑体发射光谱度量的光谱吸 收系数的平均值。在考虑容积发射和一
些有限的辐射传递情况下，都证明它是
很有用的。 2021/5/4
20
包括散射的辐射流散度
, i
1 4
4 , , i d
i
1 4
i i 4 , i di
qr
40
a eb a s i
s
4
i
4
i
,
i
,
i
di
d
对各向同性散射，,,i 1，因而， ,i 1

1

1
( u ) du
then 2 n 1 {C 4
1

1
u
B
dL dz LP
n
P n ( u ) du } ( u ) du
n
1

1
so
: [cos(
)] (cos
2 n 1 {C 4
1
g )
P n ( u ) du } * Pn (u )
jk
k ] j d
Li p i jk ( p ) k j d
i 0 j 0 k 0



Li p [ jk ( p ) j i k d ]
j 0 i 0

方程中c与方向无关，只与p有关。
将L与
的展开式代入RTE 中得： right c( p) L j p j [ Li p ji p ] j
j 0 j 0 i 0
{c( p) L j p Li p ji p } j
Bn

m
n m
n 0
!
cosdL / dz cL Bn L( z, )
n 0
The RTE can be written:
(n m)! m Pn (u ) Pnm (u) exp jm( )d m n ( n m)!
所以：
cos
4
2 4 nk nk 2k 1 2k 1
dL 4 Pk (u )d C LPk (u )d LBk Pk (u)d dz 2k 1 4 4

1 2π
2π
∫ g l (Ω l ) Ω l ⋅ Ω' f (Ω' → Ω, Ω l )dΩ l
−
−
如果再假定 g l (Ω l ) = 1 （取球面型）
58
则 Γ ( Ω' → Ω ) =
t ω [sin β − β cos β ] + l cos β 3π π
其中
β = cos −1 (Ω, Ω' ) ω = rl + t l
其中 θ s = sin
−1
sin θ ' n
尔镜面反射公式
n 为叶子的光学折射系数，F 为菲
∫ f (Ω φπ
'
→ Ω , Ω l )dΩ = rl+ + rl− + t l+ + t l− + K ( k , µ ' ) F ( n, µ ' )
2.2.4．连续植被的辐射传输方程 一般水平均匀，垂直分层介质中的辐射传输方程可表达为
其中 τ = u l ( z ) dz ，即 dτ ( z ) = ul ( z )dz
∂
∫
z
如果单片叶子的单次散射反照率是一个常数，那么辐射传输方程可变换为另一种形式。
Q
1
π
1
Γ ( Ω' → Ω ) =
1 2π
2π
∫ g l (Ω l ) | Ω l ⋅ Ω' | f (Ω' → Ω, Ω l )dΩ l
与一般辐射传输方程等式右边项相比，则
σ s ( z , Ω' → Ω ) =
− − ul ( z ) g l ( z, Ω l ) | Ω l ⋅ Ω | f (Ω' → Ω, Ω l )dΩ l ∫ 2π 2π

辐射转移方程的数值解法及其应用辐射转移是指光线、热量等物理量在介质之间传播的过程。

辐射转移方程是对辐射转移过程的数学描述，其解法对于许多领域的研究有着重要的意义。

本文将重点探讨辐射转移方程的数值解法及其应用。

一、辐射转移方程的基本形式辐射转移方程描述了辐射能量在介质中传播的规律。

以一维情况为例，辐射转移方程的基本形式可以表示为：dI(x)/dx + σ(x)I(x) = σ(x)E(x)其中，x为介质的位置，I(x)为位置x处的辐射强度，σ(x)为介质的吸收系数，E(x)为位置x处的辐射源函数。

二、数值解法辐射转移方程是一种偏微分方程，通常无法通过解析方法求得精确解。

因此，人们采用数值解法来近似求解。

1.离散化首先，将介质空间划分为若干离散的位置点，即网格。

然后，将辐射转移方程的微分形式转化为差分形式。

通过将导数项近似为差分项，可以得到离散的辐射转移方程。

2.迭代求解离散化后的辐射转移方程通常形式为：I(x) - I(x+Δx) = Δx [σ(x)I(x) - σ(x)E(x)]需要注意的是，I(x)和E(x)都是未知的。

为了求解这个方程，我们可以采用迭代的方法。

首先，给出一个初始解I0(x)，然后根据上述方程更新I1(x)，依此类推，直到解收敛为止。

三、应用领域辐射转移方程的数值解法在许多领域都有重要的应用。

以下介绍其中几个典型的应用领域。

1.大气科学辐射转移方程的数值解法在大气科学中有着广泛的应用。

通过求解辐射转移方程，可以模拟大气中太阳辐射、地球辐射等的传播过程，从而了解大气中的能量平衡、热量分布等重要参数。

2.医学影像学在医学影像学中，辐射转移方程的数值解法被用于模拟射线在人体组织中的传播过程。

通过对辐射转移方程的求解，可以得到图像中各个位置的射线强度分布，从而实现医学图像的重建和分析。

3.能源工程辐射转移方程的数值解法在能源工程中也有重要的应用。

以太阳能为例，通过求解辐射转移方程，可以模拟太阳辐射在太阳能电池板中的传播过程，从而分析电池板的效率、热量分布等关键参数，为太阳能的利用提供指导和优化。

第一章 基本概念 第二节 辐射传输 (radiance transfer) √ §1.2.1 传输方程 §1.2.2 源函数中散射的表达 §1.2.3 辐射传输方程的解
1/14
Maxwell方程组与辐射传输方程
麦克斯韦方程组描述了电磁场的基本规律。一般而 言，波长较长的电磁波波动性较为突出。在微波遥 感领域，更常看到用麦克斯韦方程组解释电磁波与 介质的相互作用。
5834离散纵标方法利用离散纵标方法可以将辐射传输方程中的散射相函数用勒让德多项式展开即用求和式代替方程中的积分式进而将原有的积分微分方程转化为微分方程组最终通过边界条件的代入求解辐射在几个特定方向由高斯点决定上的解析解
遥感物理
第五章 辐射传输方程
邓孺孺副教授
中山大学地理学院 遥感与地理信息工程系
遥感物理
II0e0/
请注意指数形式在辐射传输中的作用。
14/14
总结
两个概念：光学厚度、平面平行介质 一组不同表达形式的传输方程：
dI I J kds dI IJ
d dI IJ
d 传输方程的简单解（比尔定律）：e的指数形式
遥感物理
第一章 基本概念 第二节 辐射传输 (radiance transfer) §1.2.1 传输方程 √ §1.2.2 源函数中散射的表达 §1.2.3 辐射传输方程的解
如果在s=0处的入射强度为Iλ(0)，则在s1处， 其射出强度可以通过对上式的积分获得：
s1
I(s1)I(0)ex pk ( d)s 0
8/14
假定介质消光截面均一不变，即kλ不依赖于距离s， 并定义路径长度：
s1
u 0 ds
则此时出射强度为：
I(s1)I(0)eku
1/8

第41卷第1期东北电力大学学报Vol.41，H 2021年2月Journal Of Northeast Electric Power University Feb,2021D O I： 10. 19718/j.issn. 1005-2992.2021-01-0048-08求解辐射传输方程的多松弛格子-Boltzmann模型刘晓川\王存海2,黄勇、朱克勇1(1.北京航空航天大学航空科学与工程学院，北京100191,2.北京科技大学能源与环境工程学院，北京100083)摘 要：针对福射传输方程，文中提出了一种多松弛格子-B o l t z m a n n模型（m u l t i p l e-r e l a x a t i o n-t i m el a t t i c e B o ltz m a n n m o d e l).基于扩散尺度下的M a x w e l l迭代，辐射传输方程可以严格地从格子B o l t z m a n n方程推导得出•一些数值案例用来验证本文提出的多松弛（M R T)格子-B o l t z m a n n模型的数值特性.结果表明本文提出的多松弛格子-B o l t z m a n n模型可以稳定及精确地求解参与性介质中的瞬态及稳态辐射传输问题.并且，该模型具有二阶精度.关键词：辐射传输方程;格子-B o l t z m a n n方法;多松弛模型中图分类号：T K121 文献标识码：A辐射传输方程描述了辐射能量在介质中的传递，在许多科学和工程领域具有重要作用，例如大气辐 射传输[1]、光学层析成像[2]、天体物理学[3]及核工程[4]等.辐射传输方程是一个高维、复杂的积分微分 方程，辐射强度涉及波长、时间、空间和角度等,求其解析解十分困难.学者们提出发展了很多种数值方 法来求解辐射传输方程,如蒙特卡洛法[5]，离散坐标法[6]，有限体积法[7]，有限元法[8]等.近年来，利用格子-Boltzmann方法（L B M)来求解辐射传输方程吸引了许多学者的兴趣.L B M起源 于格子气自动机，已经发展成为了一种计算流体力学的有力数值工具[9].并且,L B M已经被拓展到求解 许多线性和非线性系统问题，例如声子输运[1°],波传播[11]，反应扩散,对流扩散等.相比于其他的求解辐射传输方程的数值方法，L B M不需要计算大量的光线轨迹，也不需要离散复杂的偏微分方程. L B M具有容易实现，高并行效率等优点•目前，对于利用L B M来求解辐射方程还不完善，发展完善的 L B M用于求解辐射传输方程是必要的.1^3111^等[14]假定了可调节的虚拟光速和辐射平衡条件，将L B M推广到分析参与性介质中的辐射 问题.M a等[~基于辐射流体力学，提出了一维辐射的格子-Boltzmann模型.Zhang等[16]通过采用全隐 式后项差分格式处理辐射方程中的瞬态项，将L B M扩展到求解参与性介质中的一维瞬态辐射传输. Mink等[171在将P1近似应用辐射传输方程的基础上提出了一种三维的格子-Boltzmann模型，然而此模 型仅适用于光学厚介质.Y i等[18]通过引入虚拟的扩散项，将辐射传输方程视为一种特殊的对流扩散方 程，从而提出了一种二维稳态射传输方程的格子-Boltzmann模型.W a n g等[19_将瞬态辐射传输方程处 理为双曲守恒方程，然后提出了 一■种求解瞬态辅射和中子输运的格子-Boltzm ann模型.目前，求解辐射方程的多松他的格子-Boltzmann模型还未见报道.本文提出了一种多松她格子-Bo­ltzmann 模型 （multiple-relaxation-time l a t t i c e Boltzmann mode丨）■基于扩散尺度下的 Maxwell 迭代，福射传收稿日期：2020-11-09基金项目：国家自然科学基金（s is M o o w g c te o i4)第一作者：刘晓川（1992-),男，在读博士研究生，主要研究方向：航空科学与工程通讯作者：黄勇（1974-),男，博士，教授，主要研究方向：航空科学与工程电子邮箱：liuxiaochuan@(刘晓川），wangcunhai@ustb_(王存海），huangy@(黄勇），zhukeyong@buaa.edu_cn(朱克勇）第1期 刘晓川等：求解辐射传输方程的多松弛格子-Bohmiami模型 49输方程可以严格地从格子Boltzmann方程推导得出，并且不引人任何限制和近似.本文发展的多松弛格 子-Boltzmann模型可以精确地求解参与性介质内的多维瞬态及稳态辐射传输问题.数值结果表明该模 型具有二阶精度和收敛速率.并且，相比于单松弛模型，多松弛模型具有更好的稳定性.该模型可以进一 步推广到求解参与性介质内的辐射传输问题.1福射传输方程的多松她格子-B o l t z m a n n模型1-1辐射传输方程考虑吸收、发射和散射介质内的辐射传输方程，其离散坐标形式可以写为[2°]dl(r,rr,t) +f f, v/(r；i T^)+/3(r)l(r,f r,t)=S(r,n r,t)，(1)cLdt公式中心为介质内的光速;/为辐射强度;r为位置坐标冶+屹为衰减系数;/r + V")+ 为离散方向，源项S可以表示为s(r,n r,t)^kaib(r,{T,t)+^J j i{r,i r')(p(n r\n n)w m',(2)47T m，=1公式中A为总的离散方向，=1,2,…，八^' = 1,2,…，yv;M；m'为对应方向的权重.考虑漫发射和反射壁面，边界条件可以写为I(rw,{r,t)^e wIb(rw,t)+^-^I(rw ,f T') \nw -HT'\w m' + (\ - p j r\rw J F",t) ,(3)17 <Q公式中:&为发射率;Pu，为反射率;广‘为外部人射辐射强度.1.2 多松弛格子-Boltzm ann模型瞬态辐射常常发生于极短的时间内，在瞬态辐射的模拟中，通常引入无量纲时间来避免过小的时间 步长.将无量纲时间T心代人方程（1)中，得到时间无量纲形式的辐射传递方程[21]di(r,n",t) +L f f. v/(r)/2m,r)= F(r,/r,f*),⑷dt'公式中F{r,n r,r) = i R s{r,{r,r)-L^i(r,f r,r),(5)公式中:心为介质的参考长度.本文提出的时间无量纲形式的辐射传输方程的格子Boltzmann方程如下/(r + c^*,t*+A t')-/(r,t*)=--^(r.t*)] + A t'X),j(6)公式中:/(r，〇为分布函数;M为变换矩阵;S = 士叫U a,…人）为松弛参数矩阵，平衡函数的表达 式为r i(r,n r x)•跑-改2c?辐射强度可以由平衡函数给出，关系如下/(r，/r，〇=-(7)(8)50东北电力大学学报第41卷L B M方法中采用D m Q n格子模型，对于一维和二维问题，本文分别采用D1Q3和D2Q9模型.对于 D1Q3模型,其格子信息为[c〇,c, ,c2] =e；c = [0 1 -l]c,c [2/3,i =0c s=—,(0： = \ll/6,i = 1,2(9)(10)M0 12 一：-1对于D2Q9模型，其格子信息为(11)M C6,C7，C8] y =.0100 1-1-111-11-1l-l- i.c, (12)「4/9，/ =:0ccs = — 〇jt=,1/9,i =]，2,3,4(13).1/36,i=5,6,7,8111111111)-4-1- 1--122224-2- 2-2- 21111010-01—1一 110-20201-1-11•(14) 00 10-111-1-100 - 20211-1-101- 11-1000000 0001-11-h13从格子Boltzm ann方程到辐射传输方程本节基于扩散尺度=7(4幻2下的Maxwell迭代，不引入任何限制和假设，从多松弛格子- Boltzmarm模型严格推导得出辐射传输方程.这种扩散尺度是针对模型中的无量纲时间步长和空间步长 的尺度.首先，令/8U，r))f，w =(叫,叫，…，叫）'时间无量纲形式的 辐射传递方程(6)可以写成矢量形式f(r + ciA t，,r+A t f)-/e9(r,〇] + A t'wF(r,t m),(15)方程（15)左边应用Taylor展开，其中微分算子D' 矩阵/(r + e,A%,«* +y{Ax)2)~^ (A x)*£)s/(r,t*),s = 1〇,=y(E,dx+ E yd y)p(ydt'),P*^s p\q\’A s d i a M e o y e m…，e8,J…，e^)，(16)(17)(18)第1期刘晓川等：求解辐射传输方程的多松弛格子-Bohzmami 模型51公式中和g 均为非负整数.令m = M •/>〃 = M •/%，将Taylor展开形式代入方程（15)并整理得到00工(A x )sDsm =- S [m - me tf] + y (A x )2FMco ,s= 1其中D ^-M D ^-y CE,SX +E ,3y V (y s r yI'*^sp \ q \E t =ME M'E y =M E M '1 .(21)**jJ、’c o定义算子A = X (4幻]5\方程（19)可重新写为s= 1m =m e " -S 'Lm + y (A x )2FS~'Mu , (22)基于扩散尺度下的Maxwell1221迭代，从m° = m 〜开始，方程（19)经过三次迭代得到：m = m" -S ~'[A x D ' + (A x )2D2 + (A x )3D ,]me ，1 + [A x S ^D 1 + (A x )2S ^,Lf2]2ma ,-+7(4.«)2厂5_|财如 + 0((4x )4) ,(23)根据矢量方程(23)的第零项及各算子作用结果，可以得到辐射传递方程a /(r ,/7",f } +L r H" • V /(r ,/T ,<*) = F (r ,/T ,t *) +0((A x )2) ,(24)dt *至此，我们从多松弛格子-Boltzmann模型出发，基于扩散尺度下的Maxwell迭代，严格推导得出了辐射 传输方程，并且可以从方程(24)理论上得出该模型具有二阶的精度.一般而言，对于对流扩散问题，计 算流体力学等问题的L B 模型，其中的松弛系数与宏观方程中的扩散系数，流体黏性系数等有定量关系. 需要指出的是，根据从多松弛格子-Boltzmann模型严格推导得出辐射传输方程可知，本文提出的多松 弛格子-Boltzmann模型中的松弛参数均是自由的，与其他参数无关.对于一维和二维L B 模型，我们取 如下的松弛参数矩阵(19)(20)S = diag( 1 ,ir,l ) ,(25)S = diag(l,l,l,ir,l,5r,1,1,1) ,(26)对流扩散方程的多松弛L B 模型也采用了同样的处理方法，其中一维模型中的松弛参数，二维模型中 的松弛参数h 和s 5与扩散系数有关，而其他的松弛参数均取1.由于松弛矩阵中的松弛参数有无限种组 合方式，因此出于通用性考虑，我们选择了这种处理方法.同时需要指出的是当松弛参数矩阵中的松弛 系数相同时，多松弛模型退化到单松弛模型，即松弛矩阵中的松弛参数均为V2 结果及分析2.1 具有高斯型发射场的一维无限大平板考虑一充满吸收发射性介质的一维无限大平板内的辐射传递问题，平板内具有一高斯型发射场，该 问题由如下方程控制^+/3l ^e -u -b )2/a 2,z,b e [0,1] ,(27)考虑如下边界条件52东北电力大学学报第41卷I (〇,〇 =f 3-]e -b 2/a \ ^>0,(28)该问题存在解析解形式，其表达式如下/(2〇 =/(0，f )e x p ( —,)| 2 - (^ + 6)}X [erf (|+^)-erf (f +a )l ^>0, (29)考虑方向f = 1. 〇，a =〇• 02,6 = 0. 5，采用L B M 来模拟衰减系数为/3 = 1，10和50 时介质内辐射强度的分布，取1〇〇个格子，无量纲时间步长取土‘ =0.000 1，单松弛模型得到的结果和解析解对比，如图1所示,L B M 得到的辐射强度分布和解析解得到的辐射强度分布吻合地很好.接下来，我们进一步研究一维多松弛模型的稳 〇.〇4定性和精度.为了研究稳定性，我们考虑衰减系数为 10 nT1的情况，取100个格子，研究不同松弛参数下|所允许的最大时间步长.数值解和解析解的相对误| 〇.〇2 差定乂为Er = ^-------------(30)丨稳定性标准为数值解和解析解的相对误差小于 10'表1给出了不同松弛参数下所允许的最大时间 步长，不同参数的最大时间步长得到是根据我们定 义的稳定性标准，然后通过数值实验得到的，可以发现多松弛模型允许的最大时间步长可以随松弛参数调整，尤其当松弛参数小于1时，所允许的时间步长 大于单松弛模型，结果表明相比单松弛模型，多松弛模型可以在更大的时间步长内保持稳定，具有更好 的稳定性•多松弛模型的碰撞过程发生在矩空间，与多个速度分布函数相关联，相比单松弛模型发生在 速度空间的碰撞，多松弛模型本身在稳定性方面展现了很大的优势,数值结果证明了多松弛模型在稳定 性上的优势.此外，表2给出了不同格子数下单松弛和多松弛模型的相对误差，可以看出多松弛模型相 比单松弛模型具有更高的精度.图1衰减系数为卢=1,丨〇和501^时LBM 得到的辐射强度分布和解析解对比表1衰减系数/3 = 1〇 n T 1,100个格子下，单松弛（B G K )和多松她（M R T )模型允许的最大时间步长sr =0• 6sr =0. 8 sr =l.O(BGK)sr = 1 • 2sr = l • 4y m a x 22.618.413.28.2 4.1W ax2.26 e-31 • 84 e-31.32 e-30. 82 e-30.41e-3表2衰减系数/3 = 10n不同格子数下,单松弛（B G K )和多松弛（M R T )模型的相对误差格子数sr =0. 65r =0. 85r = l .2sr = 1.4BGK MRT BGKMRT BGK MRT BGK MRT 100 4.24 e-27.72 e-3 1. 14 e -2 3.09 e-3 4. 14 e-3 1.71 e-3 5.79 e-3 3.02 e-3150 1.82 e-2 3.24 e-3 4.84 e-3 1.25 e-3 1.85 e-37.77 e-4 2.54 e-3 1.35 e-32001.01 e-21.79 e-32.68 e~36.83 e—41.04 e-34. 40 e-41.42 e-37.57 e -42.2受高斯型脉冲照射的一维纯散射介质考虑厚度为1 m 的一维半透明平板介质内的瞬态辐射传输问题.介质为各向同性散射，壁面和 介质温度均为〇K，无发射.介质边界为透明边界，环境为真空.平板介质的衰减系数为1 nT1，右侧边界 无照射，左侧边界受到如下法向平行光人射辐射的照射：第1期刘晓川等：求解辐射传输方程的多松弛格子-Boltoimmi 模型53lp(0,t ) = /〇exp [//(〇 ,(31)公式中:/。

一维波动方程的解法波动现象是自然界和人类生活中广泛存在的一种现象，它具有许多重要的物理意义，例如声波、光波等。

一维波动方程是描述波动的重要方程之一，本文将介绍一维波动方程的解法。

一、一维波动方程的基本形式和意义一维波动方程的基本形式为：$$\frac{\partial^2 u}{\partial t^2}-c^2\frac{\partial^2 u}{\partialx^2}=0$$其中，$u(x,t)$表示波动的幅度，$c$表示波速。

这个方程描述了介质中的一种波动现象：波动传播速度为$c$，波动在媒质中沿$x$轴方向的传播，波动的幅度随时间$t$的变化而变化。

在声波和电磁波中，$u$分别是空气压力和电场强度，$c$分别是声速和光速。

二、1. 分离变量法分离变量法是一种基本的解法，其思想是将波动方程中的未知函数$u(x,t)$表示成仅包含$x$的函数和仅包含$t$的函数的乘积形式：$$u(x,t)=X(x)T(t)$$将$u(x,t)$代入一维波动方程中，得到：$$\frac{T''(t)}{c^2T(t)}=\frac{X''(x)}{X(x)}=-\lambda^2$$其中，$\lambda$是一个常数。

由此可得到两个关于未知函数的简单微分方程：$$T''(t)+\lambda^2c^2T(t)=0$$和$$X''(x)+\lambda^2X(x)=0$$其中，第一个微分方程的解为：$$T(t)=A\cos(\lambda ct)+B\sin(\lambda ct)$$其中，A、B是常数。

第二个微分方程的解为：$$X(x)=C\cos(\lambda x)+D\sin(\lambda x)$$其中，C、D是常数。

因此，一维波动方程的通解为：$$u(x,t)=\sum_{n=1}^\infty a_n\cos(\lambda_n x)\cos(\lambda_n ct)+b_n\sin(\lambda_n x)\sin(\lambda_n ct)$$其中，$\lambda_n=n\pi/L$，$L$为介质的长度，$a_n$和$b_n$是待定常数。

中南林业科技大学偏微分方程数值解法学生姓名：***学号：********学院：理学院专业年级：08信计1班设计题目：一维热传导方程的Richardson格式2011年06月一． 问题介绍考虑一维热传导方程：（1） ,0),(22T t x f xu a t u ≤<+∂∂=∂∂ 其中a 是正常数，)(x f 是给定的连续函数。

按照定解条件的不同给法，可将方程（1）的定解问题分为两类：第一类、初值问题（也称Cauthy 问题）：求具有所需次数偏微商的函数),(t x u ，满足方程（1）（∞<<∞-x ）和初始条件：（2）),()0,(x x u ϕ= ∞<<∞-x 第二类、初边值问题（也称混合问题）：求具有所需次数偏微商的函数),(t x u ，满足方程（1）（l x <<0）和初始条件：（3） ),()0,(x x u ϕ= l x <<0及边值条件(4).0),(),0(==t l u t u T t ≤≤0 假定)(x ϕ在相应区域光滑，并且在l x ,0=满足相容条件，使上述问题有唯一充分光滑的解。

二． 区域剖分考虑边值问题（1），（4）的差分逼近。

去空间步长N l h /=和时间步长M T /=τ，其中N,M 都是正整数。

用两族平行直线：),,1,0(N j jh x x j ===),,1,0(M k k t t k ===τ 将矩形域}0;0{T t l x G ≤≤≤≤=分割成矩形网格，网格节点为),(k j t x 。

以h G 表示网格内点集合，即位于开矩形G 的网点集合；h G 表示所有位于闭矩形G 的网点集合；h Γ=h G --h G 是网格界点集合。

三． 差分格式第k+1层值通过第k 层值明显表示出来，无需求解线性代数方程组，这样的格式称为显格式。

第k+1层值不能通过第k 层值明显表示出来，而由线性代数方程组确定，这样的格式称为隐格式。

用一维辐射传递方程计算二氧化碳的辐射强迫
黄晓璜;崔国民;华泽钊;徐家良
【期刊名称】《上海理工大学学报》
【年(卷),期】2013(035)006
【摘要】为得到二氧化碳浓度变化影响下的辐射强迫值,建立了简化的大气层净辐射通量模型.基于辐射传递方程,这个模型可以计算二氧化碳浓度变化引起的大气层每一层净辐射通量以及辐射强迫值.与美国大气与环境研究中心(AER)的RRTM-LW 长波辐射传输模型在同一因素下的计算结果进行比较,得到误差不超过1％.
【总页数】5页(P599-602,613)
【作者】黄晓璜;崔国民;华泽钊;徐家良
【作者单位】上海理工大学能源与动力工程学院,上海 200093;上海理工大学能源与动力工程学院,上海 200093;上海理工大学能源与动力工程学院,上海 200093;上海气象局上海气候中心,上海 200030
【正文语种】中文
【中图分类】TM344.1
【相关文献】
1.基于地面短波辐射观测资料计算淡积云的辐射强迫及水平尺度方法研究 [J], 施红蓉;陈洪滨;夏祥鳌;范学花;张金强;曹杨;朱彦良;王姝
2.一维辐射-对流模式对云辐射强迫的数值模拟研究 [J], 刘玉芝;石广玉;赵剑琦
3.具有漫反射边界一维灰性平行平板介质中辐射传递方程的求解 [J], 程强;周怀春
4.半透明介质中辐射传递方程的反演计算及数值模拟 [J], 谈和平;余其铮;阮立明;
夏新林
5.DRESOR法对含有边界的一维各向异性散射介质辐射传递方程的求解 [J], 黄志锋;周怀春;程强;何诚
因版权原因，仅展示原文概要，查看原文内容请购买。

一维散射问题一维散射问题是量子力学中的重要研究领域，它描述的是粒子在一维势场中的传播和散射行为。

本文将从波函数的表达、散射系数的计算以及散射振幅的求解等方面来讨论一维散射问题。

一、波函数的表达在一维散射问题中，粒子的波函数可以通过薛定谔方程来描述。

对于自由粒子来说，波函数的表达式可写为：ψ(x,t) = Ae^ikx - A*^ikx其中A和A*分别是入射和反射波的振幅，k是粒子的波数，x是位置，t是时间。

而当粒子遇到势场时，波函数的表达式可以分为三个区域来讨论：入射区、势场区和散射区。

分别用ψi(x,t)、ψv(x,t)和ψr(x,t)表示，其中ψi(x,t)是入射波函数，表示粒子从远处靠近势场；ψv(x,t)是势场区的波函数，描述粒子在势场中的传播；ψr(x,t)是反射波函数，描述粒子从势场反射出来。

它们的表达式分别为：ψi(x,t) = Ae^ikxψv(x,t) = Be^-qx + Ce^qxψr(x,t) = De^-ikx其中B、C和D为待定系数，q为粒子在势场中的衰减常数。

二、散射系数的计算在一维散射问题中，我们通常关注的是粒子的散射行为，其中一个重要的参数是散射系数。

散射系数反映了入射粒子被散射到反射态的概率。

散射系数的计算可以通过比较反射波与入射波的振幅来进行。

根据波函数的连续性条件，在势场边界处，入射波的振幅与反射波的振幅之比等于散射系数的平方。

即：|R/A|^2 = |D/A|^2散射系数也可以通过波函数的归一化条件来计算。

根据波函数在正负无限远处的衰减性质，我们可以得到：1 = |A|^2 + |D|^2综合以上两个方程，我们可以求解出散射系数。

三、散射振幅的求解除了散射系数，散射振幅也是一维散射问题中需要研究的重要参数。

散射振幅描述了入射波进入势场后，被势场散射出来的波的振幅和相位。

散射振幅的求解可以通过势场边界条件和薛定谔方程进行。

根据薛定谔方程，在势场边界处，波函数的导数应该是连续的。