CN103136436B - 考虑泥沙吸附的核素输移扩散分析方法 - Google Patents

Abstract

本发明属于核素输移分析方法，具体涉及一种考虑泥沙吸附的核素输移扩散分析方法。该方法采用改进的平面二维浅水环流模型进行流速场和浓度场的模拟，引入了核素在沉积物和液相中的交换的因素，考虑了泥沙吸附对分配系数不同的核素在水中的浓度及在河床的沉积量的具体影响，采用实测的分配系数值求解浓度方程。本发明与现有的核素输移扩散分析方法相比更符合核素输移扩散的实际情况，分析结果更加准确。

Description

考虑泥沙吸附的核素输移扩散分析方法

技术领域

本发明属于核素输移分析方法，具体涉及一种考虑泥沙吸附的核素输移扩散分析方法。

背景技术

核设施运行期间排放的液态流出物中含有的放射性核素，进入受纳水体后随水体输运和迁移，当水动力条件不变时，核素在水和泥沙之间达到动态平衡。核素被水体中和河床的泥沙吸附，并沉积于河床底部，降低了水体中核素浓度；当水动力条件发生改变，吸附在泥沙中的核素又可以从吸附态转移到水体中，或将沉积在床沙中的核素随着泥沙的冲刷、再悬浮重新进入水体，增加水体中核素浓度。泥沙与水流共同成为核素输移的主要载体，影响着核素在水体中的迁移转化过程。

平面二维浅水环流模型作为一远区模型，未涉及温度分层的问题，一般假定密度不随温度发生变化，特别适用大范围内的水温及流速分布情况的计算，因而被广泛应用在冷却水工程和环境工程中。

核设施选址和运行阶段，需评价液态流出物对环境产生的影响。目前，国内外基本上都是通过采用数模计算方法进行研究与分析，计算中一般只考虑水流中核素的输运和稀释，而忽略了在沉积物和液相中的交换。这种计算方法跟核素随水体输运和迁移的实际情况有所差异，需要对平面二维浅水环流模型加以改进。

发明内容

本发明的目的在于针对现有技术的缺陷，提供一种考虑泥沙吸附的核素输移扩散分析方法，使核素的输移扩散分析结果更加准确。

本发明的技术方案如下：一种考虑泥沙吸附的核素输移扩散分析方法，该方法建立的平面二维浅水环流的基本方程如下：

沿水深平均的连续方程：

$\frac{\∂\mathrm{\ξ}}{\∂t}+\frac{\∂\left(\mathrm{Hu}\right)}{\∂x}+\frac{\∂\left(\mathrm{Hv}\right)}{\∂y}=0$

沿水深平均的动量方程：

$\frac{\∂u}{\∂t}+u\frac{\∂u}{\∂x}+v\frac{\∂u}{\∂y}=-g\frac{\∂\mathrm{\ξ}}{\∂x}-\mathrm{fv}-\frac{\mathrm{gu}}{C^{2}H}\sqrt{u^2+v^2}+\frac{{\mathrm{\τ}}_{\mathrm{sx}}}{\mathrm{\ρH}}+$

$\frac{1}{H}\frac{\∂}{\∂x}\left(H{\mathrm{\ϵ}}_x\frac{\∂u}{\∂x}\right)+\frac{1}{H}\frac{\∂}{\∂y}\left(H{\mathrm{\ϵ}}_y\frac{\∂u}{\∂y}\right)$

$\frac{\∂v}{\∂t}+u\frac{\∂v}{\∂x}+v\frac{\∂v}{\∂y}=-g\frac{\∂\mathrm{\ξ}}{\∂y}+\mathrm{fu}-\frac{\mathrm{gv}}{C^{2}H}\sqrt{u^2+v^2}+\frac{{\mathrm{\τ}}_{\mathrm{sy}}}{\mathrm{\ρH}}+$

$\frac{1}{H}\frac{\∂}{\∂x}\left(H{\mathrm{\ϵ}}_x\frac{\∂v}{\∂x}\right)+\frac{1}{H}\frac{\∂}{\∂y}\left(H{\mathrm{\ϵ}}_y\frac{\∂v}{\∂y}\right)$

浓度方程：

$\frac{\∂C}{\∂t}+u\frac{\∂C}{\∂x}+v\frac{\∂C}{\∂y}=\frac{1}{H}\frac{\∂}{\∂x}\left({\mathrm{HD}}_x\frac{\∂C}{\∂x}\right)+\frac{1}{H}\frac{\∂}{\∂y}\left({\mathrm{HD}}_y\frac{\∂C}{\∂y}\right)-\mathrm{\λC}+{\mathrm{SK}}_d+S_{\mathrm{in}}$

在以上方程中，

ζ--相对基准面水位；

H--水深；

t--时间；

u，v--水深平均流速；

x，y--坐标；

g--重力加速度；

f--柯氏力系数；

C--谢才系数；

τ_sx，τ_sy--表面风应力在x，y方向的分力；

ρ--水的密度；

ε_x，ε_y--x，y方向的广义涡粘性系数；

D_x、Dy--x、y方向核素的综合扩散系数；

S--含沙量；

K_d--核素分配系数；

S_in--挟沙力；

λ--衰变常数；

利用分步杂交法对上述深度平均的二维浅水环流模型进行离散化，采用实测的核素分配系数值求解浓度方程。

进一步，如上所述的考虑泥沙吸附的核素输移扩散分析方法，其中，所述的核素分配系数的实测方法是：采用放射性液态流出物受纳水体的水样和河床沉积物固体样进行实验，将样品与已知浓度的含有一种放射性核素的水溶液在振荡下长时间混合，直到水溶液中核素浓度稳定为止；然后进行离心使固液相分离；最后分别测定固相和液相中的核素浓度，通过计算固、液相中的核素浓度之比即可得出核素分配系数K_d。

本发明的有益效果如下：本发明采用改进的平面二维浅水环流模型进行流速场和浓度场的模拟，引入了核素在沉积物和液相中的交换的因素，即SK_d项，考虑了泥沙吸附对分配系数不同的核素在水中的浓度及在河床的沉积量的具体影响，采用实测的分配系数值求解浓度方程。本发明与现有的核素输移扩散分析方法相比更符合核素输移扩散的实际情况，分析结果更加准确。

附图说明

图1为深度平均的坐标系统示意图；

图2为计算对流部分的特征线示意图；

图3为迎风单元示意图；

图4为p_i点的集中质量区域示意图；

具体实施方式

下面结合附图和实施例对本发明进行详细的描述。

平面二维浅水环流模型作为一远区模型，未涉及温度分层的问题，一般假定密度不随温度发生变化，特别适用大范围内的水温及流速分布情况的计算，因而被广泛应用在冷却水工程和环境工程中。

平面二维浅水环流模型从不可压缩流体运动的基本方程--N-S方程出发，对于河口或大型水库，往往垂向加速度与重力加速度相比很小，可以略去，假定压强沿水深的分布为静压分布，同时考虑地球自转引起的哥氏力的作用，因而可得到垂向静压分布假定的三维流动基本方程为：

$\frac{\∂u}{\∂x}+\frac{\∂v}{\∂y}+\frac{\∂w}{\∂z}=0---(1-1)$

$\frac{\∂u}{\∂t}+u\frac{\∂u}{\∂x}+v\frac{\∂u}{\∂y}+w\frac{\∂u}{\∂z}=-\frac{1}{\mathrm{\ρ}}\frac{\∂p}{\∂x}-\mathrm{fv}+\frac{\mathrm{\μ}}{\mathrm{\ρ}}{\▿}^{2}u---(1-2)$

$\frac{\∂v}{\∂t}+u\frac{\∂v}{\∂x}+v\frac{\∂v}{\∂y}+w\frac{\∂v}{\∂z}=-\frac{1}{\mathrm{\ρ}}\frac{\∂p}{\∂x}+\mathrm{fu}+\frac{\mathrm{\μ}}{\mathrm{\ρ}}{\▿}^{2}v---(1-3)$

$\frac{\∂p}{\∂z}=-\mathrm{\ρg}---(1-4)$

式中f为哥氏力系数，f＝2ωsinφ，ω为地球自转角速度，φ为当地纬度。

本发明的模型是建立在假设水体为不可压缩，压强沿水深为静压分布的情况下，模拟水域的水平尺度远大于垂直尺度，且水平流速远大于垂向流速，这些物理量沿水深方向的变化相对沿水平方向的变化要小的多，可略去这些量沿水深方向的变化，将三维流动的基本方程组沿水深方向积分，再沿水深取平均，应用莱布尼兹公式进行变换，得到沿水深平均的二维流动基本方程。

1、沿深度积分

应用莱布尼兹公式变换

$\frac{\∂}{\∂x_k}{\∫}_{-h\left(x_k\right)}^{\mathrm{\ξ}\left(x_k\right)}f(x_1,x_2,x_3){\mathrm{dx}}_3={\∫}_{-h}^{\mathrm{\ξ}}\frac{\∂f}{\∂x_k}{\mathrm{dx}}_3+f{|}_{\mathrm{\ξ}}\frac{\∂\mathrm{\ξ}}{\∂x_k}-f{|}_{-h}\frac{\∂(-h)}{\∂x_k}---(1-5)$

2、定义下列诸量进行深度平均(坐标系统见图1)：

水体总体深度H＝h+ξ

$\stackrel{\‾}{\mathrm{\Φ}}=\frac{1}{H}{\∫}_{-h}^{\mathrm{\ξ}}\mathrm{\Φ}{\mathrm{dx}}_3---(1-6)$

(Φ表示U、V、P、T等变量，“-”表示平均)

其中：

h-为平均水面以下的水深。

ξ-相对于平均水面的水位或潮位。

x方向的平均速度分量： $u=\frac{1}{h+\mathrm{\ξ}}{\∫}_{-h}^{\mathrm{\ξ}}\mathrm{udz}---(1-7)$

y方向的平均速度分量： $v=\frac{1}{h+\mathrm{\ξ}}{\∫}_{-h}^{\mathrm{\ξ}}\mathrm{vdz}---(1-8)$

对某一确定点的任意深度流速值可以表示为垂向平均流速加上一脉动值。即：

u＝u+u′；v＝v+v′

且这些脉动值满足：

${\∫}_{-h}^{\mathrm{\ξ}}{u}^{\′}\mathrm{dz}=0;$ ${\∫}_{-h}^{\mathrm{\ξ}}{v}^{\′}\mathrm{dz}=0$

3、满足底部及水面运动条件

底部阻力沿x方向分量： ${\mathrm{\τ}}_x^b=\frac{g}{C^2}\mathrm{\ρu}{(u^2+v^2)}^{1/2};---(1-9)$

底部阻力沿y方向分量： ${\mathrm{\τ}}_y^b=\frac{g}{C^2}\mathrm{\ρv}{(u^2+v^2)}^{1/2};---(1-10)$

v^b＝0；u^b＝0；

其中，谢才系数n为底部糙率系数。

$W{|}_{\mathrm{\ξ}}=\frac{\∂\mathrm{\ξ}}{\∂t}+u{|}_{\mathrm{\ξ}}\frac{\∂\mathrm{\ξ}}{\∂x}+v{|}_{\mathrm{\ξ}}\frac{\∂\mathrm{\ξ}}{\∂y}---(1-11)$

其中，u、v分别为x、y方向流速。

将式(1-1)沿水深积分并取平均得：

$\frac{1}{H}{\∫}_{-h}^{\mathrm{\ξ}}\frac{\∂u}{\∂x}\mathrm{dz}+\frac{1}{H}{\∫}_{-h}^{\mathrm{\ξ}}\frac{\∂v}{\∂y}\mathrm{dz}+\frac{1}{H}{\∫}_{-h}^{\mathrm{\ξ}}\frac{\∂w}{\∂z}\mathrm{dz}=0---(1-12)$

根据前面的假定并应用莱布尼兹公式对方程(1-12)中的各项进行变换得：

${\∫}_{-h}^{\mathrm{\ξ}}\frac{\∂u}{\∂x}\mathrm{dz}=\frac{\∂\left(\mathrm{hu}\right)}{\∂x}-u{|}_{\mathrm{\ξ}}\frac{\∂\mathrm{\ξ}}{\∂x}+u{|}_{-h}\frac{\∂(-h)}{\∂x}---(1-13)$

${\∫}_{-h}^{\mathrm{\ξ}}\frac{\∂v}{\∂y}\mathrm{dz}=\frac{\∂\left(\mathrm{hv}\right)}{\∂y}-v{|}_{\mathrm{\ξ}}\frac{\∂\mathrm{\ξ}}{\∂y}+v{|}_{-h}\frac{\∂(-h)}{\∂y}---(1-14)$

${\∫}_{-h}^{\mathrm{\ξ}}\frac{\∂w}{\∂z}\mathrm{dz}=w{|}_{\mathrm{\ξ}}-w{|}_{-h}---(1-15)$

由底部及水面运动条件知：u|_-h＝0、v|_-h＝0、w|_-h＝0

$w{|}_{\mathrm{\ξ}}=\frac{\mathrm{dH}}{\mathrm{dt}}=\frac{\∂H}{\∂t}+u{|}_{\mathrm{\ξ}}\frac{\∂H}{\∂x}+v{|}_{\mathrm{\ξ}}\frac{\∂H}{\∂y}---(1-16)$

将(1-13)(1-14)(1-15)(1-16)代入式(1-12)，则得到沿水深平均的连续方程

$\frac{\∂\mathrm{\ξ}}{\∂t}+\frac{\∂\left(\mathrm{Hu}\right)}{\∂x}+\frac{\∂\left(\mathrm{Hv}\right)}{\∂y}=0---(1-17)$

用同样的方法可将式(1-2)、(1-3)沿水深积分并取平均，同时考虑底部阻力和水面风应力的作用，得到沿水深平均的动量方程：

$\frac{\∂u}{\∂t}+u\frac{\∂u}{\∂x}+v\frac{\∂u}{\∂y}=-g\frac{\∂\mathrm{\ξ}}{\∂x}-\mathrm{fv}-\frac{\mathrm{gu}}{C^{2}H}\sqrt{u^2+v^2}+\frac{{\mathrm{\τ}}_{\mathrm{sx}}}{\mathrm{\ρH}}+---(1-18)$

$\frac{1}{H}\frac{\∂}{\∂x}\left(H{\mathrm{\ϵ}}_x\frac{\∂u}{\∂x}\right)+\frac{1}{H}\frac{\∂}{\∂y}\left(H{\mathrm{\ϵ}}_y\frac{\∂u}{\∂y}\right)$

$\frac{\∂v}{\∂t}+u\frac{\∂v}{\∂x}+v\frac{\∂v}{\∂y}=-g\frac{\∂\mathrm{\ξ}}{\∂y}+\mathrm{fu}-\frac{\mathrm{gv}}{C^{2}H}\sqrt{u^2+v^2}+\frac{{\mathrm{\τ}}_{\mathrm{sy}}}{\mathrm{\ρH}}+---(1-19)$

$\frac{1}{H}\frac{\∂}{\∂x}\left(H{\mathrm{\ϵ}}_x\frac{\∂v}{\∂x}\right)+\frac{1}{H}\frac{\∂}{\∂y}\left(H{\mathrm{\ϵ}}_y\frac{\∂v}{\∂y}\right)$

浓度方程：

$\frac{\∂C}{\∂t}+u\frac{\∂C}{\∂x}+v\frac{\∂C}{\∂y}=\frac{1}{H}\frac{\∂}{\∂x}\left({\mathrm{HD}}_x\frac{\∂C}{\∂x}\right)+\frac{1}{H}\frac{\∂}{\∂y}\left({\mathrm{HD}}_y\frac{\∂C}{\∂y}\right)-\mathrm{\λC}+{\mathrm{SK}}_d+S_{\mathrm{in}}---(1-20)$

变量及参数说明：

ζ为相对基准面水位(m)；

h为基准面以下水深(m)；

H为水深(m)；

t为时间(s)；

ρ为水的密度(kg/m³)；

ρ_a为空气的密度(kg/m³)；

g为重力加速度(m/s²)；

T为水体超温(℃)；

T’为水体温度(℃)；

T_∞为自然水温(℃)；

K_s为水面综合散热系数(W/m²·℃)；

C为谢才系数(m^1/2/s)；

n为糙率系数(s/m^1/3)；

R为水力半径(m)；

V_W为风速(m/s)；

θ为风向角度(rad)；

f为柯氏力系数(s^-1)；

ω为地球自转角速度(rad/s)；

φ为当地的纬度(rad)；

C_D为与风速有关的无量纲系数；

C_P为定压比热(W·s/kg·℃)；

τ_s为表面风应力(N/m²)；

τ_sxτ_sy分别为表面风应力在x，y方向的分力(N/m²)；

ε_xε_y分别为x，y方向的广义涡粘性系数(m²/s)；

K_xK_y分别为x，y方向的广义热扩散系数(m²/s)；

D_x、Dy分别为x、y方向核素的综合扩散系数；

S为含沙量；

K_d为核素分配系数；

S_*为挟沙力；

x，y为坐标(m)；

u，v为水深平均流速(m/s)；

λ--衰变常数。

至此，式(1-17)、(1-18)、(1-19)、(1-20)为平面二维浅水环流的基本方程。

定解条件如下：

(1)初始条件

u(x、y、0)＝u₀(x、y)；v(x、y、0)＝v₀(x、y)；

T(x、Y、0)＝T₀(x、y)；

(2)边界条件

陆地边界：采用滑移条件。即

水边界：由实测资料推求，给出潮位变化过程线。

为了克服在求解大数下的对流扩散问题，大数下的不可压缩粘性流所遇到的困难，以及消除在计算传质传热问题时所经常出现的负度与温负温升，吴江航提出了分步杂交法，分步杂交法采用三角形网格系统，将计算的每一时间步长分成两步进行。其特点是在不规则的三角形网格上，对于对流和扩散算子利用分步的技巧，各自采用最适合它们的方法进行计算。在前半分步，对于对流算子采用特征线法，在后半分步，对于扩散算子，采用集中质量的有限元法进行计算，计算潮流时，采用有限体积法求解水位方程。

鉴于分步杂交法算法简单、稳定性好、能适应不规则的几何形状和各种边界条件等优点，本发明采用分步杂交法对方程进行离散。

在对流作用较强的近岸海域，对流效应远大于紊动扩散效应，涡粘性系数和扩散系数的取值对计算结果的敏感性并不大，因而在计算过程中经常对涡粘性系数进行适当的处理，使ε_x＝ε_y＝ε，其变化与水深流速有关。

由于因而可以将方程(1-18)(1-19)(1-20)简化成：

$\frac{\∂u}{\∂t}+u\frac{\∂u}{\∂x}+v\frac{\∂u}{\∂y}=-g\frac{\∂\mathrm{\ξ}}{\∂x}-\mathrm{fv}-\frac{\mathrm{gu}}{C^{2}H}\sqrt{u^2+v^2}+\frac{{\mathrm{\τ}}_{\mathrm{sx}}}{\mathrm{\ρH}}+\mathrm{\ϵ}(\frac{{\∂}^{2}u}{\∂x^2}+\frac{{\∂}^{2}u}{\∂y^2})---(2-1)$

$\frac{\∂v}{\∂t}+u\frac{\∂v}{\∂x}+v\frac{\∂v}{\∂y}=-g\frac{\∂\mathrm{\ξ}}{\∂y}+\mathrm{fu}-\frac{\mathrm{gv}}{C^{2}H}\sqrt{u^2+v^2}+\frac{{\mathrm{\τ}}_{\mathrm{sy}}}{\mathrm{\ρH}}+\mathrm{\ϵ}(\frac{{\∂}^{2}v}{\∂x^2}+\frac{{\∂}^{2}v}{\∂y^2})---(2-2)$

$\frac{\∂C}{\∂t}+u\frac{\∂C}{\∂x}+v\frac{\∂C}{\∂y}=\frac{1}{H}\frac{\∂}{\∂x}\left({\mathrm{HD}}_x\frac{\∂C}{\∂x}\right)+\frac{1}{H}\frac{\∂}{\∂y}\left({\mathrm{HD}}_y\frac{\∂C}{\∂y}\right)-\mathrm{\λC}+{\mathrm{SK}}_d+S_{\mathrm{in}}---(2-3)$

利用分步方法将方程式(2-1)、(2-2)、(2-3)分解为：

(1)在前半分步用特征线方法求解方程：

$\left\{\begin{array}{c}\frac{1}{2}\frac{\∂u^{\left(1\right)}}{\∂t}+u^{\left(1\right)}\frac{\∂u^{\left(1\right)}}{\∂x}+v^{\left(1\right)}\frac{\∂u^{\left(1\right)}}{\∂y}=0\\ \frac{1}{2}\frac{{\∂v}^{\left(1\right)}}{\∂t}+u^{\left(1\right)}\frac{\∂v^{\left(1\right)}}{\∂x}+v^{\left(1\right)}\frac{\∂v^{\left(1\right)}}{\∂y}=0\\ \frac{1}{2}\frac{\∂C^{\left(1\right)}}{\∂t}+u^{\left(1\right)}\frac{\∂C^{\left(1\right)}}{\∂x}+v^{\left(1\right)}\frac{\∂C^{\left(1\right)}}{\∂y}=0\end{array}\right.---(2-4)$

用改型特征线方法计算对流部分，式(2-4)的离散格式为：

$\left\{\begin{array}{c}{u}_i^{n+\frac{1}{2}}=\underset{\mathrm{\α}=1}{\overset{3}{\mathrm{\Σ}}}{\tilde{L}}_{\mathrm{\α}}^{\left(e_i\right)}{u}_{\mathrm{\α}}^{\left(e_i\right)}\left(\mathrm{n\Δt}\right)\\ v_i^{n+\frac{1}{2}}=\underset{\mathrm{\α}=1}{\overset{3}{\mathrm{\Σ}}}{\tilde{L}}_{\mathrm{\α}}^{\left(e_i\right)}{v}_{\mathrm{\α}}^{\left(e_i\right)}\left(\mathrm{n\Δt}\right)\\ C_i^{n+\frac{1}{2}}=\underset{\mathrm{\α}=1}{\overset{3}{\mathrm{\Σ}}}{\tilde{L}}_{\mathrm{\α}}^{\left(e_i\right)}{C}_{\mathrm{\α}}^{\left(e_i\right)}\left(\mathrm{n\Δt}\right)\end{array}\right.---(2-5)$

其中：α＝1，2，3为在时刻通过p_i点的特征线

$\left\{\begin{array}{c}\frac{\mathrm{dx}}{\mathrm{dt}}=2u\\ \frac{\mathrm{dy}}{\mathrm{dt}}=2v\end{array}\right.---(2-6)$

在nΔt时间层上的端点的面积坐标，如图2所示， (α＝1，2，3)分别为nΔt时刻u，v，C在三角形单元(e_i)的三个结点上的值，(e_i)为点所在单元，它与p_i点可以不在一个单元内，这些面积坐标由特征线方程(2-6)解出：

$\left.\begin{array}{c}{\tilde{L}}_1^{\left(e_i\right)}=[(a_{22}{x}_i-a_{12}{y}_i)-f]/d\\ {\tilde{L}}_2^{\left(e_i\right)}=[(a_{11}{y}_i-a_{21}{x}_i)-g]/d\\ {\tilde{L}}_3^{\left(e_i\right)}=1-{\tilde{L}}_1^{\left(e_i\right)}-{\tilde{L}}_3^{\left(e_i\right)}\end{array}\right\}---(2-7)$

其中

$a_{11}=[x_1^{\left(e_i\right)}-x_3^{\left(e_i\right)}+(u_1^{\left(e_i\right)}\left(\mathrm{n\Δt}\right)-u_3^{\left(e_i\right)}\left(\mathrm{n\Δt}\right))\mathrm{\Δt}]$

$a_{12}=[x_2^{\left(e_i\right)}-x_3^{\left(e_i\right)}+(u_2^{\left(e_i\right)}\left(\mathrm{n\Δt}\right)-u_3^{\left(e_i\right)}\left(\mathrm{n\Δt}\right))\mathrm{\Δt}]$

$a_{21}=[y_1^{\left(e_i\right)}-y_3^{\left(e_i\right)}+(v_1^{\left(e_i\right)}\left(\mathrm{n\Δt}\right)-v_3^{\left(e_i\right)}\left(\mathrm{n\Δt}\right))\mathrm{\Δt}]$

$a_{22}=[y_2^{\left(e_i\right)}-y_3^{\left(e_i\right)}+(v_2^{\left(e_i\right)}\left(\mathrm{n\Δt}\right)-v_3^{\left(e_i\right)}\left(\mathrm{n\Δt}\right))\mathrm{\Δt}]$

$b_1=x_3^{\left(e_i\right)}+u_3^{\left(e_i\right)}\left(\mathrm{n\Δt}\right)\mathrm{\Δt}$

$b_2=y_3^{\left(e_i\right)}+v_3^{\left(e_i\right)}\left(\mathrm{n\Δt}\right)\mathrm{\Δt}$

f＝a₂₂b₁-a₁₂b₂

g＝a₁₁b₂-a₁₂b₂

d＝a₁₁a₂₂-a₁₂a₂₁

在解上列各式时，利用了各物理量沿特征线为常数的性质与线形插值近似。要求α＝1，2，3。

(2)在后半分步用集中质量有限元方法求解方程：

$\left\{\begin{array}{c}\frac{1}{2}\frac{\∂u^{\left(2\right)}}{\∂t}=-g\frac{\∂\mathrm{\ξ}}{\∂x}-f{v}^{\left(2\right)}-\frac{g{u}^{\left(2\right)}}{C^{2}H}\sqrt{u^2+v^2}+\frac{{\mathrm{\τ}}_{\mathrm{sx}}}{\mathrm{\ρH}}+\mathrm{\ϵ}(\frac{{\∂}^2{u}^{\left(2\right)}}{\∂x^2}+\frac{{\∂}^2{u}^{\left(2\right)}}{\∂y^2})\\ \frac{1}{2}\frac{\∂v^{\left(2\right)}}{\∂t}=-g\frac{\∂\mathrm{\ξ}}{\∂y}+f{u}^{\left(2\right)}-\frac{g{v}^{\left(2\right)}}{C^{2}H}\sqrt{u^2+v^2}+\frac{{\mathrm{\τ}}_{\mathrm{sy}}}{\mathrm{\ρH}}+\mathrm{\ϵ}(\frac{{\∂}^2{v}^{\left(2\right)}}{\∂x^2}+\frac{{\∂}^2{v}^{\left(2\right)}}{\∂y^2})\\ \frac{1}{2}\frac{\∂C^{\left(2\right)}}{\∂t}=\frac{\∂}{\∂x}\left(K_x\frac{\∂C^{\left(2\right)}}{\∂x}\right)+\frac{\∂}{\∂y}\left(K_y\frac{\∂C^{\left(2\right)}}{\∂y}\right)-\frac{K_s{C}^{\left(2\right)}}{\mathrm{\ρ}{C}_{P}H}\end{array}\right.---(2-8)$

对(2-8)先取时间半隐式差分：

$\left\{\begin{array}{c}\frac{u^{n+1}-u^{n+\frac{1}{2}}}{2\frac{\mathrm{\Δt}}{2}}=-g\frac{\∂{\mathrm{\ξ}}^{n+\frac{1}{2}}}{\∂t}-f{v}^{n+1}-g{\left(\frac{\sqrt{u^2+v^2}}{C^{2}H}\right)}^{n+\frac{1}{2}}{u}^{n+1}+{\left(\frac{{\mathrm{\τ}}_{\mathrm{sx}}}{\mathrm{\ρH}}\right)}^{n+\frac{1}{2}}+\\ \mathrm{\ϵ}{(\frac{{\∂}^{2}u}{\∂x^2}+\frac{{\∂}^{2}u}{\∂y^2})}^{n+\frac{1}{2}}\\ \frac{v^{n+1}-v^{n+\frac{1}{2}}}{2\frac{\mathrm{\Δt}}{2}}=-g\frac{\∂{\mathrm{\ξ}}^{n+\frac{1}{2}}}{\∂t}+f{u}^{n+1}-g{\left(\frac{\sqrt{u^2+v^2}}{C^{2}H}\right)}^{n+\frac{1}{2}}{v}^{n+1}+{\left(\frac{{\mathrm{\τ}}_{\mathrm{sy}}}{\mathrm{\ρH}}\right)}^{n+\frac{1}{2}}+\\ \mathrm{\ϵ}{(\frac{{\∂}^{2}v}{\∂x^2}+\frac{{\∂}^{2}v}{\∂y^2})}^{n+\frac{1}{2}}\\ \frac{C^{n+1}-C^{n+\frac{1}{2}}}{2\frac{\mathrm{\Δt}}{2}}-\frac{\∂}{\∂x}\left(K_x\frac{\∂C^{n+\frac{1}{2}}}{\∂x}\right)+\frac{\∂}{\∂y}\left(K_y\frac{\∂C^{n+\frac{1}{2}}}{\∂y}\right)-{\left(\frac{K_s}{\mathrm{\ρ}{C}_{P}H}\right)}^{n+\frac{1}{2}}{C}^{n+\frac{1}{2}}\end{array}\right.---(2-9)$

再利用集中质量的有限元方法计算(2-9)式：得

$\left\{\begin{array}{c}{a}_{11}{u}_i^{n+1}+a_{12}{v}_i^{n+1}=u_i^{n+\frac{1}{2}}-\mathrm{\ϵ}\frac{\mathrm{\Δt}}{A_i}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{b}_{\mathrm{ij}}{u}_j^{n+\frac{1}{2}}-\frac{\mathrm{g\Δt}}{A_i}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{c}_{\mathrm{ij}}{\mathrm{\ξ}}_j^{n+\frac{1}{2}}+f_{\mathrm{xi}}^{n+\frac{1}{2}}\mathrm{\Δt}\\ a_{21}{u}_i^{n+1}+a_{22}{v}_i^{n+1}=v_i^{n+\frac{1}{2}}-\mathrm{\ϵ}\frac{\mathrm{\Δt}}{A_i}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{b}_{\mathrm{ij}}{v}_j^{n+\frac{1}{2}}-\frac{\mathrm{g\Δt}}{A_i}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{d}_{\mathrm{ij}}{\mathrm{\ξ}}_j^{n+\frac{1}{2}}+f_{\mathrm{yi}}^{n+\frac{1}{2}}\mathrm{\Δt}\\ (1+{\left(\frac{K_S}{\mathrm{\ρ}{c}_{p}H}\right)}_i^{n+\frac{1}{2}}\mathrm{\Δt})T_i^{n+1}=T_i^{n+\frac{1}{2}}-\frac{\mathrm{\Δt}}{A_i}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}(k_x{l}_{\mathrm{ij}}+k_y{m}_{\mathrm{ij}})T_j^{n+\frac{1}{2}}\end{array}\right.---(2-10)$

其中

$a_{11}=1+g{\left(\frac{\sqrt{u^2+v^2}}{C^{2}H}\right)}_i^{n+\frac{1}{2}}\mathrm{\Δt}$ a₁₂＝fΔt a₂₁＝-fΔt a₂₂＝a₁₁

$l_{\mathrm{ij}}={\∫}_{\mathrm{\Ω}}\frac{\∂N_i}{\∂x}\frac{\∂N_j}{\∂x}\mathrm{d\Ω}=\frac{1}{4}\underset{e=1}{\overset{E}{\mathrm{\Σ}}}\frac{1}{{\mathrm{\Δ}}^{\left(e\right)}}\underset{\mathrm{\α}=1}{\overset{3}{\mathrm{\Σ}}}\underset{\mathrm{\λ}=1}{\overset{3}{\mathrm{\Σ}}}\left[(y_{\mathrm{\β}}^{\left(e\right)}-y_{\mathrm{\γ}}^{\left(e\right)})(y_{\mathrm{\μ}}^{\left(e\right)}-y_v^{\left(e\right)})\right]{\mathrm{\Δ}}_{\mathrm{\αi}}^{\left(e\right)}{\mathrm{\Δ}}_{\mathrm{\λj}}^{\left(e\right)}$

$m_{\mathrm{ij}}={\∫}_{\mathrm{\Ω}}\frac{\∂N_i}{\∂y}\frac{\∂N_j}{\∂y}\mathrm{d\Ω}=\frac{1}{4}\underset{e=1}{\overset{E}{\mathrm{\Σ}}}\frac{1}{{\mathrm{\Δ}}^{\left(e\right)}}\underset{\mathrm{\α}=1}{\overset{3}{\mathrm{\Σ}}}\underset{\mathrm{\λ}=1}{\overset{3}{\mathrm{\Σ}}}\left[(x_{\mathrm{\β}}^{\left(e\right)}-x_{\mathrm{\γ}}^{\left(e\right)})(x_{\mathrm{\μ}}^{\left(e\right)}-x_v^{\left(e\right)})\right]{\mathrm{\Δ}}_{\mathrm{\αi}}^{\left(e\right)}{\mathrm{\Δ}}_{\mathrm{\λj}}^{\left(e\right)}$

b_ij＝l_ij+m_ij

$c_{\mathrm{ij}}={\∫}_{\mathrm{\Ω}}{N}_i\frac{\∂N_j}{\∂x}\mathrm{d\Ω}=\frac{1}{6}\underset{e=1}{\overset{E}{\mathrm{\Σ}}}\underset{\mathrm{\α}=1}{\overset{3}{\mathrm{\Σ}}}\underset{\mathrm{\λ}=1}{\overset{3}{\mathrm{\Σ}}}(y_{\mathrm{\μ}}^{\left(e\right)}-y_v^{\left(e\right)}){\mathrm{\Δ}}_{\mathrm{\αi}}^{\left(e\right)}{\mathrm{\Δ}}_{\mathrm{\λj}}^{\left(e\right)}$

$d_{\mathrm{ij}}={\∫}_{\mathrm{\Ω}}{N}_i\frac{\∂N_j}{\∂x}\mathrm{d\Ω}=\frac{1}{6}\underset{e=1}{\overset{E}{\mathrm{\Σ}}}\underset{\mathrm{\α}=1}{\overset{3}{\mathrm{\Σ}}}\underset{\mathrm{\λ}=1}{\overset{3}{\mathrm{\Σ}}}(x_{\mathrm{\μ}}^{\left(e\right)}-x_v^{\left(e\right)}){\mathrm{\Δ}}_{\mathrm{\αi}}^{\left(e\right)}{\mathrm{\Δ}}_{\mathrm{\λj}}^{\left(e\right)}$

$f_{\mathrm{xi}}^{n+\frac{1}{2}}={\left(\frac{{\mathrm{\τ}}_{\mathrm{sx}}}{\mathrm{\ρH}}\right)}_i^{n+\frac{1}{2}}$

$f_{\mathrm{yi}}^{n+\frac{1}{2}}={\left(\frac{{\mathrm{\τ}}_{\mathrm{sy}}}{\mathrm{\ρH}}\right)}_i^{n+\frac{1}{2}}$

其中为p_i点的迎风单元(e_i)中三个结点的坐标，如图3，Δ(e)为迎风单元的面积。A_i为三角形网格上结点p_i的集中质量区域的面积，这个区域取为相邻三角形单元内的中线所围成的闭合区域，如图4。N_i为相应于结点p_i的插值函数。整体的插值函数N_i可以写为各个单元的分块表达式：

$N_i=\underset{e=1}{\overset{E}{\mathrm{\Σ}}}\underset{\mathrm{\α}=1}{\overset{3}{\mathrm{\Σ}}}{N}_{\mathrm{\α}}^{\left(e\right)}{\mathrm{\Δ}}_{\mathrm{\αi}}^{\left(e\right)}$

式中为Boole矩阵，

$N_{\mathrm{\α}}^{\left(e\right)}=\frac{1}{2{\mathrm{\Δ}}^{\left(e\right)}}\{(y_{\mathrm{\β}}^{\left(e\right)}-y_{\mathrm{\γ}}^{\left(e\right)})x+(x_{\mathrm{\γ}}^{\left(e\right)}-x_{\mathrm{\β}}^{\left(e\right)})y+x_{\mathrm{\β}}^{\left(e\right)}{y}_{\mathrm{\γ}}^{\left(e\right)}-x_{\mathrm{\γ}}^{\left(e\right)}{y}_{\mathrm{\β}}^{\left(e\right)}\}$

α，β，γ＝1，2，3进行循环置换。

令 $e_1=u_i^{n+\frac{1}{2}}-\mathrm{\ϵ}\frac{\mathrm{\Δt}}{A_i}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{b}_{\mathrm{ij}}{u}_j^{n+\frac{1}{2}}-\frac{\mathrm{g\Δt}}{A_i}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{c}_{\mathrm{ij}}{\mathrm{\ξ}}_j^{n+\frac{1}{2}}+f_{\mathrm{xi}}^{n+\frac{1}{2}}\mathrm{\Δt}$

$e_2=v_i^{n+\frac{1}{2}}-\mathrm{\ϵ}\frac{\mathrm{\Δt}}{A_i}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{b}_{\mathrm{ij}}{v}_j^{n+\frac{1}{2}}-\frac{\mathrm{g\Δt}}{A_i}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{d}_{\mathrm{ij}}{\mathrm{\ξ}}_j^{n+\frac{1}{2}}+f_{\mathrm{xi}}^{n+\frac{1}{2}}\mathrm{\Δt}$

则解式(2-10)得：

$\left\{\begin{array}{c}{u}_i^{n+1}=(a_{22}{e}_1-a_{12}{e}_2)/(a_{11}{a}_{22}-a_{21}{a}_{12})\\ v_i^{n+1}=(a_{11}{e}_2-a_{21}{e}_1)/(a_{11}{a}_{22}-a_{21}{a}_{12})\\ C_i^{n+1}=(T_i^{n+\frac{1}{2}}-\frac{\mathrm{\Δt}}{A_i}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}(K_x{l}_{\mathrm{ij}}+K_y{m}_{\mathrm{ij}})C_i^{n+\frac{1}{2}})/\left[{(1+\frac{K_s}{\mathrm{\ρ}{c}_{p}H})}_i^{n+\frac{1}{2}}\mathrm{\Δt}\right]\end{array}\right.---(2-11)$

在前半分步求解的是双曲型方程，在特征线通过的边界点上用式(2-5)计算在特征线不能到达的边界点上要用其他的方式确定这些点的中间辅助值

方程(1-13)表示水体的每个体积元中质量守恒。在一整步中采用与u，v，T时间交错的有限体积(集中质量区域)守恒格式，则有

其中A_i为结点p_i的集中质量区域的面积，Γ_i为结点p_i的集中质量区域的边界。

式(2-12)右边的积分采用梯形公式计算，得

$H_i^{n+\frac{1}{2}}=H_i^{n-\frac{1}{2}}-\frac{\mathrm{\Δt}}{A_i}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}(p_{\mathrm{ij}}{u}_j^n+q_{\mathrm{ij}}{v}_j^n)H_j^{n-\frac{1}{2}}---(2-13)$

其中 $p_{\mathrm{ij}}=\underset{e=1}{\overset{E}{\mathrm{\Σ}}}\underset{\mathrm{\α}=1}{\overset{3}{\mathrm{\Σ}}}\underset{\mathrm{\α}=1}{\overset{3}{\mathrm{\Σ}}}{r}_{\mathrm{\α\λ}}^{\left(e\right)}{\mathrm{\Δ}}_{\mathrm{\αi}}^{\left(e\right)}{\mathrm{\Δ}}_{\mathrm{\αj}}^{\left(e\right)}$

$q_{\mathrm{ij}}=\underset{e=1}{\overset{E}{\mathrm{\Σ}}}\underset{\mathrm{\α}=1}{\overset{3}{\mathrm{\Σ}}}\underset{\mathrm{\α}=1}{\overset{3}{\mathrm{\Σ}}}{s}_{\mathrm{\α\λ}}^{\left(e\right)}{\mathrm{\Δ}}_{\mathrm{\αi}}^{\left(e\right)}{\mathrm{\Δ}}_{\mathrm{\αj}}^{\left(e\right)}$

至此，利用分步杂交法对深度平均的二维浅水环流数学模型已经全部离散化完毕。

应用本发明所提供的模式对某核电厂排放的核素进行了计算。数模模拟水域上游取自取水口上游约3km(排水口上游约4km)，下游至水库大坝，距排水口约9km。网格大小由取排水口处向外逐步扩大，即从3m到6m。

分别按连续排放和间断排放计算，计算选用的Kd值为直接采用受纳水体的水样和河床沉积物固体样进行实验所得实际测量结果。实验采用静态批式法，该方法的实验条件比较简单，可同时进行大量的平行实验，可控制和改变实验条件，且实验周期较短。根据核电厂排放源项中各核素所占比例、液态照射途径所致剂量的贡献以及便与文献资料进行对比，选取了三种核素¹⁰⁶Ru、¹³⁷Cs、⁹⁰Sr进行实验。采用放射性液态流出物受纳水体的水样和河床沉积物固体样进行实验，固液比为1∶10g/ml。其操作过程是先将样品进行预处理和分选，再将其与含有某种核素的水溶液(已知浓度)在振荡下长时间混合，直到水溶液中核素浓度稳定为止，然后进行离心使固液相分离，最后分别测定固相和液相中的核素浓度。通过计算固、液相中的核素浓度之比即可得出分配系数K_d。

计算条件和排水口浓度如下表：

表1计算条件和排水口浓度

与IAEA safety report series No.19计算的结果进行对比，结果见下表：

表2两种计算结果对比

可以看出：这两种模式所计算得到的液相浓度和有效沉积量均差别不大，均在1个数量级内，差别最大的是¹⁰⁶Ru在排放口下游1km处的有效沉积量，比IAEA的结果仅大了34％。结果表明用此模式计算核素在水中的迁移扩散是可行的。

显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其同等技术的范围之内，则本发明也意图包含这些改动和变型在内。

Claims (1)

1.一种考虑泥沙吸附的核素输移扩散分析方法，其特征在于：该方法建立的平面二维浅水环流的基本方程如下：

沿水深平均的连续方程：

$\frac{\∂\mathrm{\ζ}}{\∂t}+\frac{\∂\left(Hu\right)}{\∂x}+\frac{\∂\left(Hv\right)}{\∂y}=0$

沿水深平均的动量方程：

$\begin{array}{c}\frac{\∂u}{\∂t}+u\frac{\∂u}{\∂x}+v\frac{\∂u}{\∂y}=-g\frac{\∂\mathrm{\ζ}}{\∂x}-fv-\frac{gu}{C^{2}H}\sqrt{u^2+v^2}+\frac{{\mathrm{\τ}}_{sx}}{\mathrm{\ρ}H}+\\ \frac{1}{H}\frac{\∂}{\∂x}\left({\mathrm{H\ϵ}}_x\frac{\∂u}{\∂x}\right)+\frac{1}{H}\frac{\∂}{\∂y}\left({\mathrm{H\ϵ}}_y\frac{\∂u}{\∂y}\right)\\ \frac{\∂v}{\∂t}+u\frac{\∂v}{\∂x}+v\frac{\∂v}{\∂y}=-g\frac{\∂\mathrm{\ζ}}{\∂y}+fu-\frac{gv}{C^{2}H}\sqrt{u^2+v^2}+\frac{{\mathrm{\τ}}_{sy}}{\mathrm{\ρ}H}+\\ \frac{1}{H}\frac{\∂}{\∂x}\left({\mathrm{H\ϵ}}_x\frac{\∂v}{\∂x}\right)+\frac{1}{H}\frac{\∂}{\∂y}\left({\mathrm{H\ϵ}}_y\frac{\∂v}{\∂y}\right)\end{array}$

浓度方程：

$\frac{\∂C}{\∂t}+u\frac{\∂C}{\∂x}+v\frac{\∂C}{\∂y}=\frac{1}{H}\frac{\∂}{\∂x}\left({\mathrm{HD}}_x\frac{\∂C}{\∂x}\right)+\frac{1}{H}\frac{\∂}{\∂y}\left({\mathrm{HD}}_y\frac{\∂C}{\∂y}\right)-\mathrm{\λ}C+{\mathrm{SK}}_d+S_{in}$

在以上方程中，

ζ--相对基准面水位；

H--水深；

t--时间；

u，v--水深平均流速；

x，y--坐标；

g--重力加速度；

f--柯氏力系数；

C--谢才系数；

τ_sx，τ_sy--表面风应力在x，y方向的分力；

ρ--水的密度；

ε_x，ε_y--x，y方向的广义涡粘性系数；

D_x、D_y--x、y方向核素的综合扩散系数；

S--含沙量；

K_d--核素分配系数；

S_in--挟沙力；

λ--衰变常数；

利用分步杂交法对上述深度平均的二维浅水环流方程进行离散化，采用实测的核素分配系数值求解浓度方程，所述的核素分配系数的实测方法是：采用放射性液态流出物受纳水体的水样和河床沉积物固体样作为样品进行实验，将样品与已知浓度的含有一种放射性核素的水溶液在振荡下长时间混合，直到水溶液中核素浓度稳定为止；然后进行离心使固液相分离；最后分别测定固相和液相中的核素浓度，通过计算固、液相中的核素浓度之比即可得出核素分配系数K_d。