A Simulation Method of Incompressible Viscous Fluid Based on Signed Distance Function Representation

Abstract

The present invention discloses a method for simulating incompressible viscous fluids based on signed distance function representation, and relates to the technical fields of numerical calculation and numerical simulation. The method comprises spatially dividing a calculation area into Cartesian grids; using signed distance functions to represent virtual points in the grid at the solid wall boundary of the calculation area, calculating the perpendicular point from the virtual point to the boundary of the geometric area using a perpendicular calculation algorithm based on signed distance function representation, solving the differential information of the equation at the perpendicular point using an explicit Inverse Lax-Wendroff technique at the perpendicular point, obtaining the grid function value at the virtual point using the differential information, and simultaneously solving the algebraic relationship between the virtual point value and the internal point value, and auxiliary matrix assembly; and using the Crank-Nicolson format to advance the time step forward to realize the numerical simulation of incompressible viscous fluids based on SDF. The present invention realizes complex boundary processing while ensuring the accuracy of the algorithm.

Description

Technical Field

[0001] The present invention belongs to the technical field of numerical calculation and numerical simulation, and in particular to an incompressible viscous fluid simulation method based on signed distance function (SDF) representation. Background Art

[0002] Computer-aided engineering analysis (CAD) is a crucial component of modern industrial digitalization, with extensive and in-depth applications in industries such as aircraft, shipbuilding, and automobiles. Numerical simulation is a core technology within this field. Simulation of incompressible viscous fluids is a key area of ​​CAD, primarily used in low-speed motion scenarios, such as automobiles and submarines. Numerical simulation can help design and verify these areas more efficiently and cost-effectively.

[0003] In the existing technology, most numerical simulations in complex calculation areas rely on the generation of body-fitting grids. This method is inefficient, unstable, and difficult to guarantee calculation quality. Summary of the Invention

[0004] In order to overcome the defects in the above-mentioned prior art, the present invention provides an incompressible viscous fluid simulation method based on signed distance function (SDF) representation, which uses signed distance function to represent the calculation area, and uses explicit and implicit Inverse Lax-Wendroff methods to handle complex boundaries. It can handle complex boundary problems on Cartesian grids, avoids the generation of body-fitting grids, and ensures the computational stability and efficiency of the simulation.

[0005] To achieve the above object, the present invention adopts the following technical solutions, including:

[0006] An incompressible viscous fluid simulation algorithm based on SDF representation includes the following steps:

[0007] Step S1, performing staggered grid Cartesian discretization on the calculation area;

[0008] Step S2: locating a virtual point at the boundary of the calculation area according to the information represented by the signed distance function; the virtual point is a grid point that participates in the calculation but falls outside the calculation area;

[0009] Step S3, using a perpendicular foot calculation algorithm based on a signed distance function to calculate the perpendicular foot point where the virtual point falls on the boundary of the calculation area;

[0010] Step S4, for each perpendicular foot point falling on the boundary of the calculation area, use the internal points to calculate the extrapolated differential information of each order at the perpendicular foot point;

[0011] Step S5: for each foot point falling on the boundary of the calculation area, use the explicit Inverse Lax-Wendroff algorithm combined with the extrapolated differential information of each order to calculate the normal differential information of each physical quantity at the foot point;

[0012] Step S6, calculating the grid function value of the virtual point through the normal differential information of each physical quantity at the perpendicular point and Taylor expansion;

[0013] Step S7, combining the normal differential information of each physical quantity at the foot point and the calculation process of the grid function value of the virtual point, using the implicit Inverse Lax-Wendroff algorithm to calculate the algebraic relationship between the virtual point value and the internal point value;

[0014] Step S8, using the Crank-Nicolson format and the second-order Adams format to advance the time step, using the grid function values ​​of the virtual points and the internal point values, through the central difference format to discretize the convection term, diffusion term, pressure term and incompressible term, based on the auxiliary matrix assembly of the algebraic relationship calculated based on the implicit Inverse Lax-Wendroff algorithm, to realize the numerical simulation of incompressible viscous fluid based on the signed distance function.

[0015] Preferably, the specific process of step S1 is as follows:

[0016] Perform uniform Cartesian discretization of the calculation area horizontally and vertically, and give the calculation rectangular area ,in is the left boundary of the calculation area, is the right boundary of the calculation area, is the lower boundary of the calculation area, To calculate the upper boundary of the region, the rectangular region is spatially discretized to obtain:

[0017] ,

[0018] ,

[0019] And there are spatially discrete grid coordinates in is the subscript of the grid division, and is the number of grid divisions. Here, three groups of interlaced grids are generated. The first group consists of The second group consists of grid points of the form The third group is composed of grid points of The grid is composed of as well as , the grid division form is called staggered grid.

[0020] For the governing equations:

[0021]

[0022] in, are spatial coordinates, and is a physical variable, for Direction speed, for Direction speed, It's time, is the convection term;

[0023] ,

[0024] is the pressure, is the source term, which includes two components, namely , is the gradient operator, , is the Reynolds number, is the diffusion term;

[0025] ,

[0026] Calculation area , defining its boundaries as , and deal with the following boundary conditions: ;

[0027] Define three grid point sets as , and ,satisfy:

[0028] ,

[0029] Get the corresponding discrete physical variables . Only Discrete on, Only Discrete on, Only Sometimes we also use bold Indicates a subscript, for example ,at this time ,akin .

[0030] Preferably, in step S2,

[0031] Given a calculation area, the signed distance function is expressed as , the set of positioning virtual points is:

[0032] ,

[0033] in, ; Among them, the signed distance function The specific definition of is: for a given calculation area , the boundary of the calculation area is , for any point in space ,have:

[0034] ,

[0035] in, boundary distance, for The complement of .

[0036] Preferably, in step S3, the foot of the virtual point falling on the boundary of the calculation area is calculated by using the foot of the virtual point calculation algorithm. The specific calculation method is: for any virtual point , calculate its The foot point on , ,in, is the gradient operator.

[0037] Preferably, in step S4, for each foot point that falls on the boundary of the calculation area, the physical quantity at the foot point is calculated using the internal points The specific calculation method of each order extrapolated differential information is as follows:

[0038] S41, at the foot of the vertical Perform local coordinate rotation so that the new coordinate system Axis and foot point The normal direction is the same as that of , and the velocity field is also rotated; at this time:

[0039]

[0040] in, is the coordinate after rotation, is the coordinate of the original coordinate axis, After the rotation Direction speed, Direction speed, and is the source term after rotation, is the rotation angle, unchanged after rotation;

[0041] S42, Finding the Foot of the Perpendicular The least squares system is constructed at the nearby grid points. The specific calculation is as follows:

[0042]

[0043] in, is an indicator, for example ,have . , and Represents The set of nearby internal grid points, Represents the grid and the foot point At this time, for An estimate of for The partial derivative order of the direction, for Directional partial derivative order. , specifically:

[0044]

[0045] S43, solve according to step S42 , and ,remember , and , rewrite the three equations of step S42 as In the form of , remember the matrices corresponding to these three equations for ,have:

[0046]

[0047] in, 、 、 is a matrix No. item, and is the subscript, 、 and is a collection 、 and The number of elements in ;

[0048] S44, combined with the coordinate rotation in step S41, obtains the following equation:

[0049] ,

[0050] in, yes , or Middle elements, for Approximation, and Similarly, there are:

[0051] ,

[0052] in, , is the matrix that affects the derivatives of the rotation, , and That is, solving the intermediate parameter matrix generated in step S42 under the new coordinate system.

[0053] Preferably, in step S5, for each foot point falling on the boundary of the calculation area, the explicit Inverse Lax-Wendroff algorithm is used to calculate the normal differential information of each physical quantity at the foot point. This algorithm needs to calculate the 0th, 1st, and 2nd order normal differential information, as follows:

[0054] S51, calculate the 0th order differential information according to the boundary conditions, which is:

[0055] ,

[0056] Superscript here It is the abbreviation of Inverse Lax-Wendroff. The 0 in the superscript represents the 0th order.

[0057] S52, based on incompressibility, the first-order normal differential information is calculated as:

[0058] ,

[0059] S53, the second-order normal differential information is:

[0060]

[0061] Preferably, in step S6, the grid function value of the virtual point is calculated by using the normal differential information of each physical quantity at the perpendicular point and the Taylor expansion algorithm, including: , the value of its physical quantity is:

[0062] ,

[0063] Where, For virtual points To the foot point The distance; considering the coordinate rotation, there is:

[0064] .

[0065] Preferably, in step S7, an implicit Inverse Lax-Wendroff algorithm is used to construct an algebraic relationship between the virtual point value and the internal point value. The specific calculation formula is:

[0066]

[0067] and

[0068] ,

[0069] In the preferred step S8, the Crank-Nicolson format and the second-order Adams format are used to advance the time step. The grid function values ​​and internal function values ​​of the virtual points are used to discretize the convection term, diffusion term, pressure term and incompressible term through the central difference format. The auxiliary matrix of the algebraic relationship calculated based on the implicit Inverse Lax-Wendroff algorithm is assembled. The specific calculation is as follows:

[0070] S81, for the governing equations, the time difference format is:

[0071] ,

[0072] in, is the time interval, , for Moment , Is a natural number, indicating the time layers;

[0073] ,

[0074] Where, The right-hand term.

[0075] S82, will Grid function replace, Grid function replace, Grid function Instead, rewrite the equation in step S81 as:

[0076] ,

[0077] in, , for , The first-order difference quotient of the direction, , for , The second-order difference quotient of direction;

[0078] ,

[0079] in, , then for any variable , grid subscript , ,but:

[0080]

[0081] remember and for The first and second components of :

[0082] ,

[0083] in, for The first component is at the grid point The algebraic discretization of is:

[0084] ,

[0085] Similar to definable and .

[0086] S83, by implementing the above discretization at each grid point, a series of linear equations can be obtained. The specific equations are:

[0087]

[0088] in, is composed of all grid point functions and The vector composed of It is composed of all forms and Vector composed of. is the identity matrix, and for and The corresponding matrix, and for The corresponding matrix, for The corresponding matrix, for The corresponding matrix. Note that the calculations in the matrix and the right-hand side of the equation system may involve the values ​​of virtual points, which can be entered using the results of S6 and S7.

[0089] The beneficial effects of the present invention are:

[0090] This paper uses the Inverse Lax-Wendroff method to calculate the differential information of various orders of physical quantities at the foot point, ensuring the physical consistency of the foot point differential information with the internal computational domain and the governing equations. This allows for high-precision boundary processing on a Cartesian grid, thereby achieving high precision for the entire simulation algorithm and enabling the processing of complex boundary problems on a Cartesian grid. Furthermore, the paper explores the algebraic relationship between internal grid point values ​​and virtual point values, making this process applicable to implicit formats.

[0091] The spatial calculation format and boundary processing of the present invention both have theoretical third-order accuracy, wherein the convection term, diffusion term and pressure term are discretized using fourth-order central difference, and the boundary processing uses quadratic polynomial interpolation, which ensures third-order accuracy and can ensure the accuracy of the final simulation.

[0092] The present invention uses SDF to represent the calculation area, avoids the generation of body-fitting grids, can quickly locate virtual points and quickly calculate the perpendicular feet of virtual points, and ensures the calculation stability and efficiency of the simulation. BRIEF DESCRIPTION OF THE DRAWINGS

[0093] Figure 1 This is a flow chart of an incompressible viscous fluid simulation method based on signed distance function representation according to the present invention;

[0094] Figure 2 A two-dimensional cylindrical display diagram in an embodiment;

[0095] Figure 3 is a grid division diagram in the embodiment;

[0096] Figure 4 is a virtual point positioning diagram in the embodiment;

[0097] Figure 5 Schematic diagram of simulation results in the embodiment. DETAILED DESCRIPTION

[0098] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts are within the scope of protection of the present invention.

[0099] like Figure 1As shown in FIG, the flow chart of the incompressible viscous fluid simulation method based on signed distance function (SDF) representation of the present invention specifically includes the following steps:

[0100] Step S1, such as Figure 2 、 Figure 3 As shown in the figure, taking the flow around a cylinder as an example, the specific steps are as follows:

[0101] Perform uniform Cartesian discretization of the calculation area horizontally and vertically, and give the calculation rectangular area ,in is the left boundary of the calculation area, is the right boundary of the calculation area, is the lower boundary of the calculation area, To calculate the upper boundary of the region, the rectangular region is spatially discretized to obtain:

[0102] ,

[0103] ,

[0104] And there are spatially discrete grid coordinates in is the subscript of the grid division, and is the number of grid divisions.

[0105] For the governing equations:

[0106]

[0107] in, are spatial coordinates, and is a physical variable, for Direction speed, for Direction speed, It's time, is the convection term;

[0108] ,

[0109] is the pressure, is the source term, which includes two components, namely , is the gradient operator, , is the Reynolds number, is the diffusion term;

[0110] ,

[0111] Calculation area , defining its boundaries as , and deal with the following boundary conditions:

[0112] ,

[0113] Define three grid point sets as , and ,satisfy:

[0114] ,

[0115] Get the corresponding discrete physical variables .

[0116] S2, such as Figure 4 As shown, near the boundary of the calculation area, that is, the geometric boundary of the cylinder, locate the virtual point (the virtual point is the grid point that participates in the calculation but falls outside the calculation area). Figure 4 In the example, for a cylinder, it is divided into three layers from the inside to the outside, namely internal points, virtual points, and external points. The signed distance function of a given cylinder is expressed as , the virtual point positioning algorithm is: the positioning virtual point set is:

[0117] ,

[0118] in, ; Among them, the signed distance function The specific definition of is: for a given calculation area , the boundary of the calculation area is , for any point in space ,have:

[0119] ,

[0120] in, boundary distance, for The complement of .

[0121] S3, calculate the foot point of the virtual point on the cylinder boundary through the foot calculation algorithm. The specific calculation method is: for any virtual point , calculate its The foot point on , ,in, is the gradient operator.

[0122] S4, for each foot point that falls on the cylinder boundary, use the internal points to calculate the physical quantity at the foot point The extrapolated differential information of each order. The specific calculation is as follows:

[0123] S41, at the foot of the vertical Perform local coordinate rotation so that the new coordinate system Axis and foot point The normal direction is the same as that of , and the velocity field is also rotated; at this time:

[0124]

[0125] in, is the coordinate after rotation, is the coordinate of the original coordinate axis, After the rotation Direction speed, Direction speed, and is the source term after rotation, is the rotation angle, unchanged after rotation;

[0126] S42, Finding the Foot of the Perpendicular The least squares system is constructed at the nearby grid points. The specific calculation is as follows:

[0127]

[0128] in, is an indicator, for example ,have . , and Represents The set of nearby internal grid points, Represents the grid and the foot point At this time, for An estimate of for The partial derivative order of the direction, for Directional partial derivative order. , specifically:

[0129]

[0130] S43, solve according to step S42 , and ,remember , and , rewrite the three equations of step S42 as In the form of , remember the matrices corresponding to these three equations for ,have:

[0131]

[0132] in, 、 、 is a matrix No. item, and is the subscript, 、 and is a collection 、 and The number of elements in ;

[0133] S44, combined with the coordinate rotation in step S41, obtains the following equation:

[0134] ,

[0135] in, yes , or Middle elements, for Approximation, and Similarly, there are:

[0136] ,

[0137] in, , is the matrix that affects the derivatives due to the rotation, , and That is, solving the intermediate parameter matrix generated in step S42 under the new coordinate system.

[0138] S5, calculate the grid function value of the virtual point using the above differential information and Taylor expansion algorithm. The specific calculation method is:

[0139] S51, calculate the 0th order differential information according to the boundary conditions, which is:

[0140] ,

[0141] Superscript here It is the abbreviation of Inverse Lax-Wendroff. The 0 in the superscript represents the 0th order.

[0142] S52, based on incompressibility, the first-order normal differential information is calculated as:

[0143] ,

[0144] S53, the second-order normal differential information is:

[0145]

[0146] S6, calculate the grid function value of the virtual point through the normal differential information of each physical quantity at the perpendicular point and the Taylor expansion algorithm. , the value of its physical quantity is:

[0147] ,

[0148] Taking into account the coordinate rotation, we have:

[0149] .

[0150] S7, use the implicit Inverse Lax-Wendroff algorithm to construct the algebraic relationship between the virtual point value and the internal point value. The specific calculation formula is:

[0151]

[0152] and

[0153] ,

[0154] S8, using the Crank-Nicolson scheme and the second-order Adams scheme to advance the time step, using the grid function values ​​and internal function values ​​of virtual points, discretizing the convection term, diffusion term, pressure term, and incompressible term through the central difference scheme, and assembling the auxiliary matrix based on the algebraic relationship calculated by the implicit Inverse Lax-Wendroff algorithm. The specific calculation is as follows:

[0155] S81, for the governing equations, the time difference format is:

[0156] ,

[0157] in, is the time interval, , for Moment , Is a natural number, indicating the time layers;

[0158] ,

[0159] Where, The right-hand term.

[0160] S82, will Grid function replace, Grid function replace, Grid function Instead, rewrite the equation in step S81 as:

[0161] ,

[0162] in, , for , The first-order difference quotient of the direction, , for , The second-order difference quotient of direction;

[0163] ,

[0164] in, , then for any variable , grid subscript , ,but:

[0165]

[0166] remember and for The first and second components of :

[0167] ,

[0168] in, for The first component is at the grid point The algebraic discretization of is:

[0169] ,

[0170] Similar to definable and .

[0171] S83, by implementing the above discretization at each grid point, a series of linear equations can be obtained. The specific equations are:

[0172] ,

[0173] in, is composed of all grid point functions and The vector composed of It is composed of all forms and Vector composed of. is the identity matrix, and for and The corresponding matrix, and for The corresponding matrix, for The corresponding matrix, for The corresponding matrix. Note that the calculations in the matrix and the right-hand side of the equation group may involve the value of the virtual point, and the results of S6 and S7 can be used to enter. The simulation results are as follows Figure 5 shown.

[0174] The specific embodiments described above further illustrate the objectives, technical solutions and beneficial effects of the present invention in detail. It should be understood that the above are only specific embodiments of the present invention and are not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention should be included in the scope of protection of the present invention.

Claims

1. A method for simulating incompressible viscous fluid based on signed distance function representation, characterized in that: The following steps are involved: Step S1, performing staggered grid Cartesian discretization on the calculation area; Step S2: locating a virtual point at the boundary of the calculation area according to the information represented by the signed distance function; the virtual point is a grid point that participates in the calculation but falls outside the calculation area; Step S3, using a perpendicular foot calculation algorithm based on a signed distance function to calculate the perpendicular foot point where the virtual point falls on the boundary of the calculation area; Step S4, for each perpendicular foot point falling on the boundary of the calculation area, use the internal points to calculate the extrapolated differential information of each order at the perpendicular foot point; Step S5: for each foot point falling on the boundary of the calculation area, use the explicit Inverse Lax-Wendroff algorithm combined with the extrapolated differential information of each order to calculate the normal differential information of each physical quantity at the foot point; Step S6, calculating the grid function value of the virtual point by using the normal differential information of each physical quantity at the perpendicular point and the Taylor expansion algorithm; Step S7, combining the normal differential information of each physical quantity at the foot point and the calculation process of the grid function value of the virtual point, using the implicit Inverse Lax-Wendroff algorithm to calculate the algebraic relationship between the virtual point value and the internal point value; Step S8, using the Crank-Nicolson format and the second-order Adams format to advance the time step, using the grid function values ​​of the virtual points and the internal point values, through the central difference format to discretize the convection term, diffusion term, pressure term and incompressible term, based on the auxiliary matrix assembly of the algebraic relationship calculated based on the implicit Inverse Lax-Wendroff algorithm, to realize the numerical simulation of incompressible viscous fluid based on the signed distance function.

2. The incompressible viscous fluid simulation method based on signed distance function representation according to claim 1, characterized in that: The specific process of step S1 is as follows: Perform uniform Cartesian discretization of the calculation area horizontally and vertically, and give the calculation rectangular area ,in is the left boundary of the calculation area, is the right boundary of the calculation area, is the lower boundary of the calculation area, To calculate the upper boundary of the region, the rectangular region is spatially discretized to obtain: ， ， And there are spatially discrete grid coordinates in is the subscript of the grid division; For the governing equations: , in, are spatial coordinates, and is a physical variable, for Direction speed, for Direction speed, It's time, is the convection term; , is the pressure, is the source term, which includes two components, namely , is the gradient operator, , is the Reynolds number, is the diffusion term; , Calculation area , defining its boundaries as , and deal with the following boundary conditions: ， Define three grid point sets as , and ,satisfy: ， Get the corresponding discrete physical variables .

3. The incompressible viscous fluid simulation method based on signed distance function representation according to claim 2, characterized in that: In step S2, the given calculation area signed distance function is expressed as , the set of positioning virtual points is: ， in, is the largest grid size near the edge; where the signed distance function The specific definition of is: for a given calculation area , the boundary of the calculation area is , for any point in space ,have: , in, for point To the border distance, for The complement of .

4. The incompressible viscous fluid simulation method based on signed distance function representation according to claim 3 is characterized in that: Step S3 includes, for any virtual point , calculate its boundary in the calculation area The foot point on , ,in, is the gradient operator.

5. The incompressible viscous fluid simulation method based on signed distance function representation according to claim 4, characterized in that: In step S4, for each foot point that falls on the boundary of the calculation area, the physical quantity at the foot point is calculated using the internal points. The extrapolated differential information of each order includes: S41, at the foot of the vertical Perform local coordinate rotation so that the new coordinate system Axis and foot point The normal direction is the same as that of , and the velocity field is also rotated; at this time: , in, is the coordinate after rotation, is the coordinate of the original coordinate axis, After the rotation Direction speed, Direction speed, and is the source term after rotation, is the rotation angle, unchanged after rotation; S42, Finding the Foot of the Perpendicular The least squares system is constructed at the nearby grid points. The specific calculation is as follows: in, is the grid index, , and Represents the foot point The set of nearby internal grid points, Represents the grid and the foot point The deviation, for An estimate of for The partial derivative order of the direction, for The partial derivative order of the direction, Yes and 、 The relevant parameters reflect the and order; S43, solve according to step S42 , and ,remember , and , rewrite the three equations of step S42 as In the form of , remember the matrices corresponding to these three equations for ,have: ， in, 、 、 is a matrix No. item, and is the subscript, 、 and is a collection 、 and The number of elements in ; S44, combined with the coordinate rotation in step S41, obtains the following equation: ， in, yes , or Middle elements, for Approximation, and is the corresponding approximation, and we have: ， in, , is the matrix that affects the derivatives due to the rotation, , and That is, solving the intermediate parameter matrix generated in step S42 under the new coordinate system.

6. The incompressible viscous fluid simulation method based on signed distance function representation according to claim 5, characterized in that: The step S5 comprises: S51, calculate the 0th order differential information according to the boundary conditions, which is: ， Superscript here It is the abbreviation of Inverse Lax-Wendroff, and the 0 in the superscript represents the 0th order; S52, based on incompressibility, the first-order normal differential information is calculated as: ， S53, the second-order normal differential information is: 。 7. The incompressible viscous fluid simulation method based on signed distance function representation according to claim 6, characterized in that: In step S6, the grid function value of the virtual point is calculated by using the normal differential information of each physical quantity at the perpendicular point and the Taylor expansion algorithm, including: , the value of its physical quantity is: ， Where, For virtual points To the foot point distance; Taking into account the coordinate rotation, we have: 。 8. The incompressible viscous fluid simulation method based on signed distance function representation according to claim 7, characterized in that: Step S7 combines the calculation process of step S5 and step S6, and uses the implicit Inverse Lax-Wendroff algorithm to construct the algebraic relationship between the virtual point value and the internal point value, including: ， and 。 9. The incompressible viscous fluid simulation method based on signed distance function representation according to claim 8, characterized in that: The step S8 comprises: S81, for the governing equations, the time difference format is: ， in, is the time interval, , for Moment , Is a natural number, indicating the time layers; ， Where, is the right-hand term; S82, will Grid function replace, Grid function replace, Grid function Instead, rewrite the equation in step S81 as: ， in, , for , The first-order difference quotient of the direction, , for , The second-order difference quotient of direction; ， in, , then for any variable , grid subscript , ,but: , remember and for The first and second components of : ， in, for The first component is at the grid point The algebraic discretization of is: ； Similar can be defined , S83, performing the above discretization at each grid point to obtain a series of linear equations, thereby realizing the numerical simulation of incompressible viscous fluid based on SDF.

Images