Incompressible viscous fluid simulation method based on sign distance function representation

Abstract

The invention discloses an incompressible viscous fluid simulation method based on symbolic distance function representation, and relates to the technical field of numerical calculation and numerical simulation. The method comprises the following steps: carrying out Cartesian grid division on a calculation area in space; at the solid wall boundary of the calculation area, using a sign distance function to represent a virtual point in the positioning grid, calculating a perpendicular foot point from the virtual point to the geometric area boundary through a perpendicular foot calculation algorithm represented based on the sign distance function, and solving differential information of an equation at the perpendicular foot point through an explicit Inverse Lax-Wendroff technology on the perpendicular foot point, obtaining a lattice point function value at a virtual point through differential information, solving an algebraic relational expression between a virtual point value and an internal point value, and assisting matrix assembly; a Crank-Nicolson format is used for advancing a time step, SDF-based numerical simulation of the incompressible viscous fluid is realized, complex boundary processing is realized, and meanwhile, the precision of an algorithm is ensured.

Description

Technical Field

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

[0002] Computer-aided engineering analysis is an important part of the field of modern industrial digitalization. It has extremely wide and in-depth applications in industries such as aircraft, ships, and automobiles, and numerical simulation is a very core technology. Among them, incompressible viscous fluid simulation is a direction in computer-aided engineering. This simulation is mainly used in low-speed motion situations, such as automobiles, submarines, etc. Numerical simulation can assist relevant staff in designing and verifying related fields more efficiently and economically.

[0003] In the prior art, most numerical simulations in complex calculation areas rely on the generation of body-fitting grids. This method is inefficient, has low stability, and is difficult to ensure 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), which uses signed distance function to represent the calculation area, uses explicit and implicit InverseLax-Wendroff methods to process complex boundaries, and can process complex boundary problems on Cartesian grids, avoiding the generation of body-fitting grids, and ensuring the calculation 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 specifically 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, using the internal points to calculate the extrapolated differential information of each order at the perpendicular foot point;

[0011] Step S5, for each perpendicular foot point falling on the boundary of the calculation area, using the explicit Inverse Lax-Wendroff algorithm combined with each order of extrapolated differential information to calculate each order of normal differential information of each physical quantity at the perpendicular 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 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 algebraic relationship auxiliary matrix assembly calculated by the implicit Inverse Lax-Wendroff algorithm, to achieve 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] The calculation area is uniformly discretized horizontally and vertically in Cartesian. Given a calculation rectangular area ,in is the left boundary of the calculation area, is the right boundary of the computational region, 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 points are composed of as well as , the grid division form is called staggered grid.

[0020] For the governing equation:

[0021]

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

[0023] ,

[0024] For 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 in Discrete on, Only in Discrete on, Only in 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 The distance for The complement of .

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

[0037] Preferably, in step S4, for each perpendicular point falling on the boundary of the calculation area, the physical quantity at the perpendicular point is calculated using the internal points The specific calculation method of each order of 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, are the coordinates of the original axes, After the rotation The speed of the direction, The speed of the direction, and is the source term after rotation, is the rotation angle, It remains unchanged after rotation;

[0041] S42, Find the foot point The least squares system is constructed at 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 foot point At this time, for An estimate of for The partial derivative order of the direction, for The partial derivative order of the direction. , specifically:

[0044]

[0045] S43, solving 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 of ;

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

[0049] ,

[0050] in, yes , or Middle elements, for The approximation of and Similarly, there are:

[0051] ,

[0052] 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.

[0053] Preferably, in step S5, for each perpendicular foot point falling on the boundary of the calculation area, an explicit Inverse Lax-Wendroff algorithm is used to calculate the normal differential information of each physical quantity at the perpendicular foot point. This algorithm needs to calculate the 0th, 1st, and 2nd order normal differential information, which is 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, the first-order normal differential information is calculated based on incompressibility:

[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 foot point and the Taylor expansion algorithm, including: , the value of its physical quantity is:

[0062] ,

[0063] In the formula, For virtual point To the foot point The distance; considering the coordinate rotation, we have:

[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, and 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, and the auxiliary matrix assembly of the algebraic relationship calculated based on the implicit Inverse Lax-Wendroff algorithm is specifically calculated 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] In the formula, The right-hand term.

[0075] S82, will Using grid function replace, Using grid function replace, Using 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 group may involve the values ​​of virtual points, and the results of S6 and S7 can be used to enter them.

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

[0090] The present invention uses the Inverse Lax-Wendroff method to calculate the differential information of each order of the physical quantity at the foot point, ensures the physical consistency of the foot point differential information and the internal calculation area and the control equation, realizes high-precision boundary processing on the Cartesian grid, thereby realizing the high precision of the entire simulation algorithm, and can handle complex boundary problems under the Cartesian grid. At the same time, the present invention explores the algebraic relationship between the internal grid point value and the virtual point value, so that this process can be applied in the implicit format.

[0091] The spatial calculation format and boundary processing of the present invention 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 foot of the virtual points, and ensures the calculation stability and efficiency of the simulation. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0094] Figure 2 It is a two-dimensional cylinder display diagram in the embodiment;

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

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

[0097] Figure 5 It is a schematic diagram of the simulation results in the embodiment. DETAILED DESCRIPTION

[0098] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. 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 creative work are within the scope of protection of the present invention.

[0099] like Figure 1As shown, the flow chart of the incompressible viscous fluid simulation method based on the 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] The calculation area is uniformly discretized horizontally and vertically in Cartesian. Given a calculation rectangular area ,in is the left boundary of the calculation area, is the right boundary of the computational region, 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 equation:

[0106]

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

[0108] ,

[0109] For 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 in the figure, virtual points are located near the boundary of the calculation area, that is, the geometric boundary of the cylinder (virtual points are grid points that participate in the calculation but fall 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 set of positioning virtual points 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 The 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 of , ,in, is the gradient operator.

[0122] S4, for each foot point falling 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, are the coordinates of the original axes, After the rotation The speed of the direction, The speed of the direction, and is the source term after rotation, is the rotation angle, It remains unchanged after rotation;

[0126] S42, Find the foot point The least squares system is constructed at 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 foot point At this time, for An estimate of for The partial derivative order of the direction, for The partial derivative order of the direction. , specifically:

[0129]

[0130] S43, solving 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 of ;

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

[0134] ,

[0135] in, yes , or Middle elements, for The approximation of 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, using the above differential information and Taylor expansion algorithm to calculate the grid function value of the virtual point, 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, the first-order normal differential information is calculated based on incompressibility:

[0143] ,

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

[0145]

[0146] S6, the grid function value of the virtual point is calculated by using the normal differential information of each physical quantity at the foot 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 format and the second-order Adams format 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 format, and assembling the auxiliary matrix of the algebraic relationship calculated based on 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] In the formula, The right-hand term.

[0160] S82, will Using grid function replace, Using grid function replace, Using 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 values ​​of virtual points, which can be entered using the results of S6 and S7. The simulation results are shown in 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 description is only a specific embodiment of the present invention and is 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, using the internal points to calculate the extrapolated differential information of each order at the perpendicular foot point; Step S5, for each perpendicular foot point falling on the boundary of the calculation area, using the explicit Inverse Lax-Wendroff algorithm combined with each order of extrapolated differential information to calculate each order of normal differential information of each physical quantity at the perpendicular foot point; 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 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 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 algebraic relationship auxiliary matrix assembly calculated by the implicit Inverse Lax-Wendroff algorithm, to achieve 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: The calculation area is uniformly discretized horizontally and vertically in Cartesian. Given a calculation rectangular area ,in is the left boundary of the calculation area, is the right boundary of the computational region, 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 equation: , in, are spatial coordinates, and is a physical variable, for The speed of the direction, for The speed of the direction, It's time. is the convection term; , For 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 is 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 maximum 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 The 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 of , ,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 perpendicular point falling on the boundary of the calculation area, the physical quantity at the perpendicular 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, are the coordinates of the original axes, After the rotation The speed of the direction, The speed of the direction, and is the source term after rotation, is the rotation angle, It remains unchanged after rotation; S42, Find the foot point The least squares system is constructed at 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 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 and order; S43, solving 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 of ; S44, combined with the coordinate rotation in step S41, obtains the following equation: ， in, yes , or Middle elements, for Approximation of and is the corresponding approximation, and at the same time: ， 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. The 0 in the superscript represents the 0th order; S52, the first-order normal differential information is calculated based on incompressibility: ， 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 is 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: ， In the formula, For virtual point 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 an 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; ， In the formula, is the right-hand term; S82, will Using grid function replace, Using grid function replace, Using 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 to the definition , S83, implementing the above discretization at each grid point to obtain a series of linear equations, thereby realizing the SDF-based numerical simulation of incompressible viscous fluid.

Images