Two-phase interface area calculation method and system for simulating two-phase flow through phase field model

Abstract

The invention relates to a two-phase interface area calculation method and system for a phase field model to simulate two-phase flow, and the method comprises the following steps: S1, extracting a sequence parameter value on each grid point from a simulation result of the phase field model, and calculating a sequence parameter gradient value of each grid point; s2, calculating the body free energy density and the interface free energy density of each point of the two-phase system according to the sequence parameter value and the sequence parameter gradient value; s3, performing volume fraction on the calculated body free energy density and interface free energy density to obtain system free energy F (phi); s4, calculating a surface tension value sigma of the two-phase system according to parameter setting of the phase field model; and S5, calculating the interface surface area AF according to the system free energy F (phi) and the surface tension value sigma. According to the method disclosed by the invention, a complex interface discrete process is avoided, so that a complex grid discrete algorithm and a data structure storage mode related to the complex interface discrete process are also avoided, and the method is easy to realize parallel computing and can further improve the computing efficiency.

Description

Technical Field

[0001] The present disclosure relates to the technical field of two-phase fluid mechanics, and in particular to a two-phase interface area calculation method and system for simulating two-phase flow using a phase field model. Background Art

[0002] Two-phase flow is widely present in nature and industrial processes, such as wave breaking and chemical reactions. With the continuous improvement of computer performance and the development of algorithms, numerical simulation has become an important way to study such problems. Among them, the phase field model has a natural advantage in solving two-phase flow problems due to its physical basis.

[0003] However, the existing Sharp Interface method requires the surface to be discretized. When processing complex surfaces or for higher precision, finer face elements need to be used. Generating and processing a large number of polygonal meshes requires a lot of computing resources and time, and also places high demands on memory and storage devices. In addition, the existing method is mainly integrated in CFD post-processing software, and its code logic is complex. Its discretization process is extremely complex, involving fields such as computer graphics and data structures. In particular, it is difficult to implement parallel operations and difficult to write code to be integrated into self-developed computing programs. Summary of the invention

[0004] In order to solve the problems existing in the above-mentioned prior art, the purpose of the present disclosure is to provide a method and system for calculating the two-phase interface area for simulating two-phase flow using a phase field model. The method for calculating the two-phase interface area for simulating two-phase flow using a phase field model avoids the complex interface discretization process, and thus also avoids the complex grid discretization algorithm and data structure storage method associated therewith, and because the calculation process is based on simple mathematical operations (such as gradient calculation, integration, etc.), the code logic is clearer and the calculation process is relatively simple, making the method easy to implement parallel calculations to further improve the calculation efficiency; and the method is directly based on the thermodynamic basis of the phase field model, and compared with the commonly used Sharp Interface method, it can be more in line with the physical basis of the phase field model and can more accurately describe the interface behavior.

[0005] The method for calculating the two-phase interface area for simulating two-phase flow using a phase field model disclosed in the present disclosure comprises the following steps:

[0006] S1. Extract the order parameter value at each grid point from the phase field model simulation results, and calculate the order parameter gradient value at each grid point;

[0007] S2. Calculating the volume free energy density and the interface free energy density of each point in the two-phase system according to the order parameter value and the order parameter gradient value;

[0008] S3, performing volume integration on the calculated volume free energy density and the calculated interface free energy density to obtain system free energy F(φ);

[0009] S4. Calculate the surface tension value σ of the two-phase system according to the parameter setting of the phase field model;

[0010] S5. Calculate the interface surface area A according to the system free energy F(φ) and the surface tension value σ F .

[0011] Preferably, the step S1 specifically includes:

[0012] Read the result file of phase field model simulation and extract the order parameter value φ at each grid point in the file;

[0013] Calculate the order parameter gradient according to the order parameter value φ The order parameter gradient The calculation is as follows:

[0014] For points inside the flow field or at the periodic boundary, the central difference method is used to construct a high-order difference method:

[0015]

[0016] For the points at the boundary, the one-sided difference method is used to construct a high-order one-sided interpolation method:

[0017]

[0018] Here, i represents the x-coordinate of the point.

[0019] Preferably, the step S2 specifically includes:

[0020] The volume free energy density ψ(φ) is calculated as a function of the order parameter value, using a double potential well model:

[0021] ψ(φ)=β(φ-φ A ) 2 (φ-φ B ) 2

[0022] According to the order parameter gradient The interface free energy density is calculated, and the interface free energy density is:

[0023]

[0024] Among them, φ A and φ Bare the order parameters of the two-phase fluid; k is used to control the interface thickness. When k approaches 0, it indicates a sharp interface and also indicates the relative size of the interface free energy; β is the coefficient of the double potential well model, indicating the relative size of the body free energy.

[0025] Preferably, in step S3, the system free energy F(φ) is calculated according to the following formula:

[0026]

[0027] Wherein, ψ(φ) is the body free energy density; is the interface free energy density.

[0028] Preferably, the step S4 specifically includes:

[0029] According to the parameter setting of the phase field model, the surface tension value σ of the two-phase system is calculated according to the following formula:

[0030]

[0031] Among them, φ A and φ B are the order parameters of the two-phase fluid and satisfy φ A >φ B ; W is the parameter of the two-phase interface thickness; k is used to control the interface thickness. When k tends to 0, it indicates a sharp interface and also indicates the relative size of the interface free energy; β is the coefficient of the double potential well model, indicating the relative size of the body free energy.

[0032] Preferably, in step S5, the interface surface area A F Calculate according to the following formula:

[0033] A F =F(φ) / σ

[0034] Wherein, F(φ) is the free energy of the system; σ is the surface tension value.

[0035] The two-phase interface area calculation system used for phase field model simulation of two-phase flow includes:

[0036] at least one memory;

[0037] at least one processor;

[0038] at least one program;

[0039] The programs are stored in the memory, and the processor executes at least one of the programs to implement the two-phase interface area calculation method for simulating two-phase flow using a phase field model as described in any one of the above items.

[0040] The two-phase interface area calculation method and system for simulating two-phase flow using a phase field model disclosed in the present disclosure have the following advantages:

[0041] 1. The method for calculating the two-phase interface area for simulating two-phase flow using a phase field model disclosed in the present invention comprises the following steps: S1. extracting the order parameter value at each grid point from the simulation result of the phase field model, and calculating the order parameter gradient value at each grid point; S2. calculating the volume free energy density and the interface free energy density at each point of the two-phase system according to the order parameter value and the order parameter gradient value; S3. performing volume integration on the calculated volume free energy density and interface free energy density to obtain the system free energy F(φ); S4. calculating the surface tension value σ of the two-phase system according to the parameter setting of the phase field model; S5. calculating the interface surface area A according to the system free energy F(φ) and the surface tension value σ F The above method directly uses the order parameter and gradient of the phase field model for calculation without additional surface discretization, avoiding the need for additional discretization of the surface in the traditional discrete grid method, and reducing the requirements for memory and storage devices; the above method avoids the complex interface discretization process, and thus also avoids the related complex grid discretization algorithm and data structure storage method, and because the calculation process is based on simple mathematical operations (such as gradient calculation, integration, etc.), the code logic is clearer and the calculation process is relatively simple, making the method easy to implement parallel calculation to further improve the calculation efficiency; and the above method is directly based on the thermodynamic basis of the phase field model, compared with the commonly used Sharp Interface method, it can be more in line with the physical basis of the phase field model and can more accurately describe the interface behavior.

[0042] 2. The two-phase interface area calculation system for simulating two-phase flow with a phase field model disclosed in the present invention comprises: at least one memory; at least one processor; at least one program; the program is stored in the memory, and the processor executes at least one program to implement the two-phase interface area calculation method for simulating two-phase flow with a phase field model as described in any of the above items. The system implements the two-phase interface area calculation method for simulating two-phase flow with a phase field model. Since the calculation process mainly depends on the order parameter value on the grid point, the system can easily implement parallel calculation, which can not only further improve the calculation speed, but also better utilize the multi-core processor resources of the computer; and by using the above two-phase interface area calculation method for simulating two-phase flow with a phase field model, the complex interface discretization processing in the traditional method can be avoided, thereby reducing the demand for memory and storage devices, which is particularly important for handling large-scale two-phase flow simulation problems. BRIEF DESCRIPTION OF THE DRAWINGS

[0043] Figure 1 is a flow chart of the method for calculating the two-phase interface area for simulating two-phase flow using a phase field model according to the present disclosure;

[0044] Figure 2 It is a schematic diagram of the theoretical relationship between the dimensionless surface area of ​​a spherical droplet and the dimensionless interface thickness obtained by the two-phase interface area calculation method for simulating two-phase flow using a phase field model described in the present disclosure;

[0045] Figure 3 It is a schematic diagram of the evolution of the dimensionless surface area of ​​a stationary droplet and the correlation coefficient between the volume free energy and the interface free energy over time in the method for calculating the two-phase interface area for simulating two-phase flow with a phase field model described in the present disclosure. DETAILED DESCRIPTION

[0046] like Figure 1 - Figure 3 As shown, the method for calculating the two-phase interface area for simulating two-phase flow using a phase field model described in the present disclosure comprises the following steps:

[0047] S1. Extract the order parameter value at each grid point from the phase field model simulation results, and calculate the order parameter gradient value at each grid point;

[0048] S2. Calculate the volume free energy density and interface free energy density of each point in the two-phase system according to the order parameter value and the order parameter gradient value;

[0049] S3, perform volume integration on the calculated volume free energy density and interface free energy density to obtain the system free energy F(φ);

[0050] S4. Calculate the surface tension value σ of the two-phase system according to the parameter setting of the phase field model;

[0051] S5. Calculate the interface surface area A based on the system free energy F(φ) and the surface tension value σ F ;

[0052] The above method directly uses the order parameter and gradient of the phase field model for calculation without additional discretization, avoiding the need for additional discretization of the surface in the traditional discrete grid method, and reducing the requirements for memory and storage devices; the above method avoids the complex interface discretization process, and thus also avoids the related complex grid discretization algorithm and data structure storage method, and because the calculation process is based on simple mathematical operations (such as gradient calculation, integration, etc.), the code logic is clearer and the calculation process is relatively simple, making the method easy to implement parallel calculation to further improve the calculation efficiency; and the above method is directly based on the thermodynamic basis of the phase field model, compared with the commonly used Sharp Interface method, it can be more in line with the physical basis of the phase field model and can more accurately describe the interface behavior.

[0053] Furthermore, in this embodiment, step S1 specifically includes:

[0054] Read the result file of phase field model simulation and extract the order parameter value φ at each grid point in the file;

[0055] Calculate the order parameter gradient according to the order parameter value φ Order parameter gradient The calculation is as follows:

[0056] For points inside the flow field or at the periodic boundary, the central difference method is used to construct a high-order difference method:

[0057]

[0058] For the points at the boundary, the one-sided difference method is used to construct a high-order one-sided interpolation method:

[0059]

[0060] Where i represents the x-coordinate of the point;

[0061] The above calculation method is based on a uniform grid in the x direction.

[0062] Furthermore, in this embodiment, step S2 specifically includes:

[0063] The volume free energy density ψ(φ) is calculated as a function of the order parameter value according to the double potential well model:

[0064] ψ(φ)=β(φ-φ A ) 2 (φ-φ B ) 2

[0065] According to the order parameter gradient Calculate the interface free energy density, which is:

[0066]

[0067] Among them, φ A and φ B are the order parameters of the two-phase fluid; k is used to control the interface thickness. When k approaches 0, it indicates a sharp interface and also indicates the relative size of the interface free energy; β is the coefficient of the double potential well model, indicating the relative size of the body free energy.

[0068] Furthermore, in this embodiment, in step S3, the system free energy F(φ) is calculated according to the following formula:

[0069]

[0070] Where ψ(φ) is the volume free energy density; is the interface free energy density.

[0071] Furthermore, in this embodiment, step S4 specifically includes:

[0072] According to the parameter setting of the phase field model, the surface tension value σ of the two-phase system is calculated according to the following formula:

[0073]

[0074] Among them, φ A and φ B are the order parameters of the two-phase fluid and satisfy φ A >φ B ; W is the parameter of the two-phase interface thickness; k is used to control the interface thickness. When k tends to 0, it indicates a sharp interface and also indicates the relative size of the interface free energy; β is the coefficient of the double potential well model, indicating the relative size of the body free energy.

[0075] Furthermore, in this embodiment, in step S5, the interface surface area A F Calculate according to the following formula:

[0076]

[0077] Where F(φ) is the system free energy; σ is the surface tension value;

[0078] Combined with the above content, the principle of the two-phase interface area calculation method used for phase field model simulation of two-phase flow is further analyzed, and the above interface surface area A is explained in combination with the principle. F Derivation process of calculation method:

[0079] The phase field model is a commonly used Diffuse Interface model for solving interface flow problems. In the phase field model, the free energy functional of the two-phase flow includes the body free energy and the interface free energy, which can be expressed as:

[0080]

[0081] Among them, Ψ(φ) is a double potential well form: Ψ(φ) = β(φ-φ A ) 2 (φ-φ B ) 2 , where φ A and φ B is the order parameter of the two-phase fluid (φ A >φ B );

[0082] In the understanding of the Diffuse Interface model, φ B <φ<φ AThe area can be considered as the two-phase interface. If the sharp interface model is used, the two-phase interface can be defined as The isosurface of

[0083] Ψ(φ) is the bulk free energy density, which is the bulk free energy density of the hydrophobic effect and represents the bulk region that separates the two phases into a single phase;

[0084] is the interfacial free energy density, which, as a hydrophilic effect, tends to make the two phases mix more;

[0085] The positive constants k and β represent the relative magnitudes of the interfacial free energy and the bulk free energy, respectively. The first-order variation of the free energy F with respect to the order parameter φ gives the chemical potential μ φ ,Right now:

[0086]

[0087] For a plane interface in equilibrium, The plane solution of φ can be obtained at this time:

[0088]

[0089] Here ζ represents the distance perpendicular to the interface; The parameter representing the thickness of the interface between the two phases; the surface tension between the two phases is

[0090] Let's start with the simplest plane interface. For a plane interface perpendicular to the x-direction, the surface tension can be calculated using the following integral:

[0091]

[0092] Calculate the free energy of the flat interface system. Since the volume free energy density and the interface free energy density are equal, we can:

[0093]

[0094] Among them, A φ for The area of ​​the isosurface at ;

[0095] In gradient calculation, the central difference method is used for the internal points of the fluid to ensure its accuracy. For periodic boundaries, it is easy to construct the central difference, and for no-slip walls, a second-order difference format is constructed; and the difference method can use the first-order difference or the second-order higher-order difference method depending on the situation;

[0096] This formula shows that for a plane in equilibrium, we can use to calculate the interface area.

[0097] For a two-dimensional or three-dimensional static droplet (radius R0), its surface tension is calculated as follows:

[0098]

[0099] In this result, Ch is the Cahn number (Ch = W / R0). When the Ch number is much smaller than 1 (which is usually required in phase field simulation), this result can simplify the surface tension of the flat plate interface.

[0100] The free energy of a two-dimensional static droplet can be calculated as follows:

[0101]

[0102] When Ch is small, expand the above formula:

[0103]

[0104] Among them, α=exp(-4 / Ch)<1, so F(φ) 2D Relative to A 2D The deviation of σ=2πR0 can be obtained through the error term To make an estimate, it can be seen from Table 1 below that when Ch is small, the error term is extremely small;

[0105] Table 1 Error terms for two-dimensional and three-dimensional static droplets

[0106] Ch 0.1 0.5 1.0 1.5 2.0 <![CDATA[Err 2D ]]> 0.0000 0.0000 0.0001 0.0024 0.0104 <![CDATA[Err 3D ]]> 0.0008 0.0202 0.0806 0.0806 0.3167

[0107] For a three-dimensional stationary droplet, the volume free energy calculation method is the same as the two-dimensional case, and the result is:

[0108]

[0109] Similar to a two-dimensional droplet, when Ch is small, expansion is performed:

[0110]

[0111] Here we can use the error term To estimate F(φ) 3D With A 3D σ=4πR0 deviation; From Table 1 above, we can see that when Ch<1, the result has a good performance. Compared with the two-dimensional result, the error term of the three-dimensional result contains the Ch quadratic term without α, so a larger error will be generated in the three-dimensional case;

[0112] Figure 2The relationship between the dimensionless surface area and the interface thickness Ch in the static droplet example shown (the green dashed line is a constant of 1), where E0 represents the free energy calculated using the surface area A (since it is a static droplet, its size does not change and can be calculated accurately) and the surface tension σ, while F(φ) is the free energy calculated using its definition;

[0113] F(φ) / E0 can be expressed as That is, the dimensionless surface area. When this value approaches 1, it means that the surface area can be accurately calculated using this method. From the figure, we can see that for the two-dimensional example, the surface area calculated using this method is consistent with the theoretical surface area, and it has good accuracy even when the Ch number is large;

[0114] For the three-dimensional example, when the Ch number is small, this method can accurately calculate the interface surface area. Table 1 shows the specific error term values. When the Ch number is equal to 0.5, the relative error of the result is 2%, which is within an acceptable range. When initializing, the Ch number is generally set to 0.1. Therefore, the surface area of ​​the stationary droplet calculated using this method is sufficiently accurate.

[0115] If A φ (A φ represents the surface area calculated using the traditional method) as the standard value of the area, then the correction factor can be determined based on the Ch number, and then in A F We can also get A by making corrections based on φ ;

[0116] If the droplet undergoes complex changes, such as large deformation or breaking into multiple droplets of different sizes, then Ch will change significantly, and the whole field is not unique. Therefore, it is necessary to consider making a statistical correction based on the interface curvature, droplet size, number of droplets, etc., which is more complicated. Generally speaking, a larger Ch will appear at this time, so A F Slightly larger than A φ In fact, since the phase field model is based on the Diffuse Interface model, it can be considered that A F is the true interface area, no correction is needed;

[0117] Figure 3 The figure shows the evolution of the dimensionless surface area of ​​a stationary droplet and the correlation coefficient between the volume free energy and the interface free energy over time (the green dashed line in the figure is a constant of 1, and the pink dashed line is a constant of 1.0007, which is A F / A0 and A F / A φ ), where A0 represents the initial moment, the surface area of ​​the droplet calculated using the traditional method, and A Frepresents the surface area calculated using this method, A φ represents the surface area calculated using the traditional method, is the dimensionless time adopted;

[0118] Since the interface order parameter distribution used satisfies the hyperbolic tangent, the volume free energy density and the interface free energy density have the same value, so the theoretical values ​​of the dimensionless surface area and the correlation coefficient between the volume free energy density and the interface free energy density of the traditional method and this method should be 1;

[0119] At the initial moment, the relative surface area (A F / A0 and A φ / A0) are 1.0007 and 0.9998 respectively, and the correlation coefficient between the volume free energy density and the interface free energy density is 1, both of which are very close to their theoretical values. As time evolves, due to the shortcomings of the phase field model, the droplet will gradually become smaller due to mass loss, resulting in a decrease in the droplet volume and surface area;

[0120] When it evolves to t * =39.1, the dimensionless surface area calculated by this method is 0.9973, and the traditional method is 0.9942. This method is closer to the theoretical value;

[0121] At the same time, the correlation coefficient between the bulk free energy density and the interface free energy density is 0.9987, which is very close to its theoretical value of 1, which means that the two free energy densities always maintain a strong correlation during the static evolution of the droplet, which verifies that the bulk free energy is equal to the interface free energy in the full field angle and under time variation;

[0122] In summary, the free energy of a static droplet can be approximately expressed as (Ch<<1):

[0123] F(φ)≈A φ σ

[0124] Based on the results of flat interfaces and two-dimensional and three-dimensional static droplets, a method for calculating the two-phase interface area for phase field model simulation of two-phase flow was obtained:

[0125]

[0126] Here A F is the interfacial area calculated using this method. According to the above derivation process, it can be clearly obtained that for a flat interface in equilibrium, A F =A φ ; For static droplets, when Ch<<1, A F ≈A φ ; For a two-dimensional droplet, A F and A φVery close, if the relative error is not more than 1%, then Ch<1.96 is sufficient; but for three-dimensional droplets, it is more demanding. F and A φ If the relative error does not exceed 1%, Ch<0.36 is required, which is easy to achieve using the phase field model;

[0127] Therefore, the interface area A F It is more in line with the idea of ​​the Diffuse Interface model used in the phase field model, which essentially represents the free energy of the system and is more convenient for programming and data post-processing.

[0128] The two-phase interface area calculation system used for phase field model simulation of two-phase flow includes:

[0129] at least one memory;

[0130] at least one processor;

[0131] at least one program;

[0132] The programs are stored in the memory, and the processor executes at least one program to implement the two-phase interface area calculation method for simulating two-phase flow using a phase field model as described in any one of the above items;

[0133] The system implements a two-phase interface area calculation method for simulating two-phase flow with a phase field model. Since the calculation process mainly depends on the order parameter value on the grid, the system can easily implement parallel calculation, which can not only further improve the calculation speed, but also better utilize the computer's multi-core processor resources; and by using the above-mentioned two-phase interface area calculation method for simulating two-phase flow with a phase field model, the complex interface discretization processing in the traditional method can be avoided, thereby reducing the demand for memory and storage devices, which is particularly important for dealing with large-scale two-phase flow simulation problems.

[0134] In the description of the present disclosure, it is necessary to understand that the directions or positional relationships indicated by directional words such as "front, back, up, down, left, right", "lateral, vertical, horizontal" and "top, bottom" are usually based on the directions or positional relationships shown in the drawings, and are only for the convenience of describing the present disclosure and simplifying the description. Unless otherwise specified, these directional words do not indicate or imply that the device or element referred to must have a specific orientation or be constructed and operated in a specific orientation, and therefore cannot be understood as limiting the scope of protection of the present disclosure.

[0135] For those skilled in the art, various other corresponding changes and deformations can be made according to the technical solutions and concepts described above, and all of these changes and deformations should fall within the protection scope of the claims of the present disclosure.

Claims

1. A method for calculating the two-phase interface area for simulating two-phase flow using a phase field model, characterized in that: The following steps are involved: S1. Extract the order parameter value at each grid point from the phase field model simulation results, and calculate the order parameter gradient value at each grid point; S2. Calculating the volume free energy density and the interface free energy density of each point in the two-phase system according to the order parameter value and the order parameter gradient value; S3, performing volume integration on the calculated volume free energy density and the calculated interface free energy density to obtain system free energy F(φ); S4. Calculate the surface tension value σ of the two-phase system according to the parameter setting of the phase field model; S5. Calculate the interface surface area A according to the system free energy F(φ) and the surface tension value σ F .

2. The method for calculating the two-phase interface area for simulating two-phase flow using a phase field model according to claim 1, characterized in that: The step S1 specifically includes: Read the result file of phase field model simulation and extract the order parameter value φ at each grid point in the file; Calculate the order parameter gradient according to the order parameter value φ The order parameter gradient The calculation is as follows: For points inside the flow field or at the periodic boundary, the central difference method is used to construct a high-order difference method: For the points at the boundary, the one-sided difference method is used to construct a high-order one-sided interpolation method: Here, i represents the x-coordinate of the point.

3. The method for calculating the two-phase interface area for simulating two-phase flow using a phase field model according to claim 2, characterized in that: The step S2 specifically includes: The volume free energy density ψ(φ) is calculated as a function of the order parameter value, using a double potential well model: ψ(φ)=β(φ-φ A ) 2 (f-f) B ) 2 According to the order parameter gradient The interface free energy density is calculated, and the interface free energy density is: Among them, φ A and φ B are the order parameters of the two-phase fluid; k is used to control the interface thickness. When k approaches 0, it indicates a sharp interface and also indicates the relative size of the interface free energy; β is the coefficient of the double potential well model, indicating the relative size of the body free energy.

4. The method for calculating the two-phase interface area for simulating two-phase flow using a phase field model according to claim 3, characterized in that: In step S3, the system free energy F(φ) is calculated according to the following formula: Wherein, ψ(φ) is the body free energy density; is the interface free energy density.

5. The method for calculating the two-phase interface area for simulating two-phase flow using a phase field model according to claim 4, characterized in that: The step S4 specifically includes: According to the parameter setting of the phase field model, the surface tension value σ of the two-phase system is calculated according to the following formula: Among them, φ A and φ B are the order parameters of the two-phase fluid and satisfy φ A >φ B ; W is the parameter of the two-phase interface thickness; k is used to control the interface thickness. When k tends to 0, it indicates a sharp interface and also indicates the relative size of the interface free energy; β is the coefficient of the double potential well model, indicating the relative size of the body free energy.

6. The method for calculating the two-phase interface area for simulating two-phase flow using a phase field model according to claim 5, characterized in that: In step S5, the interface surface area A F Calculate according to the following formula: Wherein, F(φ) is the free energy of the system; σ is the surface tension value.

7. A two-phase interface area calculation system for simulating two-phase flow using a phase field model, characterized in that: include: at least one memory; at least one processor; at least one program; The programs are stored in the memory, and the processor executes at least one of the programs to implement the two-phase interface area calculation method for simulating two-phase flow using a phase field model according to any one of claims 1 to 6.

Images