CN102332040B - Three-dimensional numerical simulation method for influence of flexible net on water flow - Google Patents

Abstract

The invention discloses a three-dimensional numerical simulation method for influence of a flexible net on water flow. The method comprises the deformation calculation of the net and the calculation of a flow field at the periphery of the net. The calculation of the flow field at the periphery of the net comprises the following steps of: determining a control equation, calculating a porous medium coefficient, dividing boundary conditions and grids, and calculating numerical values. In the method, the deformation of the net in the water flow is simulated by a concentrated mass method, and a three-dimensional flow field at the periphery of the net is simulated by a Realizable k-epsilon turbulence model and a porous medium model; by the method, the deformation of the net under the action of the water flow and the three-dimensional flow field at the periphery of the net can be acquired at the same time, and the problems that the conventional calculation and measurement of the flow field at the periphery of the net is high in cost and time consumption and a three-dimensional integral flow field cannot be provided are solved; and physical experiments prove that a flow field result acquired through numerical simulation accords with an actual result (the maximal relative error is less than 5 percent), and the calculating time of the numerical simulation is less than 5 hours.

Description

A kind of Three-dimensional Numerical Simulation Method of influence of flexible net on water flow

Technical field

The invention belongs to aquacultural engineering, hydraulic engineering and field of ocean engineering, specially refer to the Three-dimensional Numerical Simulation Method of influence of flexible net on water flow in a kind of flow action process.

Background technology

In the deep-water net cage culture process, water body is mobile closely related with fish growth in the case.The flow velocity size of water body has determined the exchange velocity of the inner water body of net cage in the case, and on water quality environment impact in the case significantly, so the stream effect that subtracts of etting receives much concern; In the etting Force Calculation, the stream effect that subtracts of etting receives publicity equally, many studies show that, the stressed of etting generally is directly proportional with the quadratic power of flow velocity, even there is less flow velocity difference, also may cause larger stressed error, therefore in the Force Calculation of etting, etting be can not ignore the attenuation effect of flow velocity; At present in the culture zone, utilize that expanded metals formula seawall subtracts stream, wave absorption also is widely used, utilizing etting to reduce the culture zone flow velocity is the groundwork principle that seawall subtracts stream work.

At present, utilize the physical experiments method as the main research means of carrying out etting and subtract Flow Behavior both at home and abroad.Because etting is extraordinary ocean engineering structure, it is different from marine structure things such as being generally used for harbour, offshore oilfield fully, characteristics with large deformation, flexibility, porous, the physical experiments simulated cost is high, use duration, simultaneously because the at present restriction of observation technology is difficult to provide near the whole flow field situation in the full basin of etting.

Summary of the invention

The problems referred to above that exist for solving prior art the objective of the invention is, provide a kind of simulate etting in current distortion and etting around the Three-dimensional Numerical Simulation Method of influence of flexible net on water flow of size and Orientation of flow velocity.

A kind of Three-dimensional Numerical Simulation Method of influence of flexible net on water flow may further comprise the steps:

A, etting distortion are calculated

For the simulation of etting distortion, adopt the method for lumped mass; Suppose that etting is made of the limited lumped mass point that connects without the quality spring, by calculating the displacement of lumped mass point under current and boundary condition effect, the shape that obtains netting; Described lumped mass point is located at two ends and the centre of each mesh order pin; In order to simplify calculating, suppose that the lumped mass point at order pin two ends is shaped as circle, its hydrodynamic force coefficient is constant in direction of motion; Because the order pin can be seen cylindrical bar as, so the lumped mass point of being located in the middle of the order pin should have the hydrodynamic force character of cylinder, its hydrodynamic force coefficient is directive, and is relevant with the relative velocity direction of motion of current;

Lumped mass point represents cylinder units on the order pin, for the analysis that acts on flow force on it, adopts the method for setting up local coordinate system (τ, η, ξ) here; Set up local coordinate (τ, η, ξ) at the Lecideales pin, the τ direction of principal axis is the axis direction of cylindrical bar, and the η axle is positioned at plane that τ axle and flow velocity V form and vertical with the τ axle;

Set up whole coordinate and be about to true origin and be fixed on the water surface, x direction of principal axis and current in the same way, the z axle along the depth of water vertically upward, the y axle is positioned at plane that x axle and z axle form and vertical with the z axle; Under the whole coordinate, lumped mass point i to 1,2 unit vector are e _I1And e _I2, τ, η, ξ axle unit vector is e _τ=(x _τ, y _τ, z _τ), e _η=(x _η, y _η, z _η), e _ξ=(x _ξ, y _ξ, z _ξ)

Local coordinate and whole coordinate conversion matrix are

$\left[C\right]=\left(\begin{array}{ccc}{x}_{\mathrm{\τ}}& x_{\mathrm{\η}}& x_{\mathrm{\ξ}}\\ y_{\mathrm{\τ}}& y_{\mathrm{\η}}& y_{\mathrm{\ξ}}\\ z_{\mathrm{\τ}}& z_{\mathrm{\η}}& z_{\mathrm{\ξ}}\end{array}\right)---\left(1\right)$

The suffered bourn acting force of order pin i under local coordinate system, as follows according to morison formula:

$F_{\mathrm{D\τ}}=-\frac{1}{2}\mathrm{\ρ}{C}_{\mathrm{D\τ}}S|\stackrel{\·}{\mathrm{\τ}}-e_{\mathrm{\τ}}\·V|(\stackrel{\·}{\mathrm{\τ}}-e_{\mathrm{\τ}}\·V)$

$F_{\mathrm{D\η}}=-\frac{1}{2}\mathrm{\ρ}{C}_{\mathrm{D\η}}S|\stackrel{\·}{\mathrm{\η}}-e_{\mathrm{\η}}\·V|(\stackrel{\·}{\mathrm{\η}}-e_{\mathrm{\η}}\·V)---\left(2\right)$

$F_{\mathrm{D\ξ}}=-\frac{1}{2}\mathrm{\ρ}{C}_{\mathrm{D\ξ}}S\left|\stackrel{\·}{\mathrm{\ξ}}\right|\left(\stackrel{\·}{\mathrm{\ξ}}\right)$

In the formula: ξ is the relative velocity on three directions under the local coordinate system;

Tension force between lumped mass point can be by formula

Obtain; C ₁, C ₂Be the construction material elasticity coefficient, for coefficient C ₁, C ₂Choose, can be different according to the material of netting twine, " fibrecord " match of writing with reference to Gerhard Klust and getting; F _{D τ}, F _{D η}, F _{D ξ}Be the suffered bourn acting force of order pin on three directions under the local coordinate system; C _{D τ}, C _{D η}, C _{D ξ}Hydrodynamic force coefficient for order pin on three directions under the local coordinate system; ρ is the density of fluid; S is the etting packing; D is sealike colour; ε is the netting twine elastic elongation, ε=(l-l ₀)/l ₀, l ₀Be the netting twine original length, l is length after the netting twine distortion;

Utilize transition matrix [C] power on the local coordinate can be transferred to power under the whole coordinate, so the external force on each lumped mass point all can obtain, the lumped mass motion of point differential equation under the global coordinate system is:

$\frac{d{\stackrel{\&CenterDot;}{x}}_i}{\mathrm{dt}}=f_2(x_i,y_i,z_i,{\stackrel{\&CenterDot;}{x}}_i,{\stackrel{\&CenterDot;}{y}}_i,{\stackrel{\&CenterDot;}{z}}_i,x_1,{\mathrm{<figure-callout\; id="1"\; label="y"\; filenames="HDA0000106009800000012.png,HDA0000106009800000031.png"\; state="\{\{state\}\}">y</figure-callout>}}_1,{\mathrm{<figure-callout\; id="1"\; label="z"\; filenames="HDA0000106009800000012.png,HDA0000106009800000031.png"\; state="\{\{state\}\}">z</figure-callout>}}_1,x_2,{\mathrm{<figure-callout\; id="2"\; label="y"\; filenames="HDA0000106009800000031.png"\; state="\{\{state\}\}">y</figure-callout>}}_2,{\mathrm{<figure-callout\; id="2"\; label="z"\; filenames="HDA0000106009800000031.png"\; state="\{\{state\}\}">z</figure-callout>}}_2;t)$

$\frac{d{\stackrel{\&CenterDot;}{y}}_i}{\mathrm{dt}}=g_2(x_i,y_i,z_i,{\stackrel{\&CenterDot;}{x}}_i,{\stackrel{\&CenterDot;}{y}}_i,{\stackrel{\&CenterDot;}{z}}_i,x_1,{\mathrm{<figure-callout\; id="1"\; label="y"\; filenames="HDA0000106009800000012.png,HDA0000106009800000031.png"\; state="\{\{state\}\}">y</figure-callout>}}_1,{\mathrm{<figure-callout\; id="1"\; label="z"\; filenames="HDA0000106009800000012.png,HDA0000106009800000031.png"\; state="\{\{state\}\}">z</figure-callout>}}_1,x_2,{\mathrm{<figure-callout\; id="2"\; label="y"\; filenames="HDA0000106009800000031.png"\; state="\{\{state\}\}">y</figure-callout>}}_2,{\mathrm{<figure-callout\; id="2"\; label="z"\; filenames="HDA0000106009800000031.png"\; state="\{\{state\}\}">z</figure-callout>}}_2;t)---\left(3\right)$

$\frac{d{\stackrel{\&CenterDot;}{z}}_i}{\mathrm{dt}}=h_2(x_i,y_i,z_i,{\stackrel{\&CenterDot;}{x}}_i,{\stackrel{\&CenterDot;}{y}}_i,{\stackrel{\&CenterDot;}{z}}_i,{\mathrm{<figure-callout\; id="1"\; label="x"\; filenames="HDA0000106009800000012.png,HDA0000106009800000031.png"\; state="\{\{state\}\}">x</figure-callout>}}_1,{\mathrm{<figure-callout\; id="1"\; label="y"\; filenames="HDA0000106009800000012.png,HDA0000106009800000031.png"\; state="\{\{state\}\}">y</figure-callout>}}_1,{\mathrm{<figure-callout\; id="1"\; label="z"\; filenames="HDA0000106009800000012.png,HDA0000106009800000031.png"\; state="\{\{state\}\}">z</figure-callout>}}_1,x_2,{\mathrm{<figure-callout\; id="2"\; label="y"\; filenames="HDA0000106009800000031.png"\; state="\{\{state\}\}">y</figure-callout>}}_2,{\mathrm{<figure-callout\; id="2"\; label="z"\; filenames="HDA0000106009800000031.png"\; state="\{\{state\}\}">z</figure-callout>}}_2;t)$

Each lumped mass motion of point equation of simultaneous forms the ordinary differential equation group; The original shape of given etting utilizes fourth-order Runge-Kutta method to find the solution the ordinary differential equation group, can solve the coordinate figure of each quality point on any t moment etting, thereby determines the shape of etting;

Flow Field Calculation around B, the etting

The numerical simulation of flow field characteristic around the etting is introduced porous media model simulation plane etting, adopt finite volume method to find the solution governing equation, the flow field characteristic around the etting is carried out numerical evaluation;

B1, governing equation

The governing equation of Realizable k-ε model is:

Continuity equation:

$\frac{\&PartialD;\mathrm{\&rho;}}{\&PartialD;t}+\frac{\&PartialD;\left(\mathrm{\&rho;}{u}_i\right)}{\&PartialD;x_i}=0---\left(4\right)$

The equation of momentum:

$\frac{\&PartialD;\left({\mathrm{\&rho;u}}_i\right)}{\&PartialD;t}+\frac{\&PartialD;\left(\mathrm{\&rho;}{u}_i{u}_j\right)}{\&PartialD;x_j}=-\frac{\&PartialD;p}{\&PartialD;x_i}+\frac{\&PartialD;}{\&PartialD;x_i}(\mathrm{\&mu;}\frac{\&PartialD;u_i}{\&PartialD;x_j}-\mathrm{\&rho;}\stackrel{\&OverBar;}{u_i^{\&prime;}{u}_j^{\&prime;}})+S_i---\left(5\right)$

In the formula, t is the time, and ρ is the density of fluid, and μ is the kinetic viscosity of fluid, u _i, u _jBe the hourly value of fluid velocity component, u _i', u _j' be the pulsating quantity of speed component, p is the time average of pressure, i and j equal respectively 1,2,3 denotation coordination components; Here upper line "-" expression is measured time average to physics;

S _iEquation of momentum source item, the fluid mass S outside the porous medium border _i=0, inner on the porous medium border,

$S_i=-(D_{\mathrm{ij}}\mathrm{\&mu;u}+{\mathrm{<figure-callout\; id="1"\; label="C\; ij"\; filenames="HDA0000106009800000012.png,HDA0000106009800000031.png"\; state="\{\{state\}\}">C</figure-callout>}}_{\mathrm{ij}}\frac{1}{2}\mathrm{\&rho;}|u\left|u\right),$

$D_{\mathrm{ij}}=\left(\begin{array}{ccc}{\mathrm{<figure-callout\; id="0"\; label="D\; n"\; filenames="HDA0000106009800000031.png"\; state="\{\{state\}\}">D</figure-callout>}}_n& 0& 0\\ 0& {\mathrm{<figure-callout\; id="0"\; label="D\; t"\; filenames="HDA0000106009800000031.png"\; state="\{\{state\}\}">D</figure-callout>}}_t& 0\\ 0& 0& D_t\end{array}\right),$ $C_{\mathrm{ij}}=\left(\begin{array}{ccc}{\mathrm{<figure-callout\; id="0"\; label="C\; n"\; filenames="HDA0000106009800000031.png"\; state="\{\{state\}\}">C</figure-callout>}}_n& 0& 0\\ 0& {\mathrm{<figure-callout\; id="0"\; label="C\; t"\; filenames="HDA0000106009800000031.png"\; state="\{\{state\}\}">C</figure-callout>}}_t& 0\\ 0& 0& C_t\end{array}\right)---\left(6\right)$

In the formula, D _Ij, C _IjBe the porous medium matrix of coefficients, u is fluid velocity, D _nRepresentation is to viscosity factor, D _tExpression tangential viscous resistance coefficient, C _nRepresentation is to inertia resistance coefficient, C _tRepresent tangential inertia resistance coefficient;

The k equation:

$\frac{\&PartialD;\left(\mathrm{\&rho;k}\right)}{\&PartialD;t}+\frac{\&PartialD;\left(\mathrm{\&rho;k}{u}_i\right)}{\&PartialD;x_i}=\frac{\&PartialD;}{\&PartialD;x_j}\left[(\mathrm{\&mu;}+\frac{{\mathrm{\&mu;}}_t}{{\mathrm{\&sigma;}}_k})\frac{\&PartialD;k}{\&PartialD;x_j}\right]+G_k+G_b-\mathrm{\&rho;\&epsiv;}-Y_M+S_k---\left(7\right)$

The ε equation:

$\frac{\&PartialD;\left(\mathrm{\&rho;\&epsiv;}\right)}{\&PartialD;t}+\frac{\&PartialD;\left(\mathrm{\&rho;\&epsiv;}{u}_i\right)}{\&PartialD;x_i}=\frac{\&PartialD;}{\&PartialD;x_j}\left[(\mathrm{\&mu;}+\frac{{\mathrm{\&mu;}}_t}{{\mathrm{\&sigma;}}_{\mathrm{\&epsiv;}}})\frac{\&PartialD;\mathrm{\&epsiv;}}{\&PartialD;x_j}\right]+\mathrm{\&rho;}{\mathrm{<figure-callout\; id="1"\; label="C"\; filenames="HDA0000106009800000012.png,HDA0000106009800000031.png"\; state="\{\{state\}\}">C</figure-callout>}}_1\mathrm{S\&epsiv;}-\mathrm{\&rho;}{\mathrm{<figure-callout\; id="2"\; label="C"\; filenames="HDA0000106009800000031.png"\; state="\{\{state\}\}">C</figure-callout>}}_2\frac{{\mathrm{\&epsiv;}}^2}{k+\sqrt{\mathrm{v\&epsiv;}}}+{\mathrm{<figure-callout\; id="1"\; label="C"\; filenames="HDA0000106009800000012.png,HDA0000106009800000031.png"\; state="\{\{state\}\}">C</figure-callout>}}_{1\mathrm{\&epsiv;}}\frac{\mathrm{\&epsiv;}}{k}{C}_{3\mathrm{\&epsiv;}}{G}_b+S_{\mathrm{\&epsiv;}}---\left(8\right)$

In the formula, x _i, x _jBe the component of different directions under the cartesian coordinate system, k is the Hydrodynamic turbulence energy, and ε is turbulence dissipation rate, and μ is the kinetic viscosity of fluid, and μ t is the turbulent viscosity of fluid;

${\mathrm{<figure-callout\; id="1"\; label="C"\; filenames="HDA0000106009800000012.png,HDA0000106009800000031.png"\; state="\{\{state\}\}">C</figure-callout>}}_1=\max [0.43,\frac{\mathrm{\&eta;}}{\mathrm{\&eta;}+5}],$ $\mathrm{\&eta;}=S\frac{k}{\mathrm{\&epsiv;}},$ $S=\sqrt{2S_{\mathrm{ij}}{S}_{\mathrm{ij}}},$ $S_{\mathrm{ij}}=\frac{1}{2}(\frac{\&PartialD;u_i}{\&PartialD;x_j}+\frac{\&PartialD;u_j}{\&PartialD;x_i})---\left(9\right)$

In the equation, G _kThe Turbulent Kinetic that is produced by average velocity gradient, G _bThe Turbulent Kinetic that is produced by buoyancy, Y _MThe overall dissipation rate that represents pulsation diffusion profile in the compressible turbulent flow, C ₂, C _{1 ε}, C _{3 ε}Be constant, S is strain modulus, σ _kAnd σ _εThe turbulent flow Prandtl number of k equation and ε equation, S _kAnd S _εIt is user-defined source item; Realizable k-ε model constants value C _{1 ε}=1.44, C ₂=1.9, σ _k=1.0, σ _ε=1.2;

B2, porous medium coefficient

Porous media model is to combine one in the porous medium zone rule of thumb to be assumed to be main resistance to flow, and this Resistance Value can be calculated by following formula:

F＝S _xλA (10)

In the formula, S _xBe the equation of momentum source item of x direction, λ is porous medium thickness, and A is the porous medium area; The F direction is opposite with water (flow) direction, gets its numerical values recited herein;

Equation (6) substitution equation (10) can be obtained water resistance F _dWith lift F _lExpression formula,

$F_d=(D_n\mathrm{\&mu;u}+{\mathrm{<figure-callout\; id="1"\; label="C\; n"\; filenames="HDA0000106009800000012.png,HDA0000106009800000031.png"\; state="\{\{state\}\}">C</figure-callout>}}_n\frac{1}{2}\mathrm{\&rho;}|u\left|u\right)\mathrm{\&lambda;A}---\left(11\right)$

$F_l=(D_t\mathrm{\&mu;u}+{\mathrm{<figure-callout\; id="1"\; label="C\; t"\; filenames="HDA0000106009800000012.png,HDA0000106009800000031.png"\; state="\{\{state\}\}">C</figure-callout>}}_t\frac{1}{2}\mathrm{\&rho;}|u\left|u\right)\mathrm{\&lambda;A}---\left(12\right)$

The porous medium coefficient need to be calculated by the data of model test; If a known angle of attack is 90 ° and plane etting the Resistance Value under different in flow rate vertical with water (flow) direction, adopt least square method to function F _d(u) carry out match with the flow velocity of model test and corresponding Resistance Value, obtain optimum porous medium coefficient D _n, C _nBecause the plane etting is vertical with water (flow) direction, so the lift value is 0, D in the numerical simulation _t, C _tCan ignore;

When the plane etting becomes certain attack angle alpha with current, must be with porous medium matrix of coefficients D _IjAnd C _IjCarry out calculating after the coordinate conversion:

$D_n^{\&prime;}=\frac{D_n+{\mathrm{<figure-callout\; id="2"\; label="D\; t"\; filenames="HDA0000106009800000031.png"\; state="\{\{state\}\}">D</figure-callout>}}_t}{2}+\frac{D_n-{\mathrm{<figure-callout\; id="2"\; label="D\; t"\; filenames="HDA0000106009800000031.png"\; state="\{\{state\}\}">D</figure-callout>}}_t}{2}\cos \left(2{\mathrm{\&alpha;}}^{\&prime;}\right)$ (13)

$C_n^{\&prime;}=\frac{C_n+{\mathrm{<figure-callout\; id="2"\; label="C\; t"\; filenames="HDA0000106009800000031.png"\; state="\{\{state\}\}">C</figure-callout>}}_t}{2}+\frac{C_n-{\mathrm{<figure-callout\; id="2"\; label="C\; t"\; filenames="HDA0000106009800000031.png"\; state="\{\{state\}\}">C</figure-callout>}}_t}{2}\cos \left(2{\mathrm{\&alpha;}}^{\&prime;}\right)$

$D_t^{\&prime;}=\frac{D_n-{\mathrm{<figure-callout\; id="2"\; label="D\; t"\; filenames="HDA0000106009800000031.png"\; state="\{\{state\}\}">D</figure-callout>}}_t}{2}\sin \left(2{\mathrm{\&alpha;}}^{\&prime;}\right)$ (14)

$C_t^{\&prime;}=\frac{C_n-{\mathrm{<figure-callout\; id="2"\; label="C\; t"\; filenames="HDA0000106009800000031.png"\; state="\{\{state\}\}">C</figure-callout>}}_t}{2}\sin \left(2{\mathrm{\&alpha;}}^{\&prime;}\right)$

Wherein, α '=90 °-α, α is the angle of attack;

Adopt least square method that model test value and the theoretical value of resistance and lift are carried out match, obtain making the porous medium coefficient of total error amount minimum;

B3, boundary condition and grid are divided

Take the angle of attack as 90 ° plane etting model as example, the global coordinate system of numerical model (x, y, z) true origin is taken at the projection of etting center on the water surface, and the x axle is water (flow) direction, and the y axle is direction vertical with the x axle on the surface level, and the z axle is vertical; The numerical model left end is the speed inlet boundary; Right-hand member is defined as the free discharge border; Tank sidewall, tank bottom surface and the water surface are defined as the Gu Bi border; Carry out grid divides in the Numerical modelling zone; Trellis-type is the destructuring tetrahedral grid, porous medium zone and near refined net;

B4, numerical evaluation

Adopt finite volume method (FVM) discrete, discrete scheme adopts the Second-order Up-wind form of convective motion, and preferably SIMPLEC algorithm of stability is selected in pressure-speed coupling, and inferior relaxation factor keeps default value.

Compared with prior art, the present invention has following beneficial effect:

1, the present invention adopts the distortion of method simulation etting in current of lumped mass, utilizes simultaneously Realizable k-ε turbulence model, simulates etting three-dimensional flow field on every side in conjunction with porous media model; The distortion of adopting this method can obtain simultaneously etting under the flow action obtains the three-dimensional flow field situation around the etting simultaneously; Find that according to our physical test checking the flow field result that numerical simulation obtains coincide (maximum relative error is in 5%) with actual, the computing time of numerical simulation is in 5 hours; The present invention preferably resolve cost during Flow Field Calculation is measured around the present etting high, with duration, can't provide the problem in three-dimension integrally flow field.

2 compare with the physical experiments means, and it is fast that the present invention has computing velocity, and accuracy is high, can provide under any flow velocity the three-dimension integrally flow field situation under the full basin around the etting simultaneously.

Description of drawings

10 in the total accompanying drawing of the present invention, wherein:

Fig. 1 is that model simplification figure is calculated in the etting distortion.

Fig. 2 is local coordinate synoptic diagram on the mesh member.

Fig. 3 is that etting is reduced to the porous media model synoptic diagram.

Fig. 4 is etting angle of attack synoptic diagram.

Fig. 5 is the vertical view of model and test measuring point arrangenent diagram.

Fig. 6 is the side view of Fig. 5.

Fig. 7 is the velocity profile that passes on the horizontal section at etting center.

Fig. 8 is the velocity profile that passes on the vertical cross section at etting center.

Fig. 9 is calculated value and the trial value comparison diagram of each point velocity u.

Figure 10 is mesh shape synoptic diagram.

Embodiment

Below in conjunction with accompanying drawing the present invention is described further.Shown in Fig. 1-9, the simulation test that adopts method of the present invention to carry out is as follows:

For the simulation of etting distortion, adopt the method for lumped mass; The etting distortion is calculated model simplification figure as shown in Figure 1; The local coordinate synoptic diagram as shown in Figure 2 on the mesh member; To the numerical simulation of flow field characteristic around the etting, introduce porous media model simulation plane etting (see figure 3); The porous medium coefficient need to be calculated by the data of model test, and (see figure 4) when the plane etting becomes certain attack angle alpha with current must carry out the porous medium matrix of coefficients to calculate after the coordinate conversion;

Model test is carried out in the wave current tank of Dalian University of Technology's seashore and offshore engineering National Key Laboratory.Test water flute length 69m, wide 2m, test depth of water 0.7m.The test etting be hexagon without tubercle plane nylon (PA) net, be fixed in tank central authorities by steel bar framework.The etting detail parameters is as shown in table 1, and model and measuring point are arranged as illustrated in Figures 5 and 6.

Numerical simulation and test findings compare: set up Three-dimension Numerical Model according to test model, with rectangle porous medium regional simulation plane etting.By adopting 5mm, numerical simulation is carried out to the flow rate attenuation model of plane etting in the porous medium zone of 10mm and three kinds of thickness of 50mm, finds that within the specific limits the thickness logarithm value result of calculation in porous medium zone has no significant effect.Selecting thickness is that numerical evaluation is carried out in the porous medium zone of 50mm.By the flow velocity of model test and the data of etting water resistance (seeing Table 2), calculate porous medium coefficient D _n=75990m ^-2, C _n=8.087m ^-1

Table 1 test etting detail parameters table

Table 2 model test flow velocity and corresponding etting Resistance Value

Each point velocity value (cm/s) of table 3 model test

The flow speed data of each measuring point of model test is as shown in table 3.Set up Three-dimension Numerical Model according to the test model size, the long 3m of numerical simulation section tank, the speed entrance is arranged on 0.5m place before the etting, and with the average u=15.90cm/s of measuring point A1 ~ A3 flow velocity in the model test as the entrance flow velocity, while given entrance tubulence energy k=4.3 * 10 ^-5m ²/ s ²With turbulence dissipation rate ε=2.81 * 10 ^-8m ²/ s ³Carry out numerical evaluation.

(such as Fig. 7 and 8) can find out intuitively by the velocity flow profile situation, flow rate attenuation zone is among a small circle arranged before the etting, the flow rate attenuation regional extent is larger behind the etting, the attenuation region width is identical with etting, along with the increase with the etting distance has the trend that narrows down, etting both sides and the corresponding increase of bottom flow rate of water flow.Each point velocity u calculated value and model test value contrast situation, as shown in Figure 9.Found out that by comparing result the calculated value of each point velocity and trial value are coincide better, maximum relative error is 2.5%.Flow velocity has obvious decay behind the etting, and flow rate attenuation trend is calculated consistent with test findings.The D in etting downstream point flow velocity is less than the C point in the model test, and mean difference is 3.1%, and this difference mainly is because along journey energy loss and the impact of streaming the generation of flow velocity size.Numerical simulation is consistent with trial value, illustrates that this numerical model can take into full account these factors, and analog result is realistic.

Claims (1)

1. the Three-dimensional Numerical Simulation Method of an influence of flexible net on water flow is characterized in that: may further comprise the steps:

A, etting distortion are calculated

For the simulation of etting distortion, adopt the method for lumped mass; Suppose that etting is made of the limited lumped mass point that connects without the quality spring, by calculating the displacement of lumped mass point under current and boundary condition effect, the shape that obtains netting; Described lumped mass point is located at two ends and the centre of each mesh order pin; In order to simplify calculating, suppose that the lumped mass point at order pin two ends is shaped as circle, its hydrodynamic force coefficient is constant in direction of motion; Because the order pin can be seen cylindrical bar as, so the lumped mass point of being located in the middle of the order pin should have the hydrodynamic force character of cylinder, its hydrodynamic force coefficient is directive, and is relevant with the relative velocity direction of motion of current;

Lumped mass point represents cylinder units on the order pin, for the analysis that acts on flow force on it, adopts the method for setting up local coordinate system (τ, η, ξ) here; Set up local coordinate (τ, η, ξ) at the Lecideales pin, the τ direction of principal axis is the axis direction of cylindrical bar, and the η axle is positioned at plane that τ axle and flow velocity V form and vertical with the τ axle;

Set up whole coordinate and be about to true origin and be fixed on the water surface, x direction of principal axis and current in the same way, the z axle along the depth of water vertically upward, the y axle is positioned at plane that x axle and z axle form and vertical with the z axle; Under the whole coordinate, lumped mass point i to 1,2 unit vector are e _I1And e _I2, τ, η, ξ axle unit vector is e _τ=(x _τ, y _τ, z _τ), e _η=(x _η, y _η, z _η), e _ξ=(x _ξ, y _ξ, z _ξ)

Local coordinate and whole coordinate conversion matrix are

$\left[C\right]=\left(\begin{array}{ccc}{x}_{\mathrm{\&tau;}}& x_{\mathrm{\&eta;}}& x_{\mathrm{\&xi;}}\\ y_{\mathrm{\&tau;}}& y_{\mathrm{\&eta;}}& y_{\mathrm{\&xi;}}\\ z_{\mathrm{\&tau;}}& z_{\mathrm{\&eta;}}& z_{\mathrm{\&xi;}}\end{array}\right)---\left(1\right)$

The suffered bourn acting force of order pin i under local coordinate system, as follows according to morison formula:

$F_{\mathrm{D\&tau;}}=-\frac{1}{2}\mathrm{\&rho;}{C}_{\mathrm{D\&tau;}}S|\stackrel{\&CenterDot;}{\mathrm{\&tau;}}-e_{\mathrm{\&tau;}}\&CenterDot;V|(\stackrel{\&CenterDot;}{\mathrm{\&tau;}}-e_{\mathrm{\&tau;}}\&CenterDot;V)$

$F_{\mathrm{D\&eta;}}=-\frac{1}{2}\mathrm{\&rho;}{C}_{\mathrm{D\&eta;}}S|\stackrel{\&CenterDot;}{\mathrm{\&eta;}}-e_{\mathrm{\&eta;}}\&CenterDot;V|(\stackrel{\&CenterDot;}{\mathrm{\&eta;}}-e_{\mathrm{\&eta;}}\&CenterDot;V)---\left(2\right)$

$F_{\mathrm{D\&xi;}}=-\frac{1}{2}\mathrm{\&rho;}{C}_{\mathrm{D\&xi;}}S\left|\stackrel{\&CenterDot;}{\mathrm{\&xi;}}\right|\left(\stackrel{\&CenterDot;}{\mathrm{\&xi;}}\right)$

In the formula: $\stackrel{\&CenterDot;}{\mathrm{\&tau;}}-e_{\mathrm{\&tau;}}\&CenterDot;V,\stackrel{\&CenterDot;}{\mathrm{\&eta;}}-e_{\mathrm{\&eta;}}\&CenterDot;V,$ ξ is the relative velocity on three directions under the local coordinate system; Tension force between lumped mass point can be by formula

Obtain; C ₁, C ₂Be the construction material elasticity coefficient, for coefficient C ₁, C ₂Choose, different according to the material of netting twine, " fibrecord " match of writing with reference to Gerhard Klust and getting; F _{D τ}, F _{D η}, F _{D ξ}Be the suffered bourn acting force of order pin on three directions under the local coordinate system; C _{D τ}, C _{D η}, C _{D ξ}Hydrodynamic force coefficient for order pin on three directions under the local coordinate system; ρ is the density of fluid; S is the etting packing; D is sealike colour; ε is the netting twine elastic elongation, ε=(l-l ₀)/l ₀, l ₀Be the netting twine original length, l is length after the netting twine distortion;

Utilize transition matrix [C] to transfer the power on the local coordinate under the whole coordinate power, so the external force on each lumped mass point all can obtain, the lumped mass motion of point differential equation under the global coordinate system is:

$\frac{d{\stackrel{\&CenterDot;}{x}}_i}{\mathrm{dt}}=f_2(x_i,y_i,z_i,{\stackrel{\&CenterDot;}{x}}_i,{\stackrel{\&CenterDot;}{y}}_i,{\stackrel{\&CenterDot;}{z}}_i,x_1,y_1,z_1,x_2,y_2,z_2;t)$

$\frac{d{\stackrel{\&CenterDot;}{y}}_i}{\mathrm{dt}}=g_2(x_i,y_i,z_i,{\stackrel{\&CenterDot;}{x}}_i,{\stackrel{\&CenterDot;}{y}}_i,{\stackrel{\&CenterDot;}{z}}_i,x_1,y_1,z_1,x_2,y_2,z_2;t)---\left(3\right)$

$\frac{d{\stackrel{\&CenterDot;}{z}}_i}{\mathrm{dt}}=h_2(x_i,y_i,z_i,{\stackrel{\&CenterDot;}{x}}_i,{\stackrel{\&CenterDot;}{y}}_i,{\stackrel{\&CenterDot;}{z}}_i,x_1,y_1,z_1,x_2,y_2,z_2;t)$

Each lumped mass motion of point equation of simultaneous forms the ordinary differential equation group; The original shape of given etting utilizes fourth-order Runge-Kutta method to find the solution the ordinary differential equation group, can solve the coordinate figure of each quality point on any t moment etting, thereby determines the shape of etting;

Flow Field Calculation around B, the etting

The numerical simulation of flow field characteristic around the etting is introduced porous media model simulation plane etting, adopt finite volume method to find the solution governing equation, the flow field characteristic around the etting is carried out numerical evaluation;

B1, governing equation

The governing equation of Realizable k-ε model is:

Continuity equation:

$\frac{\&PartialD;\mathrm{\&rho;}}{\&PartialD;t}+\frac{\&PartialD;\left(\mathrm{\&rho;}{u}_i\right)}{\&PartialD;x_i}=0---\left(4\right)$

The equation of momentum:

$\frac{\&PartialD;\left({\mathrm{\&rho;u}}_i\right)}{\&PartialD;t}+\frac{\&PartialD;\left(\mathrm{\&rho;}{u}_i{u}_j\right)}{\&PartialD;x_j}=-\frac{\&PartialD;p}{\&PartialD;x_i}+\frac{\&PartialD;}{\&PartialD;x_i}(\mathrm{\&mu;}\frac{\&PartialD;u_i}{\&PartialD;x_j}-\mathrm{\&rho;}\stackrel{\&OverBar;}{u_i^{\&prime;}{u}_j^{\&prime;}})+S_i---\left(5\right)$

In the formula, t is the time, and ρ is the density of fluid, and μ is the kinetic viscosity of fluid, u _i, u _jBe the hourly value of fluid velocity component, u _i', u _j' be the pulsating quantity of speed component, p is the time average of pressure, i and j equal respectively 1,2,3 denotation coordination components; Here upper line "-" expression is measured time average to physics;

S _iEquation of momentum source item, the fluid mass S outside the porous medium border _i=0, inner on the porous medium border,

$S_i=-(D_{\mathrm{ij}}\mathrm{\&mu;u}+C_{\mathrm{ij}}\frac{1}{2}\mathrm{\&rho;}|u\left|u\right),$

$D_{\mathrm{ij}}=\left(\begin{array}{ccc}{D}_n& 0& 0\\ 0& D_t& 0\\ 0& 0& D_t\end{array}\right),$ $C_{\mathrm{ij}}=\left(\begin{array}{ccc}{C}_n& 0& 0\\ 0& C_t& 0\\ 0& 0& C_t\end{array}\right)---\left(6\right)$

In the formula, D _Ij, C _IjBe the porous medium matrix of coefficients, u is fluid velocity, D _nRepresentation is to viscosity factor, D _tExpression tangential viscous resistance coefficient, C _nRepresentation is to inertia resistance coefficient, C _tRepresent tangential inertia resistance coefficient;

The k equation:

$\frac{\&PartialD;\left(\mathrm{\&rho;k}\right)}{\&PartialD;t}+\frac{\&PartialD;\left(\mathrm{\&rho;k}{u}_i\right)}{\&PartialD;x_i}=\frac{\&PartialD;}{\&PartialD;x_j}\left[(\mathrm{\&mu;}+\frac{{\mathrm{\&mu;}}_t}{{\mathrm{\&sigma;}}_k})\frac{\&PartialD;k}{\&PartialD;x_j}\right]+G_k+G_b-\mathrm{\&rho;\&epsiv;}-Y_M+S_k---\left(7\right)$

The ε equation:

$\frac{\&PartialD;\left(\mathrm{\&rho;\&epsiv;}\right)}{\&PartialD;t}+\frac{\&PartialD;\left(\mathrm{\&rho;\&epsiv;}{u}_i\right)}{\&PartialD;x_i}=\frac{\&PartialD;}{\&PartialD;x_j}\left[(\mathrm{\&mu;}+\frac{{\mathrm{\&mu;}}_t}{{\mathrm{\&sigma;}}_{\mathrm{\&epsiv;}}})\frac{\&PartialD;\mathrm{\&epsiv;}}{\&PartialD;x_j}\right]+\mathrm{\&rho;}{C}_1\mathrm{S\&epsiv;}-\mathrm{\&rho;}{C}_2\frac{{\mathrm{\&epsiv;}}^2}{k+\sqrt{\mathrm{v\&epsiv;}}}+C_{1\mathrm{\&epsiv;}}\frac{\mathrm{\&epsiv;}}{k}{C}_{3\mathrm{\&epsiv;}}{G}_b+S_{\mathrm{\&epsiv;}}---\left(8\right)$

In the formula, x _i, x _jBe the component of different directions under the cartesian coordinate system, k is the Hydrodynamic turbulence energy, and ε is turbulence dissipation rate, and μ is the kinetic viscosity of fluid, μ _tTurbulent viscosity for fluid;

$C_1=\max [0.43,\frac{\mathrm{\&eta;}}{\mathrm{\&eta;}+5}],\mathrm{\&eta;}=S\frac{k}{\mathrm{\&epsiv;}},S=\sqrt{2S_{\mathrm{ij}}{S}_{\mathrm{ij}}},S_{\mathrm{ij}}=\frac{1}{2}(\frac{\&PartialD;u_i}{\&PartialD;x_j}+\frac{\&PartialD;u_j}{\&PartialD;x_i})---\left(9\right)$

In the equation, G _kThe Turbulent Kinetic that is produced by average velocity gradient, G _bThe Turbulent Kinetic that is produced by buoyancy, Y _MThe overall dissipation rate that represents pulsation diffusion profile in the compressible turbulent flow, C ₂, C _{1 ε}, C _{3 ε}Be constant, S is strain modulus, σ _kAnd σ _εThe turbulent flow Prandtl number of k equation and ε equation, S _kAnd S _εIt is user-defined source item; Realizable k-ε model constants value C _{1 ε}=1.44, C ₂=1.9, σ _k=1.0, σ _ε=1.2;

B2, porous medium coefficient

Porous media model is to combine one in the porous medium zone rule of thumb to be assumed to be main resistance to flow, and this Resistance Value can be calculated by following formula:

F＝S _xλA (10)

In the formula, S _xBe the equation of momentum source item of x direction, λ is porous medium thickness, and A is the porous medium area; The F direction is opposite with water (flow) direction, gets its numerical values recited herein;

Equation (6) substitution equation (10) can be obtained water resistance F _dWith lift F _lExpression formula,

$F_d=(D_n\mathrm{\&mu;u}+C_n\frac{1}{2}\mathrm{\&rho;}|u\left|u\right)\mathrm{\&lambda;A}---\left(11\right)$

$F_l=(D_t\mathrm{\&mu;u}+C_t\frac{1}{2}\mathrm{\&rho;}|u\left|u\right)\mathrm{\&lambda;A}---\left(12\right)$

The porous medium coefficient need to be calculated by the data of model test; If a known angle of attack is 90 ° and plane etting the Resistance Value under different in flow rate vertical with water (flow) direction, adopt least square method to function F _d(u) carry out match with the flow velocity of model test and corresponding Resistance Value, obtain optimum porous medium coefficient D _n, C _nBecause the plane etting is vertical with water (flow) direction, so the lift value is 0, D in the numerical simulation _t, C _tCan ignore;

When the plane etting becomes certain attack angle alpha with current, must be with porous medium matrix of coefficients D _IjAnd C _IjCarry out calculating after the coordinate conversion:

$D_n^{\&prime;}=\frac{D_n+D_t}{2}+\frac{D_n-D_t}{2}\cos \left(2{\mathrm{\&alpha;}}^{\&prime;}\right)$ (13)

$C_n^{\&prime;}=\frac{C_n+C_t}{2}+\frac{C_n-C_t}{2}\cos \left(2{\mathrm{\&alpha;}}^{\&prime;}\right)$

$D_t^{\&prime;}=\frac{D_n-D_t}{2}\sin \left(2{\mathrm{\&alpha;}}^{\&prime;}\right)$ (14)

$C_t^{\&prime;}=\frac{C_n-C_t}{2}\sin \left(2{\mathrm{\&alpha;}}^{\&prime;}\right)$

Wherein, α '=90 °-α, α is the angle of attack;

Adopt least square method that model test value and the theoretical value of resistance and lift are carried out match, obtain making the porous medium coefficient of total error amount minimum;

B3, boundary condition and grid are divided

Take the angle of attack as 90 ° plane etting model as example, the global coordinate system of numerical model (x, y, z) true origin is taken at the projection of etting center on the water surface, and the x axle is water (flow) direction, and the y axle is direction vertical with the x axle on the surface level, and the z axle is vertical; The numerical model left end is the speed inlet boundary; Right-hand member is defined as the free discharge border; Tank sidewall, tank bottom surface and the water surface are defined as the Gu Bi border; Carry out grid divides in the Numerical modelling zone; Trellis-type is the destructuring tetrahedral grid, porous medium zone and near refined net;

B4, numerical evaluation

Adopt finite volume method (FVM) discrete, discrete scheme adopts the Second-order Up-wind form of convective motion, and preferably SIMPLEC algorithm of stability is selected in pressure-speed coupling, and inferior relaxation factor keeps default value.

Images