CN105760575A - Building method of Bohai sea spilled-oil transporting and extension value forecasting system - Google Patents

Abstract

The invention discloses a Bohai sea spilled-oil transporting and extension value forecasting system.The system can conduct forecasting on a transporting trajectory of an oil film, an extended range, the oil film thickness, the oil-passing-through area, the time for spilled-oil for reaching a shore or a sensitive area and the like after oil spilling occurs on a point source and a continuous source in Bohai sea area, and a forecasting result can be displayed on a computer screen visually and dynamically.Oil spill emergency action modes, types and number of vessels and situations of laying oil fences or using oil dispersant or the like are selected according to the forecasting result, oil spill damage is relieved or eliminated, and disaster prevention and reduction are achieved.Back tracking is conducted on the source of spilled oil in the Bohai sea area by combining with finger prints and satellite remote sensing, the source of the spilled oil on the sea is found out, and therefore a troublemaker is determined.

Description

The defeated shifting of Bohai Sea marine oil overflow, extension Numerical Prediction System method for building up

Technical field

The present invention relates to marine environment science, relate to marine oil overflow disperse and transport Fast Prediction technology, be specifically related to the method for building up of the defeated shifting of marine oil overflow of a kind of Bohai Sea, extension Numerical Prediction System.

Background technology

According to estimates, the current whole world nearly 10,000,000 tons of oil every year enter ocean, account for the 5% of whole world oil total output, incredible amount.Although the reason causing offshore oil pollution is various, but topmost reason be due to oil carrier hit a submerged reef, collision accident etc. enters the oil of ocean, and the oil depot blast of the exploitation of offshore oil field, oil well accident and oil plant etc. flows into the oil of ocean, this type of accident oil spilling is the most serious, and marine eco-environment destructive power is maximum, and living marine resources infringement is the most serious, endanger human health, preventing hindering Activities of Ocean, particularly at Nearshore Shallow Sea oil spilling, this impact is particularly acute.China Bohai Sea belongs to sanidal shallow sea area, in recent years owing to the development and utilization in the Bohai Sea is day by day increased, particularly the developing rapidly of oil field development, in Bohai Sea, Bo Nan, Bo Xi and Liaodong Bay have offshore oil platform, submarine transport oil pipeline spreads all over seabed especially, add marine transportation frequently, make platform oil spilling, submerged pipeline break and the accidents such as collision of hitting a submerged reef of oil carrier constantly occur.Oil spilling is once enter ocean, experiencing the physics of complexity, chemistry and bioprocess, the defeated shifting of oil spilling, diffusion are interacted by the comprehensive of many factors, the waters of its process, marine organisms and the marine eco-environment are had a direct impact by greasy dirt, intertidal organism group can be had a strong impact on again by oil spilling to land, drastically influence the sustainable development of Bohai Sea resource.Therefore; the research afloat defeated shifting of oil spilling, flooding mechanism; development oil spilling Predicting Technique quickly and precisely; exploitation Bohai Offshore marine oil overflow Numerical Prediction Models; the afloat change procedure of Exact Forecast oil spilling, to adopting an effective measure in time, alleviates marine oil overflow and causes environment and economic loss; seawater quality and littoral aquaculture to the protection Bohai Sea marine eco-environment and the Bohai Sea, the Implement of sustainable development strategy is significant.For forecast, the prediction Bohai Sea and the disaster of littoral sea thereof and take precautions against natural calamities, mitigation provides scientific theory foundation and research method.To marine eco-environment impact assessment, the formulation etc. of administration section's contingency response system has important practical and is worth.This software is mainly used in the fields such as anti-spill response action and the maritime affairs judicial expertise at spill response center.

The analysis spread for oil spilling and prediction are basis and the foundations of oil spilling process and decision-making.American-European countries started to begin to marine oil overflow is predicted from the sixties in last century.Wherein the extension formula of the three-stage theory of Fay is still widely used, but hypothesis theoretical for Fay highlights the oil extension on static sea, and ignore sea water internal turbulence and the impact of other factors (evaporation, emulsifying etc.), being applied to the modeling of actual bay and have certain limitation, approximation is too big.Some researcheres utilize finite difference method to solve convective-diffusive equation, oil spilling is considered as source item, establish oil spilling numerical model, by solving convective-diffusive equation, numerical simulation oil spilling disperse and transport, but the method be limited in that the impact that can not consider turbulent flow well, can not the detailed process of the defeated shifting of simulated oil and extension, particularly can not simulate the process at oil spilling initial stage, well physics and chemical process can not be combined, also can only consider the stream of single partial tide for trend, its degree of accuracy is poor, it is impossible to form Forecast Mode.The theory of current particle STOCHASTIC DIFFUSION is applied to the research of marine oil spill Forecast Mode by increasing researcher, set up marine oil spill forecasting model, but in these researchs, often only consider certain single process or defeated shifting is separated with diffusion process, the oil spilling forecasting model set up, thereafter report result and actual observation result are meet better in trend, on value, gap is too big, this is the particle owing to elaioleucite being used as size constancy of volume, does not consider chemistry, the bioprocesss such as evaporation, emulsifying, dissolving." oil spilling affects model application system (SIMAP and OiIMAP) " that the states such as the U.S., Canada, Australia are commonly using, its oil spilling affects the model part of numerical forecast, though the method that make use of Monte Carlo particles track, but it is inadequate to the oil spilling defeated shifting in ocean space, diffusion mechanism research, dynamic process is not suitable for the Bohai Sea, and will be applied to the Bohai Sea need to develop again.

The research of offshore spilled oil Forecast Mode is mainly full power oil film motion value pattern by Present Domestic, though embodying the process of marine instantaneous oil spilling preferably, but the deficiency of this pattern has been the assumption that oil film nothing in change and motor process tears phenomenon, this does not meet the practical situation at sea, furthermore, more difficult simulation enters the oil (this is extremely important for the marine eco-environment) of water body and marine oil condition of overflowing continuously, and this is the limitation of method；Big mutually responsible institute of National Bureau of Oceanography and the Belgian North Sea and Escaut mouth mathematical model administrative office cooperate, " elaioleucite " concept is utilized to establish the Three-dimension Numerical Model of simulation marine oil overflow behavior, avoid the pseudo-spreading effect that advection-diffusion conventional numeric method causes, but this pattern fails to solve the evaporation of oil, oil is extended by emulsifying and solute Physicochemical interaction process, defeated divide a word with a hyphen at the end of a line into combined influence, can not quantitatively give fuel-displaced residual quantity, the key parameter such as oil concentration in oil film thickness and water body, Problems existing remains the particle that elaioleucite has been used as oil mass constancy of volume.

Particle is represented certain oil mass, physics, the chemical processes such as the afloat defeated shifting of oil spilling, extension and evaporation, emulsifying, dissolving, binding kinetics equation, with following the trail of the movement of random particles and carrying out each process after numerical simulation marine oil overflow with the loss of random particles and the change of particle oil content (bulk density), multi-mechanism after research marine oil overflow, set up oil spill for Bohai Sea comprehensive forecasting pattern, this similar research outer with Present Domestic, similar technique Integrated comparative be leading in field, the Bohai Sea.

Summary of the invention

For solving the problems referred to above, the invention provides the method for building up of the defeated shifting of marine oil overflow of a kind of Bohai Sea, extension Numerical Prediction System.

For achieving the above object, the technical scheme that the present invention takes is:

The defeated shifting of Bohai Sea marine oil overflow, extension Numerical Prediction System method for building up, comprise the steps:

S1, structure include the three-dimensional flow field of trend and wind current；

S11, set up tidal current field model:

Setting up Bohai Sea tidal current field Forecast Mode based on FVCOM model, its governing equation includes continuity equation and the equation of momentum:

Continuity equation:

$\frac{\∂\mathrm{\ζ}}{\∂t}+\frac{\∂\mathrm{Du}}{\∂x}+\frac{\∂\mathrm{Dv}}{\∂y}+\frac{\∂\mathrm{\ω}}{\∂\mathrm{\σ}}=0$

The X-direction equation of momentum:

$\begin{array}{c}\frac{\∂\mathrm{uD}}{\∂t}+\frac{\∂u^{2}D}{\∂x}+\frac{\∂\mathrm{uvD}}{\∂y}+\frac{\∂\mathrm{u\ω}}{\∂\mathrm{\σ}}-\mathrm{fvD}\\ =-\mathrm{gD}\frac{\∂\mathrm{\ζ}}{\∂x}-\frac{\mathrm{gD}}{{\mathrm{\ρ}}_o}[\frac{\∂}{\∂x}\left(D\underset{\mathrm{\σ}}{\overset{0}{\∫}}\mathrm{\ρd}{\mathrm{\σ}}^{\′}\right)+\mathrm{\σ\ρ}\frac{\∂D}{\∂x}]+\frac{1}{D}\frac{\∂}{\∂\mathrm{\σ}}\left(K_m\frac{\∂u}{\∂\mathrm{\σ}}\right)+{\mathrm{DF}}_x\end{array}$

The Y-direction equation of momentum:

$\begin{array}{c}\frac{\∂\mathrm{vD}}{\∂t}+\frac{\∂\mathrm{uvD}}{\∂x}+\frac{\∂v^{2}D}{\∂y}+\frac{\∂\mathrm{v\ω}}{\∂\mathrm{\σ}}+\mathrm{fuD}\\ =-\mathrm{gD}\frac{\∂\mathrm{\ζ}}{\∂y}-\frac{\mathrm{gD}}{{\mathrm{\ρ}}_o}[\frac{\∂}{\∂y}\left(D\underset{\mathrm{\σ}}{\overset{0}{\∫}}\mathrm{\ρd}{\mathrm{\σ}}^{\′}\right)+\mathrm{\σ\ρ}\frac{\∂D}{\∂y}]+\frac{1}{D}\frac{\∂}{\∂\mathrm{\σ}}\left(K_m\frac{\∂v}{\∂\mathrm{\σ}}\right)+{\mathrm{DF}}_y\end{array}$

Wherein: g is acceleration of gravity；F=2 Ω sin φ is Ke Shi parameter；Ω is ground tarnsition velocity；φ is geographic latitude；ζ is the height of water level counted from still water；H is the depth of water counted from still water；D=h+ ζ；U, v, w are flow velocity, and p is pressure, and ρ is density of sea water；

S12, set up wind current model:

Utilize MM5 calculate for wind current provide by time wind field, utilize FVCOM to forecast wind current, MM5 is nested through three layers, it is achieved 60h computational fields wind field data by time forecast；By cressman interpolation technique, obtaining calculating 10m wind force vector W10 above the sea in grid element, when calculating wind current, pattern input wind field is the east of 10m place wind speed W10, north component above sea, changes into sea surface wind stress τ by following equation_s:

${\mathrm{\τ}}_s=\sqrt{{\mathrm{\τ}}_{sx}^2+{\mathrm{\τ}}_{sy}^2}$

Wherein:

ρ_aFor atmospheric density, C_DFor wind-stress coefficient, it determines that formula is:

θ is wind speed northeast component angle；

S13, trend and wind current prediction result are intercoupled, constitute the three-dimensional flow field including trend and wind current；.

The theoretical model that S2, the emulsifying built under the evaporation of oil, wave action and dissolving physical and chemical process combine with Monte Carlo random theory；

S3, structure Monte Carlo particles track three dimensional non-linear coupling dynamic model；

S4, structure Bohai Sea three-dimensional oil spilling random motion numerical forecast model；

The actual oil experimental data of throwing in S5, the existing oil spill for Bohai Sea data of utilization and sea, checking model；

S6, under GIS GIS-Geographic Information System, set up software of forecasting.

The method have the advantages that

Method based on Monte Carlo particles track, structure Nonlinear Dynamical Model, the impact of oil behavior is combined by emulsifying, evaporation and the dissolving of oil and Three Dimensional Tidal Forecast Mode and wind, wave, the prediction Bohai Sea defeated shifting of marine oil overflow, extension and home to return to, and ultimately form software of forecasting, sudden marine oil spill is forecast accurately and fast.This software of forecasting can to after sea, Bohai Sea point source and continuous source generation oil spilling, the defeated shifting track of oil film, spreading range, oil film thickness, sweep that sea is long-pending, oil spilling to land or arrives at time etc. of sensitizing range and is predicted, it was predicted that result can directly perceived, Dynamic Announce on the computer screen.Select the mode of spill response action, ship type, quantity according to predicting the outcome and lay the oil fence or use the situations such as detergent, thus alleviating or eliminating oil spilling harm, reaching purpose of preventing and reducing natural disasters.By the anti-tracking that Bohai Offshore oil spilling is traced to the source, combine with finger printing and satellite remote sensing, find out the source of marine oil overflow, so that it is determined that troublemaker.By oil film is swept the prediction that sea is amassed etc., living resources loss is assessed in the investigation in conjunction with living resources.The grid resolution 100m of this software × 100m, in trajectory predictions 10h, error is less than 2km.

Accompanying drawing explanation

Fig. 1 is the flow chart of the method for building up of the defeated shifting of marine oil overflow of the embodiment of the present invention Bohai Sea, extension Numerical Prediction System.

Fig. 2 is the heat balance schematic diagram of oil film.

Oil film drift absolute fix curve and predicted position curve when Fig. 3 is for throwing oily

Fig. 4 is that in the embodiment of the present invention, the shipwreck schematic diagram that predicts the outcome in 2004.7.8. day 7:30-23:30 is taken turns in gold Jiangxi 6.

Fig. 5 is gold Jiangxi 6 in the embodiment of the present invention " take turns the shipwreck schematic diagram that predicts the outcome in 2004.7.10 06:00-11 day, 23:00 day.

Detailed description of the invention

In order to make objects and advantages of the present invention clearly understand, below in conjunction with embodiment, the present invention is further elaborated.Should be appreciated that specific embodiment described herein is only in order to explain the present invention, is not intended to limit the present invention.

As it is shown in figure 1, embodiments provide the method for building up of the defeated shifting of marine oil overflow of a kind of Bohai Sea, extension Numerical Prediction System, comprise the steps:

S1, structure include the three-dimensional flow field of trend and wind current；

S11, set up tidal current field model:

Setting up Bohai Sea tidal current field Forecast Mode based on FVCOM model, its governing equation includes continuity equation and the equation of momentum:

Continuity equation:

$\frac{\∂\mathrm{\ζ}}{\∂t}+\frac{\∂\mathrm{Du}}{\∂x}+\frac{\∂\mathrm{Dv}}{\∂y}+\frac{\∂\mathrm{\ω}}{\∂\mathrm{\σ}}=0$

The X-direction equation of momentum:

$\begin{array}{c}\frac{\∂\mathrm{uD}}{\∂t}+\frac{\∂u^{2}D}{\∂x}+\frac{\∂\mathrm{uvD}}{\∂y}+\frac{\∂\mathrm{u\ω}}{\∂\mathrm{\σ}}-\mathrm{fvD}\\ =-\mathrm{gD}\frac{\∂\mathrm{\ζ}}{\∂x}-\frac{\mathrm{gD}}{{\mathrm{\ρ}}_o}[\frac{\∂}{\∂x}\left(D\underset{\mathrm{\σ}}{\overset{0}{\∫}}\mathrm{\ρd}{\mathrm{\σ}}^{\′}\right)+\mathrm{\σ\ρ}\frac{\∂D}{\∂x}]+\frac{1}{D}\frac{\∂}{\∂\mathrm{\σ}}\left(K_m\frac{\∂u}{\∂\mathrm{\σ}}\right)+{\mathrm{DF}}_x\end{array}$

The Y-direction equation of momentum:

$\begin{array}{c}\frac{\∂\mathrm{vD}}{\∂t}+\frac{\∂\mathrm{uvD}}{\∂x}+\frac{\∂v^{2}D}{\∂y}+\frac{\∂\mathrm{v\ω}}{\∂\mathrm{\σ}}+\mathrm{fuD}\\ =-\mathrm{gD}\frac{\∂\mathrm{\ζ}}{\∂y}-\frac{\mathrm{gD}}{{\mathrm{\ρ}}_o}[\frac{\∂}{\∂y}\left(D\underset{\mathrm{\σ}}{\overset{0}{\∫}}\mathrm{\ρd}{\mathrm{\σ}}^{\′}\right)+\mathrm{\σ\ρ}\frac{\∂D}{\∂y}]+\frac{1}{D}\frac{\∂}{\∂\mathrm{\σ}}\left(K_m\frac{\∂v}{\∂\mathrm{\σ}}\right)+{\mathrm{DF}}_y\end{array}$

Wherein: g is acceleration of gravity；F=2 Ω sin φ is Ke Shi parameter；Ω is ground tarnsition velocity；φ is geographic latitude；ζ is the height of water level counted from still water；H is the depth of water counted from still water；D=h+ ζ；U, v, w are flow velocity, and p is pressure, and ρ is density of sea water；

S12, set up wind current model:

Utilize MM5 calculate for wind current provide by time wind field, utilize FVCOM to forecast wind current, MM5 is nested through three layers, it is achieved 60h computational fields wind field data by time forecast；By cressman interpolation technique, obtaining calculating 10m wind force vector W10 above the sea in grid element, when calculating wind current, pattern input wind field is the east of 10m place wind speed W10, north component above sea, changes into sea surface wind stress τ by following equation_s:

${\mathrm{\τ}}_s=\sqrt{{\mathrm{\τ}}_{sx}^2+{\mathrm{\τ}}_{sy}^2}$

Wherein:

ρ_aFor atmospheric density, C_DFor wind-stress coefficient, it determines that formula is:

θ is wind speed northeast component angle；

S13, trend and wind current prediction result are intercoupled, constitute the three-dimensional flow field including trend and wind current；.

The theoretical model that S2, the emulsifying built under the evaporation of oil, wave action and dissolving physical and chemical process combine with Monte Carlo random theory；

1) evaporation

Fuel film vaporization is by the impact of the factors such as oil point, temperature and water temperature, oil spill area, wind speed, solar radiation and oil film thickness.Assuming that:

At oil film diffusion inside unrestricted (temperature higher than 0 DEG C and oil film thickness lower than 5-10cm time substantially such)；

Oil film is thoroughly mixed；

Oil ingredient dividing potential drop in an atmosphere is negligible compared with vapour pressure.

Evaporation rate can be expressed from the next:

$N_i^e=k_{ei}\·P_i^{SAT}/RT\·\frac{M_i}{{\mathrm{\ρ}}_i}\·X\·\[m^3/m^{2}s\]$

Wherein N is evaporation rate；Ke is the defeated shifting coefficient of material；PSAT is vapour pressure；R is gas constant；T is temperature；M is molecular weight；ρ is the density of oil ingredient；I is various oil ingredient.Kei is estimated by following formula:

$k_{ei}=k\·A_{oil}^{0.045}\·{\mathrm{Sc}}_i^{-2/3}\·U_w^{0.78}$

Wherein k is evaporation coefficient；Sci is the steam Schmidts number of component i.

2) emulsifying

1. emulsion oil-in-water process is formed

Oil includes dissolving, diffusion, precipitation etc. to the movement mechanism in water body.Diffusion is most important process in initial several weeks after oil spilling occurs.Diffusion is a kind of mechanical process, and oil film is torn into oil droplet by the turbulent energy of current, forms oil-in-water emulsifying.These emulsions can be stable by surfactant, it is prevented that oil droplet returns to oil film.Under severe weather conditions, topmost diffusion power is wave breaking, and topmost diffusion is the stretching, extension compression motion of oil film as stress under tranquil weather conditions.The loss of oil composition gauge in from oil slick's pervasion to water body is calculated:

D=D_a·D_b

Wherein Da is into the component of water body；Db does not have the component returned after being into water body:

$D_a=\frac{0.11{(1+U_w)}^2}{3600}$

$D_b=\frac{1}{1+50{\mathrm{\μ}}_{oil}\·k_s\·{\mathrm{\γ}}_{ow}}$

Wherein μ_oilFor oil viscosity；γ_owFor oil-water interfacial tension.

Oil droplet returns the speed of oil film；

$\frac{{\mathrm{dV}}_{oil}}{dt}=D_a\·(1-D_b)$

2. water-in-oil emulsification process is formed

The change of Water in oil rate can be represented by following formula equilibrium equation:

$\frac{{\mathrm{dy}}_w}{dt}=R_1-R_2$

The absorption rate of R1 and R2 respectively water and release rate, be given by.

$R_1=K_1\·\frac{{(1+U_w)}^2}{{\mathrm{\μ}}_{oil}}\·(y_w^{\max }-y_w)$

$R_2=K_2\·\frac{1}{As\·Wax\·{\mathrm{\μ}}_{oil}}\·y_w$

WhereinFor maximum moisture content；y_wFor actual moisture content；As is oil Asphalt Content (weight ratio)；Wax is paraffin content (weight ratio) in oil；K1, K2 be absorptance respectively, disengages coefficient.

3) dissolve

Dissolution rate following formula represents:

$\frac{{\mathrm{dV}}_{{\mathrm{ds}}_i}}{dt}={\mathrm{Ks}}_i\·C_i^{sai}\·X_{{\mathrm{mol}}_i}\·\frac{M_i}{{\mathrm{\ρ}}_i}\·A_{oil}$

WhereinDissolubility for component i；Molar fraction for component i；Mi is the molal weight of component i；Ksi, for dissolving mass tranfer coefficient, is estimated by following formula:

Ks_l=2.36 10^-6e_i

Wherein:

4) Heat Migration in Unsaturated

Vapour pressure and viscosity temperature influence, and observe the air and the water body that find that the temperature of usual oil film is higher than surrounding, Fig. 2 is the thermal balance schematic diagram of oil film.Wherein:

Diabatic process between 1=air and oil film；

Heat radiation process between 2=air and oil film；

3=solar radiation；

4=evaporative heat loss；

Heat Migration in Unsaturated between 5=oil film and water body；

The heat radiation distributed between 6=oil film and water body and accept.

1. the Heat Migration in Unsaturated between oil film and air

Heat Migration in Unsaturated between oil film and air can be expressed as:

$H_T^{oil-air}=A_{oil}\·k_H^{oil-air}\·(T_{air}-T_{oil})$

$k_H^{oil-air}=k_m\·{\mathrm{\ρ}}_a\·C_{pa}\·{\left(\frac{S_c}{P_r}\right)}_{air}^{0.67}$

Wherein Toil is oil film temperature；Tair is atmospheric temperature；ρ_aFor atmospheric density；Cpa is the thermal capacity of air；Pr is air Prandtl number:

$P_r=\frac{C_{pa}\·{\mathrm{\ρ}}_a}{0.0241(0.18055+0.003T_{air})}$

When evaporating negligible,Can calculate by simple following formula:

$k_H^{oil-air}=5.7+3.8U_w$

2. solar radiation

Many factors are depended in the solar radiation that oil film accepts, and most important of which is the water in oil spilling position, date, moment, cloud thickness and air, dust, ozone content.Solar radiation variations in one day can be assumed that as sine curve:

$H\left(t\right)=\left(\begin{array}{cc}{K}_t\&CenterDot;H_0^{\max }\&CenterDot;\sin \left(\mathrm{\&pi;}\frac{t-t^{sunrise}}{t^{sunset}-t^{sunrise}}\right),& t^{sunrise}<t<t^{sunset}\\ 0,& otherwise\end{array}\right.$

Wherein tsunrise is sunrise moment (late night to morning number of seconds)；Tsunset is sunset moment (late night to morning number of seconds)；Td is that day is long, it may be assumed that

t^sunset=t^sunrise+T_d

Td is calculated by following formula:

Wherein φ is latitude；ζ is sun angle of inclination (sun at noon time with the angle of equatorial plane)

Interstellar radiation for high noon:

Wherein Isc is solar constant (1.353W/m)；N is number of days in a year.ω s is hour angle of sunrise, is 0 during high noon, per hour equal to 15 (are just in the morning)；Kt is coefficient, Kt=0.75 during fine day, increases along with cloud thickness and reduces.Substantial portion of solar radiation is reflected when arriving ground, and therefore net heat input is:

(1-a)·H(t)

Wherein a is diffusion coefficient (albedo).

3. evaporative heat loss

Evaporation will cause oil film thermal loss:

$H^{vapor}=\underset{i}{\&Sigma;}{N}_i\&CenterDot;{\mathrm{\&Delta;H}}_{vi}\&CenterDot;\&lsqb;W/m^2\&rsqb;$

Wherein Δ H_viHeat of vaporization for component i.The total dynamic thermal balance of oil film has considered above-mentioned various factors:

$\begin{array}{c}\frac{{\mathrm{dT}}_{oil}}{dt}=\frac{1}{\mathrm{\&zeta;}\&CenterDot;C_p\&CenterDot;h}\&CenterDot;\&lsqb;(1-a)\&CenterDot;H+(l_{air}\&CenterDot;T_{air}^4+l_{water}\&CenterDot;T_{water}^4-2l_{oil}\&CenterDot;T_{oil}^4)\&rsqb;\\ +h_{ow}\&CenterDot;(T_{water}-T_{oil})+h_{oa}\&CenterDot;(T_{air}-T_{oil})-\&Sigma;N_i\&CenterDot;{\mathrm{\&Delta;H}}_{vi}\\ +(\frac{{\mathrm{dV}}_{water}}{dt}\&CenterDot;{\mathrm{\&zeta;}}_w\&CenterDot;C_{pw}+\frac{{\mathrm{dV}}_{oil}}{dt}\&CenterDot;{\mathrm{\&zeta;}}_{oil}\&CenterDot;C_{poil})\&CenterDot;(T_{water}-T_{oil})\&CenterDot;A_{oil}\end{array}$

4. the Heat Migration in Unsaturated between oil film and water body

Heat Migration in Unsaturated between oil film and air can be expressed as:

$H_H^{oil-air}=A_{oil}\&CenterDot;k_H^{oil-air}\&CenterDot;(T_{water}-T_{oil})$

$k_H^{oil-air}=0.332+r_w\&CenterDot;C_{pw}\&CenterDot;{\mathrm{Re}}^{-0.5}\&CenterDot;{\mathrm{Pr}}_w^{-2/3}$

Wherein Cpw is the thermal capacity of water.Prw is the Prandtl number of water:

${\mathrm{Pr}}_w=C_{pw}\&CenterDot;v_w\&CenterDot;{\mathrm{\&rho;}}_w\&CenterDot;\left(\frac{1}{0.330+0.000848(T_n-273.15)}\right)$

Re is characteristic Reynold's number:

$\mathrm{Re}=\frac{v_{rel}\&CenterDot;\sqrt{\frac{4A_{oil}}{\mathrm{\&pi;}}}}{{\mathrm{\&eta;}}_u}$

Wherein vrel is the kinematic coefficient of viscosity of oil film.

(5) reflect and accept radiation

Oil film is by loss and accepts long-wave radiation.Clean receiving amount is calculated by Stefan-Boltzman formula:

$H_{total}^{rad}=\mathrm{\&sigma;}\&CenterDot;(l_{air}\&CenterDot;T_{air}^4+l_{water}{T}_{water}^4-2l_{oil}\&CenterDot;T_{oil}^4)$

Wherein, σ is Stefan-Boltzman constant [5.72 108W/ (m²K)]；The radiance of lair, Iwater, Ioil respectively big gas and water and oil.

S3, structure Monte Carlo particles track three dimensional non-linear coupling dynamic model；

The afloat diffusion process of oil spilling is actually a disperse process of turbulent flow, and the randomness of turbulent dispersion, available DSMC describes rightly.When the problem of required solution is that certain event occurs

Probability, or during certain expectation of a random variable, the method that they can pass through certain experiment, obtain the probability that this event occurs or the meansigma methods of this random variable, and by they solutions as problem.

When solving diffusion problem, utilize the randomness of diffusing phenomenon, by given turbulence intensity, time scale and population, while giving random number, try to achieve the diffusion of particle.The random number that diffusion is used, has two kinds of methods to produce, i.e. the method for the method of uniform random number and normal random number.

1) assume that a, b, c are the uniform random number between (-0.5～0.5), then:

A=a (a²+b²+c²)^1/2；

B=b (a²+b²+c²)^1/2；

C=c (a²+b²+c²)^1/2

The average represented is zero, because having discrete 1/3 distribution.Dispersion variance σ 2 and diffusion coefficient it

Between relation be: k=σ 2/ (2 Δ t).Then x, y, z direction particle displacement lx, ly, lz with the equivalence of particle Particle diffusion are:

l_x=A (6 Δ t K_x)^1/2

l_y=B (6 Δ t K_y)^1/2

l_z=C (6 Δ t K_z)^1/2

2) assume that a, b, c are the normal random number between (0～1), a, b, c average is zero, takes (0,1) normal distribution, then:

l_x=a (2 Δ t K_x)^1/2

l_y=b (2 Δ t K_y)^1/2

l_z=c (2 Δ t K_z)^1/2

Oil concentration calculates, and N number of particle when t=0, if x=y=z=0, when t=t, is converted into concentration with entering particle number in grid.The concentration tried to achieve in this way, when number of particles is few time, concentration value will be at sixes and sevens, need to do suitable smoothing, for instance around 27 grid are averaged (during two dimension 9 grid of application), so can suppress to be distributed irregular extremum.

GS4, structure Bohai Sea three-dimensional oil spilling random motion numerical forecast model；

1) oil film trajectory predictions

On the basis of the environmental forces parameter of environmental forces model offer, adopt Euler--Lagrange method for tracing, carry out the prediction of oil film centrode.Oil film center excursD speed, depend on ocean surface wind speed and subsurface runoff, being the function of room and time, its value obtains by the speed interpolation on three mesh points of oil film central point place grid, and around the speed mesh point on mesh point, the speed interpolation of Order Triangular Elements central point gets.The component velocity on xy direction on each grid node in space in certain moment is:

$\left\{\begin{array}{c}{V}_x=V_{rx}+\mathrm{\&alpha;}{V}_{wind}\sin \&lsqb;(180+{\mathrm{\&theta;}}_0+\mathrm{\&theta;})\mathrm{\&pi;}/180\&rsqb;\\ V_y=V_{ry}+\mathrm{\&alpha;}{V}_{wind}cos\&lsqb;(180+{\mathrm{\&theta;}}_0+\mathrm{\&theta;})\mathrm{\&pi;}/180\&rsqb;\end{array}\right.$

Wherein V_rx、V_ryFor the xy durection component of mesh point upper epidermis flow velocity, it is tidal current field and the superposing of wind current；Forecast wind speed on Vwind mesh point, a is wind factor, θ₀For wind direction, θ is the drift drift angle that elaioleucite is subject to wind effect.

The value of θ is relevant with the size of wind speed, and formula is:

$\mathrm{\&theta;}=\left\{\begin{array}{cc}40-8\sqrt{V_{wind}}& 0\&le;V_{wind}\&le;25m/s\\ 0& V_{wind}\&GreaterEqual;25m/s\end{array}\right.$

Elaioleucite drift orbit computing formula is:

Wherein: S₀For initial time, S is oil film central point position, V_l(x (t), y (t), be t) that Lagrange follows the trail of speed, $V_l=\sqrt{V_x^2+V_y^2}.$

2) oil slick spread is defeated moves prediction

The diffusion process that shearing flow and turbulent flow cause belongs to random motion, available random walk about method realize simulation.Whole particle cloud cluster diffusion process in water body is caused due to the random motion of each particle.Can describe with following formula for water surface Stochastic Model:

r_a'=R (6k_aΔt)^1/2

Wherein: r_a' for α=(x, y, z) turbulence diffusion length on direction；R is uniform random number between [-1,1].k_aFor the turbulent diffusivity on α direction, Δ t is time step.

The drift of oil spilling is advection process, diffusion process, the coefficient result of wind.

The displacement within the Δ t period of the i-th particle is represented by:

x_i=u_iΔt+r_x′

y_i=v_iΔt+r_y′

Wherein: u_i=u_Stream+u_Wind；v_i=v_Stream+v_Wind；

r_x', r_y' it is the random displacement on x, y direction

u_Stream, u_Wind, v_Stream, v_Wind, all obtained by environmental dynamics model.

Owing to each particle represents certain oil mass, expanding area and the oil film thickness of oil spilling can be calculated according to the position at mark particle place and representative oil mass.

The actual oil experimental data of throwing in S5, the existing oil spill for Bohai Sea data of utilization and sea, checking model

During 17 days 4 May of calendar year 2001,55 points of Chinese Marine University and spill response center, Yantai are combined at off-lying sea, Yantai point (122 ° 14 ' 48 " E, 38 ° 00 ' 35 " N) place and have been carried out throwing oil test, and are tracked location.Oil film drift absolute fix curve and predicted position curve when Fig. 3 is for throwing oily, it is known that measured value is substantially identical with predictive value.

S6, under GIS GIS-Geographic Information System, set up software of forecasting；

Oil spilling prediction and warning visualization system includes three parts task: Functional mould design, and map datum processes and database establishment, and GIS visualization system secondary development designs.

1) Functional mould design

By Fortran language, tidal current field forecasting model, wind current prediction model and oil spilling model are programmed

And generate dynamic link library (DynamicLinkLibrary, DLL) file, the program module debugged is generated DLL.Dll file encapsulates all of functions of modules, when system development, by GIS platform, allow generalized information system part reserve corresponding open interface, the dll file of the Fortran language generation of each module is embedded, carries out secondary development in conjunction with ArcEngine assembly, some the order control on visualization interface is improved in exploitation as requested, in order to call DLL model to complete correlation computations, providing the function support of core for visualization system, result is shown on interface the most at last.

2) map datum processes and database establishment

Building storehouse according to the various data that the module design phase needs, the GIS map data of foundation can be divided into the attribute data of spatial data and its correspondence content.Spatial data is divided into again raster data and vector data, including various gradient maps and thematic map.

3) GIS visualization system secondary development design

GIS visualization system secondary development design is after keying in corresponding functional module, on the basis of GIS database system, adopt GIS-Geographic Information System tool software and visual programming technique to carry out the various predictions of integrated project development, finally realize the displaying of relevant information and provide useful information with aid decision.Wherein, system designs and develops selection VisualBasic.net2005, ArcEngine9.2 of language platform.

Originally be embodied as first with have fourth-order accuracy imperial lattice--storehouse tower method solves elaioleucite advection process, and the method is higher than the Euler's method precision of single order and second order, make the prediction of oil spilling more meet reality；By Monte Carlo particles track method, achieve dynamic process first in the Bohai Sea and combine the efflorescences such as evaporation, emulsifying, comprehensively solve by the motion of elaioleucite and change in volume and loss.

Wherein, the Runge-Kutta methods of fourth-order accuracy solves elaioleucite advection process and is:

After there is oil spilling, oil film issues life shifting movement in wind, subsurface runoff effect, and this is the oil film Lagrangian Drift Process under the environmental forces effect that the shearing stress of wind, subsurface runoff synthesize substantially.The traveling locus of coastal waters oil spilling can be described by Lagrangian particle tracing model, thus problem is converted to solving this Lagrangian model equation.Lagrangian particle tracing path is based on following vector equation:

$\frac{{\mathrm{dx}}_i}{dt}=v_a(x_i,t)+v_d(x_i,t)$

X in formula_iFor particle coordinate, v_aFor the advection speed of coordinate i, the value in its discrete time is determined by hydrodynamic model, v_dFor random velocity.

Runge-Kutta methods solves a kind of technology substantially indirectly using Taylor series, and each step needs to calculate functional value 4 times, and its truncated error is 0 (h⁵), this differential equation Runge-Kutta methods solves and can be expressed as:

$a_i=\mathrm{\&Delta;}t\&lsqb;v_a(x_i^n,t^n)+v_d(x_i^n,t^n)\&rsqb;$

$b_i=\mathrm{\&Delta;}t\{v_a\&lsqb;x_i^n+\frac{1}{2}{a}_i,t^{n+1/2}\&rsqb;+v_d\&lsqb;x_i^n+\frac{1}{2}{a}_i,t^{n+1/2}\&rsqb;\}$

$c_i=\mathrm{\&Delta;}t\{v_a\&lsqb;x_i^n+\frac{1}{2}{b}_i,t^{n+1/2}\&rsqb;+v_d\&lsqb;x_i^n+\frac{1}{2}{b}_i,t^{n+1/2}\&rsqb;\}$

$d_i=\mathrm{\&Delta;}t\&lsqb;v_a(x_i^n+c_i,t^{n+1})+v_d(x_i^n+c_i,t^{n+1})\&rsqb;$

$x_i^{n+1}\&cong;x_i^n+\frac{1}{6}(a_i+2b_i+2c_i+d_i)$

Above-mentioned Runge-Kutta methods needs to ask the velocity amplitude of intermediate time, and intermediate time speed interpolation method is tried to achieve, and quadravalence interpolation speed is:

$v_i^{n+1/2}\left(x\right)=\frac{5}{16}{v}_i^{n+1}\left(x\right)+\frac{5}{16}{v}_i^n\left(x\right)-\frac{5}{16}{v}_i^{n-1}\left(x\right)+\frac{1}{16}{v}_i^{n-2}\left(x\right)$

Formula has been used the flow speed value in n+1, n, n-1, n-2 moment.

Embodiment

Morning on July 8th, 2004, one S.Korean nationality merchant ship " SAEHANACAXY " wheel and a Chinese merchant ship " gold Jiangxi 6 " wheel 121 ° 23 ' 17 of waters in the middle part of Bohai Haixia "; 38 ° 21 ' 51 " place crashes, " gold Jiangxi 6 " wheel sinks immediately,, there is oil spilling spill and leakage immediately in total fuel oil 90 tons on ship.Utilizing this software, this oil spill accident has carried out spilled oil drift track and DIFFUSION PREDICTION (prediction figure is as shown in Figure 4, Figure 5), Fig. 1 is: " gold Jiangxi 6 " wheel shipwreck oil spilling prediction place: 121 ° 23 ' 17 ", 38 ° 21 ' 51 "；Predicted time: 2004.7.8. day 7:30-23:30 predicted time interval: 1h, it was predicted that wind speed, wind direction: northeaster 2m/s；Fig. 2 is " gold Jiangxi 6 " wheel shipwreck oil spilling prediction place: 121 ° 23 ' 17 "; 38 ° 21 ' 51 ", predicted time: 2004.7.10 06:00-11 day, predicted time day 23:00 interval: 3h, it was predicted that wind speed, wind direction: southeaster turns northeaster 8m/s, continuous oil spilling 0.1 ton hour

Owing to lacking after oil spilling the Monitoring Data data to oil film, it is difficult to the accuracy and precision judging to predict the outcome, but the description according to the operating personnel that rushes to the scene, trajectory predictions had doped the particular location of oil exactly on 8th, and model DIFFUSION PREDICTION result on the 9th at sea finds that with Field Force the place of oil spilling band coincide with the time；The result of prediction on the 12nd is also substantially identical with the position finding oil film.

The above is only the preferred embodiment of the present invention; it should be pointed out that, for those skilled in the art, under the premise without departing from the principles of the invention; can also making some improvements and modifications, these improvements and modifications also should be regarded as protection scope of the present invention.

Claims (1)

1. the defeated shifting of Bohai Sea marine oil overflow, extension Numerical Prediction System method for building up, it is characterised in that comprise the steps:

S1, structure include the three-dimensional flow field of trend and wind current；

S11, set up tidal current field model:

Setting up Bohai Sea tidal current field Forecast Mode based on FVCOM model, its governing equation includes continuity equation and the equation of momentum:

Continuity equation:

$\frac{\&part;\mathrm{\&zeta;}}{\&part;t}+\frac{\&part;Du}{\&part;x}+\frac{\&part;Dv}{\&part;y}+\frac{\&part;\mathrm{\&omega;}}{\&part;\mathrm{\&sigma;}}=0$

The X-direction equation of momentum:

$\begin{array}{c}\frac{\&part;uD}{\&part;t}+\frac{\&part;u^{2}D}{\&part;x}+\frac{\&part;uvD}{\&part;y}+\frac{\&part;u\mathrm{\&omega;}}{\&part;\mathrm{\&sigma;}}-fvD\\ =-gD\frac{\&part;\mathrm{\&zeta;}}{\&part;x}-\frac{gD}{{\mathrm{\&rho;}}_o}\&lsqb;\frac{\&part;}{\&part;x}\left(D\underset{\mathrm{\&sigma;}}{\overset{0}{\&Integral;}}{\mathrm{\&rho;d\&sigma;}}^{\&prime;}\right)+\mathrm{\&sigma;}\mathrm{\&rho;}\frac{\&part;D}{\&part;x}\&rsqb;+\frac{1}{D}\frac{\&part;}{\&part;\mathrm{\&sigma;}}\left(K_m\frac{\&part;u}{\&part;\mathrm{\&sigma;}}\right)+{\mathrm{DF}}_x\end{array}$

The Y-direction equation of momentum:

$\begin{array}{c}\frac{\&part;vD}{\&part;t}+\frac{\&part;uvD}{\&part;x}+\frac{\&part;v^{2}D}{\&part;y}+\frac{\&part;v\mathrm{\&omega;}}{\&part;\mathrm{\&sigma;}}+fuD\\ =-gD\frac{\&part;\mathrm{\&zeta;}}{\&part;y}-\frac{gD}{{\mathrm{\&rho;}}_o}\&lsqb;\frac{\&part;}{\&part;y}\left(D\underset{\mathrm{\&sigma;}}{\overset{0}{\&Integral;}}{\mathrm{\&rho;d\&sigma;}}^{\&prime;}\right)+\mathrm{\&sigma;}\mathrm{\&rho;}\frac{\&part;D}{\&part;y}\&rsqb;+\frac{1}{D}\frac{\&part;}{\&part;\mathrm{\&sigma;}}\left(K_m\frac{\&part;v}{\&part;\mathrm{\&sigma;}}\right)+{\mathrm{DF}}_y\end{array}$

Wherein: g is acceleration of gravity；F=2 Ω sin φ is Ke Shi parameter；Ω is ground tarnsition velocity；φ is geographic latitude；ζ is the height of water level counted from still water；H is the depth of water counted from still water；D=h+ ζ；U, v, w are flow velocity, and p is pressure, and p is density of sea water；

S12, set up wind current model:

Utilize MM5 calculate for wind current provide by time wind field, utilize FVCOM to forecast wind current, MM5 is nested through three layers, it is achieved 60h computational fields wind field data by time forecast；By cressman interpolation technique, obtaining calculating 10m wind force vector W10 above the sea in grid element, when calculating wind current, pattern input wind field is the east of 10m place wind speed W10, north component above sea, changes into sea surface wind stress τ by following equation_s:

${\mathrm{\&tau;}}_s=\sqrt{{\mathrm{\&tau;}}_{sx}^2+{\mathrm{\&tau;}}_{sy}^2}$

Wherein:

ρ_aFor atmospheric density, C_DFor wind-stress coefficient, it determines that formula is:

θ is wind speed northeast component angle；

S13, trend and wind current prediction result are intercoupled, constitute the three-dimensional flow field including trend and wind current；.

The theoretical model that S2, the emulsifying built under the evaporation of oil, wave action and dissolving physical and chemical process combine with Monte Carlo random theory；

S3, structure Monte Carlo particles track three dimensional non-linear coupling dynamic model；

S4, structure Bohai Sea three-dimensional oil spilling random motion numerical forecast model；

The actual oil experimental data of throwing in S5, the existing oil spill for Bohai Sea data of utilization and sea, checking model；

S6, under GIS GIS-Geographic Information System, set up software of forecasting.