深水环境下海面溢油轨迹及其归宿的数值模拟方法

技术领域

本发明涉及一种深水环境下海面溢油轨迹及其归宿的数值模拟方法，属于海洋污染分析及防治技术领域。

背景技术

在海洋中，蕴含着丰富的油气资源、水产资源及药物资源等。海面溢油一旦发生，就会在海面上形成油膜。油膜不仅会阻碍大气与海水之间的氧气交换，减少阳光进入海水，同时也会附着在海洋生物的表面上，从而导致海洋生物因缺氧而死亡。海洋溢油污染具有影响海域范围广、持续时间长、对海洋生物和生态环境破坏大的特点。

关于海面溢油轨迹的数值模拟方法，主要有对流扩散方程法和油粒子模型法两种。由于前者在模拟长时间的溢油运动时，可能会产生数值扩散，因此其模拟的溢油轨迹与实际结果偏差较大。而油粒子模型在模拟溢油轨迹时分为两步，第一步是考虑在初期阶段，溢油在自身重力、表面张力及惯性力作用下的扩展过程；第二步是油粒子的漂移过程，这是模拟溢油轨迹的核心。在深水环境下，油粒子在表面波浪的作用下产生漂移运动。因此，水质点的运动是深水环境下模拟溢油轨迹的关键，而风场是影响水质点运动的关键因素之一。

在空气和水的环境下，油粒子除了做漂移运动外，同时还发生着蒸发、乳化、溶解及分散等过程，这些会引起油污物理化学性质的变化，如成份、密度、粘度。这属于溢油归宿的问题，通常采用实验测定法或者经验模型法来计算。

发明内容

发明目的：为了克服现有海面溢油轨迹及其归宿的数值模拟方法的不足，本发明提供了一种深水环境下海面溢油轨迹及其归宿的数值模拟方法。

技术方案：一种深水环境下海面溢油轨迹及其归宿的数值模拟方法，包括以下步骤：

1.设置海面、溢油及其环境参数，并计算溢油中所含油粒子总数；

溢油的参数为{t₀,t_end,T_o,ρ_o,μ_o,V_o,d_o,}，其中，t₀为模拟溢油的初始时间，t_end为模拟结束时间，T_o为溢油的温度，ρ_o为溢油的密度，μ_o为溢油的粘度，V_o为溢油的体积，d_o为油粒子的直径；海面参数为其中，U为海面上10m处的风速，为风区长度；水的参数{ρ_w,μ_w}，其中，ρ_w为水的密度，μ_w为水的粘度；

计算油粒子总数其中，floor(·)表示下取整，初始化油粒子的位置，$>X_i^0=0,Y_i^0=0,Z_i^0=0,$>i＝1,2,...,L。

2.根据海洋环境参数和海浪的色散关系，计算海面的频谱S(ω)，并对海面的频谱转换为波数谱S(k_m,k_n)；

$>S\left(\omega \right)=\frac{a{g}^2}{{\omega }^5}\exp [-\frac{5}{4}{\left(\frac{{\omega }_p}{\omega }\right)}^4]{\gamma }^{\beta },$>其中，$>\beta =\exp [-\frac{{(\omega -{\omega }_p)}^2}{2{\delta }^2{\omega }_p^2}],$>$>\delta =\left(\begin{array}{cc}0.07& {\omega }_p\le \omega \\ 0.09& {\omega }_p>\omega \end{array},,\right)$>$>{\omega }_p=7\pi {\left(\frac{g^2\tilde{F}}{U^3}\right)}^{-0.33},$>ω_p为谱峰值角频率，ω为海浪的角频率，γ为峰升因子，δ为峰形参量，g为引力常数，α为尺度系数；

海面波数进行网格离散化处理(k_m,k_n)，k_m为海面在x方向上的波数，m＝0,1,2,...,M-1，即在x方向上海面波数离散化M个点；k_n为海面在y方向上的波数，n＝0,1,2,...,N-1，即在y方向上海面波数离散化N个点，其中，L₁为海面在x方向的长度，L₂为海面在y方向的长度；

计算每个网格点上的海浪波数$>k_{\mathrm{mn}}=\sqrt{k_m^2+k_n^2}$>和角频率$>{\omega }_{\mathrm{mn}}=\sqrt{g{k}_{\mathrm{mn}}},$>则海面的波数谱$>S(k_m,k_n)=\frac{\mathrm{gS}\left({\omega }_{\mathrm{mn}}\right)}{2\sqrt{g{k}_{\mathrm{mn}}}}.$>

3.根据Longuet-Higgins海浪模型，计算海面的铅直位移和水质点的速度，分别记为η和{u,v,w}；

根据Longuet-Higgins海浪模型，深水海面在z方向上的铅直位移$>\eta =\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}{a}_{\mathrm{mn}}\cos (k_{m}x+k_{n}y-{\omega }_{\mathrm{mn}}t+{\phi }_{\mathrm{mn}}),$>其中，$>a_{\mathrm{mn}}=\sqrt{2S(k_m,k_n)\Delta k_x{k}_y},$>t为时间，φ_mn为[-π,π]均匀分布的相位噪声，Δk_x为x方向上相邻波数之差，Δk_y为y方向上相邻波数之差；则深水海面上水质点在{x,y,z}方向上的速度{u,v,w}为，

$>\left(\begin{array}{c}u=\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}{a}_{\mathrm{mn}}{\omega }_{\mathrm{mn}}{e}^{k_{\mathrm{mn}}\eta }\cos {\theta }_{\mathrm{mn}}\cos (k_{m}x+k_{n}y-{\omega }_{\mathrm{mn}}t+{\phi }_{\mathrm{mn}})\\ v=\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}{a}_{\mathrm{mn}}{\omega }_{\mathrm{mn}{}^{}}{e}^{k_{\mathrm{mn}}\eta }\sin {\theta }_{\mathrm{mn}}\cos (k_{m}x+k_{n}y-{\omega }_{\mathrm{mn}}t+{\phi }_{\mathrm{mn}})\left(a\right),\\ w=\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}{a}_{\mathrm{mn}}{\omega }_{\mathrm{mn}}{e}^{k_{\mathrm{mn}}\eta }\sin (k_{m}x+k_{n}y-{\omega }_{\mathrm{mn}}t+{\phi }_{\mathrm{mn}})\end{array}\right)$>其中，$>{\theta }_{\mathrm{mn}}=\arctan \left(\frac{k_n}{k_m}\right).$>

4.根据拉格朗日法和漂移过程，计算溢油轨迹；

(1).令l＝1，m＝0；

(2).令式(a)中的$>x=X_l^{\left(m\right)},$>$>y=Y_l^{\left(m\right)},$>$>z=Z_l^{\left(m\right)},$>$>t=\tilde{t},$>求得第l个油粒子在时刻的速度为则第l个油粒子在时刻的位置为$>\left(\begin{array}{c}{X}_l^{(m+1)}=X_l^{\left(m\right)}+\tilde{u}\mathrm{\Delta t}+\stackrel{}{}[2R-1]\sqrt{6K_x\mathrm{\Delta t}}\\ Y_l^{(m+1)}=Y_l^{\left(m\right)}+\tilde{v}\mathrm{\Delta t}+[2R-1]\sqrt{6K_y\mathrm{\Delta t}}\\ Z_l^{(m+1)}=Z_l^{\left(m\right)}+\tilde{w}\mathrm{\Delta t}+\frac{g({\rho }_w-{\rho }_o)d_o^2}{18{\mu }_w}+[2R-1]\sqrt{6k_z\mathrm{\Delta t}}\end{array}\right),$>其中，R为[0,1]均匀分布的随机数，Δt为时间步长，K_x为溢油在x方向上的扩散系数，K_y为溢油在y方向上的扩散系数，K_z为溢油在z方向上的扩散系数，若l＜L，则l＝l+1，并重复该步骤；

(3).若则m＝m+1，重复步骤(2)，否则得到所有油粒子在t_end时刻的位置。

5.根据经验模型，计算溢油在归宿时刻的的蒸发量、分散量、溶解量、含水量、密度及粘度；

溢油扩展后的面积$>A_o=2270{\left(\frac{{\rho }_w-{\rho }_o}{{\rho }_o}{V}_o\right)}^{\frac{2}{3}}\sqrt{t^{\prime }}+40{\left(\frac{{\rho }_w-{\rho }_o}{{\rho }_o}\right)}^{\frac{1}{3}}{U}^{\frac{4}{3}}{t}^{\prime },$>其中，t′为扩展时间；计算溢油在t_end时刻的蒸发量$>F_{\mathrm{ev}}=\frac{T_o}{10.3T_G}\ln [1+0.0026\left(\frac{T_G{A}_o{U}^{0.78}}{T_o{V}_o}\right)(t_{\mathrm{end}}-t_0)\exp (6.3-\frac{10.3T_0}{T_o})],$>其中，T₀为F_ev＝0时溢油初始沸点的温度，T_G为溢油沸点曲线的梯度；分散量$>F_d=\frac{0.11V_o{(1+U)}^2}{1+50\tau \frac{V_o\sqrt{{\mu }_o}}{A_o}},$>其中，τ为油水界面张力；溶解量$>F_{\mathrm{dis}}=K_d{A}_o{S}_0{e}^{-\upsilon (t_{\mathrm{end}}-t_0)},$>其中，υ为衰减常数，K_d溶解常数，S₀纯水中油的溶解量；乳化过程后油的含水量$>B_w=\frac{4}{5}[1-\exp (-\frac{4.5}{8}\times {10}^{-5}{(1+U)}^2(t_{\mathrm{end}}-t_0))],$>密度$>\tilde{\rho }=(1-Y_w)[F_{\mathrm{ev}}(0.6{\rho }_o-0.34)+{\rho }_o]+Y_w{\rho }_w,$>粘度$>\tilde{\mu }={\mu }_o{10}^{4F_{\mathrm{ev}}}\exp \left[\frac{2.5B_w}{1-0.654B_w}\right].$>

有益效果：与现有技术相比，本发明所提供的深水环境下海面溢油轨迹及其归宿的数值模拟方法，通过JONSWAP经验海浪谱来体现风场对水质点运动的影响。采用JONSWAP经验海浪谱和Longuet-Higgins海浪模型，计算深水环境下海面的铅直位移和水质点的速度，具有较高的运算效率；根据拉格朗日法和经验模型，计算溢油轨迹及其归宿。

附图说明

图1为本发明实施例的流程图。

具体实施方式

下面结合具体实施例，进一步阐明本发明，应理解这些实施例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。

图1为本发明的流程图。如图1所示，深水环境下海面溢油轨迹及其归宿的数值模拟方法，包括以下步骤：

1.设置海面、溢油及其环境参数，并计算溢油中所含油粒子总数；

溢油的参数为{t₀,t_end,T_o,ρ_o,μ_o,V_o,d_o,}，其中，t₀为模拟溢油的初始时间，t_end为模拟结束时间，T_o为溢油的温度，ρ_o为溢油的密度，μ_o为溢油的粘度，V_o为溢油的体积，d_o为油粒子的直径；海面参数为其中，U为海面上10m处的风速，为风区长度；水的参数为{ρ_w,μ_w}，其中，ρ_w为水的密度，μ_w为水的粘度；

计算油粒子总数其中，floor(·)表示下取整，初始化油粒子的位置，$>X_i^0=0,Y_i^0=0,Z_i^0=0,$>i＝1,2,...,L。

2.根据海洋环境参数和海浪的色散关系，计算海面的频谱S(ω)，并对海面的频谱转换为波数谱S(k_m,k_n)；

$>S\left(\omega \right)=\frac{\alpha g^2}{{\omega }^5}\exp [-\frac{5}{4}{\left(\frac{{\omega }_p}{\omega }\right)}^4]{\gamma }^{\beta },$>其中，$>\beta =\exp [-\frac{{(\omega -{\omega }_p)}^2}{2{\delta }^2{\omega }_p^2}],$>$>\delta =\left(\begin{array}{cc}0.07& {\omega }_p\le \omega \\ 0.09& {\omega }_p>\omega \end{array}\right),$>$>{\omega }_p=7\pi {\left(\frac{g^2\tilde{F}}{U^3}\right)}^{-0.33},$>ω_p为谱峰值角频率，ω为海浪的角频率，γ为峰升因子，δ为峰形参量，g为引力常数，α为尺度系数；

海面波数进行网格离散化处理(k_m,k_n)，k_m为海面在x方向上的波数，m＝0,1,2,...,M-1，即在x方向上海面波数离散化M个点；k_n为海面在y方向上的波数，n＝0,1,2,...,N-1，即在y方向上海面波数离散化N个点，其中，L₁为海面在x方向的长度，L₂为海面在y方向的长度；

计算每个网格点上的海浪波数$>k_{\mathrm{mn}}=\sqrt{k_m^2+k_n^2}$>和角频率$>{\omega }_{\mathrm{mn}}=\sqrt{g{k}_{\mathrm{mn}}},$>则海面的波数谱$>S(k_m,k_n)=\frac{\mathrm{gS}\left({\omega }_{\mathrm{mn}}\right)}{2\sqrt{g{k}_{\mathrm{mn}}}}.$>

3.根据Longuet-Higgins海浪模型，计算海面的铅直位移和水质点的速度，分别记为η和{u,v,w}；

根据Longuet-Higgins海浪模型，深水海面在z方向上的铅直位移$>\eta =\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}{a}_{\mathrm{mn}}\cos (k_{m}x+k_{n}y-{\omega }_{\mathrm{mn}}t+{\phi }_{\mathrm{mn}}),$>其中，$>a_{\mathrm{mn}}=\sqrt{2S(k_m,k_n)\Delta k_x\Delta k_y},$>t为时间，φ_mn为[-π,π]均匀分布的相位噪声，Δk_x为x方向上相邻波数之差，Δk_y为y方向上相邻波数之差；则深水海面上水质点在{x,y,z}方向上的速度{u,v,w}为，

$>\left(\begin{array}{c}u=\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}{a}_{\mathrm{mn}}{\omega }_{\mathrm{mn}}{e}^{k_{\mathrm{mn}}\eta }\cos {\theta }_{\mathrm{mn}}\cos (k_{m}x+k_{n}y-{\omega }_{\mathrm{mn}}t+{\phi }_{\mathrm{mn}})\\ v=\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}{a}_{\mathrm{mn}}{\omega }_{\mathrm{mn}}{e}^{k_{\mathrm{mn}}\eta }\sin {\theta }_{\mathrm{mn}}\cos (k_{m}x+k_{n}y-{\omega }_{\mathrm{mn}}t+{\phi }_{\mathrm{mn}})\\ w=\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}{a}_{\mathrm{mn}}{\omega }_{\mathrm{mn}}{e}^{k_{\mathrm{mn}}\eta }\sin (k_{m}x+k_{n}y-{\omega }_{\mathrm{mn}}t+{\phi }_{\mathrm{mn}})\end{array}\right)---\left(a\right),$>其中，$>{\theta }_{\mathrm{mn}}=\arctan \left(\frac{k_n}{k_m}\right).$>

4.根据拉格朗日法和漂移过程，计算溢油轨迹；

(1).令l＝1，m＝0；

(2).令式(a)中的$>x=X_l^{\left(m\right)},$>$>y=Y_l^{\left(m\right)},$>$>z=Z_l^{\left(m\right)},$>$>t=\tilde{t},$>求得第l个油粒子在时刻的速度为则第l个油粒子在时刻的位置为其中，R为[0,1]均匀分布的随机数，Δt为时间步长，K_x为溢油在x方向上的扩散系数，K_y为溢油在y方向上的扩散系数，K_z为溢油在z方向上的扩散系数，若l＜L，则l＝l+1，并重复该步骤；

(3).若则m＝m+1，重复步骤(2)，否则得到所有油粒子在t_end时刻的位置。

5.根据经验模型，计算溢油在归宿时刻的的蒸发量、分散量、溶解量、含水量、密度及粘度；

溢油扩展后的面积$>A_o=2270{\left(\frac{{\rho }_w-{\rho }_o}{{\rho }_o}{V}_o\right)}^{\frac{2}{3}}\sqrt{t^{\prime }}+40{\left(\frac{{\rho }_w-{\rho }_o}{{\rho }_o}{V}_o\right)}^{\frac{1}{3}}{U}^{\frac{4}{3}}{t}^{\prime },$>其中，t′为扩展时间；计算溢油在t_end时刻的蒸发量$>F_{\mathrm{ev}}=\frac{T_o}{10.3T_G}\ln [1+0.0026\left(\frac{T_G{A}_o{U}^{0.78}}{T_o{V}_o}\right)(t_{\mathrm{end}}-t_0)\exp (6.3-\frac{10.3T_0}{T_o})],$>其中，T₀为F_ev＝0时溢油初始沸点的温度，T_G为溢油沸点曲线的梯度；分散量$>F_d=\frac{0.11V_o{(1+U)}^2}{1+50\tau \frac{V_o\sqrt{{\mu }_o}}{A_o}},$>其中，τ为油水界面张力；溶解量$>F_{\mathrm{dis}}=K_d{A}_o{S}_0{e}^{-\upsilon (t_{\mathrm{end}}-t_0)},$>其中，υ为衰减常数，K_d溶解常数，S₀纯水中油的溶解量；乳化过程后油的含水量$>B_w=\frac{4}{5}[1-\exp (-\frac{4.5}{8}\times {10}^{-5}{(1+U)}^2(t_{\mathrm{end}}-t_0))],$>密度$>\tilde{\rho }=(1-Y_w)[F_{\mathrm{ev}}(0.6{\rho }_o-0.34)+{\rho }_o]+Y_w{\rho }_w,$>粘度$>\tilde{\mu }={\mu }_o{10}^{4F_{\mathrm{ev}}}\exp \left[\frac{2.5B_w}{1-0.654B_w}\right].$>