Bridge scour simulation and boundary setting method for SPH water-sand two-phase flow

Abstract

The application discloses a bridge scour simulation and boundary setting method for SPH water-sand two-phase flow, and comprises the following steps: model construction, including setting a bed sand layer area range and defining an inlet and outlet boundary area range; particle construction, including setting fluid particles and bed sand particles of a flow field area based on a multiphase flow theory, setting soil yield strength based on a HBP model and introducing a DP yield criterion, initializing inlet and outlet boundary particle parameters, setting particle type identification, and setting particle position identification; N-S control equation solving, obtaining corresponding flow force parameter results; entering a sediment solving module, including sediment particle conversion and viscosity updating calculation; entering an inlet and outlet boundary updating module, including inlet and outlet particle position retrieval, particle type conversion, and particle supplement and deletion; completing a calculation cycle and entering the next cycle; through secondary development based on the SPH multiphase flow theory, the application can realize numerical simulation of foundation local scour by applying the inlet and outlet boundaries to the multiphase flow model and introducing various structure and state models.

Description

TECHNICAL FIELD

[0001] The application belongs to the field of bridge scour simulation, and particularly relates to a bridge scour simulation method for SPH water-sand two-phase flow and a boundary setting method. BACKGROUND

[0002] Bridge foundation scour disease is one of the main reasons for the failure of the function of the bridge structure and the loss of safety performance, which has attracted widespread attention from many scholars. Numerical simulation is still one of the most effective methods for studying bridge scour, which has many advantages such as low cost, high efficiency, and short cycle. The mainstream algorithm used today is to obtain fluid dynamics based on the Euler form solution of the N-S fluid control equation, and to obtain the numerical solution of the bed elevation change by introducing a numerical bed sand transport model. However, the numerical simulation method based on the Euler form often has difficulty in converging when solving the breaking of complex fluid surfaces and wave, tsunami overflows, etc., or requires a large amount of solving time and resources, and there is no real soil model, which makes it difficult to track the trajectory of soil evolution. The Smoothed Particle Hydrodynamics (SPH) method is a numerical solution method based on the Lagrangian form. Compared with the Euler method, the particle itself has mass, which can ensure mass conservation without additional calculation; when simulating complex free surface flow, it is not necessary to track the fluid boundary and different fluid interfaces. When simulating scour problems, a soil model can be constructed without introducing a numerical sand transport model, which can further improve efficiency and accuracy. However, at present, there are few studies and application cases of scour simulation based on SPH, the simulation process is not clear, there is a lack of modular design scheme, and there is a lack of sediment calculation module, and the multiphase flow model cannot apply import and export boundaries, which need to be further solved. SUMMARY

[0003] The purpose of the present application is to provide a bridge scour simulation method for SPH water-sand two-phase flow and a boundary setting method, to apply import and export boundaries to the multiphase flow model, and to introduce a plurality of theoretical or experimental criteria for sediment particle transformation and evolution judgment, and to simulate the local scour of the foundation in detail.

[0004] To achieve the above-mentioned purpose, the present application provides the following technical scheme: a bridge scour simulation method for SPH water-sand two-phase flow and a boundary setting method, comprising the following steps:

[0005] S1, based on the range of the bed sand layer region, the range of the bed sand layer region is divided into an import boundary region range, an export boundary region range and a region range within a drainage basin;

[0006] S2, based on the basic principle of multiphase flow, a fluid particle model of Newton fluid and a soil particle model of non-Newton fluid are constructed respectively; based on the physical properties of the shear stress and the strain rate nonlinearity of the non-Newton fluid, a HBP model is introduced to obtain a primary soil viscosity model μ ori ; a DP yield criterion is introduced to calculate the material yield stress τy the material yield stress τ y replacing the material yield stress τ ori in the primary soil viscosity model μ c , a soil viscosity calculation model μ HBP is constructed;

[0007] Different types of labels are then given to different fluid phases: the preset water body is a Newtonian fluid, type label mk1, the preset soil body is a non-Newtonian fluid, type label mk2, the preset bed load particles are non-Newtonian fluids, type label mk3, and the preset suspended load particles are Newtonian fluids, type label mk4;

[0008] Then, based on the fluid particle model of the Newtonian fluid, the soil particle model of the non-Newtonian fluid, and the inlet boundary region range, the outlet boundary region range, and the region range within the flow basin, the particle position label of the preset inlet boundary region range is code2, the particle position label of the outlet boundary region range is code3, and the particle position label of the region range within the flow basin is code1; according to the water depth H of the position where the outermost particles of the inlet and outlet boundaries are located, the inlet and outlet particle flow velocity initialization functions V(H) are respectively preset;

[0009] S3, particle coupling calculation is performed on the inlet and outlet boundary particles and the particles within the flow basin, based on all selected particles, all particles within the nuclear function radius of a single particle are retrieved in turn, according to the type label of the retrieved particles, the corresponding particle viscosity is called to perform viscous stress coupling calculation with the retrieved particles, and the speed and position pos information of all selected particles are updated; S4, according to the speed of the particles obtained in step S3, the fluid shear stress τ b experienced by all soil particles with the identification number mk2 is calculated in turn using the Einstein logarithmic flow velocity distribution formula, based on the soil viscosity calculation model μ HBP , and the Shield criterion of sediment particles is used to calculate the corresponding sediment incipient critical stress τ cr,0 ;

[0010] According to the sediment incipient critical stress τ cr,0 , the modified sediment incipient critical stress τ cr considering the slope effect is calculated; according to the fluid shear stress τ b and the modified sediment incipient critical stress τ cr , the starting state of all soil particles with the type label mk2 is judged in turn, if τ cr ≥ τ b , τ cr is substituted into the soil viscosity calculation model μ HBP , the yield stress τ y is replaced, and the soil viscosity calculation model μ HBPand the soil particle is updated to be a bed load particle, and the type label mk2 is updated to be mk3, otherwise the soil viscosity calculation model μ is not updated HBP ;

[0011] According to the updated bed load particle and its type label mk3, the critical flow velocity u of the bed load particle is calculated by using the Mastbergen formula lift If the actual flow velocity u of the particle satisfies u lift , the bed load particle is updated to be a suspended load particle, and the type label mk3 is updated to be mk4, and the equivalent viscosity μ of the suspended load particle is calculated lift The soil viscosity calculation model μ is updated HBP to be the model μ lift If the condition is not satisfied, the bed load particle is not updated

[0012] According to the updated suspended load particle and its type label mk4, the critical flow velocity u of the suspended load particle is calculated by using the Mastbergen formula set If the actual flow velocity u of the particle satisfies u set , the suspended load particle is updated to be a bed load particle, and the type label mk4 is updated to be mk3, and τ cr,0 is substituted into the soil viscosity calculation model μ HBP , the yield stress τ y is replaced, and the soil viscosity calculation model μ is updated HBP Otherwise, the suspended load particle and the model μ are not updated lift ;

[0013] S5, according to the particle position pos information obtained in step S3, the inlet particle is updated, and the particle position pos information with the position label code2 is searched in sequence, and compared with the inlet boundary area range In, if there is a particle exceeding the corresponding range of In, the particle is converted to be a flow domain particle, and the position label is converted to be code1; then new inlet particles are supplemented in the outermost layer of the inlet boundary, and the position label code2 and the type label mk1 are preset, and the particle flow velocity is initialized; otherwise, according to the particle position pos information obtained in step S3, the outlet particle is updated, and the particle position information with the position label code3 is searched in sequence, and compared with the outlet boundary area range Out, if there is a particle exceeding the corresponding range of Out, the particle is deleted

[0014] According to the particle position pos information obtained in step S3, the outlet particle is updated, and the particle position pos information with the position label code1 is searched in sequence, and compared with the outlet boundary area range Out, if there is a particle entering the corresponding range of Out, the particle is converted to be an outlet boundary particle, and the position label is converted to be code3 and the type label is mk1

[0015] According to the particle position pos information obtained in step S3, inlet particle updating is performed, and the position code1 particle position pos information is sequentially searched, compared with the inlet boundary region range In, and if there is a particle entering the range corresponding to In, the particle is deleted;

[0016] S6, return to perform steps S3 to S5 until all the particles in the inlet boundary region range, the outlet boundary region range, and the region range in the drainage basin are traversed.

[0017] Further, the aforementioned step S1 includes the following sub-steps:

[0018] S1.1, based on the bed sand layer region range, preset the bed sand layer outermost position distance from the inlet and outlet boundary to exceed the kernel function radius h adopted by the SPH algorithm by one time;

[0019] S1.2, define the inlet boundary region range set as In[x min ,x max ,y min ,y max ,z min ,z max ], wherein x min ,x max are the minimum X-axis coordinate value and the maximum X-axis coordinate value of the inlet boundary region range, y min ,y max are the minimum Y-axis coordinate value and the maximum Y-axis coordinate value of the inlet boundary region range, and z min ,z max are the minimum Z-axis coordinate value and the maximum Z-axis coordinate value of the inlet boundary region range;

[0020] S1.3, define the outlet boundary region range set as Out[x’ min ,x’ max ,y’ min ,y’ max ,z’ min ,z’ max ], wherein x’ min ,x’ max are the minimum X-axis coordinate value and the maximum X-axis coordinate value of the inlet boundary region range, y’ min ,y’ max are the minimum Y-axis coordinate value and the maximum Y-axis coordinate value of the inlet boundary region range, and z’ min ,z’ max are the minimum Z-axis coordinate value and the maximum Z-axis coordinate value of the inlet boundary region range.

[0021] Further, in the aforementioned step S2, the primary soil body viscosity model μ ori is constructed according to the following formula:

[0022]

[0023] where τ c is the material yield stress, II D is the second invariant of the fluid strain rate tensor, m is the stress exponent growth coefficient, μ is the water viscosity, and n is the power related to the shear stress.

[0024] Further, in the aforementioned step S2, the material yield stress τ y is calculated according to the following formula:

[0025] |τ y | = αp + β

[0026] where p is the hydrostatic pressure acting on the saturated sediment particles, and α and β are given by the Mohr-Coulomb yield criterion parameters as follows:

[0027]

[0028] where θ is the internal friction angle and c is the soil cohesion.

[0029] Further, in the aforementioned step S2, a soil viscosity calculation model μ HBP is constructed according to the following formula:

[0030]

[0031] Further, in the aforementioned step S2, according to the water depth H of the outermost particles at the inlet and outlet boundaries, the inlet and outlet particle velocity initialization functions V(H) are preset as follows:

[0032] For the inlet particle velocity initialization function V(H), it includes:

[0033] (1) V = K, K is a constant, and the velocity is constant;

[0034] (2) V = aH + b, a and b are undetermined coefficients, and the velocity is linearly distributed with depth;

[0035] (3) V = mH 2 +nH + k, m, n, and k are undetermined coefficients, and the velocity is second-order distributed with depth;

[0036] For the outlet particle velocity initialization function V(H), it includes: V = K, K is a constant, and the velocity is constant.

[0037] Further, in the aforementioned step S4, the fluid shear stress τ b is calculated according to the following formula:

[0038]

[0039] Where d is the characteristic particle size; Δu is the velocity difference between the sediment particle and the nearby water particle; κ is the von Kármán constant; and ρ is the density of the sediment particle.

[0040] Calculate the critical stress τ for sediment initiation using the following formula. cr,0 :

[0041] τ cr，0 =θ cr ·(ρ s -ρ)gd,

[0042] Where, θ cr ρ is the critical Shields number. s ρ is the density of saturated sediment, ρ is the density of water, g is the acceleration due to gravity, and d is the particle size of soil.

[0043] The modified critical stress τ for sediment initiation, taking into account slope effects, is calculated using the following formula. cr :

[0044]

[0045] Where η is a constant, α, β, and γ are the angles between the slope normal vector and the X, Y, and Z axes, respectively, and θ is the angle between the slope normal vector and the X, Y, and Z axes. x θ y θ z These are the angles between the unit vector of the flow direction and the positive directions of the X, Y, and Z axes, respectively. Indicates the internal friction angle of mud and sand;

[0046] Calculate the critical velocity u using the following formula. lift :

[0047]

[0048] Where, α i n is the sediment transport coefficient. s Let d be the normal vector of the bed surface. * The particle size coefficient of sediment. ρ s ρ is the saturated sediment density, ρ is the water density, g is the gravitational acceleration, d is the soil particle size, μ is the water viscosity, and θ is the saturated sediment density. b θ is the actual Shields number of soil particles. cr The critical Shields number;

[0049] Calculate the equivalent viscosity μ using the following formula. lift :

[0050]

[0051] Where μ is the viscosity of water, C vThe concentration of the sediment particles within the radius of the kernel function.

[0052] The critical flow velocity u is calculated according to the following formula set :

[0053]

[0054] Wherein, mu is the viscosity of the water body, d is the particle size of the soil body, d * is the sediment particle size coefficient.

[0055] Further, in the step S5, the particle exceeding the range of In should satisfy the following conditions:

[0056] pos x ≥x max ,

[0057] Wherein, pos x is the particle position information X axis coordinate;

[0058] The particle exceeding the range designated by Out should satisfy the following conditions:

[0059] pos x ≥x' max ||pos y ≥y' max ||pos z ≥z' max ,

[0060] Wherein, pos z ,pos z is the particle position information Y, Z axis coordinate;

[0061] The particle entering the range designated by Out should satisfy the following conditions:

[0062] pos x ≥x' min ,

[0063] The particle entering the range designated by In should satisfy the following conditions:

[0064] pos x ≤x max .

[0065] Further, the aforementioned import flow velocity is directed from the X axis negative direction to the X axis positive direction, and the gravity direction is the Z axis negative direction.

[0066] Compared with the prior art, the present application has the following beneficial effects:

[0067] 1. The present application is based on the secondary development of SPH multiphase flow algorithm, the fine simulation of the whole process of sediment particle scouring and evolution is realized by setting the import and export boundaries and introducing multiple theoretical or experimental criteria, so as to realize the modular design of the numerical simulation of local scouring of bridge foundation, and the present application has the characteristics of easy programming, high accuracy and strong operability.

[0068] 2. The present application realizes the numerical simulation of scouring based on the secondary development of SPH multiphase flow, the particles themselves have mass, and the mass conservation can be ensured without additional calculation; when simulating the scouring problem, the soil body model can be constructed, and the numerical sand transport model does not need to be introduced, which can further speed up the efficiency and improve the accuracy. BRIEF DESCRIPTION OF DRAWINGS

[0069] Figure 1 It is the bridge scouring simulation and boundary setting method flow chart of the present application facing SPH water-sand two-phase flow.

[0070] Figure 2 It is the bed sand layer and import and export boundary area range diagram defined by the present application.

[0071] Figure 3 It is the import boundary particle distribution and conversion diagram simulated by the present application.

[0072] Figure 4 It is the export boundary particle distribution and conversion diagram simulated by the present application.

[0073] Figure 5 It is the modified sediment incipient critical stress τ cr force diagram considering the effect of slope modified by the present application;

[0074] Figure 6 It is the sediment distribution diagram based on the present application.

[0075] Figure 7 It is an enlarged view of the interface between bedload and suspended load. DETAILED DESCRIPTION

[0076] In order to better understand the technical content of the present application, specific embodiments are described below with reference to the accompanying drawings.

[0077] In the present application, aspects of the present application are described with reference to the accompanying drawings, which show many illustrative embodiments. The embodiments of the present application are not limited to the drawings described. It should be understood that the present application is realized by any one of the above-mentioned various concepts and embodiments, and the concepts and embodiments described in detail below, because the concepts and embodiments disclosed in the present application are not limited to any embodiment. In addition, some aspects disclosed in the present application can be used alone, or in any suitable combination with other aspects disclosed in the present application.

[0078] The distribution of import / export channels and sediment layers upon which this invention is based is as follows: Figure 2 As shown, the inlet / outlet boundary 1 is located on the left, and the sediment layer 2 is located on the right, with a distance of L between their edges.

[0079] The inlet boundary particle distribution and transformation upon which this invention is based are as follows: Figure 3 As shown, inlet region particles 5 and watershed particles 4 are distributed above the boundary particle 3; particles that exceed the boundary of the inlet region are designated as particles to be converted 6, and are converted into watershed particles 4 in subsequent steps; at the same time, new particles 7 will be added to the inlet region to supplement it.

[0080] The inlet boundary particle distribution and transformation upon which this invention is based are as follows: Figure 4 As shown, exit region particles 11 and watershed particles 10 are distributed above the boundary particle 9; particles that exceed the exit region boundary 14 are designated as particles to be excluded 13 and deleted in subsequent steps; at the same time, watershed particles 10 that newly enter the exit region boundary 14 are designated as particles to be transformed 12 and transformed into exit region particles 11 in subsequent steps.

[0081] The modified critical stress τ for sediment initiation that takes into account slope effects, as described in this invention. cr Sediment particles are subjected to forces as follows Figure 5 As shown, the sediment particle is a point mass on the sloping bed surface, and the forces acting on it include underwater gravity W (including gravity and buoyancy) and water flow drag F. D Lift force F of water flow L Resistance to starting force F C (In this article, this refers to the Coulomb friction of the soil.) It is a unit vector in an orthogonal rectangular coordinate system; The plane in question is a sloping bed surface, in which Located in the XOZ plane, making an angle α with the X-axis, It lies in the YOZ plane, making an angle β with the Y-axis; The unit normal vector of the sloping surface makes an angle γ with the positive Z-axis, and the drag force of the water flow is... Direction and flow velocity They are in the same direction.

[0082] The sediment distribution on which this invention is based is as follows: Figure 6 As shown, the riverbed is composed of ordinary soil 15, and the scour pit is composed of yielding soil 16. Inside the scour pit, there are bedload 17 and suspended soil 18. Bedload 17 is adjacent to the yielding soil 16, and suspended soil 18 is distributed on top of bedload 17. Figure 7 This is an enlarged view of the interface between bedload and suspended mass.

[0083] The specific implementation method of the application takes the SPH open source calculation software Dualsphysics as an example, the solving of the multiphase flow N-S control equation is processed by the embedded algorithm, the conversion of different flow phases of particles can be realized based on the change of particle type label, and the program can automatically bring in the corresponding solving module for solving by identifying the particle type label and the position label.

[0084] As shown in the flowchart of the application, the bridge scour simulation and boundary setting method for the SPH water-sand two-phase flow includes the following steps: Figure 1

[0085] S1, based on the bed sand layer area range, the bed sand layer area range is divided into an import boundary area range, an export boundary area range and an area range in a drainage basin, including the following sub-steps:

[0086] S1.1, based on the bed sand layer area range, the outermost position of the bed sand layer is preset to be more than one time of the radius h of the kernel function adopted by the SPH algorithm from the import and export boundary.

[0087] S1.2, the import boundary area range set is defined as In[x min ,x max ,y min ,y max ,z min ,z max ], wherein x min ,x max are the minimum X-axis coordinate value and the maximum X-axis coordinate value of the import boundary area range respectively, y min ,y max are the minimum Y-axis coordinate value and the maximum Y-axis coordinate value of the import boundary area range respectively, and z min ,z max are the minimum Z-axis coordinate value and the maximum Z-axis coordinate value of the import boundary area range respectively.

[0088] S1.3, the export boundary area range set is defined as Out[x’ min ,x’ max ,y’ min ,y’ max ,z’ min ,z’ max ], wherein x’ min ,x’ max are the minimum X-axis coordinate value and the maximum X-axis coordinate value of the import boundary area range respectively, y’ min ,y’ max are the minimum Y-axis coordinate value and the maximum Y-axis coordinate value of the import boundary area range respectively, and z’ min ,z’ max ​are the minimum and maximum Z-axis coordinate values of the inlet boundary region range respectively. The inlet flow velocity is directed from the X-axis negative direction to the X-axis positive direction, and the gravity direction is the Z-axis negative direction.

[0089] S2, a fluid particle model of Newtonian fluid and a soil particle model of non-Newtonian fluid are respectively constructed based on the principle of multiphase flow; based on the physical property of non-Newtonian fluid shear stress and strain rate nonlinearity, the HBP model is introduced to obtain a primary soil viscosity model μ ori , as follows:

[0090]

[0091] wherein τ c is the material yield stress, II D is the second invariant of the fluid strain rate tensor, m is the stress exponent growth coefficient, μ is the water viscosity, and n is the power related to the shear stress.

[0092] Then the DP yield criterion is introduced to calculate the material yield stress τ y , as follows:

[0093] |τ y | = αp + β,

[0094] wherein p is the static water pressure acting on the saturated sediment particles, and α and β are both given by the Mohr-Coulomb yield criterion parameters, as follows:

[0095]

[0096] wherein θ is the internal friction angle, and c is the soil cohesion.

[0097] The material yield stress τ y is replaced in the primary soil viscosity model μ ori , the material yield stress τ c is replaced, and a soil viscosity calculation model μ HBP is constructed, as follows:

[0098]

[0099] Then different type labels are given for different fluid phases: the preset water body is Newtonian fluid, the type label is mk1, the preset soil body is non-Newtonian fluid, the type label is mk2, the preset bed load particle is non-Newtonian fluid, the type label is mk3, and the preset suspended load particle is Newtonian fluid, the type label is mk4.

[0100] Then based on the fluid particle model of Newtonian fluid, the soil particle model of non-Newtonian fluid, and the bed sand layer import boundary area range, the outlet boundary area range, the area range in the drainage basin, the particle position index of the import boundary area range is preset as code2, the particle position index of the outlet boundary area range is preset as code3, and the particle position index of the area range in the drainage basin is preset as code1; according to the water depth H of the position where the outermost particles of the import and outlet boundaries are located, the import and outlet particle flow velocity initialization functions V(H) are respectively preset, and the specific contents are as follows:

[0101] For the import particle flow velocity initialization function V(H), the following contents are included:

[0102] (1) V = K, K is a constant, and the flow velocity is constant;

[0103] (2) V = aH + b, a and b are undetermined coefficients, and the flow velocity is linearly distributed with the depth;

[0104] (3) V = mH 2 +nH+k, m, n and k are undetermined coefficients, and the flow velocity is distributed in the second order with the depth;

[0105] For the outlet particle flow velocity initialization function V(H), the following contents are included: V = K, K is a constant, and the flow velocity is constant.

[0106] S3, particle coupling calculation is performed on the import and outlet boundary particles and the particles in the drainage basin, based on all selected particles, all particles in the radius of the core function of a single particle are searched in turn, according to the searched particle type index, the corresponding particle viscosity is called to perform viscous stress coupling calculation with the searched particle, and the speed and position pos information of all selected particles are updated; S4, according to the speed of the particle obtained in step S3, the fluid shear stress τ b experienced by the soil particles with the identification number mk2 is calculated in turn by using the Einstein logarithmic flow velocity distribution formula, as follows:

[0107]

[0108] Wherein, d is the characteristic particle size of the particle; Δu is the velocity difference between the sediment particle and the nearby water particle, κ is the von Karman constant, which can be taken as 0.41 generally, and ρ is the density of the sediment particle.

[0109] Based on the soil viscosity calculation model μ HBP , the Shield criterion of the sediment particle is used to calculate the corresponding sediment starting critical stress τ cr,0 , as follows:

[0110] τ cr，0 =θ cr ·(ρ s -ρ)gd,

[0111] wherein θ cr is the critical Shields number, only related to the sediment parameters, ρ s is the saturated sediment density, ρ is the water density, g is the gravity acceleration, and d is the soil particle size.

[0112] According to the sediment incipient critical stress τ cr,0 , the modified sediment incipient critical stress τ cr considering the slope effect is calculated as follows:

[0113]

[0114] wherein η = 0.7, α, β, and γ are the angles between the slope normal vector and the X, Y, and Z axis directions, respectively, θ x , θ y , and θ z are the angles between the flow direction unit vector and the X, Y, and Z axis positive directions, respectively, represents the sediment internal friction angle;

[0115] According to the fluid shear stress τ b and the modified sediment incipient critical stress τ cr , the incipient state of the soil particle with the particle type label mk2is determined in sequence. If τ cr ≥ τ b , τ cr is substituted into the soil viscosity calculation model μ HBP , the yield stress τ y is replaced, the soil viscosity calculation model μ HBP is updated, and the soil particle is updated to a bed load particle, and the corresponding type label mk2is updated to mk3, otherwise the soil viscosity calculation model μ HBP is not updated.

[0116] According to the updated bed load particle and its type label mk3, the critical flow velocity u lift is calculated by using the Mastbergen formula as follows:

[0117]

[0118] wherein α i is the sediment transport coefficient, n s is the bed normal vector, d * is the sediment particle size coefficient, ρ s is the saturated sediment density, ρ is the water density, g is the gravity acceleration, d is the soil particle size, μ is the water viscosity, θ b is the actual Shields number of the soil particle, and θ cr is the critical Shields number.

[0119] If the actual flow velocity of the particle u is greater than or equal to the threshold value u lift , the bed load particle is updated to the suspended load particle, and the type code mk3 is updated to mk4, and the equivalent viscosity μ lift of the suspended load particle is calculated as follows:

[0120]

[0121] where μ is the viscosity of the water body, and C v is the concentration of the sediment particles in the kernel function radius.

[0122] The soil viscosity calculation model μ HBP is updated to μ lift , and if the condition is not met, the bed load particle is not updated.

[0123] According to the updated suspended load particle and the type code mk4, the threshold flow velocity u set of the suspended load particle is calculated by using the Mastbergen formula as follows:

[0124]

[0125] where μ is the viscosity of the water body, d is the particle size of the soil, and d * is the sediment particle size coefficient.

[0126] If the actual flow velocity of the particle u is less than or equal to the threshold value u set , the suspended load particle is updated to the bed load particle, the type code mk4 is updated to mk3, τ cr,0 is substituted into the soil viscosity calculation model μ HBP , the yield stress τ y is replaced, the soil viscosity calculation model μ HBP is updated, and otherwise the suspended load particle and the model μ lift are not updated.

[0127] S5, according to the particle position pos information obtained in step S3, the inlet particle is updated, the particle position pos information with the position code code2 is searched in sequence, and is compared with the inlet boundary region range In, if there is a particle exceeding the range corresponding to In, the particle is converted to a watershed particle, the position code is converted to code1, and then new inlet particles are supplemented in the outermost layer of the inlet boundary, and the position code code2 and the type code mk1 are preset, and the particle flow velocity is initialized.

[0128] The particle exceeding the range of In should meet the following condition:

[0129] pos x ≥ x max ,

[0130] where posx is the X-axis coordinate of the particle position information X;

[0131] Otherwise, the outlet particle is updated according to the particle position pos information obtained in step S3, the particle position information with the position code code3 is searched in sequence, and compared with the outlet boundary region range Out. If there is a particle exceeding the range corresponding to Out, the particle is deleted. The particle exceeding the range specified by Out should satisfy the following conditions:

[0132] pos x ≥ x' max || pos y ≥ y' max || pos z ≥ z' max ,

[0133] pos z pos z are the Y-axis and Z-axis coordinates of the particle position information X;

[0134] The outlet particle is updated according to the particle position pos information obtained in step S3, the particle position pos information with the position code code1 is searched in sequence, and compared with the outlet boundary region range Out. If there is a particle entering the range corresponding to Out, the particle is converted into an outlet boundary particle, the position code is converted into code3, and the type code is mk1. The particle entering the range specified by Out should satisfy the following conditions:

[0135] pos x ≥ x' min ,

[0136] The inlet particle is updated according to the particle position pos information obtained in step S3, the particle position pos information with the position code code1 is searched in sequence, and compared with the inlet boundary region range In. If there is a particle entering the range corresponding to In, the particle is deleted. The particle entering the range specified by In should satisfy the following conditions:

[0137] pos x ≤ x max .

[0138] S6, return to execute steps S3 to S5 until all the particles in the inlet boundary region range, the outlet boundary region range, and the region range in the drainage basin are traversed.

[0139] Although the present application has been described above with reference to a preferred embodiment, it is not intended to limit the present application. Those skilled in the art, without departing from the spirit and scope of the present application, can make various modifications and improvements. Therefore, the scope of protection of the present application should be defined by the claims.

Claims

1. A bridge scour simulation method and boundary setting method for SPH water-sand two-phase flow, characterized in that, Comprising the following steps: S1, based on the bed sand layer area range, the bed sand layer area range is divided into import boundary area range, export boundary area range and flow area range; S2. Based on the fundamental principles of multiphase flow, construct fluid particle models for Newtonian fluids and soil particle models for non-Newtonian fluids respectively; based on the nonlinear physical properties of shear stress and strain rate in non-Newtonian fluids, introduce the HBP model to obtain the primary soil viscosity model μ. ori Introducing the DP yield criterion, the yield stress τ of the material is calculated. y The yield stress τ of the material y Replace the primary soil viscosity model μ ori The material yield stress τ c Construct a soil viscosity calculation model μ HBP ; Then different type labels are given for different fluid phases: preset water body is Newton fluid, type label mk1, preset soil body is non-Newton fluid, type label mk2, preset bed load particle is non-Newton fluid, type label mk3, and preset suspended load particle is Newton fluid, type label mk4; Then, based on the fluid particle model of Newton fluid, the soil body particle model of non-Newton fluid, and the bed sand layer import boundary area range, export boundary area range and flow area range, the particle position label of the preset import boundary area range is code2, the particle position label of the export boundary area range is code3, and the particle position label of the flow area range is code1; according to the water depth H of the position where the outermost particles of the import and export boundaries are located, the import and export particle flow velocity initialization functions V(H) are preset respectively; S3, particle coupling calculation is performed on the import, export boundary particles and flow area particles, based on all selected particles, all particles within the kernel function radius of a single particle are searched in turn, according to the type label of the searched particles, the corresponding particle viscosity is called to perform viscous stress coupling calculation with the searched particles, and the speed and position pos information of all selected particles are updated; S4, the velocity of the particle obtained according to step S3, the fluid shear stress τ experienced by the soil particle with the particle identification number mk2 is calculated in turn using the Einstein logarithmic flow velocity distribution formula b ; based on the soil viscosity calculation model μ HBP , the corresponding sediment incipient critical stress τ is calculated using the Shield criterion for sediment particles cr,0 ; According to the sediment incipient critical stress τ cr,0 , a modified sediment incipient critical stress τ cr considering the slope effect is calculated; according to the fluid shear stress τ b and the modified sediment incipient critical stress τ cr , the starting state of the soil particle with the particle type label mk2 is sequentially judged, if τ cr ≥τ b , τ cr is substituted into the soil viscosity calculation model μ HBP , the yield stress τ y is replaced, the soil viscosity calculation model μ HBP is updated, and the soil particle is updated to the bed load particle, while the type label mk2 is updated to mk3, otherwise the soil viscosity calculation model μ HBP is not updated. Based on the updated bedload particles and their type designation mk3, the critical velocity u was calculated using the Mastbergen formula. lift If the actual particle velocity u ≥ u lift Then, the bedload particles are updated to suspended particles, and the type label mk3 is updated to mk4 accordingly. The equivalent viscosity μ of the suspended particles is then calculated. lift Update the soil viscosity calculation model μ HBP For model μ lift If the condition is not met, the transport particle will not be updated. According to the updated suspended particle and its type label mk4, the critical flow velocity u set , if the actual flow velocity of the particle u set , the suspended particle is updated to the pushed particle, the type label mk4 is updated to mk3, τ cr,0 is substituted into the soil viscosity calculation model μ HBP , the yield stress τ y is replaced, the soil viscosity calculation model μ HBP is updated, otherwise the suspended particle and the model μ lift are not updated; S5, according to the particle position pos information obtained in step S3, the import particle is updated, the particle position pos information with position label code2 is searched in turn and compared with the import boundary area range In, if there is a particle exceeding the range corresponding to In, the particle is converted into a flow area particle, and the position label is converted into code1; then new import particles are supplemented in the outermost layer of the import boundary, and the position label code2 and the type label mk1 are preset and the particle flow velocity is initialized; otherwise, according to the particle position pos information obtained in step S3, the export particle is updated, the particle position information with position label code3 is searched in turn and compared with the export boundary area range Out, if there is a particle exceeding the range corresponding to Out, the particle is deleted; According to the particle position pos information obtained in step S3, the export particle is updated, the particle position pos information with position label code1 is searched in turn and compared with the export boundary area range Out, if there is a particle entering the range corresponding to Out, the particle is converted into an export boundary particle, the position label is converted into code3, and the type label is mk1; According to the particle position pos information obtained in step S3, the import particle is updated, the particle position pos information with position label code1 is searched in turn and compared with the import boundary area range In, if there is a particle entering the range corresponding to In, the particle is deleted; S6, return to execute step S3 to step S5 until all particles in the import boundary area range, export boundary area range and flow area range are traversed.

2. The SPH water-sand two-phase flow oriented bridge scour simulation and boundary setting method according to claim 1, characterized in that, Step S1 comprises the following substeps: S1.1, based on the bed sand layer area range, the bed sand layer outermost position distance from the import and export boundary exceeds the nuclear function radius h adopted by the SPH algorithm one time; S1.2, define the import boundary area range set as In[x min ,x max ,y min ,y max ,z min ,z max ], wherein x min ,x max are respectively the minimum X-axis coordinate value and the maximum X-axis coordinate value of the import boundary area range, y min ,y max are respectively the minimum Y-axis coordinate value and the maximum Y-axis coordinate value of the import boundary area range, and z min ,z max are respectively the minimum Z-axis coordinate value and the maximum Z-axis coordinate value of the import boundary area range; S1.3, define the export boundary region range set as Out[x min ,y max ,z min ], wherein x max ,y min ,z max are the minimum and maximum X, Y, Z axis coordinate values of the import boundary region range respectively. min ,x max are the minimum and maximum X axis coordinate values of the import boundary region range respectively, y min ,y max are the minimum and maximum Y axis coordinate values of the import boundary region range respectively, and z min ,z max are the minimum and maximum Z axis coordinate values of the import boundary region range respectively.

3. The SPH water-sand two-phase flow oriented bridge scour simulation and boundary setting method according to claim 1, characterized in that, In step S2, a primary soil body viscosity model μ is constructed as follows ori : where τ c is the yield stress of the material, II D is the second invariant of the fluid strain rate tensor, m is the stress exponent growth coefficient, μ is the water body viscosity, and n is the power related to the shear stress.

4. The bridge scour simulation and boundary setting method for SPH water-sand two-phase flow according to claim 3, characterized in that, In step S2, the material yield stress τ is calculated as follows y : |τ y |=αp+β Wherein, p is the hydrostatic pressure acting on the saturated sediment particles, and both a and β are given by Mohr-Coulomb yield criterion parameters, as follows: Wherein, θ is the internal friction angle, and c is the soil cohesion.

5. The SPH water-sand two-phase flow oriented bridge scour simulation and boundary setting method according to claim 4, characterized in that, In step S2, a soil viscosity calculation model μ is constructed according to the following equation HBP :

6. The SPH water-sand two-phase flow oriented bridge scour simulation and boundary setting method according to claim 5, characterized in that, In step S2, according to the water depth H of the position where the outermost particles of the import and export boundary are located, the import and export particle flow rate initialization function V(H) is preset as follows: For the import particle flow rate initialization function V(H), it includes: (1) V=K, K is a constant, and the flow rate is constant; (2) V=aH+b, a and b are undetermined coefficients, and the flow rate is linearly distributed with the depth; (3) V = mH 2 + nH + k, m, n, k are undetermined coefficients, the flow rate is distributed with depth to the second order; For the export particle flow rate initialization function V(H), it includes: V=K, K is a constant, and the flow rate is constant.

7. The SPH water-sand two-phase flow oriented bridge scour simulation and boundary setting method according to claim 6, characterized in that, In step S4, the fluid shear stress τ is calculated as follows b : Wherein, d is the characteristic particle size of the particle; Δu is the velocity difference between the sediment particles and the nearby water particles, κ is the von Karman constant, and ρ is the density of the sediment particles; The incipient critical stress of sediment τ is calculated by the following equation cr,0 : τ cr，0 = θ cr · (ρ s - ρ) gd, where θ cr is the critical Shields number, p s is the density of the saturated sediment, p is the density of the water, g is the acceleration of gravity, and d is the particle size of the soil. The modified sediment incipient critical stress τ considering the slope effect is calculated according to the following formula cr : wherein η is a constant value, α, β, and γ are the angles between the normal vector of the slope and the X, Y, and Z axis directions, respectively, θ x , θ y , and θ z are the angles between the unit vector of the flow direction and the X, Y, and Z axis positive directions, respectively, denotes the internal friction angle of the sediment. The critical flow velocity u is calculated according to the following formula lift : wherein α i is the sediment transport coefficient, n s is the bed normal vector, d * is the sediment particle size coefficient, p s is the saturated sediment density, p is the water density, g is the gravitational acceleration, d is the soil particle size, m is the water viscosity, q b is the actual Hillz number of the soil particle, q cr is the critical Hillz number; The equivalent viscosity μ is calculated according to the following formula lift : where μ is the viscosity of the water body, C v is the concentration of sediment particles within the kernel radius. The critical flow velocity u is calculated according to the following formula set : wherein μ is the viscosity of the water body, d is the particle size of the soil body, d * is the silt particle size coefficient.

8. The SPH water-sand two-phase flow oriented bridge scour simulation and boundary setting method according to claim 7, characterized in that, In step S5, the particles exceeding the range of In should satisfy the following conditions: pos x ≥ x max , wherein pos x is the particle position information X-axis coordinate; The particles exceeding the range specified by Out should satisfy the following conditions: pos x ≥ x' max || pos y ≥ y' max || pos z ≥ z' max , wherein pos z pos z is the particle position information Y, Z-axis coordinate; The particles entering the range specified by Out should satisfy the following conditions: pos x ≥ x' min , The particles entering the range specified by In should satisfy the following conditions: pos x ≤ x max . 9.The SPH water-sand two-phase flow oriented bridge scour simulation and boundary setting method according to claim 1, wherein, The import flow rate is directed from the negative direction of the X-axis to the positive direction of the X-axis, and the gravity direction is the negative direction of the Z-axis.

Images