A method for optimizing the simulation efficiency of multiphase flows based on the multiphase lattice Boltzmann flux method

Abstract

The present invention discloses a method for optimizing the simulation efficiency of multiphase flow based on the multiphase lattice Boltzmann flux method, and the steps are as follows: establish a physical model corresponding to the multiphase flow problem and perform grid division on the entire computational domain; input the initial conditions and boundary conditions required for calculation according to the problem type; perform spatial discretization based on the finite volume method through the Cahn–Hilliard equation controlled by the phase field and the lattice Boltzmann flux solver controlled by the flow field, set the global time step, use the dual-time stepping method during the flow solution process, and use the second-order TVD Runge-Kutta time format to perform the advancement calculation on the inner iteration pseudo-time step; continuously calculate the local time step during the solution process and apply it to the iteration process of the dual-time stepping method until the inner iteration residual meets the convergence requirement, so as to obtain the flow field phase field numerical values of the entire computational domain at the next moment. The present invention can effectively reduce the calculation time of multiphase flow problems and improve the calculation stability.

Description

Technical Field

[0001] The present invention relates to the field of multiphase fluid calculations, and more specifically to a method for optimizing the simulation efficiency of multiphase flows based on the multiphase lattice Boltzmann flux method. Background Art

[0002] The problem of multiphase flows has always been a matter of great concern in computational fluid dynamics and other physical disciplines. It widely exists in fields such as aerospace, chemical production, and energy development. For example, it can be found in aircraft icing, oil extraction, etc. Therefore, the research on multiphase flow problems has always been the focus of the scientific community, and its related theoretical algorithms have broad development prospects.

[0003] To study multiphase flow problems and effectively simulate multiphase flows such as incompressible fluids, many numerical methods have been proposed, such as the volume of fluid method (VOF), marker and cell method (MAC), diffuse interface method (DI), level set method, and so on. Among them, in recent years, the diffuse interface method has stronger adaptability when dealing with complex multiphase flow problems and has been significantly popularized. As a method with strong adaptability, the diffuse interface method has been widely improved and applied at home and abroad. The existing diffuse interface methods are mainly divided into two categories. The first category is the continuum method based on macroscopic conservation laws, and the second category is the kinetic method based on the mesoscopic lattice Boltzmann equation. The main idea of the continuum method is to directly solve the Cahn–Hilliard equation of the phase field and the Navier–Stokes equation of the flow field, and it can relatively easily apply high-order upwind and other formats to obtain relatively accurate simulation results. The kinetic method is based on the lattice Boltzmann equation. Since it considers the collision motion of microscopic particles, the kinetic method can better simulate fluid motion from a microscopic perspective.

[0004] In recent years, based on these two methods, a multiphase lattice Boltzmann flux solver (MLBFS) has emerged. This method can combine the advantages of the two current mainstream diffuse interface methods and can relatively accurately calculate the flow field of multiphase flows. However, due to the large density and high viscosity difference of many multiphase flow problems, there are problems such as low stability and long calculation time when calculating multiphase flow problems. For the calculation of multiphase flows, a new optimization algorithm is needed to improve the calculation efficiency and stability. Summary of the Invention

[0005] The object of the present invention is to provide a method for optimizing the simulation efficiency of multiphase flows based on the multiphase lattice Boltzmann flux method in view of the current problems of the multiphase lattice Boltzmann flux algorithm. This efficiency optimization method can reduce the calculation time of multiphase flow problems and is applicable to multiphase flow problems with large density ratios, high viscosity ratios, etc.

[0006] The object of the present invention is solved by the following technical solutions:

[0007] A method for optimizing the simulation efficiency of multiphase flow based on the multiphase lattice Boltzmann flux method, characterized in that the steps of the efficiency optimization method are as follows:

[0008] S1: Establish a physical model corresponding to the multiphase flow problem and divide the entire computational domain into grids;

[0009] S2: Input the initial conditions and boundary conditions required for the calculation according to the problem type;

[0010] S3: Based on the finite volume method, perform spatial discretization through the Cahn–Hilliard equation controlled by the phase field and the lattice Boltzmann flux solver controlled by the flow field. Set the global time step, use the dual-time stepping method during the flow solution process, and use the second-order TVD Runge-Kutta time format to perform the advancement calculation on the inner iteration pseudo-time step;

[0011] S4: Continuously calculate the local time step during the solution process of step S3 and apply it to the iteration process of the dual-time stepping method until the inner iteration residual meets the preset convergence requirements, and obtain the flow field phase field numerical values of the entire computational domain at the next moment.

[0012] In step S4, in order to accelerate the convergence of the inner iteration residual, an implicit residual smoothing method is adopted to increase the local time step to accelerate the convergence speed.

[0013] The formula of the implicit residual smoothing method is:

[0014]

[0015] In formula (7): ε I and ε J are the relaxation coefficients in the I and J directions of the computational domain respectively; is the residual value recorded at the grid center (I-1,J) after the first smoothing in the I direction; is the residual value recorded at the grid center (I,J) after the first smoothing in the I direction; is the residual value recorded at the grid center (I+1,J) after the first smoothing in the I direction; is the initial residual value recorded at the grid center (I,J) before smoothing; is the final residual value recorded at the grid center (I-1,J) after the second smoothing in the J direction; is the final residual value recorded at the grid center (I,J) after the second smoothing in the J direction; is the final residual value recorded at the grid center (I+1,J) after the second smoothing in the J direction; Record the residual value at the center of the (I, J) grid after the first smoothing in the I direction.

[0016] The physical model in step S1 is a two-phase flow model.

[0017] The initial distribution of the phase field in the initial conditions in step S2 is obtained from the phase field distance sign function, and the selected phase field distance sign function is:

[0018]

[0019] In formula (1): C represents the order parameter in the phase field, that is, the volume fraction of the heavier mass fluid in the corresponding control volume; d is the vertical distance from this control volume to the interface of the multiphase fluid, and its sign is determined by the initial control volume density. For two-phase flow, the distance of the control volume at the high-density fluid is positive, and the distance of the control volume at the low-density fluid is negative; ξ is the defined interface thickness.

[0020] The boundary conditions in step S2 are correctly selected according to the research object.

[0021] The Cahn–Hilliard equation for phase field control in step S3 is:

[0022]

[0023] In formula (2): C represents the order parameter in the phase field, that is, the volume fraction of the heavier mass fluid in the corresponding control volume; t is time; is the divergence operator; u is the flow field velocity vector; Γ is the mass flow rate; μ C represents the chemical potential, and μ C is:

[0024]

[0025] In formula (3): σ is the surface tension coefficient; ξ is the interface thickness.

[0026] The lattice Boltzmann flux solver for flow field control in step S3 is:

[0027]

[0028] In formula (4): W is the macroscopic conserved quantity; t is time; is the divergence operator; F is the flow flux; S is the source term;

[0029]

[0030] In formula (5): p is the macroscopic pressure at the grid center of the grid where it is located; ρ is the macroscopic pressure at the grid center of the grid where it is located; u is the macroscopic horizontal velocity; v is the macroscopic vertical velocity; c sis the speed of sound, given by the lattice velocity model; T is the transpose of the matrix; α is the direction of particle motion represented; e α represents the particle velocity in the corresponding direction, given by the lattice velocity model; is the particle distribution function in the corresponding direction after collision; is the particle distribution function in the corresponding direction before collision; is the comprehensive distribution function; is the gradient operator; F s is the surface tension; r is the displacement in the corresponding direction; t is the time; w α is the weight coefficient in the corresponding direction, given by the lattice velocity model; η = min(λ, 0.025), λ = tanh(|C L -C R | / (2(C L +C R +0.2))), C L is the order parameter of the control volume on the left side of a certain side of the corresponding grid, C R is the order parameter of the control volume on the right side of a certain side of the corresponding grid; δ t is the given particle motion time; τ is the relaxation time.

[0031] The original formula for time marching of the unsteady flow problem is: In the formula, Ω is the control volume, and R is the residual.

[0032] The formula for the dual-time stepping method in steps S3 and S4 is:

[0033]

[0034] In formula (6): Ω is the control volume; is the macroscopic conserved quantity at the pseudo-time step; t * is the pseudo-time; R is the residual; is the macroscopic conserved quantity at the nth real-time step; is the macroscopic conserved quantity at the (n - 1)th real-time step; Δt is the global time step; R * is the final residual recorded at the grid center after implicit residual smoothing.

[0035] The present invention has the following advantages compared with the prior art:

[0036] The method for optimizing the efficiency of multiphase flow simulation provided by the present invention solves the Cahn–Hilliard equation controlled by the phase field and calculates the macroscopic physical quantities of the flow field by using the lattice Boltzmann flux solver. During the solving process, based on the dual-time-step advancement equation, local time steps and implicit residue smoothing techniques are introduced in the inner iteration process, improving the multiphase lattice Boltzmann flux method (MLBFS). This not only retains the good mesoscopic characteristics of the original method in solving multiphase flow problems but also enables the reduction of calculation time and the improvement of calculation stability when calculating multiphase flow problems with large density ratios and high viscosity ratios, enhancing the adaptability of MLBFS in the field of multiphase fluid calculations. Description of the Drawings

[0037] In order to more clearly illustrate the technical solutions in the embodiments of the present invention or in the prior art, the following will briefly introduce the drawings required for use in the description of the embodiments or the prior art. Obviously, the drawings in the following description are some embodiments of the present invention. For those of ordinary skill in the art, without creative efforts, other drawings can also be obtained based on these drawings.

[0038] Att Figure 1 It is a flowchart of a method for optimizing the efficiency of a multiphase lattice Boltzmann flux method (MLBFS) provided by an embodiment of the present invention;

[0039] Att Figure 2 It is a schematic diagram of discrete velocities of a D2Q9 lattice velocity model provided by an embodiment of the present invention;

[0040] Att Figure 3 It is a schematic diagram of a geometric model of a single-bubble rising problem provided by an embodiment of the present invention;

[0041] Att Figure 4 It is a comparison diagram of the bubble interface position and the original MLBFS results for a single-bubble rising problem based on a low reference velocity provided by an embodiment of the present invention;

[0042] Att Figure 5 It is a result diagram of the number of inner iteration steps for a single-bubble rising problem based on a low reference velocity provided by an embodiment of the present invention;

[0043] Att Figure 6 It is a comparison diagram of the bubble interface position and the original MLBFS results for a single-bubble rising problem based on a high reference velocity provided by an embodiment of the present invention;

[0044] Att Figure 7 It is a result diagram of the number of inner iteration steps for a single-bubble rising problem based on a high reference velocity provided by an embodiment of the present invention. Detailed Embodiments

[0045] In order to enable those skilled in the art to better understand the technical solutions in the present invention, the following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. Based on the embodiments in the present invention, all other embodiments obtained by those of ordinary skill in the art without creative efforts shall fall within the protection scope of the present invention.

[0046] Although the present invention provides method operation steps or device structures as shown in the following embodiments or drawings, more or fewer operation steps or module units may be included in the method or device based on routine or non-creative labor. In steps or structures where there is no necessary causal relationship logically, the execution order of these steps or the module structure of the device is not limited to the execution order or module structure shown in the embodiments or drawings of the present invention. When the described method or module structure is applied to an actual device or terminal product, it can be executed sequentially or in parallel according to the method or module structure shown in the embodiments or drawings (for example, in an environment of parallel processors or multi-threaded processing, and even including an implementation environment of distributed processing).

[0047] Figure 1 The method flow chart of an efficiency optimization method based on the multiphase lattice Boltzmann flux method (MLBFS) provided for an embodiment of the present invention. The idea of the present invention is based on the dual-time-step advancement equation, introducing local time steps and implicit residual smoothing techniques in the inner iteration process to improve the computational efficiency of MLBFS for multiphase flow problems.

[0048] The flow of a multiphase flow simulation efficiency optimization method based on the multiphase lattice Boltzmann flux method is as Figure 1 shown and includes:

[0049] S1: Establish a physical model corresponding to the multiphase flow problem and perform grid division on the entire computational domain.

[0050] S2: Input the initial conditions and boundary conditions required for the calculation according to the problem type; the initial distribution of the phase field in the initial conditions is obtained by the phase field distance sign function, and the selected phase field distance sign function is:

[0051]

[0052] In Equation (1): C represents the order parameter in the phase field, that is, the volume fraction of the heavier mass fluid corresponding to the control volume; d is the vertical distance from the control volume to the interface of the multiphase fluid, and its sign is determined by the density of the initial control volume. For two-phase flow, the distance of the control volume at the high-density fluid is positive, and the distance of the control volume at the low-density fluid is negative; ξ is the defined interface thickness.

[0053] S3: Based on the finite volume method, perform spatial discretization through the Cahn–Hilliard equation controlled by the phase field and the lattice Boltzmann flux solver controlled by the flow field. Set the global time step. During the flow solution process, use the dual-time stepping method and the second-order TVD Runge-Kutta time format to perform the advancement calculation on the pseudo-time step of the inner iteration.

[0054] When solving the phase field change using the Cahn–Hilliard equation controlled by the phase field, its equation is:

[0055]

[0056] In Equation (2): C represents the order parameter in the phase field, that is, the volume fraction of the heavier mass fluid corresponding to the control volume; t is time; is the divergence operator; u is the flow field velocity vector; Γ is the mass flow rate; μ C represents the chemical potential, and μ C is:

[0057]

[0058] In Equation (3): σ is the surface tension coefficient; ξ is the interface thickness.

[0059] When using the lattice Boltzmann flux solver, its equation is:

[0060]

[0061] In Equation (4): W is the macroscopic conserved quantity; t is time; is the divergence operator; F is the flow flux; S is the source term; The expansions of the terms in Equation (4) are:

[0062]

[0063] In Equation (5): p is the macroscopic pressure at the grid center; ρ is the macroscopic pressure at the grid center; u is the macroscopic horizontal velocity; v is the macroscopic vertical velocity; c s is the sound speed, given by the lattice velocity model; T is the transpose of this matrix; α is the direction of particle motion represented; e α represents the particle velocity in the corresponding direction, given by the lattice velocity model; is the particle distribution function in the corresponding direction after collision; is the particle distribution function in the corresponding direction before collision; is the comprehensive distribution function; is the gradient operator; F s is the surface tension; r is the displacement in the corresponding direction; t is the time; w α is the weight coefficient in the corresponding direction, given by the lattice velocity model; η = min(λ, 0.025), λ = tanh(|C L - C R | / (2(C L + C R + 0.2))), C L is the order parameter of the control volume on the left side of a certain side of the corresponding grid, C R is the order parameter of the control volume on the right side of a certain side of the corresponding grid; δ t is the given particle motion time; τ is the relaxation time; Since the present invention mainly uses the D2Q9 lattice velocity model as shown in Figure 2 : Therefore, e α represents the particle velocity in the corresponding direction, F s is the surface tension, and the corresponding coefficients are w0 = 4 / 9, w1 = w2 = w3 = w4 = 1 / 9, w5 = w6 = w7 = w8 = 1 / 36.

[0064] In the calculation process, the formula of the dual-time stepping method is adopted to utilize the advantages of inner iteration calculation. To accelerate the calculation speed, the second-order TVD Runge-Kutta time format is adopted. The formula of the dual-time stepping method is:

[0065]

[0066] In Equation (6): Ω is the control volume; is the macroscopic conserved quantity at the pseudo-time step; t * is the pseudo-time; R is the residual; is the macroscopic conserved quantity at the nth real-time step; is the macroscopic conserved quantity at the (n - 1)th real-time step; Δt is the global time step; R * is the final residual recorded at the grid center after implicit residual smoothing.

[0067] S4: During the solution process of step S3, continuously calculate the local time step and apply it to the iteration process of the dual-time stepping method until the inner iteration residual meets the preset convergence requirement, and obtain the flow field phase field numerical values of the entire computational domain at the next moment; and in order to accelerate the convergence of the inner iteration residual, the implicit residual smoothing method is adopted to increase the local time step to accelerate the convergence speed. The formula of the implicit residual smoothing method is:

[0068]

[0069] In Equation (7): ε I and ε J are the relaxation coefficients in the directions of computational domain I and computational domain J, respectively; is the residual value recorded at the center of the (I-1, J) grid after the first smoothing in the I direction; is the residual value recorded at the center of the (I, J) grid after the first smoothing in the I direction; is the residual value recorded at the center of the (I+1, J) grid after the first smoothing in the I direction; is the initial residual value recorded at the center of the (I, J) grid before smoothing; is the final residual value recorded at the center of the (I-1, J) grid after the second smoothing in the J direction; is the final residual value recorded at the center of the (I, J) grid after the second smoothing in the J direction; is the final residual value recorded at the center of the (I+1, J) grid after the second smoothing in the J direction; is the residual value recorded at the center of the (I, J) grid after the first smoothing in the I direction. In the embodiment given in the present invention: ε I = ε J = 0.75.

[0070] Embodiment

[0071] The present invention further elaborates on the multi-phase flow simulation efficiency optimization method and its superiority provided by the present invention based on a specific embodiment of the multi-phase lattice Boltzmann flux method.

[0072] This embodiment simulates the problem of single bubble rising under a large density ratio. The description of the geometric model is as Figure 3 shown. The computational domain grid is divided into 241×481, the reference length is the bubble diameter D = 120, the entire dimensionless computational domain is [0, 2D]×[0, 4D], the initial bubble is placed at the position (D, D), the upper and lower wall surfaces adopt no-slip boundary conditions, while the left and right wall surfaces adopt periodic boundary conditions.

[0073] This single bubble rising problem is determined by two dimensionless numbers. The definition methods of the Reynolds number and the Eötvös number are:

[0074]

[0075] In Equation (8): The characteristic velocity g is the gravitational acceleration; σ is the surface tension coefficient; the density ρ H of the heavier fluid = 1, the density ρ L of the lighter fluid = 0.001; U = 0.0012, Re = 35, Eo = 125, the interface thickness ξ = 4, the viscosity ratio μ H / μL = 100, reference time T = D / U.

[0076] Correspondingly, to ensure that the calculation does not diverge, a global time step of Δt = 1 is used when calculating with the original MLBFS, and to illustrate the impact of the present invention on the selectable global time, a global time step of Δt = 10 is used when using the efficiency optimization algorithm of the present invention, and the calculation time range is t = 0 - 7.

[0077] For this embodiment, the relaxation time τ can take values as:

[0078]

[0079] In Equation (9): δ t = δ x = 0.498, c = δ x / δ t .

[0080] Correspondingly, under the condition that the reference speed U = 0.0012 calculated by the original MLBFS method and the method of the present invention after efficiency optimization, the comparison diagram of the bubble interface at t = 5 is as Figure 4 shown. Among them, present represents the interface position of the bubble after using the efficiency optimization method, while original is the result of the bubble interface position calculated by the original MLBFS method.

[0081] As shown, the interface calculation result of the efficiency optimization method for the multiphase flow problem is basically consistent with the interface calculation result of the original MLBFS method for this embodiment, indicating that the efficiency optimization method proposed by the present invention does not affect the calculation accuracy of the multiphase flow problem, and proving the feasibility of the efficiency optimization method.

[0082] Furthermore, when calculating using the efficiency optimization algorithm, the number of inner iteration steps required for each global time step advancement is extracted, and the extracted results are arranged and distributed in chronological order, as Figure 5 shown.

[0083] According to Figure 5 shown, the distribution of the number of inner iteration steps is relatively scattered. To clarify the effectiveness of the present invention for calculating efficiency optimization, the number of inner iteration steps of the original MLBFS is extracted and compared with the total number of steps of the efficiency optimization algorithm under the condition of low reference speed. The comparison of the total number of inner iteration steps and the efficiency improvement of the single bubble rising problem based on the low reference speed using the method of the present invention and the total number of steps of the original MLBFS is shown in Table 1.

[0084] Example Method type Characteristic velocity U Total number of inner iteration steps Improve efficiency 1 Original method 0.0012 700000 0 2 Efficiency optimization method 0.0012 600516 14.2％

[0085] Table 1 Comparison of solutions for single-bubble rising problems at low reference speeds

[0086] As shown in Table 1, compared with the total number of steps required for the original MLBFS calculation, the total number of inner iterations required for the efficiency optimization method has decreased from 700,000 steps to 600,516 steps, and the calculation efficiency has increased by 14.2%. It can be seen that for the calculation of multiphase flow problems, the efficiency optimization algorithm proposed by the present invention is significantly effective.

[0087] Based on an embodiment of the present application, while keeping other settings unchanged, the magnitude of the characteristic velocity is corrected to clarify the influence of the present invention on the stability of multiphase flow calculations and the reliability of efficiency optimization.

[0088] For modifying the characteristic velocity, under the condition of keeping other initial calculation settings unchanged, with the characteristic velocity U = 0.006, the efficiency optimization method of the present invention is used to calculate multiphase flow problems. Under this setting condition, for the original MLBFS method, when the characteristic velocity exceeds 0.0012, the calculation is prone to instability problems and diverges, resulting in the final calculation failure. Therefore, a larger characteristic velocity is selected to verify that the method proposed by the present invention can improve stability while reducing the calculation period, thereby achieving a greater breakthrough in calculation efficiency.

[0089] Under the setting condition of U = 0.006, the global time step of Δt = 10 is used with the efficiency optimization algorithm of the present invention, and the calculation time range is t = 0 - 7.

[0090] Correspondingly, the comparison diagram of the bubble interface at t = 5 under the condition of the reference speed U = 0.0012 for the original MLBFS method and the reference speed U = 0.006 calculated by the present invention method after efficiency optimization is as Figure 6 shown. Among them, present represents the interface position of the bubble after using the efficiency optimization method, while original is the result of the bubble interface position calculated by the original MLBFS method.

[0091] The interface calculation results of the efficiency optimization method for multiphase flow problems shown are basically consistent with the interface calculation results of the original MLBFS method for this embodiment, indicating that the efficiency optimization method proposed by the present invention does not affect the calculation accuracy of multiphase flow problems, proving the feasibility of the efficiency optimization method.

[0092] Furthermore, when using the efficiency optimization algorithm for calculation, the number of inner iterations required for each global time step advancement is extracted, and the extracted results are arranged in chronological order, as Figure 7 shown.

[0093] According to Figure 7As shown, after improving the characteristic speed using the efficiency optimization method, this type of high density ratio multiphase flow problem can be effectively calculated, and the variation range of the number of inner iteration steps is significantly reduced compared to Figure 5 , indicating that the efficiency optimization method proposed in the present invention can significantly improve the stability of the calculation for multiphase flow problems.

[0094] To clarify the effectiveness of the present invention in optimizing the calculation efficiency, the number of inner iteration steps of the original MLBFS was extracted and compared with the total number of steps of the efficiency optimization algorithm under the condition of high reference speed.

[0095] The comparison of the total number of inner iterations of the single bubble rising problem based on high reference speed using the method of the present invention with the total number of steps of the original MLBFS and the comparison results of the efficiency improvement are shown in Table 2.

[0096] Example Method type Characteristic velocity U Total number of inner iteration steps Improve efficiency 1 Original method 0.0012 700000 0 2 Efficiency optimization method 0.0060 185845 73.5％

[0097] Table 2 Comparison of solving the single bubble rising problem based on high reference speed

[0098] According to Table 2, compared with the total number of steps required for the original MLBFS calculation, the total number of inner iterations required for the efficiency optimization method calculation has decreased from 700,000 steps to 185,845 steps, and the calculation efficiency has increased by 73.5%. It can be seen that for the calculation of multiphase flow problems, the efficiency optimization algorithm proposed in the present invention is significantly effective.

[0099] Combined with Figures 1 - 7 Tables 1 - 2, compared with the original MLBFS method, the efficiency optimization method proposed in the present invention, based on the dual time method, can significantly increase the global time step and reduce the calculation time by introducing local time step and implicit residual smoothing technology. For multiphase flow problems with limited reference speed, the reference speed can be increased to a large extent, the calculation period can be reduced, the calculation efficiency can be improved, and the calculation stability has also been improved to a certain extent.

[0100] The above embodiments are only used to illustrate the technical idea of the present invention, and the protection scope of the present invention cannot be limited thereby. Any changes made on the basis of the technical solution according to the technical idea proposed by the present invention fall within the protection scope of the present invention; technologies not covered by the present invention can be realized through existing technologies.

Claims

1. A method for optimizing the simulation efficiency of multiphase flow based on the multiphase lattice Boltzmann flux method, characterized in that: The steps of the efficiency optimization method are as follows: S1: Establish a physical model corresponding to the multiphase flow problem and perform grid division on the entire computational domain; S2: Input the initial conditions and boundary conditions required for the calculation according to the problem type; S3: Based on the finite volume method, perform spatial discretization through the Cahn-Hilliard equation controlled by the phase field and the lattice Boltzmann flux solver controlled by the flow field. Set the global time step, use the dual-time stepping method in the flow solution process, and use the second-order TVD Runge-Kutta time format to perform the advancement calculation on the pseudo-time step of the inner iteration; S4: Continuously calculate the local time step during the solution process of step S3 and apply it to the iteration process of the dual-time stepping method until the inner iteration residual meets the preset convergence requirement, and obtain the flow field phase field numerical values of the entire computational domain at the next moment; In step S4, in order to accelerate the convergence of the inner iteration residual, an implicit residual smoothing method is used to increase the local time step to accelerate the convergence speed; The formula of the implicit residual smoothing method is: In Equation (7): ε I and ε J are the relaxation coefficients in the directions of computational domain I and computational domain J, respectively; is the residual value recorded at the grid center of (I-1, J) after the first smoothing in the I direction; is the residual value recorded at the grid center of (I, J) after the first smoothing in the I direction; is the residual value recorded at the grid center of (I+1, J) after the first smoothing in the I direction; is the initial residual value recorded at the grid center of (I, J) before smoothing; is the final residual value recorded at the grid center of (I-1, J) after the second smoothing in the J direction; is the final residual value recorded at the grid center of (I, J) after the second smoothing in the J direction; is the final residual value recorded at the grid center of (I+1, J) after the second smoothing in the J direction; is the residual value recorded at the grid center of (I, J) after the first smoothing in the I direction.

2. The method for optimizing the simulation efficiency of multiphase flow based on the multiphase lattice Boltzmann flux method according to claim 1, wherein: The physical model in step S1 is a two-phase flow model.

3. The method for optimizing the multiphase flow simulation efficiency based on the multiphase lattice Boltzmann flux method according to claim 1, characterized in that: The initial distribution of the phase field in the initial conditions in step S2 is obtained from the phase field distance sign function, and the selected phase field distance sign function is: In formula (1): C represents the order parameter in the phase field, that is, the volume fraction of the heavier mass fluid in the corresponding control volume; d is the perpendicular distance from the control volume to the interface of the multiphase fluid, and its sign is determined by the initial control volume density. For two-phase flow, the distance of the control volume at the high-density fluid is positive, and the distance of the control volume at the low-density fluid is negative; ξ is the defined interface thickness.

4. The method for optimizing the multiphase flow simulation efficiency based on the multiphase lattice Boltzmann flux method according to claim 1, characterized in that: The Cahn-Hilliard equation controlled by the phase field in step S3 is: In Equation (2): C represents the order parameter in the phase field, i.e., the volume fraction of the heavier mass fluid in the control volume; t is time; is the divergence operator; u is the flow field velocity vector; Γ is the mass flow rate; μ C represents the chemical potential, and μ C is: In formula (3): σ is the surface tension coefficient; ξ is the interface thickness.

5. The method for optimizing the multiphase flow simulation efficiency based on the multiphase lattice Boltzmann flux method according to claim 1, wherein: The lattice Boltzmann flux solver controlled by the flow field in step S3 is: In Equation (4): W is the macroscopic conserved quantity; t is time; is the divergence operator; F is the flow flux; S is the source term; In Equation (5): p is the macroscopic pressure at the grid cell center; ρ is the macroscopic pressure at the grid cell center; u is the macroscopic horizontal velocity; v is the macroscopic vertical velocity; c s is the speed of sound, given by the lattice velocity model; T is the transpose of the matrix; α is the particle motion direction represented; e α represents the particle velocity in the corresponding direction, given by the lattice velocity model; is the particle distribution function in the corresponding direction after collision; is the particle distribution function in the corresponding direction before collision; is the comprehensive distribution function; is the gradient operator; F s is the surface tension; r is the displacement in the corresponding direction; t is the time; w α is the weight coefficient in the corresponding direction, given by the lattice velocity model; η = min(λ, 0.025), λ = tanh(|C L - C R | / (2(C L + C R + 0.2))), C L is the order parameter of the control volume on the left side of a certain side of the corresponding grid, and C R is the order parameter of the control volume on the right side of a certain side of the corresponding grid; δ t is the given particle motion time; τ is the relaxation time.

6. The method for optimizing the efficiency of multiphase flow simulation based on the multiphase lattice Boltzmann flux method according to claim 1, characterized in that: The formula of the dual-time stepping method in steps S3 and S4 is: In Equation (6): Ω is the volume of the control volume; is the macroscopic conserved quantity at the pseudo time step; t * is the pseudo time; R is the residual; is the macroscopic conserved quantity at the nth real time step; is the macroscopic conserved quantity at the (n - 1)th real time step; Δt is the global time step; R * is the final residual recorded at the grid center after implicit residual smoothing.

Images