Lattice Boltzmann flux method suitable for super-large-density-ratio axisymmetric multiphase flow

Abstract

The invention relates to a lattice Boltzmann flux method suitable for super-large-density-ratio axisymmetric multiphase flow, which comprises the following steps of: S1, establishing a plurality of continuous control body unit grids under a polar coordinate system for a calculation region of a multiphase fluid mechanics problem with an axisymmetric characteristic; s2, setting initial conditions, and determining an overall time step length of simulation; s3, the time step length of the lattice-Boltzmann flow step and the relaxation time of the lattice-Boltzmann evolution equation and the relaxation time of the Cahn-Hilliard equation are determined; s4, calculating a time iteration macroscopic physical quantity; s5, obtaining macroscopic quantities on all unit interfaces; s6, calculating the flux of adjacent interfaces of each control body unit; s7, obtaining a macroscopic physical quantity at the center of the cell of the control body at the next time step; and S8, repeating the steps S3-S7 until a preset termination condition is met. The method has stability and precision, is high in calculation efficiency, and can stably simulate the axial symmetry multiphase flow problem with the super-large density ratio and the large viscosity ratio.

Description

Technical Field

[0001] The present invention relates to the field of fluid mechanics, and more specifically, to a lattice Boltzmann flux method applicable to axisymmetric multiphase flows with extremely large density ratios. Background Art

[0002] Multiphase flow problems with axisymmetric characteristics widely exist in the engineering field, such as the collision between droplets, the impact of droplets on liquid films, and the formation of droplets in microfluidic devices. The study of these phenomena is of great significance for solving applications in specific fields. Numerical simulation methods are powerful tools for studying multiphase flow problems. Fluid systems with axisymmetric characteristics can be simplified to two-dimensional problems. Therefore, establishing an axisymmetric multiphase flow model specifically can greatly reduce the computational amount and data storage amount of axisymmetric multiphase flow simulations. In an engineering environment, there are also many cases where the density ratio and viscosity ratio between different phases of axisymmetric multiphase flows are relatively large. Such problems require numerical simulation methods to have high accuracy and stability. The accuracy and stability limit the application of many axisymmetric multiphase flow numerical simulation methods in this important field.

[0003] The axisymmetric multiphase flow model based on the lattice Boltzmann method has received much attention in many aspects due to the natural parallelism of the collision and flow steps during its evolution and the advantage of dealing with complex convection-diffusion problems compared with the traditional model based on the macroscopic Navier-Stokes equations. However, the standard lattice Boltzmann method is built on a uniform orthogonal grid, which is not very flexible and efficient for handling hydrodynamic problems that require local grid refinement. In the literature (L. Yang, C. Shu, Y. Yu, Y. Wang, G. Hou, A mass-conserved fractional step axisymmetric lattice Boltzmann flux solver for incompressible multiphase flows with large density ratio, Physics of Fluids 32 (10) (2020) 103308), an axisymmetric multiphase lattice Boltzmann flux solver was developed. This method is based on the finite volume method, can use non-uniform grids, is relatively flexible to apply, and locally applies the lattice Boltzmann method to the calculation of interface fluxes, greatly simplifying the complexity of flux calculation. However, this method has the following three problems: 1) To ensure the calculation accuracy at the multiphase interface, this method introduces a mass correction term similar to other axisymmetric lattice Boltzmann multiphase flow simulation methods. The correction steps are rather cumbersome, reducing the calculation efficiency and simulation stability. 2) This method uses a two-step prediction-correction method when simulating the flow field and multiphase interface, and the axisymmetric effect is reflected in the correction step. The two-step method is rather cumbersome to implement, affecting the calculation efficiency and introducing additional errors. 3) When calculating the interface flux, the axisymmetric effect is not taken into account at the control volume interface, resulting in inconsistent control equations at the control volume center and the interface, affecting the calculation accuracy and stability. The literature (L. Yang, X. Yang, Y. Yang, G. Hou, Y. Wang, An improved axisymmetric interfacial lattice Boltzmann flux solver for large-density-ratio multiphase flows, Physics of Fluids 36 (2024) 023338) corrected the second problem above, but the first and third problems still exist. Due to these problems, the above methods cannot stably simulate axisymmetric multiphase flows with extremely large density ratios and viscosity ratios. Summary of the Invention

[0004] The technical problem to be solved by the present invention is to provide a lattice Boltzmann flux method applicable to axisymmetric multiphase flows with an extremely large density ratio. This method takes into account the axisymmetric effect simultaneously at the center of the finite volume control volume and the control volume interface, ensuring the physical consistency of the simulation results, laying a solid theoretical foundation for the high precision and stability of the model, and removing the complex mass correction term. It has high stability and precision, high computational efficiency, and can stably simulate axisymmetric multiphase flow problems with extremely large density ratios and viscosity ratios.

[0005] The technical solution adopted by the present invention to solve its technical problems is: constructing a lattice Boltzmann flux method applicable to axisymmetric multiphase flows with an extremely large density ratio, including the following steps: S1. For the computational domain of the multiphase hydrodynamics problem with axisymmetric characteristics, establish a continuous plurality of control volume unit grids in the polar coordinate system; S2. Set the initial conditions and determine the overall time step of the simulation; S3. Determine the time step of the lattice-Boltzmann flow step and the relaxation times of the lattice-Boltzmann evolution equation and the Cahn-Hilliard equation; S4. For each node at the control volume unit interface, calculate the time-iterative macroscopic physical quantities; S5. Obtain the macroscopic quantities on all unit interfaces; S6. Calculate the fluxes at the adjacent interfaces of each control volume unit; S7. Use the three-step Runge-Kutta method to update the evolution equation of the finite volume lattice Boltzmann flux solver to obtain the macroscopic physical quantities at the center of the control volume cell at the next time step; S8. Repeat steps S3 - S7 until the preset termination conditions are met.

[0006] According to the above solution, in step S1, finite volume method is used for discretization to establish a continuous plurality of control volume unit grids in the polar coordinate system where , respectively represent the radial and axial coordinates, is the angular coordinate; The central node of each unit stores the macroscopic physical quantity parameters of the fluid, ( ); among them, is the density of the fluid, is the velocity of the fluid, and p is the pressure of the fluid point; is the order parameter, and the order parameter is the control parameter for simulating multiphase flow based on the Cahn-Hilliard equation, used to indicate whether the target unit is at the interface of multiple phases or in a certain phase.

[0007] According to the above scheme, in the step S2, the initial conditions at include the density at the center of each control cell, pressure , velocity , order parameter at the center of each control cell, or within the cell, calculate the relaxation times of the lattice-Boltzmann evolution equation and the Cahn-Hilliard equation, where is the position coordinate of the cell boundary point. For the two-dimensional axisymmetric lattice-Boltzmann method, it corresponds to the discrete velocity :

[0008] The weights in each direction are:

[0009]

[0010] where is the mobility in the Cahn-Hilliard equation, is the adjustment parameter, is the dynamic viscosity coefficient of the fluid, .

[0011] According to the above scheme, in the step S4, interpolate and calculate the macroscopic physical quantity ( ) at the position and time

[0012] using the following formula: where represents the macroscopic physical quantity to be interpolated; and represent the center point coordinates of two adjacent cells; represents the coordinate at the interface of two adjacent cells; calculate the equilibrium distribution function of the Cahn-Hilliard equation at the position and time , and the equilibrium distribution function of the evolution equation ;

[0013]

[0014] Among them, , is the chemical potential in the Cahn-Hilliard equation and is solved by the following formula:

[0015]

[0016] Among them, is the surface tension coefficient, is the interface thickness.

[0017] According to the above scheme, in the step S5, the macroscopic quantities on all unit interfaces at the moment are obtained through the following formula; .

[0018] According to the above scheme, in the step S6, the flux at the control volume unit interface is calculated , ;

[0019] .

[0020] According to the above scheme, in the step S7, the evolution equation of the finite volume lattice Boltzmann flux solver is updated using the three-step Runge-Kutta method to obtain the macroscopic physical quantities at the center of the control volume cell at the next time step t + △t ( ), and the time step of the update iteration is only restricted by the CFL condition:

[0021] Among them, , , is the volume of the control volume unit , is the area of the -th unit interface surrounding , is the outward unit normal vector at the -th unit interface.

[0022] Applying the lattice Boltzmann flux method of the present invention for axisymmetric multiphase flows with ultra-high density ratios has the following beneficial effects: 1. The present invention simultaneously considers the axisymmetric effect at the center and interface of the finite volume control volume, both of which can accurately revert to the macroscopic control equations, ensuring the physical consistency of the simulation results and laying a solid theoretical foundation for the high precision and stability of axisymmetric multiphase flow simulations; 2. The present invention removes the complex mass correction term, is simple and efficient to implement, reduces the computational complexity, and improves the stability of the simulation; 3. The present invention can stably simulate axisymmetric multiphase flow problems with ultra-high density ratios and large viscosity ratios. BRIEF DESCRIPTION OF THE DRAWINGS

[0023] The present invention will be further described below in conjunction with the drawings and embodiments. In the drawings: Figure 1 is the computational flow chart of the lattice Boltzmann flux method of the present invention for axisymmetric multiphase flows with ultra-high density ratios; Figure 2 is a schematic diagram of the local application of the axisymmetric lattice-Boltzmann method on the interface of the finite volume control volume unit of the present invention; Figure 3 is the initial shape diagram of the droplet simulated by the present invention; Figure 4 is the shape diagram of the droplet oscillation evolution simulated by the present invention at time steps of (a) t = 0.0, (b) t = 12000, (c) t = 23000, and (d) t = 190000. DETAILED DESCRIPTION OF THE INVENTION

[0024] In order to have a clearer understanding of the technical features, objectives, and effects of the present invention, the specific embodiments of the present invention will now be described in detail with reference to the drawings.

[0025] As Figures 1-4 shown, the lattice Boltzmann flux method of the present invention for axisymmetric multiphase flows with ultra-high density ratios includes the following steps: S1. For the computational domain of the multiphase hydrodynamics problem with axisymmetric characteristics, use the finite volume method for discretization to establish a continuous grid of multiple control volume units in the polar coordinate system under which . respectively represent the radial and axial coordinates, is the angular coordinate. The central node of each unit stores the macroscopic physical quantity parameters of the fluid ( ), where is the density of the fluid, is the velocity of the fluid, and p is the pressure at the fluid point; is the order parameter, which is the control parameter for simulating multiphase flow based on the Cahn-Hilliard equation and is used to represent whether the target cell is at the interface of multiple phases or in a certain phase.

[0026] As Figure 2 shown, and respectively represent the central nodes of the control units and . The evolution of the fluid will be achieved through the flux exchange between the cell boundaries. Applying the axisymmetric lattice-Boltzmann method at the cell boundaries can greatly simplify the calculation of the flux at the interface and is also beneficial to the improvement of calculation efficiency, accuracy, and stability.

[0027] S2. Set the initial conditions, that is, the density , pressure , velocity , order parameter at the center of each control cell. Determine the overall time step t of the simulation. S3. Determine the time step of the lattice-Boltzmann flow step, generally taking . is the spatial step of the lattice-Boltzmann. The step size should be selected such that the position is within the cell or . The units in the lattice-Boltzmann method are generally taken as dimensionless units. Accordingly, the relaxation times and corresponding to the evolution equation of the lattice-Boltzmann and the Cahn-Hilliard equation can be calculated. Here, refers to the position coordinates of the cell boundary points. For the two-dimensional axisymmetric lattice-Boltzmann method, the D2Q9 model is selected, and the corresponding discrete velocities are:

[0028] The weights in each direction are:

[0029] and are expressed as:

[0030] Among them, is the mobility in the Cahn-Hilliard equation, is the adjustment parameter, is the dynamic viscosity coefficient of the fluid, .

[0031] S4. Interpolate and calculate the macroscopic physical quantity at position at time using the following formula ( );

[0032] where, represents the macroscopic physical quantity to be interpolated; and represent the central point coordinates of two adjacent cells; represents the coordinate at the interface of two adjacent cells.

[0033] Calculate the equilibrium distribution function of the Cahn-Hilliard equation, and the equilibrium distribution function of the evolution equation at position and time using the following formula;

[0034]

[0035] where, , is the chemical potential in the Cahn-Hilliard equation, which is solved using the following formula:

[0036]

[0037] where, is the surface tension coefficient, is the interface thickness.

[0038] S5. Obtain the macroscopic quantities on all cell interfaces at time using the following formula;

[0039] S6. Calculate the fluxes , at the interfaces of the control volume cells;

[0040]

[0041] where

[0042] where,

[0043] is the component of the surface tension; is the external force; is the Kronecker function, .

[0044] S7. Update the evolution equation of the finite volume lattice Boltzmann flux solver using the three-step Runge-Kutta method to obtain the macroscopic physical quantities at the cell center at the next time step . The time step for the update iteration is only restricted by the CFL condition.

[0045]

[0046] where , , is the volume of the control volume cell , is the area of the -th cell interface surrounding , is the outward unit normal vector at the -th cell interface. As can be seen from the above steps, the macroscopic physical quantities of the flow field at the control cell interface can be reconstructed from the physical quantities at the control cell center without storage, saving memory space.

[0047] S8. Repeat steps S3 - S7 until the preset termination condition is met.

[0048] Example The droplet oscillation problem is a dynamic basic example for verifying the multiphase flow numerical model.

[0049] S1. In the initial state, set the droplet to be elliptical. The shape of the ellipse is described by the order parameter as follows:

[0050] where and are the semi-major axis length along the radial direction and the semi-major axis length along the axial direction of the elliptical droplet, respectively. is the radius of the droplet when it finally reaches the equilibrium state. According to the mass conservation, it can be known that . The elliptical droplet represented by the above order parameter is as shown in Figure 3 .

[0051] The initial elliptical droplet is in an unstable state and will eventually become a stable circular state under the action of surface tension, viscous force, etc.

[0052] S2. Establish a uniform orthogonal grid, , with the center of the droplet at . Set the left boundary to a symmetric boundary condition and the remaining boundaries to no-slip wall boundaries. The size of the droplet is , , the mobility , the surface tension coefficient , the kinematic viscosity coefficient of the liquid phase = the kinematic viscosity coefficient of the gas phase , the thickness of the two-phase interface , and the density ratio of the liquid phase to the gas phase is taken as a relatively large value . The density of the liquid phase , and the density of the gas phase .

[0053] S3. Set at the center of each control cell, ), where the initial velocity , the density of the fluid point , the initial pressure is 0, and the order parameter is: .

[0054] Determine the overall time step t of the simulation, which is taken as 1.0 here. S4. Determine the time step of the lattice-Boltzmann flow step, and the step size is selected such that the position is within the cell or , which is taken as 0.45 here. Select the D2Q9 model and calculate the relaxation times and of the corresponding lattice-Boltzmann evolution equation and the Cahn-Hilliard equation;

[0055] Among them, , is the adjustment parameter, .

[0056] S5. Interpolate and calculate the macroscopic physical quantities ( ) at the position and time ) according to the following formula.

[0057]

[0058] Then calculate at the position and time , the equilibrium distribution function of the Cahn-Hilliard equation , and the equilibrium distribution function of the evolution equation

[0059]

[0060]

[0061] S6. The macroscopic quantities at all unit interfaces at a certain moment can be obtained through the following formula ;

[0062] S7. Next, calculate the flux at the control volume unit interface , ;

[0063]

[0064] S8. Use the three-step Runge-Kutta method to update the evolution equation of the finite volume lattice Boltzmann flux solver, and obtain the macroscopic physical quantities at the cell center at the next time step . The time step of the update iteration is only restricted by the CFL condition.

[0065]

[0066] S9. Repeat steps S3 - S8 until the preset termination condition is met.

[0067] For example Figure 4 In the evolution process of the simulated oscillating droplet problem shown, the oscillating droplet will come to rest and its shape will become a standard circle, with a theoretical radius of 79.7. The radius obtained in this embodiment is 79.3, and the relative error is only 0.5%. For the case where the two-phase density ratio is 10000 (ultra-large density ratio), the accuracy is relatively high.

[0068] The embodiments of the present invention have been described above in conjunction with the accompanying drawings. However, the present invention is not limited to the above specific embodiments. The above specific embodiments are merely illustrative and not restrictive. Under the inspiration of the present invention, those of ordinary skill in the art can also make many forms without departing from the purpose of the present invention and the scope protected by the claims. These all fall within the protection scope of the present invention.

Claims

1. A lattice Boltzmann flux method suitable for axisymmetric multiphase flow with ultra-large density ratio, characterized in that: The following steps are involved: S1. For the calculation area of ​​the multiphase fluid mechanics problem with axisymmetric characteristics, establish multiple continuous control volume unit grids in the polar coordinate system; S2, set initial conditions and determine the overall time step of the simulation; S3, determine the time step of the lattice-Boltzmann flow step and the relaxation time of the lattice-Boltzmann evolution equation and the Cahn-Hilliard equation; S4, for each node at the interface of the control volume unit, calculate the time iteration macroscopic physical quantity; S5. Obtain the macroscopic quantities on all unit interfaces; S6, calculating the flux at the adjacent interface of each control volume unit; S7, using the three-step Runge-Kutta method to update the evolution equations of the finite volume lattice Boltzmann flux solver to obtain the macroscopic physical quantities at the center of the control volume unit cell at the next time step; S8. Repeat steps S3-S7 until a preset termination condition is met.

2. The lattice Boltzmann flux method for ultra-large density ratio axisymmetric multiphase flow according to claim 1, characterized in that: In step S1, the finite volume method is used for discretization to establish a polar coordinate system. A continuous grid of multiple control volume elements , Represent the radial and axial coordinates respectively, is the angle coordinate; The central node of each unit stores the macroscopic physical quantity parameters of the fluid. );in, is the density of the fluid, is the velocity of the fluid, p is the pressure at the fluid point; It is the order parameter, which is the control parameter for simulating multiphase flow based on the Cahn-Hilliard equation. It is used to indicate whether the target unit is at the interface of multiple phases or in a certain phase.

3. The lattice Boltzmann flux method for ultra-large density ratio axisymmetric multiphase flow according to claim 2, characterized in that: In the step S2, The initial conditions include the density at the center of each control cell ,pressure ,speed , Sequence Parameters .

4. The lattice Boltzmann flux method for ultra-large density ratio axisymmetric multiphase flow according to claim 3, characterized in that: In step S3, the time step of the lattice-Boltzmann flow step is determined ,Pick , the step size selection position In the cell or Calculate the relaxation time of the lattice-Boltzmann evolution equation and the Cahn-Hilliard equation within and , is the position coordinate of the unit boundary point, which corresponds to the discrete velocity for the two-dimensional axisymmetric lattice-Boltzmann method : The weight of each direction for: and It is expressed as: in, is the mobility in the Cahn-Hilliard equation, To adjust the parameters, is the dynamic viscosity coefficient of the fluid, .

5. The lattice Boltzmann flux method for ultra-large density ratio axisymmetric multiphase flow according to claim 4, characterized in that: In step S4, the position is calculated by interpolation according to the following formula: Place, time The macroscopic physical quantity when ): in, Represents the macroscopic physical quantity to be interpolated; and Represents the coordinates of the center points of two adjacent cells; Represents the coordinates at the interface between two adjacent units; The position is calculated by the following formula Place, time , the equilibrium distribution function of the Cahn-Hilliard equation , and the equilibrium distribution function of the evolution equation ; in, , is the chemical potential in the Cahn-Hilliard equation, which can be solved using the following formula: in, is the surface tension coefficient, is the interface thickness.

6. The lattice Boltzmann flux method for ultra-large density ratio axisymmetric multiphase flow according to claim 5, characterized in that: In step S5, the following formula is used to obtain Macroscopic quantities on all unit interfaces at all times; 。 7. The lattice Boltzmann flux method for ultra-large density ratio axisymmetric multiphase flow according to claim 6, characterized in that: In step S6, the flux at the interface of the control volume unit is calculated. , ; 。 8. The lattice Boltzmann flux method for ultra-large density ratio axisymmetric multiphase flow according to claim 7, characterized in that: In step S7, the three-step Runge-Kutta method is used to update the evolution equation of the finite volume lattice Boltzmann flux solver to obtain the macroscopic physical quantity ( ), update the time step of the iteration Only subject to CFL conditions: in, , , It is the control volume unit The volume of It is surrounded No. The area of ​​the unit interface, It is The outward-pointing unit normal vector at the interface of each element.

Images