CN104750943A - Gas-liquid two-phase interface area calculation method based on numerical simulation - Google Patents

Abstract

The invention provides a method for calculating gas-liquid two-phase flowing phase interface area according to numerical simulation results. The method includes the steps: extracting data including position coordinates, grid side lengths, target fluid volume fraction and volume fraction gradient vectors corresponding to grids from the two-phase flowing VOF (volume of fluid) numerical simulation results; classifying phase interface morphology in the grids into five classes according to the target fluid volume fraction and the gradient vectors thereof; determining phase interface plane equation key parameters according to the grid side lengths, the target fluid volume fraction and the volume fraction gradient vectors; respectively calculating areas of phase interfaces of the five classes of the morphology according to phase interface plane equations, the grid side lengths and the volume fraction gradient vectors. Phase interface area data within specified time and a space range in a two-phase fluid system can be effectively acquired, and basic data can be provided for heat transfer, mass transfer and quantitative analysis of physical-chemical reaction states in various two-phase fluid systems.

Description

A kind of gas-liquid two-phase interfacial area computing method based on numerical simulation

Art

The present invention relates to two-phase flow quantitative test field, particularly relate to a kind of gas-liquid two-phase interfacial area computing method based on numerical simulation.

Background technology

Gas-liquid two-phase interfacial area calculates the surface of contact area typically referring to and obtain gas-liquid two-phase in gas-liquid two-phase system.

The test determination method that existing gas liquid film amasss mainly contains the optical instrument direct method of measurement and the chemical indirect method of measurement two class.The optical instrument direct method of measurement adopts electronic optical instrument to measure the size distribution of bubble, further appliance computer image processing techniques obtains phase contact area, but the method is only applicable to the dispersed system of continuous phase light-permeable, application surface is narrower, and is difficult to obtain flow field picture and accurate quantitative result clearly.Chemical method is the method that application is more at present, the method, by measuring solute absorption rate and solvents concentration, adopts Danckwerts drawing method to obtain phase contact area, but, the value of chemical assay is with the solute selected---solvent series is relevant, and the accuracy of experiment test means is poor.

Current numerical simulation technology can realize the simulation of the distribution sign in system to gas-liquid two-phase and interface thereof preferably, each mesh definition in VOF (Volume of Fluid) method stream field fluid volume function (i.e. the volume of target fluid and the ratio of mesh volume), the distribution characteristics at and interface mutually each according to the known polyphasic flow of this functional value.PLIC (Piecewise Linear Interface Calculation) is the important free surface tracking scheme of one for VOF method analog result, in each computing grid, adopting the method approximate reconstruction phase interface of " replace curve by straight line ", being widely used because having the high and feature that calculated amount is little of precision.The employing such as Zhang Qin PLIC-VOF method is reconstructed the phase interface in float glass process between glass metal and gas, Liu etc. utilize PLIC-VOF method to follow the trail of the motion of wave, Maurya etc. apply the Evolution History that PLIC-VOF method simulates phase interface in the upper gas-liquid film evaporation device process of tiltedly flat board, Yang Junying etc. adopt PLIC-VOF method to simulate position and the shape of flow front phase interface in resin transfer molding (RTM) process mold filling process, and these researchs show that the phase interface of gas-liquid two-phase dispersed system can more adequately be followed the trail of and reconstruct to PLIC-VOF method.But the method that application numerical simulation obtains interfacial area there is not yet bibliographical information.

From analyzing above, nowadays method for numerical simulation can realize the simulation of the distribution sign in system to gas-liquid two-phase and interface thereof preferably, traditional method has a lot of requirement to the physics of two-phase material, chemical property, and degree of accuracy is not high, has significant limitation in actual applications.

Summary of the invention

The invention provides a kind of gas-liquid two-phase interfacial area computing method based on numerical simulation, in order to calculate VOF numerical simulation gas-liquid two-phase system in the contact area of gas-liquid.

In order to calculate the gas-liquid two-phase interfacial area based on VOF numerical simulation, its step is as follows:

A. the two-phase flow VOF numerical simulation result adopting rectangular parallelepiped grid to divide for all or part of zoning calculates the phase contact area of Zone Full or appointed area, read each grid values analog result data in rectangular parallelepiped grid region to be analyzed, comprise grid element center coordinate (x _i, y _i, z _i), the grid length of side (δ _{i_x}, δ _{i_y}, δ _{i_z}), gas, liquid volume of fluid (c in grid _{i_g}, c _{i_l}) and the gradient vector of arbitrary volume of fluid i in aforementioned each variable represents i-th grid, and x, y, z represent three change in coordinate axis direction of three-dimensional cartesian coordinate system, and g, l represent gas phase and liquid phase respectively;

B. the order (if there is equal situation, then maintaining coherent element original relative rank) from small to large by element absolute value is to vector (n _{i_x}δ _{i_x}, n _{i_y}δ _{i_y}, n _{i_z}δ _{i_z}) element sequence, the new vector identification of reordering is (n _{i_a}δ _{i_a}, n _{i_b}δ _{i_b}, n _{i_c}δ _{i_c}) (in this expression formula, a, b, c are different and all belong to set { x, y, z}, such as, if abs (n _{i_x}δ _{i_x})≤abs (n _{i_y}δ _{i_y})≤abs (n _{i_z}δ _{i_z}), then a=x, b=y, c=z);

C. former coordinate system x-y-z is transformed to by rotation, symmetry operation namely represent that following matrix is taken advantage of on the equal right side of the row vector in position, direction in former coordinate system:

$A=\left(\begin{array}{ccc}\frac{\mathrm{abs}\left(n_{i\_x}\right)}{n_{i\_x}}(a==x)& \frac{\mathrm{abs}\left(n_{i\_x}\right)}{n_{i\_x}}(a==y)& \frac{\mathrm{abs}\left(n_{i\_x}\right)}{n_{i\_x}}(a==z)\\ \frac{\mathrm{abs}\left(n_{i\_y}\right)}{n_{i\_y}}(b==x)& \frac{\mathrm{zbs}\left(n_{i\_y}\right)}{n_{i\_y}}(b==y)& \frac{\mathrm{abs}\left(n_{i\_y}\right)}{n_{i\_y}}(b==z)\\ \frac{\mathrm{abs}\left(n_{i\_z}\right)}{n_{i\_z}}(c==x)& \frac{\mathrm{abs}\left(n_{i\_y}\right)}{n_{i\_z}}(c==y)& \frac{\mathrm{abs}\left(n_{i\_z}\right)}{n_{i\_z}}(c==z)\end{array}\right),$

In each element calculating formula of matrix, p==q (p=a, b, c, q=x, y, z) is logical operation, and it is 1 that equation sets up operation result, otherwise operation result is 0, and regulation works as n _{i_q}when=0 after coordinate transform, the normal vector at interface is computing grid is followed successively by δ along the length of side of new coordinate three change in coordinate axis direction _{i_a}, δ _{i_b}, δ _{i_c};

D. classify according to table 1 pair phase interface pattern, in table, c=min (c _{i_g}, c _{i_l}), v=δ _{i_a}δ _{i_b}δ _{i_c}for rectangular parallelepiped grid volume, n _a, n _b, n _crepresent abs (n respectively _{i_a}), abs (n _{i_b}), abs (n _{i_c}), in table, each symbol all refers in particular to variable corresponding to i-th computing grid in addition, in order to express easily omits subscript i_, such as δ _arepresent δ _{i_a};

E. the type belonging to phase interface and table 2 method determination phase interface plane equation n _{i_a}x+n _{i_b}y+n _{i_c}undetermined constant d in z=d, in table, each symbol all refers in particular to variable corresponding to i-th computing grid in addition, in order to express easily omits subscript i_, such as δ _arepresent δ _{i_a};

F. the calculating formula of phase contact area determined by the type belonging to phase interface, phase plane equation and table 3, calculates the interfacial area A of this grid _i, in table, each symbol all refers in particular to variable corresponding to i-th computing grid in addition, in order to express easily omits subscript i_, such as δ _arepresent δ _{i_a};

G. total phase contact area is by cumulative for the phase contact area of grids all in region to be analyzed,

The criterion of five class basic interfaces in table 1 rectangular parallelepiped

The calculating formula of table 2. five class phase interface plane equation constant d

The calculating formula of table 3 five class phase contact area

Good effect of the present invention:

1. the optical instrument direct method of measurement that the method applicability is more traditional and the chemical indirect method of measurement are more extensive, not by the physics of two-phase material, the restriction of chemical property, can process the insurmountable problem of many classic methods.

2. the method form is simple, is easy to computing machine and realizes, effectively can obtain the phase contact area information in two-phase flow system in different time, spatial dimension, can be widely used in the industry such as metallurgy, chemical industry.

Accompanying drawing explanation

Fig. 1 is that the gas-liquid two-phase interfacial area based on numerical simulation of the embodiment of the present invention calculates five class basic interface grid schematic diagram.

Fig. 2 is the gas-liquid two-phase interfacial area calculation flow chart based on numerical simulation of the embodiment of the present invention.

Embodiment

Below in conjunction with accompanying drawing and illustrative examples, the present invention is further described.

Fig. 1 is that the five class rectangular parallelepiped grid interface topographies that the present invention is based on VOF numerical simulation divide schematic diagram.5 class patterns shown in Fig. 1 represent 55 kinds of patterns in original coordinate system, such as first, upper left rectangular parallelepiped (phase interface cut hexahedron angle number be 1 or 7) by space bit to symmetry and turning operation can reappear pattern in 8 in original coordinate system (phase interface cuts hexahedral 8 angles respectively).

Fig. 2 is computer implemented method process flow diagram of the present invention.

Below further describe concrete technical scheme of the present invention, so that those skilled in the art understands the present invention further, and do not form the restriction to its right.

Embodiment 1 calculates gas-liquid two-phase interfacial area by two-phase flow VOF numerical simulation result, and its step is as follows:

A. the two-phase flow VOF numerical simulation result adopting rectangular parallelepiped grid to divide for all or part of zoning calculates the phase contact area of Zone Full or appointed area, read each grid values analog result data in rectangular parallelepiped grid region to be analyzed, comprise grid element center coordinate (x _i, y _i, z _i), the grid length of side (δ _{i_x}, δ _{i_y}, δ _{i_z}), gas, liquid volume of fluid (c in grid _{i_g}, c _{i_l}) and the gradient vector of arbitrary volume of fluid i in aforementioned each variable represents i-th grid, and x, y, z represent three change in coordinate axis direction of three-dimensional cartesian coordinate system, and g, l represent gas phase and liquid phase respectively;

B. the order (if there is equal situation, then maintaining coherent element original relative rank) from small to large by element absolute value is to vector (n _{i_x}δ _{i_x}, n _{i_y}δ _{i_y}, n _{i_z}δ _{i_z}) element sequence, the new vector identification of reordering is (n _{i_a}δ _{i_a}, n _{i_b}δ _{i_b}, n _{i_c}δ _{i_c}) (in this expression formula, a, b, c are different and all belong to set { x, y, z}, such as, if abs (n _{i_x}δ _{i_x})≤abs (n _{i_y}δ _{i_y})≤abs (n _{i_z}δ _{i_z}), then a=x, b=y, c=z);

C. former coordinate system x-y-z is transformed to by rotation, symmetry operation namely represent that following matrix is taken advantage of on the equal right side of the row vector in position, direction in former coordinate system:

$A=\left(\begin{array}{ccc}\frac{\mathrm{abs}\left(n_{i\_x}\right)}{n_{i\_x}}(a==x)& \frac{\mathrm{abs}\left(n_{i\_x}\right)}{n_{i\_x}}(a==y)& \frac{\mathrm{abs}\left(n_{i\_x}\right)}{n_{i\_x}}(a==z)\\ \frac{\mathrm{abs}\left(n_{i\_y}\right)}{n_{i\_y}}(b==x)& \frac{\mathrm{zbs}\left(n_{i\_y}\right)}{n_{i\_y}}(b==y)& \frac{\mathrm{abs}\left(n_{i\_y}\right)}{n_{i\_y}}(b==z)\\ \frac{\mathrm{abs}\left(n_{i\_z}\right)}{n_{i\_z}}(c==x)& \frac{\mathrm{abs}\left(n_{i\_y}\right)}{n_{i\_z}}(c==y)& \frac{\mathrm{abs}\left(n_{i\_z}\right)}{n_{i\_z}}(c==z)\end{array}\right),$

In each element calculating formula of matrix, p==q (p=a, b, c, q=x, y, z) is logical operation, and it is 1 that equation sets up operation result, otherwise operation result is 0, and regulation works as n _{i_q}when=0 after coordinate transform, the normal vector at interface is computing grid is followed successively by δ along the length of side of new coordinate three change in coordinate axis direction _{i_a}, δ _{i_b}, δ _{i_c};

D. classify according to table 1 pair phase interface pattern, in table, c=min (c _{i_g}, c _{i_l}), v=δ _{i_a}δ _{i_b}δ _{i_c}for rectangular parallelepiped grid volume, n _a, n _b, n _crepresent abs (n respectively _{i_a}), abs (n _{i_b}), abs (n _{i_c}), in table, each symbol all refers in particular to variable corresponding to i-th computing grid in addition, in order to express easily omits subscript i_, such as δ _arepresent δ _{i_a};

E. the type belonging to phase interface and table 2 method determination phase interface plane equation n _{i_a}x+n _{i_b}y+n _{i_c}undetermined constant d in z=d, in table, each symbol all refers in particular to variable corresponding to i-th computing grid in addition, in order to express easily omits subscript i_, such as δ _arepresent δ _{i_a};

F. the calculating formula of phase contact area determined by the type belonging to phase interface, phase plane equation and table 3, calculates the interfacial area A of this grid _i, in table, each symbol all refers in particular to variable corresponding to i-th computing grid in addition, in order to express easily omits subscript i_, such as δ _arepresent δ _{i_a};

G. total phase contact area is by cumulative for the phase contact area of grids all in region to be analyzed,

The criterion of five class basic interfaces in table 1 rectangular parallelepiped

The calculating formula of table 2. five class phase interface plane equation constant d

The calculating formula of table 3 five class phase contact area

Adopt the test experience that the present embodiment method is carried out.

These interfacial area computing method are applied to during copper matte regulus PS converter copper making period converting process numerical simulation result analyzes, obtain gas liquid film and amass, thus calculate and to blast in gas oxygen concentration and sulfur dioxide concentration with the variation relation of melt height.Calculating oxygen utilization rate is thus 93.2%.This result and bibliographical information and actual measured results (95% ~ 98%) basically identical, this can confirm in this paperly to extract based on VOF numerical simulation the accuracy that gas liquid film amasss method to a certain extent.

Claims (1)

1. the gas-liquid two-phase interfacial area computing method based on numerical simulation, it is characterized in that: the two-phase flow VOF numerical simulation result calculating Zone Full adopting rectangular parallelepiped grid to divide according to all or part of zoning or the phase contact area of appointed area, its step is as follows:

A. read each grid values analog result data in region to be analyzed, comprise grid element center coordinate (x _i, y _i, z _i), the grid length of side (δ _{i_x}, δ _{i_y}, δ _{i_z}), gas, liquid volume of fluid (c in grid _{i_g}, c _{i_l}) and the gradient vector of arbitrary volume of fluid i in aforementioned each variable represents i-th grid, and x, y, z represent three change in coordinate axis direction of three-dimensional cartesian coordinate system, and g, l represent gas phase and liquid phase respectively;

B. element absolute value order is from small to large pressed to vector (n _{i_x}δ _{i_x}, n _{i_y}δ _{i_y}, n _{i_z}δ _{i_z}) element sequence, if there is equal situation, then maintain coherent element original relative rank, the new vector identification of reordering is (n _{i_a}δ _{i_a}, n _{i_b}δ _{i_b}, n _{i_c}δ _{i_c}), wherein, a, b, c are different and all belong to set { x, y, z};

C. former coordinate system x-y-z is transformed to by rotation, symmetry operation namely represent that following matrix is taken advantage of on the equal right side of the row vector in position, direction in former coordinate system:

$\left(\begin{array}{ccc}\frac{\mathrm{abs}\left(n_{i\_x}\right)}{n_{i\_x}}(a==x)& \frac{\mathrm{abs}\left(n_{i\_x}\right)}{n_{i\_x}}(a==y)& \frac{\mathrm{abs}\left(n_{i\_x}\right)}{n_{i\_x}}(a==z)\\ \frac{\mathrm{abs}\left(n_{i\_y}\right)}{n_{i\_y}}(b==x)& \frac{\mathrm{abs}\left(n_{i\_y}\right)}{n_{i\_y}}(b==y)& \frac{\mathrm{abs}\left(n_{i\_y}\right)}{n_{i\_y}}(b==z)\\ \frac{\mathrm{abs}\left(n_{i\_z}\right)}{n_{i\_z}}(c==x)& \frac{\mathrm{abs}\left(n_{i\_z}\right)}{n_{i\_z}}(c==y)& \frac{\mathrm{abs}\left(n_{i\_z}\right)}{n_{i\_z}}(c==z)\end{array}\right),$

In each element calculating formula of matrix, p==q is logical operation, and it is 1 that equation sets up operation result, otherwise operation result is 0, wherein p=a, b, c, q=x, y, z, and regulation works as n _{i_q}when=0 after coordinate transform, the normal vector at interface is computing grid is followed successively by δ along the length of side of new coordinate three change in coordinate axis direction _{i_a}, δ _{i_b}, δ _{i_c};

The condition of D. classifying to phase interface pattern is as follows:

type1： $c\≤\frac{{\left(n_a{\mathrm{\δ}}_a\right)}^3}{6n_a{n}_b{n}_{c}v}$

type2： $\frac{{\left(n_a{\mathrm{\&delta;}}_a\right)}^3}{6n_a{n}_b{n}_{c}v}<c\&le;\frac{{\left(n_b{\mathrm{\&delta;}}_b\right)}^3-{(n_b{\mathrm{\&delta;}}_b-n_a{\mathrm{\&delta;}}_a)}^3}{6n_a{n}_b{n}_{c}v}$

type3： $\left\{\begin{array}{ccc}{n}_a{\mathrm{\&delta;}}_a+n_b{\mathrm{\&delta;}}_b\&GreaterEqual;n_c{\mathrm{\&delta;}}_c& \mathrm{and}& \frac{{\left(n_b{\mathrm{\&delta;}}_b\right)}^3-{(n_b{\mathrm{\&delta;}}_b-n_a{\mathrm{\&delta;}}_a)}^3}{6n_a{n}_b{n}_{c}v}<c\&le;\frac{{\left(n_c{\mathrm{\&delta;}}_c\right)}^3-{(n_c{\mathrm{\&delta;}}_c-n_a{\mathrm{\&delta;}}_a)}^3-{(n_c{\mathrm{\&delta;}}_c-n_b{\mathrm{\&delta;}}_b)}^3}{6n_a{n}_b{n}_{c}v}\\ & \mathrm{or}& \\ n_a{\mathrm{\&delta;}}_a+n_b{\mathrm{\&delta;}}_b<n_c{\mathrm{\&delta;}}_c& \mathrm{and}& \frac{{\left(n_b{\mathrm{\&delta;}}_b\right)}^3-{(n_b{\mathrm{\&delta;}}_b-n_a{\mathrm{\&delta;}}_a)}^3}{6n_a{n}_b{n}_{c}v}<c\&le;\frac{{(n_a{\mathrm{\&delta;}}_a+n_b{\mathrm{\&delta;}}_b)}^3-{\left(n_a{\mathrm{\&delta;}}_a\right)}^3-{\left(n_b{\mathrm{\&delta;}}_b\right)}^3}{6n_a{n}_b{n}_{c}v}\end{array}\right.$

type4： $\left\{\begin{array}{c}{n}_a{\mathrm{\&delta;}}_a+n_b{\mathrm{\&delta;}}_b\&GreaterEqual;n_c{\mathrm{\&delta;}}_c\\ \mathrm{and}\\ \frac{{\left(n_c{\mathrm{\&delta;}}_c\right)}^3-{(n_c{\mathrm{\&delta;}}_c-n_a{\mathrm{\&delta;}}_a)}^3-{(n_c{\mathrm{\&delta;}}_c-n_b{\mathrm{\&delta;}}_b)}^3}{6n_a{n}_b{n}_{c}v}<c\&le;\frac{{(n_a{\mathrm{\&delta;}}_a+n_b{\mathrm{\&delta;}}_b)}^3-{\left(n_a{\mathrm{\&delta;}}_a\right)}^3-{\left(n_b{\mathrm{\&delta;}}_b\right)}^3-{(n_a{\mathrm{\&delta;}}_a+n_b{\mathrm{\&delta;}}_b-n_c{\mathrm{\&delta;}}_c)}^3}{6n_a{n}_b{n}_{c}v}\end{array}\right.$

type5： $\left\{\begin{array}{c}{n}_a{\mathrm{\&delta;}}_a+n_b{\mathrm{\&delta;}}_b<n_c{\mathrm{\&delta;}}_c\\ \mathrm{and}\\ \frac{{(n_a{\mathrm{\&delta;}}_a+n_b{\mathrm{\&delta;}}_b)}^3-{\left(n_a{\mathrm{\&delta;}}_a\right)}^3-{\left(n_b{\mathrm{\&delta;}}_b\right)}^3}{6n_a{n}_b{n}_{c}v}<c\&le;\frac{{\left(n_c{\mathrm{\&delta;}}_c\right)}^3-{(n_c{\mathrm{\&delta;}}_c-n_a{\mathrm{\&delta;}}_a)}^3-{(n_c{\mathrm{\&delta;}}_c-n_b{\mathrm{\&delta;}}_b)}^3+{(n_c{\mathrm{\&delta;}}_c-n_a{\mathrm{\&delta;}}_a-n_b{\mathrm{\&delta;}}_b)}^3}{6n_a{n}_b{n}_{c}v}\end{array}\right.$

Above-mentioned various in, c=min (c _{i_g}, c _{i_l}), v=δ _{i_a}δ _{i_b}δ _{i_c}for rectangular parallelepiped grid volume, n _a, n _b, n _crepresent abs (n respectively _{i_a}), abs (n _{i_b}), abs (n _{i_c}), in formula, each symbol all refers in particular to variable corresponding to i-th computing grid in addition, in order to express easily omits subscript i_;

E. the type belonging to phase interface and following method determination phase interface plane equation n _{i_a}x+n _{i_b}y+n _{i_c}undetermined constant d in z=d:

type1：d ³＝6n _an _bn _ccv and d≤n _a

type2：d ³-(d-n _aδ _a) ³＝6n _an _bn _ccv and n _a＜d≤n _b

type3：d ³-(d-n _aδ _a) ³-(d-n _bδ _b) ³＝6n _an _bn _ccv and n _b＜d and(d≤n _cor d≤n _a+n _b)

type4：d ³-(d-n _aδ _a) ³-(d-n _bδ _b) ³-(d-n _cδ _c) ³＝6n _an _bn _ccv and n _c＜d≤n _a+n _b

type5：d ³-(d-n _aδ _a) ³-(d-n _bδ _b) ³+(d-n _aδ _a-n _bδ _b) ³＝6n _an _bn _ccv and n _a+n _b＜d≤n _c

In formula, each symbol all refers in particular to variable corresponding to i-th computing grid, in order to express easily omits subscript i_;

The calculating formula A of the type F. belonging to phase interface, phase plane equation and following method determination phase contact area _i:

type1： $A_i=\frac{\sqrt{n_a^2+n_b^2+n_c^2}}{2n_a{n}_b{n}_c}\&CenterDot;d^2$

type2： $A_i=\frac{\sqrt{n_a^2+n_b^2+n_c^2}}{2n_a{n}_b{n}_c}\&CenterDot;[d^2-{(d-n_a{\mathrm{\&delta;}}_a)}^2]$

type3： $A_i=\frac{\sqrt{n_a^2+n_b^2+n_c^2}}{2n_a{n}_b{n}_c}\&CenterDot;[d^2-{(d-n_a{\mathrm{\&delta;}}_a)}^2-{(d-n_b{\mathrm{\&delta;}}_b)}^2]$

type4： $A_i=\frac{\sqrt{n_a^2+n_b^2+n_c^2}}{2n_a{n}_b{n}_c}\&CenterDot;[d^2-{(d-n_a{\mathrm{\&delta;}}_a)}^2-{(d-n_b{\mathrm{\&delta;}}_b)}^2-{(d-n_c{\mathrm{\&delta;}}_c)}^2]$

type5： $A_i=\frac{\sqrt{n_a^2+n_b^2+n_c^2}}{2n_a{n}_b{n}_c}\&CenterDot;[d^2-{(d-n_a{\mathrm{\&delta;}}_a)}^2-{(d-n_b{\mathrm{\&delta;}}_b)}^2+{(d-n_c{\mathrm{\&delta;}}_c-n_b{\mathrm{\&delta;}}_b)}^2]$

In formula, each symbol all refers in particular to variable corresponding to i-th computing grid, in order to express easily omits subscript i_;

G. total phase contact area is by cumulative for the phase contact area of grids all in region to be analyzed,