CN113177373B - Fluid distribution calculation method based on VOF principle - Google Patents

Fluid distribution calculation method based on VOF principle Download PDF

Info

Publication number
CN113177373B
CN113177373B CN202110467994.8A CN202110467994A CN113177373B CN 113177373 B CN113177373 B CN 113177373B CN 202110467994 A CN202110467994 A CN 202110467994A CN 113177373 B CN113177373 B CN 113177373B
Authority
CN
China
Prior art keywords
omega
vof
unit
interface
fluid
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202110467994.8A
Other languages
Chinese (zh)
Other versions
CN113177373A (en
Inventor
李超
段文洋
赵彬彬
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Harbin Engineering University
Original Assignee
Harbin Engineering University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Harbin Engineering University filed Critical Harbin Engineering University
Priority to CN202110467994.8A priority Critical patent/CN113177373B/en
Publication of CN113177373A publication Critical patent/CN113177373A/en
Application granted granted Critical
Publication of CN113177373B publication Critical patent/CN113177373B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00—Details relating to the application field
    • G06F2113/08—Fluids
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Fluid Mechanics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Computing Systems (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • Algebra (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

The invention belongs to the technical field of CFD (computational fluid dynamics), and particularly relates to a fluid distribution calculation method based on a VOF (Voltage over fiber) principle. The invention can be suitable for any polyhedral mesh, overcomes the defect that the initial value of the VOF of a CFD mesh unit can not be accurately calculated in the prior art, and effectively improves the calculation efficiency and saves the memory of a computer by dividing an interface unit into a plurality of sub-tetrahedral or sub-triangular units and calculating the VOF value of the interface subunit by means of an equal division ratio by adopting multi-stage 'division-judgment'. The method can analyze the error of the VOF initialization, if the error does not meet the preset precision, the maximum segmentation-judgment level needs to be improved, iterative calculation is carried out until the relevant conditions are met, and the calculation parameters are adjusted through error feedback, so that the high-precision calculation of the VOF initialization is realized.

Description

Fluid distribution calculation method based on VOF principle
Technical Field
The invention belongs to the technical field of CFD (computational fluid dynamics), and particularly relates to a fluid distribution calculation method based on a VOF (Voltage over fiber) principle.
Background
In the case of numerical simulation of a complementary miscible multiphase Fluid using cfd (computational Fluid dynamics) technology, the vof (volume of Fluid) method is widely used for capturing a motion interface; in the VOF method, the exact fluid interface position is replaced by a discrete fluid volume fraction VOF value, and therefore the fluid volume fraction VOF needs to be initialized before solving. Since the CFD technique discretizes a continuous physical computation field into a grid, which consists of a series of cells, the VOF initialization is the calculation of the fractional volume of fluid in each cell.
Consider the computational domain Ω - Ω formed by two immiscible fluids a and BA∪ΩBAs shown in fig. 1, and defining a fluid indicator function H (i.e., fluid interface expression) as follows:
Figure BDA0003044073160000011
where x represents a position vector, then Ω is the grid celliVolume fraction of inner, fluid A, VOFiIs defined as:
Figure BDA0003044073160000012
so VOFiThe physical meaning of (A) is that the fluid A is in the unit omegaiVolume ratio (area ratio in 2D); and a unit omegaiDensity of internal composite fluid ρiKinetic viscosity μiIsofluid properties are according to VOFiThe values of (a) are calculated as:
ρi＝ρAVOFi+ρB(1-VOFi)
μi＝μAVOFi+μB(1-VOFi)
the density value can influence flow field variables such as pressure and the like, so that an accurate VOF initial value has a positive promotion effect on ensuring the accuracy, initial convergence and the like of CFD solving.
For a given fluid interface (e.g. region Ω)AThe boundary surface of (b) is referred to as a full cell, the cell completely inside the boundary surface is referred to as an empty cell, and the other cells are referred to as interface cells, and therefore:
Figure BDA0003044073160000013
for a regular cell, the initialization of the VOF is easier, but for an irregular polyhedral cell, a specific method is required to accurately calculate the initial value of the VOF.
In mainstream CFD software, the VOF initialization is usually simply processed, or only assigned according to a full unit or an empty unit, or the initialization precision is not high, so that an obvious numerical error may be caused under a specific condition; as shown in fig. 2(a), when the theoretical volume fraction value of the interface element is exactly 0.5, in the commercial software StarCCM +, only the result shown in fig. 2(b) can be obtained by means of a simple user-defined function and the element centroid position, i.e., the VOF of the element has only 0 or 1 value, while in the open-source software OpenFOAM, the result shown in fig. 2(c) can be obtained by means of a setFields function, i.e., the VOF of the element also has only 0 or 1 value, and jagged jumps can occur at the interface; of course, writing advanced user-defined functions in different software can provide accurate or reasonable approximations of the VOF for the regular-shaped grid shown in FIG. 2(a), but it is difficult for any polyhedral cell to be irregular in shape.
Disclosure of Invention
The invention aims to overcome the defect that the VOF initial value of a CFD grid unit cannot be accurately calculated in the prior art, and provides a high-precision VOF initial value for CFD solution so as to better ensure the accuracy and the initial convergence of the CFD solution.
The purpose of the invention is realized by the following technical scheme: the method comprises the following steps:
step 1: acquiring a CFD grid of computational domains in which two immiscible fluids A and B are distributed; preprocessing the CFD mesh of the computational domain, perfecting the topological relations between units and nodes, between units and surfaces, between surfaces and nodes, and calculating the mass center of unit bodies and unit surfaces;
step 2: reading in an expression H of a fluid interface of two immiscible fluids A and B in a CFD grid of a computational domain;
and step 3: for each unit omega in the CFD grid of the computational domainiLabeling and calculating the fluid volume fraction thereof;
if unit omegaiAll located within the fluid interface, unit omega is formediMarking as a full cell; if unit omegaiAll located outside the fluid interface, unit omega is formediMarking as an empty cell; otherwise, the unit omega isiMarking as interface unit, obtaining total number N of interface unit in CFD grid of calculation domainjm；
Unit omegaiFluid volume fraction of (VOF)iRepresenting the fluid volume distribution ratio of immiscible fluids A and B; the fluid volume fraction of the full unit is VOF i1 denotes the unit ΩiOnly one fluid a is distributed; the fluid volume fraction of the empty unit is VOF i0 denotes the unit ΩiOnly one fluid B is distributed;
for the interface unit ΩiFluid volume fraction of VOFiThe calculation method comprises the following steps:
step 3.1: interface unit omegaiDivided into subunits omegaij；
If the CFD grid of the computational domain is a 3D grid, the 3D interface unit omega is usediDivided into tetrahedral subunits omegaij(ii) a If 3D interface unit omegaiIncluding a polygonal surface, the polygonal surface is divided into a plurality of triangular surfaces, and then the 3D interface unit omega is formediThe division into tetrahedral subunits omega according to the centroid and the respective triangular surfacesij；
If the CFD grid of the computational domain is a 2D grid, the computational domain is divided into omega units according to the 2D interfaceiThe centroid and each edge node of (2D) interface unit omegaiDivided into triangular subunits omegaij；
Step 3.2: subunit omegaijLabeling and calculating the fluid volume fraction thereof;
if subunit ΩijAll the nodes are located in the fluid interface, the subunit omega is connectedijMarking as a full subunit; if subunit ΩijAll located outside the fluid interface, the subunit omega is connectedijMarking as a null subunit; otherwise, the subunit omega is connectedijMarking as an interface subunit;
the fluid volume fraction of the full subunit is VOF ij1 is ═ 1; the fluid volume fraction of the void cell is VOF ij0; for the interface subunit ΩijFluid volume fraction of VOFijThe calculation method comprises the following steps:
step 3.2.1: setting a maximum segmentation-judgment level R, a volume equal-dividing ratio of a 3D unit or an area equal-dividing ratio m of a 2D unit; initialization r 1, subunit ΩijIs a divided target unit;
step 3.2.2: equally dividing each divided target unit into m subunits, marking all subunits, and acquiring the number N of full subunits generated by divisionijrFinding out all interface subunits;
step 3.2.3: if R is less than R, all interface subunits generated by segmentation in the segmentation-judgment level R are used as target units segmented in the next segmentation-judgment level R +1, R is made to be R +1, and the step is returned to 3.2.2;
step 3.2.4: compute interface subunit omegaijFluid volume fraction (VOF) ofij；
Figure BDA0003044073160000031
Step 3.3: computing interface Unit omegaiFluid volume fraction of (VOF)i；
Figure BDA0003044073160000032
Wherein, | ΩiL represents the volume of a 3D cell or the area of a 2D cell; n is a radical ofsubIs an interface unit omegaiDivided subunit omegaijThe number of (c);
and 4, step 4: each unit omega in CFD grid of output calculation domainiFluid volume fraction of (VOF)iThe distribution of the two immiscible fluids a and B in the calculated domain is obtained.
The present invention may further comprise:
outputting each unit omega in the CFD grid of the calculation domain in the step 4iFluid volume fraction (VOF) ofiWhether the preset precision is met or not needs to be checked in the prior art, and the specific method comprises the following steps:
step 4.1: calculating the 3D volume or 2D area V of the fluid contained in all the interface elementsjmRelative error between two adjacent segmentation-decision levels k and k-1
Figure BDA0003044073160000033
Figure BDA0003044073160000034
Figure BDA0003044073160000041
Figure BDA0003044073160000042
Wherein k is more than or equal to 2 and less than or equal to R;
Figure BDA0003044073160000043
represents the u interface unit;
Figure BDA0003044073160000044
as an interface unit
Figure BDA0003044073160000045
Fluid volume fraction in "segmentation-decision" level k;
Figure BDA0003044073160000046
as an interface unit
Figure BDA0003044073160000047
Subunit omega ofujFluid volume fraction in "segmentation-decision" level k;
Figure BDA0003044073160000048
step 4.2: calculating the fluid 3D volume or 2D area V contained by all cells in the CFD gridallRelative error between two adjacent segmentation-decision levels k and k-1
Figure BDA0003044073160000049
Figure BDA00030440731600000410
Figure BDA00030440731600000411
Wherein N iscellTo calculate the number of all cells contained in the CFD mesh of the domain;
step 4.3: judgment of
Figure BDA00030440731600000412
And
Figure BDA00030440731600000413
whether the preset precision is met or not; if it is
Figure BDA00030440731600000414
Or
Figure BDA00030440731600000415
If the preset precision is not met, increasing the value of the maximum segmentation-judgment level R, and returning to the step 3; if it is
Figure BDA00030440731600000416
And
Figure BDA00030440731600000417
all meet the preset precision, and then each unit omega in the CFD grid of the calculation domain is outputiFluid volume fraction of (VOF)iThe distribution of the two immiscible fluids a and B in the calculated domain is obtained.
The invention has the beneficial effects that:
(1) the method overcomes the defect that the VOF initial value of the CFD grid unit cannot be accurately calculated in the prior art, and can be suitable for any polyhedral grid;
(2) the invention adopts multi-level 'segmentation-judgment' to calculate the VOF value of the interface related unit by means of an equal proportion, thereby effectively improving the calculation efficiency and saving the memory of a computer;
(3) the method can analyze the error of the VOF initialization, can adjust the calculation parameters through error feedback, and realizes the high-precision calculation of the VOF initialization.
Drawings
Fig. 1 is a schematic diagram of the physical significance of a CFD mesh and its fractional volume of fluid per cell, VOF.
Fig. 2(a) is a theoretical (or target) VOF cloud when testing mainstream CFD software.
FIG. 2(b) is an initialization cloud of the commercial software Star-CCM +.
Fig. 2(c) is an initialization cloud of the open source software OpenFOAM.
Fig. 3 is an overall flow chart of the present invention.
FIG. 4 shows an interface unit ΩiDivided into several subunits omegaijSchematic representation of (a).
FIG. 5 is a schematic diagram of a plane-partitioned tetrahedron formed by connecting edges of the tetrahedron.
FIG. 6 is a schematic diagram of a triangle being divided by connecting the middle points of the triangle.
Fig. 7 is a schematic diagram of stepwise division of tetrahedrons.
Fig. 8 is a schematic diagram of stepwise division of triangles.
Fig. 9(a) is a 3D view of a CFD numerical sink grid.
Fig. 9(b) is a schematic diagram of a CFD numerical sink grid side view and a target water plane (i.e., water-air interface).
Fig. 9(c) is the VOF initial value for water at the center of a grid cell calculated by the mainstream open source software OpenFOAM.
FIG. 9(d) is the VOF initial value calculated by the present invention for water at the center of the grid.
FIG. 9(e) is a 0.5 iso-surface side view of the water VOF values calculated by the OpenFOAM software.
FIG. 9(f) is a 0.5 iso-surface side view of the water VOF values calculated by the present invention.
Fig. 9(g) shows the VOF initialization error of the present invention.
Fig. 10 is a table of VOF initialization error forms in the present invention.
Detailed Description
The invention is further described below with reference to the accompanying drawings.
The invention provides a fluid distribution calculation method based on a VOF principle, which can be used for any polyhedral mesh. Reading a CFD grid file; reading an expression of a fluid interface; judging the position of the unit and marking according to the classification of full, empty and interface, dividing the interface unit into a plurality of sub-tetrahedral or sub-triangular units, judging the position of the sub-unit and marking according to the classification of full, empty and interface; the VOF value of the interface subunit is calculated by adopting multi-stage 'segmentation-judgment' by means of an equal division ratio, so that the calculation efficiency is effectively improved, and the memory of a computer is saved; taking the volume or area weighted average of all the subunits of the interface unit as the VOF value of the interface unit; the unit VOF values of type full and empty are 1 and 0, respectively; after traversing all the units, the VOF initialization error needs to be calculated, if the error does not meet the preset precision, the maximum segmentation-judgment level needs to be improved, iterative calculation is carried out until the relevant conditions are met, and the high-precision VOF initial value can be provided on any polyhedral grid.
A fluid distribution calculation method based on a VOF principle comprises the following steps:
step 1: acquiring a CFD grid of computational domains in which two immiscible fluids A and B are distributed; preprocessing the CFD mesh of the computational domain, perfecting the topological relations between units and nodes, between units and surfaces, between surfaces and nodes, and calculating the mass center of unit bodies and unit surfaces;
step 2: reading in an expression H of a fluid interface of two immiscible fluids A and B in a CFD grid of a computational domain;
and 3, step 3: for each cell Ω in the CFD grid of the computational domainiLabeling and calculating the fluid volume fraction thereof;
if unit omegaiAll located within the fluid interface, unit ΩiMarking as a full cell; if unit omegaiAll located outside the fluid interface, unit omega is formediMarking as an empty cell; otherwise, the unit omega isiMarking as interface unit, obtaining total number N of interface unit in CFD grid of calculation domainjm；
The fluid volume fraction of the full unit is VOF i1 is ═ 1; the fluid volume fraction of the empty cell is VOF i0; unit omegaiFluid volume fraction of (VOF)iRepresenting the fluid volume distribution ratio of immiscible fluids A and B; VOF i1 denotes the unit ΩiOnly one fluid a is distributed; VOF i0 denotes the unit ΩiOnly one fluid B is distributed;
for the interface unit omegaiFluid volume fraction of VOFiThe calculation method comprises the following steps:
step 3.1: interface unit omegaiDivided into subunits omegaij；
If the CFD grid of the computational domain is a 3D grid, the 3D interface unit omega is usediDivided into tetrahedral subunits omegaij(ii) a If 3D interface unit omegaiIncluding polygonal surface, the polygonal surface is divided into several triangular surfaces, and then the 3D interface unit omega is formediThe division into tetrahedral subunits omega according to the centroid and the respective triangular surfacesij；
If the CFD grid of the computational domain is a 2D grid, the computational domain is divided into omega units according to the 2D interfaceiThe centroid and each edge node of (2D) interface unit omegaiDivided into triangular subunits omegaij；
Step 3.2: subunit omegaijLabeling and calculating the fluid volume fraction thereof;
if subunit ΩijAll nodes are located within the fluid interface, then willSubunit ΩijMarking as a full subunit; if subunit ΩijAll located outside the fluid interface, the subunit omega is connectedijMarking as a null subunit; otherwise, the subunit omega is connectedijMarking as an interface subunit;
the fluid volume fraction of the full subunit is VOF ij1; the fluid volume fraction of the void cell is VOF ij0; for the interface subunit ΩijFluid volume fraction of VOFijThe calculation method comprises the following steps:
step 3.2.1: setting a maximum segmentation-judgment level R, a volume equal-dividing ratio of a 3D unit or an area equal-dividing ratio m of a 2D unit; initialization r is 1, subunit ΩijIs a divided target unit;
step 3.2.2: equally dividing each divided target unit into m subunits, marking all subunits, and acquiring the number N of full subunits generated by divisionijrFinding out all interface subunits;
step 3.2.3: if R is less than R, all interface subunits generated by segmentation in the segmentation-judgment level R are used as target units segmented in the next segmentation-judgment level R +1, R is made to be R +1, and the step is returned to 3.2.2;
step 3.2.4: compute interface subunit omegaijFluid volume fraction (VOF) ofij；
Figure BDA0003044073160000071
Step 3.3: computing interface Unit omegaiFluid volume fraction (VOF) ofi；
Figure BDA0003044073160000072
Wherein, | ΩiL represents the volume of a 3D cell or the area of a 2D cell; n is a radical ofsubIs an interface unit omegaiThe divided sub-unit omegaijThe number of (2);
and 4, step 4: calculating the 3D volume or 2D area V of the fluid contained in all the interface elementsjmRelative error between two adjacent segmentation-decision levels k and k-1
Figure BDA0003044073160000073
Figure BDA0003044073160000074
Figure BDA0003044073160000075
Figure BDA0003044073160000076
Wherein k is more than or equal to 2 and less than or equal to R;
Figure BDA0003044073160000077
represents the u interface unit;
Figure BDA0003044073160000078
as an interface unit
Figure BDA0003044073160000079
Fluid volume fraction in "segmentation-decision" level k;
Figure BDA00030440731600000710
as an interface unit
Figure BDA00030440731600000711
Subunit omega ofujFluid volume fraction in "segmentation-decision" level k;
Figure BDA00030440731600000712
and 5: computing all cell inclusions in a CFD meshOf a fluid 3D volume or 2D area VallRelative error between two adjacent segmentation-decision levels k and k-1
Figure BDA00030440731600000713
Figure BDA00030440731600000714
Figure BDA00030440731600000715
Wherein N iscellTo calculate the number of all cells contained in the CFD mesh of the domain;
and 6: judgment of
Figure BDA0003044073160000081
And with
Figure BDA0003044073160000082
Whether the preset precision is met or not; if it is
Figure BDA0003044073160000083
Or
Figure BDA0003044073160000084
If the preset precision is not met, increasing the value of the maximum segmentation-judgment level R, and returning to the step 3; if it is
Figure BDA0003044073160000085
And with
Figure BDA0003044073160000086
All meet the preset precision, and then each unit omega in the CFD grid of the calculation domain is outputiFluid volume fraction of (VOF)iThe distribution of the two immiscible fluids a and B in the calculated domain is obtained.
The invention is not limited to solving for the distribution of two immiscible fluids a and B in a computational domain, but can also be used to solve for the distribution of each fluid in a computational domain when there are multiple immiscible fluids in the computational domain. The distribution of the fluid A in the calculation domain is firstly obtained by the invention, then one fluid is taken from the fluid B as A, and the distribution of the fluid A in the calculation domain is calculated again. And repeating the steps until the distribution of all the fluids in the calculation domain is obtained.
Compared with the prior art, the invention has the beneficial effects that:
(1) the method overcomes the defect that the VOF initial value of the CFD grid unit cannot be accurately calculated in the prior art, and can be suitable for any polyhedral grid;
(2) the invention adopts multi-level 'segmentation-judgment' to calculate the VOF value of the interface related unit by means of an equal proportion, thereby effectively improving the calculation efficiency and saving the memory of a computer;
(3) the method can carry out error analysis on the VOF initialization, can adjust the calculation parameters through error feedback, and realizes high-precision calculation of the VOF initialization.
Example 1:
step 1, reading in a CFD (computational fluid dynamics) grid file, such as an Ansys Fluent format msh grid file, wherein the geometrical topological relation contained in the grid file is usually incomplete, and data of unit bodies or plane centroids do not exist, so that the read CFD grid file needs to be preprocessed, namely, the topological relations between units and nodes, between units and surfaces, and between surfaces and nodes are perfected, and the centroids of the unit bodies and the unit surfaces are calculated;
step 2, reading an expression H of the fluid interface, wherein the expression is a step function related to a space coordinate x, namely when the space point coordinate is positioned in the fluid interface (including a boundary), a function value is 1, otherwise, the function value is 0; solving the most common fluid region shape for CFD: the cuboid can be obtained by reading in the boundary coordinates of the cuboid: x is a radical of a fluorine atommin，xmax，ymin，ymax，zmin，zmaxThe reading in of H is realized;
Figure BDA0003044073160000087
step 3, judging the ith unit omega in the gridiIf all the unit nodes are positioned in the fluid interface plane, the unit nodes are marked as full units, if all the unit nodes are positioned outside the fluid interface plane, the unit nodes are marked as empty units, and if all the unit nodes are positioned outside the fluid interface plane, the unit nodes are marked as interface units;
step 4, if omegaiVOF of full cell i1, if ΩiIs void cell then VOFi＝0；
Step 5, if omegaiDividing the interface unit into a plurality of sub-units omegaij(j＝1，2，...，Nsub) Omega in 3D gridijIs tetrahedral, omega in a 2D gridijIs triangular; as shown in FIG. 4, the polyhedral cell Ω is formed according to the centroid and the side nodes of the 2D celliDividing the triangle into a plurality of triangle subunits; if the 3D cell includes a polygonal surface, the polygonal surface is first divided into a plurality of triangular surfaces according to the method shown in fig. 4, and then the 3D cell is divided into a plurality of tetrahedrons according to the body centroid and each triangular surface.
Step 6, judging subunit omegaijPosition of (d) if ΩijVOF if the subunit is fullij1, if ΩijIs a void sub-unit VOFij＝0；
Step 7, if omegaijIs an interface subunit, for omegaijVOF calculation by means of aliquot using multi-level' segmentation-judgmentijThe values, namely:
Figure BDA0003044073160000091
where R denotes a maximum "segmentation-judgment" level, R is usually set to 5, and m is a volume equivalence ratio of a 3D unit or an area equivalence ratio of a 2D unit, and a specific segmentation manner may be, but is not limited to, that shown in fig. 5(m is 8) and fig. 6(m is 4); in each level r of the "division-judgment", each target unit is equally divided into m sub-unitsThe unit judges the positions of all newly generated subunits and marks the subunits according to the classification of full, empty and interface, and only takes all interface subunits in the level r as the divided target units in the next level r +1, and the specific operation is shown in fig. 7(3D) and fig. 8(2D), so that the calculation efficiency can be improved, and the computer memory can be saved; when r is 1, the target unit is Ωij，NijrIs omegaijThe number of full subunits generated by the new segmentation in level r.
Step 8, calculating all VOFsij(j＝1，2，...，Nsub) As a unit omega, a 3D volume (or 2D area) weighted average ofiVOF ofiValues, i.e. using the formula:
Figure BDA0003044073160000092
step 9, repeating the steps 3 to 8 until all the units omega are traversedi(i＝1，2，...，Ncel1)；
Step 10, outputting initial values VOF of fluid volume fractions of all unitsi(i＝1，2，...，Ncell)。
Calculating the initialized errors of the two VOFs, if the errors do not meet the preset precision, increasing the value of the maximum segmentation-judgment level R, and repeating the step 9 until the errors meet the requirement or the maximum number of times of R value modification is reached; wherein the first error is the volume (area in 2D) V of fluid contained by all interface elementsjmRelative error between k and k-1 in two adjacent segmentation-decision levels
Figure BDA0003044073160000093
Figure BDA0003044073160000101
Here, k is 2. ltoreq. k.ltoreq.R,
Figure BDA0003044073160000102
the following formula is used for calculation:
Figure BDA0003044073160000103
wherein the content of the first and second substances,
Figure BDA0003044073160000104
represents the unit of the i-th interface unit,
Figure BDA0003044073160000105
is composed of
Figure BDA0003044073160000106
Volume fraction value in order k:
Figure BDA0003044073160000107
Figure BDA0003044073160000108
is a subunit omegaijVolume fraction value in order k:
Figure BDA0003044073160000109
the second error is the volume (area in 2D) V of fluid contained by all cellsallRelative error between k and k-1 in two adjacent segmentation-decision levels
Figure BDA00030440731600001010
Figure BDA00030440731600001011
Figure BDA00030440731600001012
Likewise, k is more than or equal to 2 and less than or equal to R,
Figure BDA00030440731600001013
calculated by formula (5) (only need to be)
Figure BDA00030440731600001014
Replaced by omegai)；
Typically, the maximum "split-judge" level R is set to 5, with a preset accuracy of
Figure BDA00030440731600001015
If the predetermined accuracy is not satisfied when the level k is equal to R, then increase R by 1 or according to
Figure BDA00030440731600001016
And
Figure BDA00030440731600001017
increasing R straight with the change rule of k, and setting the maximum number of times of R value modification to be 3; in terms of output form, error
Figure BDA00030440731600001018
And
Figure BDA00030440731600001019
the output may be performed in a list format, which may be, but is not limited to, the format shown in fig. 10:
and the final VOF initial value of the CFD grid unit can be independently output as a binary file, and the format is as follows:
(VOFi，i＝1，2，...，Ncell)
the number i of the unit is the same as the number of the grid unit, and meanwhile, a VOF value file in a unit center type Tecplot format can be output, so that a VOF cloud map can be conveniently and visually checked.
To further illustrate the calculation effect of the present invention, consider a CFD numerical water tank as shown in fig. 9(a), wherein the water tank has a length of 5m, a height of 1m, a width of 0.1m, and a target waterline height of 0.5m (indicated by a red dotted line in fig. 9 (b)), the grid of the water tank is composed of 2492 cells (cell types include triangular prism and polyhedron), and the grid file is in Ansys Fluent format; for a water area in the water tank, the expression h (x) of the interface thereof may be read in a rectangular parallelepiped manner according to formula (9), where:
[xmin，xmax]×[ymin，ymax]×[zmin，zmax]＝[0，5]×[0，0.5]×[0，0.1]
the maximum "division-determination" level R is set to 5, and the unit division method shown in fig. 5 is adopted, i.e., the target tetrahedral unit is equally divided into 8 sub-tetrahedrons (m is 8), and the following equation (2) is provided:
Figure BDA0003044073160000111
this computational accuracy is equivalent to splitting the interface subunit ΩijUniformly dividing into m at one timeR(＝8532768) subordinate subunits, and because the calculation is performed by adopting multi-stage 'division-judgment', compared with one-time uniform division, the calculation efficiency is effectively improved, and the memory of the computer is saved; set the preset precision to
Figure BDA0003044073160000112
In order to compare the calculation effect of the present invention, the mainstream open source CFD software OpenFOAM and the present invention are used for comparison calculation, and the output VOF initial value file is imported into the visualization software ParaView, so that the OpenFOAM and the mesh unit VOF value of the present invention are obtained as shown in fig. 9(c) and fig. 9(d), respectively, and it can be seen that: the VOF value of OpenFOAM is only 0 and 1, and the fractional value of the interface unit can be calculated by the method; furthermore, in the VOF method, the fluid interface is generally represented by a 0.5 isosurface (or 2D isosurface) of the VOF value, and in this example, OpenFOAM and the calculated VOF of the present invention equal to 0.5 isosurface are shown in fig. 9(e) and 9(f), respectively, as is apparent: the equivalent surface obtained by OpenFOAM is very wrinkled and has very different true horizontal state from the water line surfaceLarge; the isosurface obtained by the invention is relatively very horizontal; therefore, the invention can provide a reasonable initial value of the VOF under the condition of poor grid quality, the error of the VOF is shown in figure 9(g), and the preset precision requirement is met; therefore, the invention can provide a high-precision VOF initial value for CFD solution, thereby more accurately calculating the variables such as density, viscosity, pressure and the like of the composite fluid in the flow field, and further better ensuring the accuracy and initial convergence of the CFD solution.
The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention, and various modifications and changes may be made by those skilled in the art. Any modification, equivalent replacement, or improvement made within the spirit and principle of the present invention should be included in the protection scope of the present invention.

Claims (2)

1. A fluid distribution calculation method based on a VOF principle is characterized by comprising the following steps:
step 1: acquiring a CFD grid of computational domains in which two immiscible fluids A and B are distributed; preprocessing the CFD mesh of the computational domain, perfecting the topological relations between units and nodes, between units and surfaces, between surfaces and nodes, and calculating the mass center of unit bodies and unit surfaces;
step 2: reading in an expression H of a fluid interface of two immiscible fluids A and B in a CFD grid of a computational domain;
and 3, step 3: for each cell Ω in the CFD grid of the computational domainiLabeling and calculating the fluid volume fraction thereof;
if unit omegaiAll located within the fluid interface, unit omega is formediMarking as a full cell; if unit omegaiAll located outside the fluid interface, unit omega is formediMarking as an empty cell; otherwise, the unit omega isiMarking as interface unit, obtaining total number N of interface unit in CFD grid of calculation domainjm；
Unit omegaiFluid volume fraction of (VOF)iIndicating the ratio of the fluid volume distributions of immiscible fluids A and BExample (c); the fluid volume fraction of the full unit is VOFi1 denotes the unit ΩiOnly one fluid a is distributed; the fluid volume fraction of the empty cell is VOFi0 denotes the unit ΩiOnly one fluid B is distributed;
for the interface unit omegaiFluid volume fraction of VOFiThe calculation method comprises the following steps:
step 3.1: interface unit omegaiDivided into subunits omegaij；
If the CFD grid of the computational domain is a 3D grid, the 3D interface unit omega is usediDivided into tetrahedral subunits omegaij(ii) a If 3D interface unit omegaiIncluding a polygonal surface, the polygonal surface is divided into a plurality of triangular surfaces, and then the 3D interface unit omega is formediThe division into tetrahedral subunits omega according to the centroid and the triangular surfacesij；
If the CFD grid of the computational domain is a 2D grid, the computational domain is divided into omega units according to the 2D interfaceiThe centroid and each edge node of (2D) interface unit omegaiDivided into triangular subunits omegaij；
Step 3.2: subunit omegaijLabeling and calculating the fluid volume fraction thereof;
if subunit ΩijAll the nodes are located in the fluid interface, the subunit omega is connectedijMarking as a full subunit; if subunit ΩijAll located outside the fluid interface, the subunit omega is connectedijMarking as a null subunit; otherwise, the subunit omega is connectedijMarking as an interface subunit;
the fluid volume fraction of the full subunit is VOFij1 is ═ 1; the fluid volume fraction of the void cell is VOFij0; for the interface subunit ΩijFluid volume fraction of VOFijThe calculation method comprises the following steps:
step 3.2.1: setting a maximum segmentation-judgment level R, a volume equal-dividing ratio of a 3D unit or an area equal-dividing ratio m of a 2D unit; initialization r 1, subunit ΩijIs a divided target unit;
step 3.2.2: equally dividing each divided target unit into m subunits, marking all subunits, and acquiring the number N of full subunits generated by divisionijrFinding out all interface subunits;
step 3.2.3: if R is less than R, all interface subunits generated by segmentation in the segmentation-judgment level R are taken as target units segmented in the next segmentation-judgment level R +1, R is made equal to R +1, and the process returns to the step 3.2.2;
step 3.2.4: compute interface subunit omegaijFluid volume fraction of (VOF)ij；
Figure FDA0003044073150000021
Step 3.3: computing interface Unit omegaiFluid volume fraction of (VOF)i；
Figure FDA0003044073150000022
Wherein, | ΩiL represents the volume of a 3D cell or the area of a 2D cell; n is a radical ofsubIs an interface unit omegaiDivided subunit omegaijThe number of (2);
and 4, step 4: each unit omega in CFD grid of output calculation domainiFluid volume fraction of (VOF)iThe distribution of the two immiscible fluids a and B in the calculated domain is obtained.
2. A method according to claim 1, wherein the fluid distribution calculation method based on the VOF principle includes: outputting each unit omega in the CFD grid of the computation domain in the step 4iFluid volume fraction of (VOF)iWhether the preset precision is met or not needs to be checked in the prior art, and the specific method comprises the following steps:
step 4.1: calculating the 3D volume or 2D area V of the fluid contained in all the interface elementsjmRelative error between two adjacent segmentation-decision levels k and k-1
Figure FDA0003044073150000023
Figure FDA0003044073150000024
Figure FDA0003044073150000025
Figure FDA0003044073150000026
Wherein k is more than or equal to 2 and less than or equal to R;
Figure FDA0003044073150000027
represents the u interface unit;
Figure FDA0003044073150000028
as an interface unit
Figure FDA0003044073150000029
Fluid volume fraction in "segmentation-decision" level k;
Figure FDA00030440731500000210
as an interface unit
Figure FDA00030440731500000211
Subunit omega ofujFluid volume fraction in "segmentation-decision" level k;
Figure FDA0003044073150000031
step 4.2: calculating the fluid 3D volume or 2D area V contained by all the cells in the CFD gridallRelative error between two adjacent segmentation-decision levels k and k-1
Figure FDA0003044073150000032
Figure FDA0003044073150000033
Figure FDA0003044073150000034
Wherein N iscellTo calculate the number of all cells contained in the CFD mesh of the domain;
step 4.3: judgment of
Figure FDA0003044073150000035
And
Figure FDA0003044073150000036
whether the preset precision is met or not; if it is
Figure FDA0003044073150000037
Or
Figure FDA0003044073150000038
If the preset precision is not met, increasing the value of the maximum segmentation-judgment level R, and returning to the step 3; if it is
Figure FDA0003044073150000039
And
Figure FDA00030440731500000310
all meet the preset precision, and then each unit omega in the CFD grid of the calculation domain is outputiFluid volume fraction of (VOF)iThe distribution of the two immiscible fluids a and B in the calculated domain is obtained.
CN202110467994.8A 2021-04-28 2021-04-28 Fluid distribution calculation method based on VOF principle Active CN113177373B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110467994.8A CN113177373B (en) 2021-04-28 2021-04-28 Fluid distribution calculation method based on VOF principle

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110467994.8A CN113177373B (en) 2021-04-28 2021-04-28 Fluid distribution calculation method based on VOF principle

Publications (2)

Publication Number Publication Date
CN113177373A CN113177373A (en) 2021-07-27
CN113177373B true CN113177373B (en) 2022-06-21

Family

ID=76926938

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110467994.8A Active CN113177373B (en) 2021-04-28 2021-04-28 Fluid distribution calculation method based on VOF principle

Country Status (1)

Country Link
CN (1) CN113177373B (en)

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102129517A (en) * 2011-03-10 2011-07-20 西安交通大学 High-precision two-phase fluid interface capturing method
CN107066731A (en) * 2017-04-13 2017-08-18 中南大学 The method that gas distributional pattern in two phase flow is recognized according to numerical simulation result
CN111911148A (en) * 2020-09-16 2020-11-10 西南石油大学 Method for simulating gas-liquid two-phase seepage characteristics in shale natural fracture

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20150242545A1 (en) * 2014-02-21 2015-08-27 Junghyun Cho Method of Simulation of Moving Interfaces using Geometry-Aware Volume of Fluid Method
US11250186B2 (en) * 2019-10-04 2022-02-15 Autodesk, Inc. Machine learning approach to piecewise linear interface construction

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102129517A (en) * 2011-03-10 2011-07-20 西安交通大学 High-precision two-phase fluid interface capturing method
CN107066731A (en) * 2017-04-13 2017-08-18 中南大学 The method that gas distributional pattern in two phase flow is recognized according to numerical simulation result
CN111911148A (en) * 2020-09-16 2020-11-10 西南石油大学 Method for simulating gas-liquid two-phase seepage characteristics in shale natural fracture

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
Verification of Volume-of-Fluid (VOF) simulation for thin liquid film applications;S. Balachandran 等;《2009 3rd International Conference on Energy and Environment (ICEE)》;20091208;449-455 *
基于CFD的船舶破舱进水时域模拟;郑宇 等;《舰船科学技术》;20171018;第39卷(第19期);29-33 *
基于自适应网格的二维孤立波生成;张运兴 等;《第二十九届全国水动力学研讨会论文集（上册）》;20180831;212-218 *
气-液-固体系中夹杂物去除的CFD-VOF-DPM数值模拟;聂海棋;《中国优秀博硕士学位论文全文数据库(硕士)工程科技Ⅰ辑》;20200615(第06期);B023-72 *

Also Published As

Publication number Publication date
CN113177373A (en) 2021-07-27

Similar Documents

Publication Publication Date Title
CN107767457B (en) STL digital-analog generating method based on point cloud rapid reconstruction
Dyadechko et al. Moment-of-fluid interface reconstruction
CN110851929B (en) Two-dimensional leaf-type optimization design method and device based on self-adaptive grid
JP2008165804A (en) Flow simulation calculating method and system
CN112507600B (en) Construction method of symmetrical boundary conditions of semi-implicit method of moving particles
CN113178011B (en) Cut grid THINC method for solving VOF convection equation
CN113177373B (en) Fluid distribution calculation method based on VOF principle
CN117217108A (en) Aerodynamic force analysis method for navigation aircraft flight simulation based on CFD
Dannelongue et al. Three‐dimensional adaptive finite element computations and applications to non‐Newtonian fluids
CN110705189A (en) Method for establishing sedimentation air flotation tank air flotation zone hydrodynamics model
Schmidt et al. Direct interface tracking of droplet deformation
CN114396892B (en) Track curvature measuring method for track traffic curve
van der Pijl Computation of bubbly flows with a mass-conserving level-set method
Loch The level set method for capturing interfaces with applications in two-phase flow problems
CN108073776B (en) Complex river network trunk and branch intersection grid drawing and river center continent grid processing method
CN114003977A (en) Method for optimizing electrode configuration of constructed wetland-microbial fuel cell system
Ong et al. Immersed boundary method with irrotational discrete delta vector for droplet simulations of large density ratio
CN113761778B (en) High-precision geometric model particle dispersion method, device and system based on gridless method
Peng et al. An approach of dynamic mesh adaptation for simulating 3‐dimensional unsteady moving‐immersed‐boundary flows
Sugaya et al. Grid metrics modification approach for flow simulation around 3D geometries on Cartesian CFD method
Schlipf et al. On the application of hybrid meshes in hydraulic machinery CFD simulations
Seshadri Sharp Interface Immersed Boundary Framework for All-Speed Flow Solver
Jarkowski et al. Comparison between sliding and chimera grids
CN117763931A (en) Complex geometric model particle discrete method, system and electronic equipment
Jemison An asymptotically preserving method for multiphase flow

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant