CN104778356A - Value simulation method of convection-diffusion mass-transfer process - Google Patents

Abstract

The invention provides a value simulation method of a convection-diffusion mass-transfer process, and belongs to the technical field of underground water pollution program value simulation. According to the method, a flow net is used as a calculation unit for carrying out underground water solute transport value analysis, according to special properties of the flow net, the flow line and the equipotential line of the flow net are used for scattering a solution domain. In a flow net unit, two flow lines form a sealed curve surface called as a flow pipe. The flow pipe inner and outer substance exchange is generated through molecular diffusion under the concentration gradient effect, so that in the convection-diffusion mass-transfer process, the substance exchange between the flow pipes cannot be influenced by convection, and only relies on the molecular diffusion under the concentration gradient, and the difficult-to-measure parameter of the dispersion coefficient is not needed. The value simulation method of the convection-diffusion mass-transfer process provided by the invention has the advantages that the flow net is used as a calculation unit for the value simulation of the convection-diffusion mass-transfer process, and the problem that the water power dispersion coefficient which is difficult to be accurately determined is required in the existing convection-diffusion mass-transfer problem value solution method is solved.

Description

A kind of method for numerical simulation of convection-diffusion effect mass transport process

Technical field

The invention belongs to groundwater contamination problem numerical simulation technology field, relate to a kind of method for numerical simulation of convection-diffusion effect mass transport process.

Background technology

Economic fast development is along with day by day serious water pollution problems.The water safety of refuse landfill filter liquor seepage, seawater intrusion, nuclear waste and productive life waste water serious threat people.Pollutant is the main path of contaminant transportation with ground water movement under convection current with diffusion acting in conjunction.Groundwater Contamination Prediction needs to set up the mathematical model describing transport of landfill pollutants with control, and wherein the determination of dispersion parameter and seepage field is that Contaminants Transport model is set up and the key solved.But due to the complicacy of seepage field and transport of landfill pollutants, hydrodynamic dispersion difficult parameters is with Accurate Measurement.Hydrodynamic dispersion process is made up of molecular diffusion and mechanical dispersion, different from the intrinsic parameter of this measurable material of molecular diffusivity, the mechanical dispersion coefficient that indoor or field survey obtains is by various factors, there is obvious size effect, the accuracy of parameter and the popularity of engineer applied very limited.Researchist had once furtherd investigate affects hydrodynamic dispersion parameter Different factor, and especially size effect is on the impact of hydrodynamic dispersion parameter, but achievement not yet has the directive significance of ubiquity.Numerical simulation carries out the main method that Convection-diffusion Equations solves, and hydrodynamic dispersion parameter is the key parameter of Convection-Diffusion Process numerical solution, in hydrodynamic dispersion difficult parameters with under the prerequisite of Accurate Determining, precision and the engineer applied of Convection-diffusion Equations numerical solution result are greatly limited.

Summary of the invention

In view of this, the object of the present invention is to provide a kind of method for numerical simulation of convection-diffusion effect mass transport process, the method is according to the special nature of drift net, take drift net as the numerical analysis that computing unit carries out groundwater solute transfer, no longer need dispersion coefficient when analyzing, thus solve hydrodynamic dispersion difficult parameters required in existing Convection-diffusion Equations method of value solving accurately to determine this difficult problem.

For achieving the above object, the invention provides following technical scheme:

A method for numerical simulation for convection-diffusion effect mass transport process, the method comprises the following steps:

S1: by following formula, solves stable state plane seepage field potential function,

$\frac{\∂\left(k_x\frac{\∂H}{\∂x}\right)}{\∂x}+\frac{\∂\left(k_y\frac{\∂H}{\∂y}\right)}{\∂y}=0,$

Wherein, k _x, k _yfor solving the infiltration coefficient in x, y direction in territory; H is head potential function;

S2: by following formula, solves stable state plane seepage field stream function,

ψ＝∫u _xdy-u _ydx

Wherein, u _x, u _yfor the seepage velocity in x, y direction;

S3: utilize method of interpolation, obtains the isoline of head potential function and stream function, and solves the intersecting point coordinate of isoline;

S4: utilize straight-line segment to be connected by adjacent isoline intersection point, will solve territory and be divided into quadrilateral mesh, described quadrilateral network is drift net cell;

S5: set up the mass-conservation equation in drift net cell;

To unit P, four unit adjacent around it are respectively W, S, E, N, and the Solute mass conservation equation in unit P is

$a_P{C}_P=a_W{C}_W+a_E{C}_E+a_S{C}_S+a_N{C}_N+a_P^0{C}_P^0+S_u,$

Wherein, C _p, C _w, C _e, C _s, C _nfor t+ Δ t controls the solution concentration in volume P, W, E, S, N; for t controls the pollutant levels in volume P; S _ufor the coefficient relevant with source sink term; a _p, a _w, a _e, a _s, a _n, for discrete equation coefficient;

S6: the coefficient a of Solute mass conservation equation in computing unit _w, a _e, a _s, a _n; During employing central difference, each coefficient is

$\left\{\begin{array}{c}{a}_W=D_w+\frac{E_w}{2};a_E=D_e-\frac{F_e}{2}\\ a_S=D_s+\frac{F_s}{2};a_N=D_n-\frac{F_n}{2}\end{array}\right.,$

Wherein, F be by control volume boundary interface to flow; D is the coefficient of diffusional resistance on interface; According to the character of drift net, n, s border be 0 to flow; And the flow of fluid all equals the flow q flowing tube fluid on w, e border; The expression formula of F and the D on each interface is:

F _w＝qC _w；F _e＝qC _e；F _n＝F _s＝0

$D_w=\frac{D_m{A}_w}{\mathrm{\δ}{x}_{\mathrm{WP}}};D_e=\frac{D_m{A}_e}{\mathrm{\δ}{x}_{\mathrm{EP}}};D_s=\frac{D_m{A}_s}{\mathrm{\δ}{y}_{\mathrm{SP}}};D_n=\frac{D_m{A}_n}{\mathrm{\δ}{y}_{\mathrm{NP}}}$

△F＝F _e-F _w+F _n-F _s＝q(C _e-C _w)

Wherein, q is the mass-flux in stream pipe; D _mfor coefficient of diffusion; A is interface perpendicular to the projected area on dispersal direction; C _w, C _efor controlling the solution concentration of w, e two on interface of volume P; δ x, δ y are the distance controlling volume node and adjacent control volume node;

S7: by following formula, the coefficient of Solute mass conservation equation in computing unit

$a_P^0=\frac{1}{2\mathrm{\Δt}}\underset{i=0}{\overset{3}{\mathrm{\Σ}}}(x_i{y}_{i+1}-x_{i+1}{y}_i),$

Wherein, x _i, y _ifor controlling the coordinate of volume P tetra-angle points; Δ t is the time step of numerical simulation;

S8: by following formula, the coefficient S of Solute mass conservation equation in computing unit _u, S _p,

$\stackrel{\‾}{S}\mathrm{\ΔV}=S_u+S_p{C}_P,$

Wherein, for source sink term is controlling the mean value in volume;

S9: by following formula, the coefficient a of Solute mass conservation equation in computing unit _p,

$a_P=a_W+a_E+a_S+a_N+a_P^0+\mathrm{\ΔF}-S_P;$

S10: Solute mass conservation equation is listed to each computing unit, and solve the linear equation line be made up of each element quality conservation equation.

Further, the u in described S2 _x, u _yvalue by solving hydraulic Head Distribution and infiltration coefficient obtains.

Further, in described S4, the density of drift net cell is controlled by the number of S3 medium value line.

Further, the C in described S6 _w, C _evalue utilize linear interpolation method to obtain.

Beneficial effect of the present invention is: the method for numerical simulation of a kind of convection-diffusion effect mass transport process provided by the invention, the method is according to the feature of drift net, be that the numerical analysis that computing unit carries out groundwater solute transfer has following advantage with drift net: 1. form a closed tubular surface by two streamlines and be called stream pipe, because streamline can not intersect, therefore the inside and outside fluid of stream pipe can not pass through tube wall, namely the fluid flowed in pipe does not have fluid system to exchange with the external world, therefore, the flowing of flowing in pipe is the One-Dimensional flows along grain direction.In solute transfer process, the material exchange between stream pipe by the impact no longer by convection current, and only relies on the molecular diffusion under concentration gradient.Therefore no longer dispersion coefficient is needed when analyzing.2. constant flow in stream pipe, therefore the position drift net unit of high flow rate is closeer, and low flow velocity place drift net unit is dredged.High flow rate place be all generally solute transfer when analyzing discrete unit need the position of encryption, and the loose unit in low flow velocity place also can save computational resource, namely drift net cell can preferably self-adaptation stratum and seepage field heterogeneous body to the requirement of unit density.

Accompanying drawing explanation

In order to make the object, technical solutions and advantages of the present invention clearly, below in conjunction with accompanying drawing, the present invention is described in further detail, wherein:

Fig. 1 is the process flow diagram of the method for the invention;

Fig. 2 is drift net cell schematics;

Fig. 3 is certain seepage field of utilizing this method to carry out stress and strain model based on drift net unit;

Fig. 4 is the contrast of numerical simulation result and the laboratory test results utilizing this method to carry out.

Embodiment

Below in conjunction with accompanying drawing, the preferred embodiments of the present invention are described in detail.

For hydrodynamic dispersion difficult parameters required in Convection-diffusion Equations numerical solution process accurately to determine this difficult problem, it is that computing unit carries out the method for Computer Numerical Simulation to convection-diffusion effect mass transport process that the present invention to propose under a kind of complex geological condition with drift net.The method, according to the special nature of drift net, utilizes the streamline of drift net and equipotential line to carry out discrete to solving territory.In drift net unit, form a closed tubular surface by two streamlines and be called stream pipe.Stream is managed inside and outside mass exchange and is occurred by the molecular diffusion under concentration gradient effect, therefore in convection-diffusion effect mass transport process, material exchange between stream pipe is by the impact no longer by convection current, and only relying on the molecular diffusion under concentration gradient, this is difficult to the parameter that accurately records namely no longer to need dispersion coefficient.

As shown in Figure 1, the method comprises the following steps:

S1: solve stable state plane seepage field potential function.Potential function is obtained by the stable percolation field equation solved under certain boundary value condition, as shown in the formula:

$\frac{\∂\left(k_x\frac{\∂H}{\∂x}\right)}{\∂x}+\frac{\∂\left(k_y\frac{\∂H}{\∂y}\right)}{\∂y}=0,$

In formula, k _x, k _yfor solving the infiltration coefficient in x, y direction in territory; H is head potential function.

S2: solve stable state plane seepage field stream function.After obtaining hydraulic Head Distribution, stream function can be tried to achieve according to following formula:

ψ＝∫u _xdy-u _ydx，

In formula, u _x, u _yfor the seepage velocity in x, y direction, can by solving the hydraulic Head Distribution that obtains and infiltration coefficient is tried to achieve.

S3: utilize method of interpolation, obtains the isoline of head potential function and stream function, and tries to achieve the intersecting point coordinate of isoline.

S4: utilize straight-line segment to be connected by adjacent isoline intersection point, will solve territory and be divided into quadrilateral mesh, described quadrilateral network is drift net cell; The density of drift net cell is controlled by the number of S3 medium value line.

S5: set up the mass-conservation equation in drift net cell, to unit P, four unit adjacent around it are respectively W, S, E, N, and the Solute mass conservation equation so in unit P is:

$a_P{C}_P=a_W{C}_W+a_E{C}_E+a_S{C}_S+a_N{C}_N+a_P^0{C}_P^0+S_u,$

In formula, C _p, C _w, C _e, C _s, C _nfor t+ Δ t controls the solution concentration in volume P, W, E, S, N; for t controls the pollutant levels in volume P; S _ufor the coefficient relevant with source sink term; a _p, a _w, a _e, a _s, a _n, for discrete equation coefficient.

S6: the coefficient a of Solute mass conservation equation in computing unit _w, a _e, a _s, a _n, during employing central difference, each coefficient is:

$\left\{\begin{array}{c}{a}_W=D_w+\frac{E_w}{2};a_E=D_e-\frac{F_e}{2}\\ a_S=D_s+\frac{F_s}{2};a_N=D_n-\frac{F_n}{2}\end{array}\right.,$

In formula, F be by control volume boundary interface to flow; D is the coefficient of diffusional resistance on interface.According to the character of drift net, n, s border be 0 to flow; And the flow of fluid all equals the flow q flowing tube fluid on w, e border.The expression formula of F and the D on each interface is:

F _w＝qC _w；F _e＝qC _e；F _n＝F _s＝0

$D_w=\frac{D_m{A}_w}{\mathrm{\δ}{x}_{\mathrm{WP}}};D_e=\frac{D_m{A}_e}{\mathrm{\δ}{x}_{\mathrm{EP}}};D_s=\frac{D_m{A}_s}{\mathrm{\δ}{y}_{\mathrm{SP}}};D_n=\frac{D_m{A}_n}{\mathrm{\δ}{y}_{\mathrm{NP}}}$

△F＝F _e-F _w+F _n-F _s＝q(C _e-C _w)

In formula, q is the mass-flux in stream pipe; D _mfor coefficient of diffusion; A is interface perpendicular to the projected area on dispersal direction; C _w, C _efor controlling the solution concentration of w, e two on interface of volume P, linear interpolation method is utilized to obtain; δ x, δ y are the distance controlling volume node and adjacent control volume node.

S7: the coefficient of Solute mass conservation equation in computing unit

$a_P^0=\frac{1}{2\mathrm{\Δt}}\underset{i=0}{\overset{3}{\mathrm{\Σ}}}(x_i{y}_{i+1}-x_{i+1}{y}_i),$

X _i, y _ifor controlling the coordinate of volume P tetra-angle points; Δ t is the time step of numerical simulation.

S8: the coefficient S of Solute mass conservation equation in computing unit _u, S _p:

$\stackrel{\‾}{S}\mathrm{\ΔV}=S_u+S_p{C}_P$

In formula, for source sink term is controlling the mean value in volume.

S9: the coefficient a of Solute mass conservation equation in computing unit _p:

$a_P=a_W+a_E+a_S+a_N+a_P^0+\mathrm{\ΔF}-S_P.$

S10: Solute mass conservation equation is listed to each computing unit, and solve the linear equation line be made up of each element quality conservation equation.

The cell schematics that Fig. 2 this method utilizes drift net to divide.In seepage field, due to unevenness and the anisotropy on stratum, streamline may orthogonal also possibility oblique angle with equipotential line.Fig. 2 is the general type of equipotential line and streamline grid division local in drift net.Wherein control module is P, and 4 control modules be adjacent are W, S, E, N.Control module is 4 angle points of P, and namely the intersection point of equipotential line and streamline, is respectively ω _p, η _p, λ _p, ε _p, the coordinate of 4 angle points in rectangular coordinate system is respectively (x ₀, y ₀), (x ₁, y ₁), (x ₂, y ₂), (x ₃, y ₃).The node coordinate controlling volume P is:

$\left\{\begin{array}{c}{c}_{\mathrm{Px}}=\frac{\underset{i=0}{\overset{3}{\mathrm{\Σ}}}(x_i^P+x_{i+1}^P)(x_i^P{y}_{i+1}^P-x_{i+1}^P{y}_i^P)}{3\underset{i=0}{\overset{n-1}{\mathrm{\Σ}}}(x_i^P{y}_{i+1}^P-x_{i+1}^P{y}_i^P)}\\ c_{\mathrm{Py}}=\frac{\underset{i=0}{\overset{3}{\mathrm{\Σ}}}(y_i^P+y_{i+1}^P)(x_i^P{y}_{i+1}^P-x_{i+1}^P{y}_i^P)}{3\underset{i=0}{\overset{n-1}{\mathrm{\Σ}}}(x_i^P{y}_{i+1}^P-x_{i+1}^P{y}_i^P)}\end{array}\right.$

Other node coordinate controlling volume can obtain by same procedure, therefore can obtain the internodal distance of adjacent control volume, as:

$\mathrm{\δ}{x}_{\mathrm{WP}}=\sqrt{(c_{\mathrm{Px}}^2-c_{\mathrm{Wx}}^2)+(c_{\mathrm{Py}}^2-c_{\mathrm{Wy}}^2)}$

In formula, (C _wx, C _wy) for controlling the node coordinate of volume W.Control volume node control adjacent with other volume is internodal to be obtained apart from available above formula similar type.

For trying to achieve the interfacial area of vertical and dispersal direction, as A _w, first make straight line slope be respectively k ₁, k ₂, that is:

$k_1=\frac{c_{\mathrm{Py}}-c_{\mathrm{Wy}}}{c_{\mathrm{Px}}-c_{\mathrm{Wx}}};k_2=\frac{y_3-y_0}{x_3-x_0}$

Then straight line between angle be:

${\mathrm{\α}}_{\mathrm{WP}}=\arctan \frac{|k_1-k_2|}{1+k_1{k}_2}$

Therefore can obtain:

$A_w=\sin {\mathrm{\α}}_{\mathrm{WP}}\·L_{\stackrel{\‾}{{\mathrm{\ϵ}}_P{\mathrm{\ω}}_P}}=\frac{|k_1-k_2|\·\sqrt{(x_0^2-x_3^2)+(y_0^2-y_3^2)}}{\sqrt{1+k_1^2{k}_2^2+k_1^2+k_2^2}}$

All the other each interface projected areas can obtain by above formula similar type.

Fig. 3 is certain seepage field of utilizing this method to carry out stress and strain model based on drift net unit.Rectangle seepage field length is 30cm, is highly 15cm.Up-and-down boundary is water proof border, and left margin enters water boundaries for determining head, and right margin is for determine head water outlet border.Determining seepage field generation steady seepage under head difference effect.Utilize seepage field drift net unit to carry out discrete to seepage field, the limited bulk length of side after discrete is equal, gets Δ x=Δ y=5.0mm, and namely the width of every root stream pipe is 5.0mm.

Fig. 4 is the calculated examples of Fig. 3 institute representation model.Apply to determine head at right boundary to make one dimension steady seepage occurs in model, seepage velocity is 4.2 × 10 ^-4cm/s.Enter water boundaries inject relative concentration be the sodium chloride solution of 1.0, get time step Δ t=1.0s, utilize method described in the unit of drift net shown in Fig. 3 and this invention to simulate the migration process of sodium chloride in seepage field.The change of 3 measuring points monitoring concentration of sodium chloride solution is evenly set along solute transfer direction, and compares with indoor model test result, can find out method described in this invention and model test result very identical.

What finally illustrate is, above preferred embodiment is only in order to illustrate technical scheme of the present invention and unrestricted, although by above preferred embodiment to invention has been detailed description, but those skilled in the art are to be understood that, various change can be made to it in the form and details, and not depart from claims of the present invention limited range.

Claims (4)

1. a method for numerical simulation for convection-diffusion effect mass transport process, is characterized in that: the method comprises the following steps:

S1: by following formula, solves stable state plane seepage field potential function,

$\frac{\∂\left(k_x\frac{\∂H}{\∂x}\right)}{\∂x}+\frac{\∂\left(k_y\frac{\∂H}{\∂y}\right)}{\∂y}=0,$

Wherein, k _x, k _yfor solving the infiltration coefficient in x, y direction in territory; H is head potential function;

S2: by following formula, solves stable state plane seepage field stream function,

ψ＝∫u _xdy-u _ydx

Wherein, u _x, u _yfor the seepage velocity in x, y direction;

S3: utilize method of interpolation, obtains the isoline of head potential function and stream function, and solves the intersecting point coordinate of isoline;

S4: utilize straight-line segment to be connected by adjacent isoline intersection point, will solve territory and be divided into quadrilateral mesh, described quadrilateral network is drift net cell;

S5: set up the mass-conservation equation in drift net cell;

To unit P, four unit adjacent around it are respectively W, S, E, N, and the Solute mass conservation equation in unit P is

$a_P{C}_P=a_W{C}_W+a_E{C}_E+a_S{C}_S+a_N{C}_N+a_P^0{C}_P^0+S_u,$

Wherein, C _p, C _w, C _e, C _s, C _nfor t+ Δ t controls the solution concentration in volume P, W, E, S, N; for t controls the pollutant levels in volume P; S _ufor the coefficient relevant with source sink term; a _p, a _w, a _e, a _s, a _n, for discrete equation coefficient;

S6: the coefficient a of Solute mass conservation equation in computing unit _w, a _e, a _s, a _n; During employing central difference, each coefficient is

$\left\{\begin{array}{c}{a}_W=D_w+\frac{F_w}{2};a_E=D_e-\frac{F_e}{2}\\ a_S=D_s+\frac{F_s}{2};a_N=D_n-\frac{F_n}{2}\end{array}\right.,$

Wherein, F be by control volume boundary interface to flow; D is the coefficient of diffusional resistance on interface; N, s border be 0 to flow; On w, e border, the flow of fluid all equals the flow q flowing tube fluid; The expression formula of F and the D on each interface is:

F _w＝qC _w；F _e＝qC _e；F _n＝F _s＝0

$D_w=\frac{D_m{A}_w}{\mathrm{\δ}{x}_{\mathrm{WP}}};D_e=\frac{D_m{A}_e}{\mathrm{\δ}{x}_{\mathrm{EP}}};D_s=\frac{D_m{A}_s}{\mathrm{\δ}{y}_{\mathrm{SP}}};D_n=\frac{D_m{A}_n}{\mathrm{\δ}{y}_{\mathrm{NP}}}$

ΔF＝F _e-F _w+F _n-F _s＝q(C _e-C _w)

Wherein, q is the mass-flux in stream pipe; D _mfor coefficient of diffusion; A is interface perpendicular to the projected area on dispersal direction; C _w, C _efor controlling the solution concentration of w, e two on interface of volume P; δ x, δ y are the distance controlling volume node and adjacent control volume node;

S7: by following formula, the coefficient of Solute mass conservation equation in computing unit

$a_P^0=\frac{1}{2\mathrm{\Δt}}\underset{i=0}{\overset{3}{\mathrm{\Σ}}}(x_i{y}_{i+1}-x_{i+1}{y}_i),$

Wherein, x _i, y _ifor controlling the coordinate of volume P tetra-angle points; Δ t is the time step of numerical simulation;

S8: by following formula, the coefficient S of Solute mass conservation equation in computing unit _u, S _p,

$\stackrel{\‾}{S}\mathrm{\ΔV}=S_u+S_p{C}_P,$

Wherein, for source sink term is controlling the mean value in volume;

S9: by following formula, the coefficient a of Solute mass conservation equation in computing unit _p,

$a_P=a_W+a_E+a_S+a_N+a_P^0+\mathrm{\ΔF}-S_P;$

S10: Solute mass conservation equation is listed to each computing unit, and solve the linear equation line be made up of each element quality conservation equation.

2. the method for numerical simulation of a kind of convection-diffusion effect mass transport process according to claim 1, is characterized in that: the u in described S2 _x, u _yvalue by solving hydraulic Head Distribution and infiltration coefficient obtains.

3. the method for numerical simulation of a kind of convection-diffusion effect mass transport process according to claim 1, is characterized in that: in described S4, the density of drift net cell is controlled by the number of S3 medium value line.

4. the method for numerical simulation of a kind of convection-diffusion effect mass transport process according to claim 1, is characterized in that: the C in described S6 _w, C _evalue utilize linear interpolation method to obtain.