A method for simulating weakly compressible fluid flow with velocity inlet using SPH

Abstract

The invention relates to a method for simulating the flow of a weakly compressible fluid with a velocity inlet by using SPH, comprising the following steps: initializing parameter information, the parameter information comprising vertex coordinates of a calculation domain, vertex coordinates of a fluid domain, vertex coordinates of a wall boundary, vertex coordinates of an inlet boundary, vertex coordinates of an outlet boundary, fluid particle spacing, wall particle spacing, inlet particle spacing, outlet particle spacing, total time, and time step; generating fluid particles, wall particles, and inlet particles; listing the discrete format of a control equation describing the flow of the weakly compressible fluid and a boundary condition application method and performing iterative calculation; and outputting results: outputting the position, velocity, pressure, and density of all fluid particles, the position, velocity, pressure, and density of all wall particles, and the position, velocity, pressure, and density of all inlet particles.

Description

Technical Field

[0001] The invention relates to the technical field of SPH fluid numerical simulation, in particular to a method for simulating the flow of a weakly compressible fluid with a velocity inlet by using SPH. Background Art

[0002] Smooth Particle Hydrodynamics (SPH) is a numerical method and an interpolation method. It was first introduced into the field of astrophysics[1][2] and is now widely used in fluid numerical simulation. SPH fluid numerical simulation is based on the SPH numerical method to solve mathematical models. The continuous mathematical model is discretized in time and space, so that the changing trend of various physical properties in the fluid motion process can be automatically solved by computer based on the initial parameter conditions.

[0003] In the SPH numerical simulation process, the area of ​​the fluid flow problem being studied is called the problem domain. In the case of a large problem domain, the critical fluid flow problem in a small area in the problem domain can be studied. This small area is called the computational domain. For example, in the problem of pipeline fluid flow, usually only the changes in the fluid flow state in the curved part of the pipeline are concerned. The entire pipeline area is the problem domain, and the curved part of the pipeline is the computational domain. The size of the computational domain directly determines the amount of calculation, so the computational domain cannot be expanded indefinitely and needs to be controlled within a certain range. The fluid flows into the computational domain from the inlet boundary and flows out of the computational domain from the outlet boundary. The inlet boundary and the outlet boundary are also the intersection of the computational domain and the problem domain.

[0004] SPH numerical simulation regards all substances as particles. For example, the fluid is composed of fluid particles, the inlet boundary is composed of inlet particles, the wall boundary is composed of wall particles, and the outlet boundary is composed of outlet particles. Therefore, the velocity inlet boundary (referred to as "velocity inlet") here refers to the velocity of the given inlet particles. For fluid particles, there are usually four physical properties that represent their flow state: position, velocity, density and pressure. As we all know, the position of fluid particles can be calculated by the velocity of fluid particles. Therefore, the physical properties that must be determined to determine the flow state of fluid particles are fluid particle velocity, fluid particle density and fluid particle pressure. Inlet particles, wall particles, and outlet particles can all be regarded as special fluid particles, and they also need to have the physical properties of fluid particles. Therefore, when realizing the velocity inlet, just like determining the flow state of fluid particles, it is necessary to determine the inlet particle pressure and inlet particle density based on the inlet particle velocity. In the weakly compressible fluid flow simulation, the inlet particle density will change slightly. There is an equation relationship between the inlet particle density and the inlet particle pressure. As long as the inlet particle pressure is determined, the inlet particle density can be determined. However, there is no direct equation relationship between the inlet particle pressure and the inlet particle velocity. Therefore, it is difficult to determine the inlet particle pressure when realizing the velocity inlet.

[0005] However, determining a reasonable inlet particle pressure based on the inlet particle velocity is the key to ensuring the correctness of the simulation results. For this problem, the methods previously used include using the pressure of the outlet particles as the pressure of the inlet particles, but this only applies to the case where the outlet boundary and the inlet boundary are completely consistent [3]; there is a velocity inlet method without processing, but non-physical pressure waves will continue to be generated inside the fluid [4]; there is a velocity inlet method based on the average value, which ignores the differences between the pressures of all inlet particles before the end of the simulation [5]; there is a SPH interpolation method that transfers all physical properties of fluid particles to the inlet particles [6], that is, the physical properties of the inlet particles are the weighted superposition of the physical properties of each fluid particle in the local area around the inlet particles. This local area is called the support domain. When the support domain is complete, that is, the support domain is full of fluid particles, the physical properties of the inlet particles obtained by this method [6] can reach second-order approximation accuracy. However, the support domain is truncated by the inlet boundary and is not filled with fluid particles, which will cause certain errors. Another method is to use a mirroring method to mirror the inlet particles into the fluid domain in the form of virtual nodes, and then use the weighted superposition value of the pressure of each fluid particle in the local area around the virtual node as the pressure of the inlet particle [7]. This method avoids the error caused by the support domain being truncated by the inlet boundary, but it is difficult to implement on complex boundaries and has low scalability. In response to the above problems, the present invention proposes a new method for simulating the flow of weakly compressible fluids with velocity inlets using SPH.

[0006] References:

[0007] [1]Lucy L BA numerical approach to the testing of the fissionhypothesis.The Astrophysical Journal,1977,82(12):1013-1024.

[0008] [2]Gingold RA, Monaghan J J. Smoothed particle hydrodynamics: theory and application to non-spherical stars. Monthly Notices of the Royal Astronomical Society, 1977, 181: 375-389.

[0009] [3] Morris JP, Patrick JF, Zhu Y. Modeling low Reynolds number incompressible flows using SPH. Journal of Computational Physics, 1997, 136 (1): 214-226.

[0010] [4]Hou Q,Kruisbrink ACH,Pearce FR,Tijsseling AS,Yue T.Smoothedparticle hydrodynamics simulations of flow separation at bends.Computers&Fluids,2014,90:138-146.

[0011] [5] Wang Zhe. Research on sound wave propagation under moving boundary conditions based on SPH numerical simulation[D]. Tianjin University.

[0012] [6]Carlos EAR,Jaime K,Leonardo DGS,JoséM D,Eduardo DL CS.Nonreflecting outlet boundary conditions for incompressible flows usingSPH.Computers&Fluids,2017,159:177-188.

[0013] [7] Negi P, Ramachandran P, Haftu A. An improved non-reflecting outletboundary condition for weakly-compressible SPH. Computer Methods in Applied Mechanics and Engineering, 2020,367. Summary of the invention

[0014] The purpose of the present invention is to propose a new method for simulating the flow of weakly compressible fluids with velocity inlets using SPH. The method of the present invention can simulate the flow of weakly compressible fluids, can use the inlet boundary and the outlet boundary to achieve continuous inflow and outflow of fluids in the calculation domain, and can obtain the wall boundary pressure distribution that conforms to the actual situation in the weakly compressible fluid with velocity inlets simulated by SPH, which can guide engineering practice, such as the material selection problem of pipeline walls. It has the characteristics of small error and wide application range. The technical solution of the present invention is as follows:

[0015] A method for simulating weakly compressible fluid flow with a velocity inlet using SPH, comprising the following steps:

[0016] Step 1, initializing parameter information, the parameter information includes vertex coordinates of the computational domain, vertex coordinates of the fluid domain, vertex coordinates of the wall boundary, vertex coordinates of the inlet boundary, vertex coordinates of the outlet boundary, fluid particle spacing, wall particle spacing, inlet particle spacing, outlet particle spacing, total time, and time step;

[0017] Step 2: Generate particles as follows:

[0018] 1) Generate fluid particles, evenly distribute the fluid particles in the fluid domain, and add initial fluid particle information, the initial fluid particle information includes the position of the fluid particles, the fluid particle pressure, the fluid particle x-direction velocity, the fluid particle y-direction velocity, the fluid particle density, the fluid particle mass, the fluid particle sound velocity, and the fluid particle type identification number;

[0019] 2) Generate wall particles, evenly distribute the wall particles on the wall boundary, and add initial wall particle information, the initial wall particle information includes wall particle position, wall particle pressure, wall particle x-direction velocity, wall particle y-direction velocity, wall particle density, wall particle mass, wall particle sound velocity, and wall particle type identification number;

[0020] 3) Generate inlet particles, evenly distribute the inlet particles at the inlet boundary, and add initial inlet particle information, which includes inlet particle position, inlet particle pressure, inlet particle x-direction velocity, inlet particle y-direction velocity, inlet particle density, inlet particle mass, inlet particle sound velocity, and inlet particle type identification number;

[0021] Step 3: List the discrete format of the control equations and the way to apply boundary conditions and perform iterative calculations:

[0022] 1) Discrete format of the control equations

[0023] The governing equations describing the flow of weakly compressible fluids are discretized using the SPH numerical method in the following format:

[0024]

[0025]

[0026] p i =c 2 (ρ i -p0) (3)

[0027] Among them, equation (1) is the discrete continuity equation, equation (2) is the discrete momentum equation, and equation (3) is the state equation. i -x j ,h) is the kernel function of SPH, h is the SPH smooth length, dx is the particle distance; the kernel function of SPH selects the cubic spline function, and R represents x i and x j The relative distance between them, R = |x i -x j | / h, then the kernel function W(x i -x j ,h)The expression is:

[0028]

[0029] i is the current fluid particle, and its adjacent 2h range is its support domain. j is one of the fluid particles in the support domain of the current fluid particle. i and x j are the positions of the current fluid particle i and the fluid particle j in its support domain; N is the number of all fluid particles in the support domain of the current particle i; m j is the mass of fluid particle j, ρ i and ρ j are the densities of the current fluid particle i and the fluid particle j in its support domain, p i and p j are the pressures of the current fluid particle i and the fluid particle j in the support domain, vi and v j are the velocities of the current fluid particle i and the fluid particle j in its support domain; t is time; c is the speed of sound; ρ0 is the standard density of the fluid; ρ i is the density of the current fluid particle i; Π ij The calculation expression is as follows:

[0030]

[0031] Among them, α=0.01~0.03, β=1~2; is the average of the sound velocity of the current fluid particle i and the sound velocity of the fluid particle j in its support domain; is the average of the density of the current fluid particle i and the density of the fluid particle j in its support domain; is the average of the smooth length of the current fluid particle i and the smooth length of the fluid particle j in its support domain;

[0032] 2) Boundary condition application method: Free slip boundary conditions are applied to the wall boundary, and the velocity of the wall particles is obtained by the following formula:

[0033]

[0034] Where n is the normal direction along the wall boundary, τ is the tangent direction along the wall boundary, and v w is the velocity of the current wall particle, and its adjacent 2h range is its support domain, v f is the velocity of the fluid particles in the current wall particle support domain, ΔV f =m f / ρ f is the volume of the fluid particles in the current wall particle support domain, where m f is the mass of the fluid particles in the current wall particle support domain, ρ f is the density of the fluid particles in the current wall particle support domain, W wf is the kernel function of SPH, W wf =W(R,h), R = |x w -x f | / h,x w is the current position of the wall particle, x f The position of the fluid particles in the current wall particle support domain, the pressure p of the wall particles w Obtained through:

[0035]

[0036] Among them, p w is the current pressure of the particle on the wall, and its adjacent 2h range is its support domain, p fis the pressure of the fluid particles in the current wall particle support domain;

[0037] The free outflow boundary condition is applied to the outlet boundary: when the fluid particles flow into the outlet domain, their fluid particle type identification number is changed to the outlet particle type identification number, and they become new outlet particles; the velocity and pressure of the outlet particles remain unchanged;

[0038] Velocity conditions are imposed on the inlet boundary: The pressure of the inlet particle is obtained as follows:

[0039]

[0040] Among them, p in is the pressure of the inlet particle, and its surrounding 2h range is its support domain, p F is the pressure of the fluid particles in the current inlet particle support domain; ΔV F =m F / ρ F is the volume of the fluid particle within the current inlet particle support domain, where m F is the mass of the fluid particles in the current inlet particle support domain, ρ F is the density of fluid particles in the current inlet particle support domain, W in F is the kernel function of SPH, W in F =W(R,h)(R=|x in -x F | / h), x in is the current position of the inlet particle, x F The position of the fluid particle within the current inlet particle support domain;

[0041] 3) Iterate the calculation, and when the total time is up, end the loop;

[0042] Step 4, output results: output the position, velocity, pressure, density of all fluid particles, the position, velocity, pressure, density of all wall particles, and the position, velocity, pressure, density of all inlet particles.

[0043] Furthermore, in the initialization parameter information described in step 1:

[0044] 1) The fluid particle spacing, inlet particle spacing and wall particle spacing are the same;

[0045] 2) The inlet boundary, outlet boundary and wall boundary are all located outside the fluid domain and are fixed areas. The vertex coordinates of the inlet boundary, outlet boundary and wall boundary will not change.

[0046] Furthermore, the particles generated in step 2:

[0047] 1) The number of layers of inlet particles and wall particles is not less than two;

[0048] 2) The sound velocity of fluid particles, the sound velocity of inlet particles and the sound velocity of wall particles are all the same.

[0049] Furthermore, the velocity condition is applied to the inlet boundary described in step 3 2):

[0050] 1) Obtain the pressure of the current inlet particle through the interpolation method described by equation (8).

[0051] 2) The pressure of the inlet particle is determined by the pressure of the fluid particles in the support domain of the current inlet particle.

[0052] Furthermore, the inlet particles are supplemented at the corresponding position of the outermost layer of the inlet boundary in the calculation process ⑦ in step three 3): when the inlet particle flows into the fluid domain, a new fluid particle is generated in the fluid domain, and the new fluid particle inherits the position, velocity, pressure, density, mass, and sound speed properties of the original inlet particle; the coordinate value of the inlet particle flowing into the fluid domain in the normal direction of the inlet boundary is changed to the coordinate value in the normal direction of the inlet boundary minus the length of the inlet boundary.

[0053] Furthermore, the iterative calculation process is as follows:

[0054] ① Determine whether the current time has reached the total time; if yes, end the calculation; if not, repeat the calculation process ② to ⑧;

[0055] ② For each wall particle, the velocity and pressure of each wall particle are updated according to equations (6) and (7) using the fluid particles searched by each wall particle in the support domain; for each inlet particle, the pressure of each inlet particle is updated according to equation (8) using the fluid particles searched by each inlet particle in the support domain;

[0056] ③ Each fluid particle takes turns to participate in the calculation as the current fluid particle. The velocity derivative of the current fluid particle is calculated according to equation (2) using other fluid particles searched in the support domain by the current fluid particle. The velocity derivative of the current fluid particle is used to update the velocity of the current fluid particle after 1 / 2 time step.

[0057] ④ Update the position of each fluid particle after 1 / 2 time step according to the speed of each fluid particle; for the fluid particles entering the outlet domain after the update, change their fluid particle type identification number to the outlet particle type identification number, and become new outlet particles, and the speed and pressure of the outlet particles remain unchanged; and repeat the calculation process ②;

[0058] ⑤ Each fluid particle takes turns to participate in the calculation as the current fluid particle. The density derivative of the current fluid particle is calculated according to equation (1) using other fluid particles searched in the support domain by the current fluid particle. The density of the current fluid particle is updated after one time step using the density derivative of the current fluid particle. The pressure of the current fluid particle is updated according to equation (3) using the density of the current fluid particle.

[0059] ⑥Continue to update the positions of all fluid particles after 1 / 2 time step according to the speed of all fluid particles. For the fluid particles entering the outlet domain after the update, change their fluid particle type identification number to the outlet particle type identification number and become new outlet particles. The speed and pressure of the outlet particles remain unchanged, and repeat the calculation process ③;

[0060] ⑦ Update the position of each inlet particle according to its velocity. If an inlet particle enters the fluid domain after the update, the inlet particle is added at the corresponding position of the outermost layer of the inlet boundary.

[0061] ⑧Current time = current time + 1 time step, repeat calculation process ①.

[0062] Furthermore, it is characterized in that h is set to 1.3dx.

[0063] Furthermore, the speed of sound c is 10 to 15 v max , v max is the maximum velocity of the fluid particles during the entire simulation.

[0064] The present invention relates to the technical field of SPH fluid numerical simulation, and proposes a new method for simulating the flow of weakly compressible fluids with velocity inlets using SPH. The present invention is mainly applied to engineering problems related to fluid flow phenomena and can guide practice. It is achieved by the following steps: Step 1, system initialization: initialization parameter information; Step 2, particle generation: including generation of fluid particles, generation of wall particles and generation of inlet particles; Step 3, listing the discrete format of the control equation and the boundary condition application method and iterative calculation: using the SPH method to discretize the control equation that describes the flow of weakly compressible fluids, using the interpolation method to determine the pressure of the inlet particles, and finally iterative calculation starting from time 0 based on the initial parameter conditions; Step 4, outputting the results: including the position, velocity, pressure and density of all types of particles at the final moment. The characteristics of the method of the present invention are: 1) no non-physical pressure wave is generated inside the fluid; 2) the error caused by the truncation of the support domain by the inlet boundary is reduced; 3) it has strong scalability for complex boundaries. BRIEF DESCRIPTION OF THE DRAWINGS

[0065] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0066] Figure 1 Program flow chart

[0067] Figure 2 Physical Model of Curved Pipe Problem

[0068] Figure 3 The results of the new method of the present invention using SPH to simulate the flow of weakly compressible fluids with velocity inlets are compared with the partial exact solution: the downstream fluid width in the exact solution is 0.5255m, while the downstream fluid width in the present invention is 0.5298m.

[0069] Figure 4 Comparison of pressure distribution results on wall A′D: including the results of the new method of simulating weakly compressible fluid flow with velocity inlet using SPH in the present invention (VI (interpolation)); the results of Fluent commercial software (Fluent); and the results of the velocity inlet method without any processing (VI (do-nothing)) [4]. DETAILED DESCRIPTION

[0070] The present invention mainly studies the inlet boundary, especially the velocity inlet boundary. The specific implementation method is illustrated by the curved pipe problem, which is described as follows: when the fluid flows through a bend or a channel where the area suddenly expands or contracts, the velocity of the fluid and the direction of movement of the fluid will change dramatically, causing flow separation. The curved pipe problem is a complex problem involving the treatment of wall boundaries, inlet boundaries and outlet boundaries. The physical model of the curved pipe problem can be found in Figure 2 . The length of the wall boundaries A′D and DB′ is 2m, and the length of the wall boundaries AC and CB is 3m. The length of the inlet boundary AA′ is 1m, and the length of the outlet boundary BB′ is 1m. The area surrounded by ACBB′A′ is the area where the fluid is located, referred to as the "fluid domain". The present invention is further described in detail below in conjunction with the accompanying drawings. Figure 1 The program flow chart of the proposed method is given. It includes:

[0071] Step 1: System initialization: Initialize parameter information. Specific parameter information includes

[0072] The vertex coordinates of the computational domain are (-0.2, -1), (-0.2, 4), (5.5, -1), (5.5, 4); the vertex coordinates of the fluid domain are (0, 0), (0, 1), (2.5, 0), (2.5, 1); the vertex coordinates of the wall boundary A′D are (-0.015, 0.999), (2.0, 0.999), (2.0, 1.015 ), (-0.015, 1.015); the vertex coordinates of the wall boundary AC are (-0.015, 0.001), (3.0, 0.001), (3.0, -0.015), (-0.015, -0.015); the vertex coordinates of the wall boundary CB are (2.999, -0.015), (2.999, 3.0), (3.015, 3. 0), (3.015,-0.015); the vertex coordinates of the wall boundary DB′ are (1.999,0.999), (1.999,3.0), (2.015,3.0), (2.015,0.999); the vertex coordinates of the inlet boundary are (0,0), (0,0.999), (-0.03,0.999), (-0.03,0); the vertex coordinates of the outlet boundary are (2.015,3), (2.015,3.015), (3.0,3), (3.0,3.015), the fluid particle spacing is 0.01m, the wall particle spacing is 0.01m, the inlet particle spacing is 0.01m, the outlet particle spacing is 0.01m, the total time is 8s, and the length of one time step (time step) is 5×10- 6 s.

[0073] Step 2: Generate particles, including:

[0074] 1) Generate fluid particles, evenly distribute fluid particles in the fluid domain, and add initial fluid particle information. Evenly distribute fluid particles in the fluid domain ACBB′A′ with the fluid particle spacing as the unit. Fluid particle pressure, fluid particle x-direction velocity, fluid particle y-direction velocity, fluid particle density, fluid particle mass, fluid particle sound velocity, and fluid particle type identification number are 0, 0, 0, and 1000 kg / m respectively. 3 , 0.1kg, 25m / s, 1;

[0075] 2) Generate wall particles, evenly distribute them on the wall boundary, and add initial wall particle information. Arrange two layers of wall particles on wall A'D, wall DB', wall AC, and wall CB, respectively, with the wall particle spacing as the unit. The wall particle pressure, wall particle x-direction velocity, wall particle y-direction velocity, wall particle density, wall particle mass, wall particle sound velocity, and wall particle type identification number are 0, 0, 0, and 1000 kg / m, respectively. 3 , 0.1kg, 25m / s, 0;

[0076] 3) Generate inlet particles, evenly distribute inlet particles at the inlet boundary, and add initial inlet particle information. Arrange three layers of inlet particles at the inlet boundary AA′ with the wall particle spacing as the unit. The inlet particle pressure, inlet particle x-direction velocity, inlet particle y-direction velocity, inlet particle density, inlet particle mass, inlet particle sound velocity, and inlet particle type identification number are 0, 1m / s, 0, 1000kg / m 3 , 0.1kg, 25m / s, 0;

[0077] Step 3: List the discrete format of the control equations and the way to apply boundary conditions and perform iterative calculations:

[0078] 1) Discrete format of the control equation. The control equation describing the flow of weakly compressible fluid is discretized using the SPH numerical method in the following format:

[0079]

[0080]

[0081] p i =c 2 (ρ i -p0) (3)

[0082] Wherein, (1) is the discrete continuity equation (hereinafter referred to as the “continuity equation”), (2) is the discrete momentum equation (hereinafter referred to as the “momentum equation”), and (3) is the state equation. i -x j ,h) is the kernel function of SPH (hereinafter referred to as "kernel function"), h is the SPH smoothing length, which is 1.3dx, dx is the particle spacing, and the kernel function usually selects the cubic spline function, and R represents x i and x j The relative distance between them, R = |x i -x j | / h, then the kernel function W(x i -x j ,h)The expression is:

[0083]

[0084] i is the current fluid particle, and its adjacent 2h range is its support domain. j is one of the fluid particles in the support domain of the current fluid particle. i and x j are the positions of the current fluid particle i and the fluid particle j in its support domain. N is the number of all fluid particles in the support domain of the current particle i. j is the mass of fluid particle j, ρ i and ρj are the densities of the current fluid particle i and the fluid particle j in its support domain, p i and p j are the pressures of the current fluid particle i and the fluid particle j in the support domain, v i and v j are the velocities of the current fluid particle i and the fluid particle j in its support domain. t is time. c is the speed of sound, which is 25 m / s. ρ0 is the standard density of the fluid, which is 1000 kg / m 3 ρ i is the density of the current fluid particle i. Π ij The calculation expression is as follows:

[0085]

[0086] Where α=0.16, β=0. It is the average of the sound velocity of the current fluid particle i and the sound velocity of the fluid particle j in its support domain. It is the average of the density of the current fluid particle i and the density of the fluid particle j in its support domain. It is the average of the smooth length of the current fluid particle i and the smooth length of the fluid particle j in its support domain.

[0087] 2) Boundary condition application method. The free slip boundary condition is applied to the wall boundary, and the velocity of the wall particle is obtained by the following formula:

[0088]

[0089] Where n is the normal direction along the wall boundary, τ is the tangent direction along the wall boundary, and v w is the velocity of the current wall particle, and its adjacent 2h range is its support domain, v f is the velocity of the fluid particles in the current wall particle support domain, ΔV f =m f / ρ f is the volume of the fluid particles in the current wall particle support domain, where m f is the mass of the fluid particles in the current wall particle support domain, ρ f is the density of the fluid particles in the current wall particle support domain, and the kernel function W wf =W(R,h)(R=|x w -x f | / h), x w is the current position of the wall particle, x f The position of the fluid particle in the current wall particle support domain is calculated by the kernel function as shown in equation (4). The pressure of the wall particle is obtained by the following method:

[0090]

[0091] Among them, p w is the current pressure of the particle on the wall, and its adjacent 2h range is its support domain, p f is the pressure of the fluid particles within the current wall particle support domain.

[0092] The free outflow boundary condition is applied to the outlet boundary. When the fluid particles flow to the outlet domain, their fluid particle type identification number is changed to the outlet particle type identification number, and they become new outlet particles. The velocity and pressure of the outlet particles remain unchanged.

[0093] The velocity condition is applied to the inlet boundary. The pressure of the inlet particle is obtained as follows:

[0094]

[0095] Among them, p in is the pressure of the inlet particle, and its surrounding 2h range is its support domain, p F is the pressure of the fluid particles in the current inlet particle support domain. ΔV F =m F / ρ F is the volume of the fluid particle within the current inlet particle support domain, where m F is the mass of the fluid particles in the current inlet particle support domain, ρ F is the density of fluid particles in the current inlet particle support domain, and the kernel function W in F =W(R,h)(R=|x in -x F | / h), x in is the current position of the inlet particle, x F The position of the fluid particle within the current inlet particle support domain, the kernel function calculation method is shown in equation (4).

[0096] 3) Iterative calculation. The current time starts from 0 and the simulation is 8s in total. The specific calculation process is as follows:

[0097] ① Determine whether the current time has reached 8s. If yes, end the calculation; if no, repeat the calculation process ② to ⑧;

[0098] ② For each wall particle, the velocity and pressure of each wall particle are updated according to equations (6) and (7) using the fluid particles searched by each wall particle in the support domain; for each inlet particle, the pressure of each inlet particle is updated according to equation (8) using the fluid particles searched by each inlet particle in the support domain.

[0099] ③ Each fluid particle takes turns to participate in the calculation as the current fluid particle. The velocity derivative of the current fluid particle is calculated according to equation (2) using other fluid particles searched in the support domain by the current fluid particle. The velocity derivative of the current fluid particle is updated according to the new velocity equal to the original velocity plus the velocity derivative multiplied by time 2.5×10 -6 The velocity of the current fluid particle after s;

[0100] ④ Update 2.5×10 with the velocity of each fluid particle according to the new position equal to the original position plus the velocity multiplied by the time -6 The position of each fluid particle after s. For the fluid particles entering the outlet domain after the update, change their fluid particle type identification number 1 to the outlet particle type identification number 5, and become new outlet particles. The outlet particle speed and outlet particle pressure remain unchanged. Repeat the calculation process ②;

[0101] ⑤ Each fluid particle takes turns to participate in the calculation as the current fluid particle. The density derivative of the current fluid particle is calculated according to equation (1) using other fluid particles searched in the support domain by the current fluid particle. The density derivative of the current fluid particle is updated according to the new density equal to the original density plus the density derivative multiplied by time 5×10 -6 The density of the current fluid particle after s is used to update the pressure of the current fluid particle according to equation (3);

[0102] ⑥ Use the velocity of all fluid particles to continue updating 2.5×10 based on the new position equal to the original position plus the velocity multiplied by the time -6 The positions of all fluid particles after s, for the fluid particles entering the outlet domain after the update, change their fluid particle type identification number 1 to outlet particle type identification number 5, becoming new outlet particles, and the outlet particle speed and outlet particle pressure remain unchanged. Then repeat the calculation process ③;

[0103] ⑦ Update 5×10 with the velocity of each entrance particle according to the new position equal to the original position plus the velocity multiplied by time -6 After s, if an inlet particle enters the fluid domain after the update, a new fluid particle is generated in the fluid domain, and the new fluid particle inherits the position, velocity, pressure, density, mass, speed of sound and other properties of the original inlet particle flowing into the fluid domain. Then, the horizontal coordinate value in the position of the inlet particle flowing into the fluid domain is changed to the horizontal coordinate value in the position of the inlet particle flowing into the fluid domain minus 0.03m, so as to supplement the inlet particle.

[0104] ⑧Current time = Current time + 5 × 10 -6 s. Repeat the calculation process ①.

[0105] Step 4: Output the results. When the total time is up, end the loop and output the position, velocity, pressure, density of all fluid particles, the position, velocity, pressure, density of all wall particles, and the position, velocity, pressure, density of all inlet particles to a file as the final result.

[0106] Further, the system initialization described in step 1:

[0107] 1) The fluid particle spacing, inlet particle spacing and wall particle spacing are all 0.01m;

[0108] 2) The inlet boundary, outlet boundary and wall boundary are all located outside the fluid domain and are fixed areas. The vertex coordinates of the inlet boundary, outlet boundary and wall boundary will not change.

[0109] Furthermore, the generated particles described in step 2:

[0110] 1) There are three layers of inlet particles and two layers of wall particles;

[0111] 2) The sound velocity of fluid particles, the sound velocity of inlet particles and the sound velocity of wall particles are all 25m / s.

[0112] Furthermore, the velocity condition is applied to the inlet boundary described in step 3 2):

[0113] 1) Obtain the pressure of the current inlet particle through the interpolation method described by equation (8).

[0114] 2) The pressure of the inlet particle is determined by the pressure of the fluid particles in the support domain of the current inlet particle.

[0115] Furthermore, the inlet particle is supplemented at the corresponding position of the outermost layer of the inlet boundary in the calculation process ⑦ in step 3): when the abscissa of the inlet particle position is greater than 0, a new fluid particle will be generated in the fluid domain, and the new fluid particle inherits the position, velocity, pressure, density, mass, sound speed and other properties of the original inlet particle with abscissa greater than 0. Then, the abscissa value of the inlet particle position with abscissa greater than 0 is changed to the abscissa value of the inlet particle position flowing into the fluid domain minus 0.03m.

[0116] Figure 3 It shows that the results of the new method of simulating the flow of weakly compressible fluid with velocity inlet by SPH of the present invention are consistent with the exact solution in position distribution. In addition, the downstream fluid width of the exact solution is 0.5255m, and the downstream fluid width of the present invention is 0.5298m, with an error of no more than 1%. Figure 4is the pressure distribution result of the wall A′D. Among them, VI (interpolation) is the result of the new method of simulating the flow of weakly compressible fluid with velocity inlet using SPH of the present invention, and the others are the result of Fluent commercial software (Fluent) and the result of the velocity inlet method without any processing (VI (do-nothing)) [4], which are used for comparison with the present invention. It can be seen that the result of the new method of simulating the flow of weakly compressible fluid with velocity inlet using SPH of the present invention eliminates the non-physical pressure waves and obtains the same result as the current commercial software result (Fluent) in terms of pressure distribution. This shows the effectiveness of the new method of simulating the flow of weakly compressible fluid with velocity inlet using SPH proposed by the present invention.

[0117] The present invention is mainly applied to engineering problems involving inlets and can guide practice. For example, when applied to curved pipe problems, suitable materials can be selected to build pipes according to the maximum pressure at the pipe wall in the simulation results to prevent pipe rupture and avoid material waste; when applied to aerodynamic acoustic simulation, it can guide sound generation; when applied to jet simulation with gravity, the corresponding relationship between the inlet pressure of the fountain and the height of the fountain can be obtained, thereby controlling the height of each fountain rising and presenting a beautiful musical fountain performance.

[0118] Although the present invention has been described above in conjunction with the drawings, the present invention is not limited to the above-mentioned specific embodiments, which are merely illustrative rather than restrictive. Under the guidance of the present invention, ordinary technicians in this field can make many modifications without departing from the purpose of the present invention, all of which are within the protection of the present invention.

Claims

1. A method for simulating weakly compressible fluid flow with a velocity inlet using SPH, comprising the following steps: Step 1, initializing parameter information, the parameter information includes vertex coordinates of the computational domain, vertex coordinates of the fluid domain, vertex coordinates of the wall boundary, vertex coordinates of the inlet boundary, vertex coordinates of the outlet boundary, fluid particle spacing, wall particle spacing, inlet particle spacing, outlet particle spacing, total time, and time step; Step 2: Generate particles as follows: 1) Generate fluid particles, evenly distribute the fluid particles in the fluid domain, and add initial fluid particle information, the initial fluid particle information includes the position of the fluid particles, the fluid particle pressure, the fluid particle x-direction velocity, the fluid particle y-direction velocity, the fluid particle density, the fluid particle mass, the fluid particle sound velocity, and the fluid particle type identification number; 2) Generate wall particles, evenly distribute the wall particles on the wall boundary, and add initial wall particle information, the initial wall particle information includes wall particle position, wall particle pressure, wall particle x-direction velocity, wall particle y-direction velocity, wall particle density, wall particle mass, wall particle sound velocity, and wall particle type identification number; 3) Generate inlet particles, evenly distribute the inlet particles at the inlet boundary, and add initial inlet particle information, which includes inlet particle position, inlet particle pressure, inlet particle x-direction velocity, inlet particle y-direction velocity, inlet particle density, inlet particle mass, inlet particle sound velocity, and inlet particle type identification number; Step 3: List the discrete format of the control equations and the way to apply boundary conditions and perform iterative calculations: 1) Discrete format of the control equations The governing equations describing the flow of weakly compressible fluids are discretized using the SPH numerical method in the following format: Among them, equation (1) is the discrete continuity equation, equation (2) is the discrete momentum equation, and equation (3) is the state equation. i -x j ,h) is the kernel function of SPH, h is the SPH smooth length, dx is the particle distance; the kernel function of SPH selects the cubic spline function, and R represents x i and x j The relative distance between them, R = |x i -x j | / h, then the kernel function W(x i -x j ,h)The expression is: i is the current fluid particle, and its adjacent 2h range is its support domain. j is one of the fluid particles in the support domain of the current fluid particle. i and x j are the positions of the current fluid particle i and the fluid particle j in its support domain; N is the number of all fluid particles in the support domain of the current particle i; m j is the mass of fluid particle j, ρ i and ρ j are the densities of the current fluid particle i and the fluid particle j in its support domain, p i and p j are the pressures of the current fluid particle i and the fluid particle j in the support domain, v i and v j are the velocities of the current fluid particle i and the fluid particle j in its support domain; t is time; c is the speed of sound; ρ0 is the standard density of the fluid; ρ i is the density of the current fluid particle i; Π ij The calculation expression is as follows: Among them, α=0.01~0.03, β=1~2; is the average of the sound velocity of the current fluid particle i and the sound velocity of the fluid particle j in its support domain; is the average of the density of the current fluid particle i and the density of the fluid particle j in its support domain; is the average of the smooth length of the current fluid particle i and the smooth length of the fluid particle j in its support domain; 2) Boundary condition application method: Free slip boundary conditions are applied to the wall boundary, and the velocity of the wall particles is obtained by the following formula: Where n is the normal direction along the wall boundary, τ is the tangent direction along the wall boundary, and v w is the velocity of the current wall particle, and its adjacent 2h range is its support domain, v f is the velocity of the fluid particles in the current wall particle support domain, ΔV f =m f / ρ f is the volume of the fluid particles in the current wall particle support domain, where m f is the mass of the fluid particles in the current wall particle support domain, ρ f is the density of the fluid particles in the current wall particle support domain, W wf is the kernel function of SPH, W wf =W(R,h), R = |x w -x f | / h,x w is the current position of the wall particle, x f The position of the fluid particles in the current wall particle support domain, the pressure p of the wall particles w Obtained through: Among them, p w is the current pressure of the particle on the wall, and its adjacent 2h range is its support domain, p f is the pressure of the fluid particles in the current wall particle support domain; The free outflow boundary condition is applied to the outlet boundary: when the fluid particles flow into the outlet domain, their fluid particle type identification number is changed to the outlet particle type identification number, and they become new outlet particles; the velocity and pressure of the outlet particles remain unchanged; Velocity conditions are imposed on the inlet boundary: The pressure of the inlet particle is obtained as follows: Among them, p in is the pressure of the inlet particle, and its surrounding 2h range is its support domain, p F is the pressure of the fluid particles in the current inlet particle support domain; ΔV F =m F / ρ F is the volume of the fluid particle within the current inlet particle support domain, where m F is the mass of the fluid particles in the current inlet particle support domain, ρ F is the density of fluid particles in the current inlet particle support domain, W inF is the kernel function of SPH, W inF =W(R,h)(R=|x in -x F | / h), x in is the current position of the inlet particle, x F The position of the fluid particle within the current inlet particle support domain; 3) Iterative calculation, when the total time is up, the loop ends. The iterative calculation process is as follows: ① Determine whether the current time has reached the total time; if yes, end the calculation; if not, repeat the calculation process ② to ⑧; ② For each wall particle, use the fluid particles searched by each wall particle in the support domain according to equations (6) and (7) Update the velocity of each wall particle and the pressure of each wall particle; for each inlet particle, update the pressure of each inlet particle according to equation (8) using the fluid particles searched by each inlet particle in the support domain; ③ Each fluid particle takes turns to participate in the calculation as the current fluid particle. The velocity derivative of the current fluid particle is calculated according to equation (2) using other fluid particles searched in the support domain by the current fluid particle. The velocity derivative of the current fluid particle is used to update the velocity of the current fluid particle after 1 / 2 time step. ④ Update the position of each fluid particle after 1 / 2 time step according to the speed of each fluid particle; for the fluid particles entering the outlet domain after the update, change their fluid particle type identification number to the outlet particle type identification number, and become new outlet particles, and the speed and pressure of the outlet particles remain unchanged; and repeat the calculation process ②; ⑤ Each fluid particle takes turns to participate in the calculation as the current fluid particle. The density derivative of the current fluid particle is calculated according to equation (1) using other fluid particles searched in the support domain by the current fluid particle. The density of the current fluid particle is updated after one time step using the density derivative of the current fluid particle. The pressure of the current fluid particle is updated according to equation (3) using the density of the current fluid particle. ⑥Continue to update the positions of all fluid particles after 1 / 2 time step according to the speed of all fluid particles. For the fluid particles entering the outlet domain after the update, change their fluid particle type identification number to the outlet particle type identification number and become new outlet particles. The speed and pressure of the outlet particles remain unchanged, and repeat the calculation process ③; ⑦ Update the position of each inlet particle according to its velocity. If an inlet particle enters the fluid domain after the update, the inlet particle is added at the corresponding position of the outermost layer of the inlet boundary. ⑧Current time = current time + 1 time step, repeat calculation process ①; Step 4, output results: output the position, velocity, pressure, density of all fluid particles, the position, velocity, pressure, density of all wall particles, and the position, velocity, pressure, density of all inlet particles.

2. The method for simulating weakly compressible fluid flow with velocity inlet using SPH as claimed in claim 1, characterized in that: In the initialization parameter information described in step 1: 1) The fluid particle spacing, inlet particle spacing and wall particle spacing are the same; 2) The inlet boundary, outlet boundary and wall boundary are all located outside the fluid domain and are fixed areas. The vertex coordinates of the inlet boundary, outlet boundary and wall boundary will not change.

3. The method for simulating the flow of a weakly compressible fluid with a velocity inlet using SPH as claimed in claim 1, characterized in that: Particles generated in step 2: 1) The number of layers of inlet particles and wall particles is not less than two; 2) The sound velocity of fluid particles, the sound velocity of inlet particles and the sound velocity of wall particles are all the same.

4. The method for simulating weakly compressible fluid flow with velocity inlet using SPH as claimed in claim 1, characterized in that: The velocity condition applied to the inlet boundary described in step 3 2) is: 1) Obtain the pressure of the current inlet particle by the interpolation method described by equation (8); 2) The pressure of the inlet particle is determined by the pressure of the fluid particles in the support domain of the current inlet particle.

5. The method for simulating weakly compressible fluid flow with velocity inlet using SPH as claimed in claim 1, characterized in that: In step three (3), the inlet particles are supplemented at the corresponding position of the outermost layer of the inlet boundary as described in calculation process (7): when the inlet particle flows into the fluid domain, a new fluid particle is generated in the fluid domain, and the new fluid particle inherits the position, velocity, pressure, density, mass, and sound speed properties of the original inlet particle; the coordinate value of the inlet particle flowing into the fluid domain in the normal direction of the inlet boundary is changed to the coordinate value in the normal direction of the inlet boundary minus the length of the inlet boundary.

6. The method for simulating weakly compressible fluid flow with velocity inlet using SPH as claimed in claim 1, characterized in that: h is set to 1.3dx.

7. The method for simulating weakly compressible fluid flow with velocity inlet using SPH as claimed in claim 1, characterized in that: The speed of sound c is 10~15v max , v max is the maximum velocity of the fluid particles during the entire simulation.

Images