A numerical simulation method and device based on high-precision immersed boundary method

Abstract

The embodiment of the present application discloses a numerical simulation method and device based on a high-precision immersed boundary method, and relates to the field of computer technology. The method obtains the left eigenvector and the right eigenvector of the flux at the target half-node; converts the positive and negative flux values ​​at the template node into feature space variables according to the left eigenvector; nonlinearly weights the feature space variables at the template node by the WENO-type nonlinear interpolation method, obtains the left state value and the right state value of the feature space variable at the half-node, and the nonlinear weighting coefficient is determined by the Z-type nonlinear weight; according to the left state value, the right state value and the right eigenvector, the numerical approximation of the flux at the target half-node is obtained. Therefore, by introducing the WENO-type nonlinear interpolation method, the feature space variables at the template node are nonlinearly weighted, and while maintaining numerical stability, numerical oscillation can be significantly reduced, thereby improving the resolution and accuracy of the numerical approximation of the flux at the target half-node.

Description

Technical Field

[0001] The present application relates to the field of computer technology, and in particular to a numerical simulation method and device based on a high-precision immersed boundary method. Background Art

[0002] In the field of computational fluid dynamics (CFD), high-order precision formats have smaller dissipation, dispersion errors and higher spatial resolution than low-order precision formats, so they can more accurately depict complex flow field structures. Among them, the third-order weighted compact nonlinear scheme (WCNS) format, as a high-order precision format, has been widely used in structured grids. However, in Cartesian grids, when the WCNS format is used to deal with flow problems with complex geometric shapes, non-physical oscillations are prone to occur, which affects the stability and accuracy of the calculation.

[0003] In order to cope with the computational challenges of complex geometric flow fields under Cartesian grids, the relevant technologies use the third-order WCNS-JS format and the third-order WCNS-Z format based on the third-order precision virtual unit immersed boundary method. These two formats enhance the stability of the format by introducing nonlinear weights and achieve high-order accuracy in smooth areas.

[0004] However, the third-order WCNS-JS format and the third-order WCNS-Z format have order reduction phenomenon at the first-order and second-order extreme points of the solution, resulting in reduced resolution and insufficient accuracy. Summary of the invention

[0005] Based on the above problems, the present application provides a numerical simulation method and device based on a high-precision immersed boundary method, which can improve the resolution and accuracy of the numerical approximation of the flux at the target half-node.

[0006] The embodiments of the present application disclose the following technical solutions:

[0007] In a first aspect, the present application discloses a numerical simulation method based on a high-precision immersed boundary method, the method comprising:

[0008] Get the left eigenvector and right eigenvector of the flux at the target half-node;

[0009] According to the left eigenvector, the positive flux value and the negative flux value at the template node are converted into feature space variables;

[0010] The feature space variables at the template node are nonlinearly weighted by a high-precision weighted essential non-oscillating WENO-type nonlinear interpolation method to obtain the left state value and the right state value of the feature space variable at the half-node, wherein the nonlinear weighting coefficient is determined by the Z-type nonlinear weight; according to the left state value, the right state value and the right eigenvector, the numerical approximation of the flux at the target half-node is obtained by converting the feature space into the physical space.

[0011] Optionally, the weighting coefficient of the nonlinear weighting is determined by a Z-type nonlinear weight and a global smoothing factor, wherein the formula of the global smoothing factor is specifically as follows:

[0012]

[0013] in, is the global smoothing factor, β n is the sub-smoothing factor of the nth sub-template.

[0014] Optionally, obtaining the left eigenvector and the right eigenvector of the flux at the target half-node includes:

[0015] Get the Jacobian matrix of the flux, which is the matrix of all partial derivatives of the flux with respect to its independent variable, and the independent variable is U=(ρ,ρu,ρv,E) T , where ρ is the fluid density, u and v are the velocities of the fluid in the x and y directions respectively, and E is the total energy per unit volume;

[0016] Performing eigendecomposition on the Jacobian matrix to obtain a left eigenvector and a right eigenvector of the flux at the target half-node, wherein the Jacobian matrix is ​​the product of a diagonal matrix, a left eigenvector and a right eigenvector of the flux at the target half-node.

[0017] Optionally, the method further comprises:

[0018] For nodes near the object surface outside the flow field, a two-variable quadratic polynomial of the objective function at the virtual unit is constructed according to the flow field information in the interpolation template, wherein the interpolation template is determined based on the projection point of the virtual unit on the object surface, the object surface normal, the inward unit normal vector and the corresponding tangent vector;

[0019] The constrained least square method is used to fit the bivariate quadratic polynomial of the objective function to obtain a third-order approximation at the virtual unit.

[0020] Optionally, the method further comprises:

[0021] When there is a large gradient or discontinuity in the flow field near the object surface, the WENO type extrapolation method is used, and equidistant interpolation points are taken in the normal direction in combination with the geometric characteristics of the object surface;

[0022] Nonlinear weighting is performed based on the equidistant interpolation points to obtain a third-order approximation at the virtual unit, wherein the weight coefficient of the nonlinear weighting is determined according to the smoothing factor at the equidistant interpolation points.

[0023] In a second aspect, the present application discloses a numerical simulation device based on a high-precision immersed boundary method, the device comprising: a vector acquisition module, a variable acquisition module, a state acquisition module and a numerical approximation module;

[0024] The vector acquisition module is used to obtain the left eigenvector and the right eigenvector of the flux at the target half node;

[0025] The variable acquisition module is used to convert the positive flux value and the negative flux value at the template node into a feature space variable according to the left eigenvector;

[0026] The state acquisition module is used to perform nonlinear weighting on the feature space variables at the template nodes by a high-precision weighted essential non-oscillating WENO type nonlinear interpolation method to obtain the left state value and the right state value of the feature space variables at the half node, wherein the nonlinear weighting coefficient is determined by a Z-type nonlinear weight;

[0027] The numerical approximation module is used to obtain the numerical approximation of the flux at the target half-node by converting the feature space into the physical space according to the left state value, the right state value and the right eigenvector.

[0028] Optionally, the weighting coefficient of the nonlinear weighting is determined by a Z-type nonlinear weight and a global smoothing factor, wherein the formula of the global smoothing factor is specifically as follows:

[0029]

[0030] in, is the global smoothing factor, and βn is the sub-smoothing factor of the nth sub-template.

[0031] Optionally, the vector acquisition module specifically includes: a first acquisition submodule and a second acquisition submodule;

[0032] The first acquisition submodule is used to obtain the Jacobian matrix of the flux, where the Jacobian matrix is ​​a matrix of all partial derivatives of the flux with respect to its independent variable, and the independent variable is U=(ρ,ρu,ρv,E) T , where ρ is the fluid density, u and v are the velocities of the fluid in the x and y directions respectively, and E is the total energy per unit volume;

[0033] The second acquisition submodule is used to perform eigendecomposition on the Jacobian matrix to obtain the left eigenvector and the right eigenvector of the flux at the target half-node, wherein the Jacobian matrix is ​​the product of a diagonal matrix and the left eigenvector and the right eigenvector of the flux at the target half-node.

[0034] Optionally, the device further comprises: a polynomial construction module and a polynomial fitting module;

[0035] The polynomial construction module is used to construct a two-variable quadratic polynomial of the objective function at the virtual unit for the nodes near the object surface outside the flow field according to the flow field information in the interpolation template, wherein the interpolation template is determined based on the projection point of the virtual unit on the object surface, the object surface normal, the inward unit normal vector and the corresponding tangent vector;

[0036] The polynomial fitting module is used to fit the binary quadratic polynomial of the objective function by using the constrained least square method to obtain a third-order approximation at the virtual unit.

[0037] Optionally, the device further comprises: an interpolation point acquisition module and an approximate value acquisition module;

[0038] The interpolation point acquisition module is used to adopt the WENO type extrapolation method and take equidistant interpolation points in the normal direction in combination with the geometric characteristics of the object surface when there is a large gradient or discontinuity in the flow field near the object surface;

[0039] The approximate value acquisition module is used to perform nonlinear weighting based on the equidistant interpolation points to obtain a third-order approximate value at the virtual unit, wherein the weight coefficient of the nonlinear weighting is determined according to the smoothness factor at the equidistant interpolation points.

[0040] Compared with the prior art, this application has the following beneficial effects:

[0041] The present application discloses a numerical simulation method and device based on a high-precision immersed boundary method, the method comprising: obtaining a left eigenvector and a right eigenvector of a flux at a target half-node; converting the positive and negative flux values ​​at a template node into a feature space variable according to the left eigenvector; performing nonlinear weighting on the feature space variable at the template node by a WENO-type nonlinear interpolation method, obtaining a left state value and a right state value of the feature space variable at the half-node, wherein the weighting coefficient of the nonlinear weighting is determined by a Z-type nonlinear weight; obtaining a numerical approximation of the flux at the target half-node according to the left state value, the right state value and the right eigenvector; for nodes near the object surface outside the flow field, a constrained weighted least squares virtual unit immersed boundary method is proposed to handle complex object surface boundaries well. Therefore, the present application nonlinearly weights the feature space variables at the template node by introducing a WENO-type nonlinear interpolation method, which can significantly reduce numerical oscillation while maintaining numerical stability, thereby improving the resolution and accuracy of the numerical approximation of the flux at the target half-node. In addition, the virtual cell immersed boundary method can solve the inviscid compressible flow problems around irregular geometries with high accuracy, efficiency and robustness. BRIEF DESCRIPTION OF THE DRAWINGS

[0042] In order to more clearly illustrate the embodiments of the present application or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present application. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative labor.

[0043] Figure 1 A flowchart of a numerical simulation method based on a high-precision immersed boundary method provided in an embodiment of the present application;

[0044] Figure 2 A schematic diagram of a high-precision GCIB method interpolation template provided in an embodiment of the present application;

[0045] Figure 3 A schematic diagram of a numerical simulation device based on a high-precision immersed boundary method provided in an embodiment of the present application. DETAILED DESCRIPTION

[0046] As described above, in order to cope with the computational challenges of complex geometric flow fields under Cartesian grids, the relevant technologies use the third-order WCNS-JS format and the third-order WCNS-Z format based on the third-order precision virtual cell immersed boundary method. These two formats enhance the stability of the format by introducing nonlinear weights and achieve high-order accuracy in smooth areas. The ghost cell immersed boundary method (GCIB) effectively handles the flow around complex fluid shapes by applying solid wall boundary conditions at the virtual cells outside the flow field.

[0047] However, the third-order WCNS-JS format and the third-order WCNS-Z format have order reduction phenomena at the first-order and second-order extreme points of the solution, resulting in reduced resolution and insufficient accuracy. At the same time, the accuracy of the GCIB method when calculating compressible flows is generally no more than the second order, which also leads to reduced resolution and insufficient accuracy, making it difficult to describe the fine flow field structure.

[0048] After research, the inventor proposed a numerical simulation method and device based on a high-precision immersed boundary method, the method comprising: obtaining a left eigenvector and a right eigenvector of the flux at a target half-node; converting the flux value at the template node into a feature space variable according to the left eigenvector; performing nonlinear weighting on the feature space variable at the template node by a WENO-type nonlinear interpolation method to obtain a left state value and a right state value of the feature space variable at the half-node, wherein the weighting coefficient of the nonlinear weighting is determined by a Z-type nonlinear weight; obtaining a numerical approximation of the flux at the target half-node according to the left state value, the right state value and the right eigenvector; and for nodes near the object surface outside the flow field, a constrained weighted least squares virtual unit immersed boundary method is proposed to handle complex object surface boundaries well.

[0049] Thus, firstly, by interpolating and approximating in the feature space, the present application can more accurately handle the complex changes in the flow field, especially those complex flow fields containing shock waves and rich scales. By obtaining the left eigenvector and right eigenvector of the flux at the target half-node, and converting the positive and negative flux values ​​at the template node into feature space variables based on these eigenvectors, the characteristic information in the flow field can be captured more effectively. Secondly, the present application adopts a high-precision WENO-type nonlinear interpolation method to perform nonlinear weighting on the feature space variables, which significantly improves the accuracy and stability of the calculation. In particular, the weighting coefficient is determined by the Z-type nonlinear weight. This weight design enables high-order accuracy to be achieved in the smooth area of ​​the flow field, and the weight can be automatically adjusted near the discontinuity to avoid the generation of non-physical oscillations. Finally, by converting the feature space to the physical space, the numerical approximation of the flux at the target half-node is obtained, and this step ensures the accuracy and reliability of the calculation results. Combined with the virtual unit immersed boundary method, the numerical simulation method can solve the problem of inviscid compressible flow around irregular geometric shapes with high accuracy, high efficiency and robustness.

[0050] In general, this application combines the advantages of high-precision WENO-type nonlinear interpolation and virtual unit immersed boundary method, and can solve complex flow field problems with high precision, high efficiency and robustness, providing a powerful numerical tool for scientific research and engineering applications.

[0051] In order to enable those skilled in the art to better understand the solution of the present application, the technical solution in the embodiments of the present application will be clearly and completely described below in conjunction with the drawings in the embodiments of the present application. Obviously, the described embodiments are only part of the embodiments of the present application, not all of the embodiments. Based on the embodiments in the present application, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of this application.

[0052] See also Figure 1 , which is a flow chart of a numerical simulation method based on a high-precision immersed boundary method provided in an embodiment of the present application. The method includes:

[0053] S101: Obtain the left eigenvector and the right eigenvector of the flux at the target half-node.

[0054] First, the Roe average of the state value at the half-node (i.e., the midpoint between two adjacent nodes) is calculated. Roe average is an averaging method used to process numerical flux calculations near discontinuities (such as shock waves). The core idea of ​​Roe average is to estimate the state value at the half-node by weighted averaging based on the state values ​​of adjacent nodes. This state value will be used for subsequent numerical flux calculations.

[0055] In some specific implementations, the formula for calculating the Roe average of the state value at the half-node can be shown as the following formula (1):

[0056] (1)

[0057] Among them, i+1 / 2 represents a half node, ρ i and ρ i＋1 Represent the density of the adjacent left and right nodes of the half node, and denote the speed of the left and right nodes of the half node and the speed of sound respectively. The speed of sound is , P is pressure, specific heat ratio γ=1.4.

[0058] After obtaining the Roe average state value, it can be used to calculate the Jacobian matrix of the flux F and the left and right eigenvectors of the Jacobian matrix. Among them, the flux F is usually a vector function that describes the transfer of mass, momentum or energy of the fluid in different directions. The Jacobian matrix is ​​a matrix of all partial derivatives of the function value with respect to its independent variables. The characteristic directions are vectors corresponding to the eigenvalues ​​of the Jacobian matrix A. They describe the characteristic directions of fluid flow, and fluid disturbances in these directions will propagate at specific speeds (ie, eigenvalues). It can be understood that the Jacobian matrix and its eigenvectors of the flux F are calculated to find the characteristic directions and characteristic speeds of the fluid.

[0059] In some specific implementations, the Jacobian matrix A of the flux F can be expressed as follows:

[0060] (2)

[0061] Where A is the Jacobian matrix, F is the flux, and U=(ρ,ρu,ρv,E) T is the unknown quantity, ρ is the fluid density, u and v are the velocities of the fluid in the x and y directions respectively, and E is the total energy per unit volume.

[0062] In some specific implementations, the left eigenvector of the flux F is L i＋1 / 2 , the right eigenvector is R i＋1 / 2 ,and (That is, the Jacobian matrix is ​​the product of the diagonal matrix, the left eigenvector and the right eigenvector of the flux at the target half-node). is a diagonal matrix whose diagonal elements are the four eigenvalues ​​of matrix A, namely .

[0063] It can be understood that the left and right eigenvectors of the Jacobian matrix are used to transform the flux function from the physical space to the feature space so as to perform interpolation calculation in the feature space.

[0064] S102: According to the left eigenvector, the flux value at the template node is converted into a feature space variable.

[0065] First, the flux F is split. Flux splitting is an important numerical technique that helps to design a stable and efficient numerical format to deal with complex phenomena in compressible flows. Specifically, the flux function is divided into two parts, one is the positive part and the other is the negative part, that is, F=F ﹢ ﹢F - Among them, F ± =1 / 2(F±α x U), is the maximum characteristic speed.

[0066] After the flux is split into positive and negative fluxes, these split fluxes will be used to calculate the feature space variables. This is because in the feature space, the splitting of flux can more easily handle complex flow phenomena such as shock waves and vortices. Therefore, it is necessary to construct a template based on the third-order format and convert the flux values ​​at the template nodes into Transformed into feature space variables, that is .

[0067] S103: Perform nonlinear weighting on the feature space variables at the template nodes by using a high-precision weighted essential non-oscillating WENO-type nonlinear interpolation method to obtain the left state value and the right state value of the feature space variables at the half-node, wherein the weighting coefficient of the nonlinear weighting is determined by the Z-type nonlinear weight.

[0068] When solving the inviscid Euler equation, in order to obtain a high-precision numerical solution, the spatial derivatives need to be discretized. After discretizing the spatial derivatives using the fourth-order central difference format, it is necessary to calculate the numerical approximation of the flux derivatives at the half-nodes. In view of the numerical stability of the format, the WENO-type nonlinear interpolation method based on eigendecomposition is chosen to calculate the eigenspace variable values ​​at these half-nodes.

[0069] Specifically, we first construct a three-point template S={x i-1 ,x i ,x i+1}, and perform a third-order linear approximation on the three-point template. At the same time, two sub-templates S0={x i-1 ,x i} and S1={x i ,x i＋1}, and perform second-order linear approximation on these sub-templates. Then, these three linear approximations are linearly weighted to obtain the following formula (3):

[0070] (3)

[0071] in, is the difference approximation at the half node, and formula (3) is the high-order interpolation on the three-point template S and sub-templates S0 and S1, two low-order interpolation and A hybrid format in which the three linear weights satisfy , such as the value .

[0072] Furthermore, in order to avoid numerical oscillation near discontinuities in high-order formats, nonlinear weights are used instead of linear weights. In this way, high-order accuracy can be achieved in smooth areas while avoiding interpolation across discontinuities. In some specific practical methods, the calculation formula of nonlinear weights is shown in the following formula (4):

[0073] (4)

[0074] The nonlinear weight is calculated based on the global smoothing factor and the smoothing factors of the two sub-templates. The global smoothing factor is defined as a function of some form, which is used to measure the smoothness of the data on the template. The smoothing factors of the two sub-templates are also defined similarly. The nonlinear weight in this formula is It can be expressed as the following formula (5):

[0075] (5)

[0076] in, is the nonlinear weight, is the intermediate variable, is the linear weight, is the global smoothing factor, β k is the sub-smoothing factor, ε is a very small positive number, and the global smoothing factor Defined as , the smoothing factors of the two sub-templates S0 and S1 are , the smoothness factor of sub-template S2 is . The above smooth factors are placed at point x i Performing Taylor expansion at and substituting it into formula (5) proves that: ,for ; ; All three cases hold. Therefore, the designed WCNS format can recover the designed third-order accuracy at the first-order and second-order extreme points, so that high-order accuracy can be achieved in the smooth region while maintaining numerical stability near the discontinuity.

[0077] S104: According to the left state value, the right state value and the right eigenvector, a numerical approximation of the flux at the target half-node is obtained by converting the eigenspace to the physical space.

[0078] The state values ​​in these feature spaces are converted back to the physical space to obtain the actual numerical approximation of the flux at the target half-node. This is achieved by utilizing the right eigenvector. The right eigenvector describes the transformation relationship from the feature space to the physical space, which allows the state value in the feature space to be mapped back to the flux value in the physical space. Therefore, according to the left state value, the right state value and the right eigenvector, the following formula (6) is calculated to return to the physical space and obtain the numerical approximation of the flux at the target half-node:

[0079] (6)

[0080] in, is the numerical approximation of the flux at the target half-node, R i+1 / 2 is the tensor at the target half-node, and denote the forward and reverse average state vectors at the target half-node, respectively.

[0081] It can be understood that the above embodiment is to give the flux F at the half node x i＋1 / 2 In practical applications, the above method can also be used to obtain the flux F at the half node x i－1 / 2 The numerical approximation at , and based on the flux F at the half node x i＋1 / 2 At and half nodes x i－1 / 2 The value at is approximated, and the node x is obtained i The flux derivative F x The numerical approximation of F i ’ . Among them, we get node x i The flux derivative F x The numerical approximation of F i ’ The formula can be shown as follows (7):

[0082] (7)

[0083] Among them, F i ’ For node x i The flux derivative F x The numerical approximation of and They represent the flux F at the half node x i±1 / 2 The values ​​at are approximate.

[0084] It is understandable that the above WCNS format is used to calculate the entire flow field information, but when the grid nodes are located near the object surface, some template points are located outside the flow field. Therefore, the physical quantities on these template points need to be given. For complex curved boundaries, the third-order accuracy GCIB method for calculating the physical quantities on these template points is given below:

[0085] S105: constructing a binary quadratic polynomial of the objective function at the virtual unit according to the flow field information in the interpolation template, wherein the interpolation template is determined based on the projection point of the virtual unit on the object plane, the object plane normal, the inward unit normal vector and the corresponding tangent vector.

[0086] All cells in the computational domain are divided into different types of cells, including solid cells, fluid cells, and virtual cells. Solid cells represent the area occupied by solid substances and usually do not participate in fluid calculations. Fluid cells represent the area occupied by fluid substances and are the main objects of fluid calculations. Virtual cells are a layer of cells located outside the boundary of the fluid domain and are used to handle boundary conditions, especially when the boundary is a complex shape.

[0087] The virtual unit is set up to provide a continuous flow field information on the boundary so that the boundary conditions can be applied more effectively. To achieve this, an interpolation template needs to be constructed for the virtual unit. The interpolation template is an elliptical area centered on the projection point and determined by the normal vector and the tangent vector. The part of this elliptical area inside the flow field is used as the interpolation template, which contains the flow field data points used for interpolation calculation. Among them, the projection point refers to the projection point of the virtual unit on the object plane (i.e., the boundary). The normal vector refers to the inward unit normal vector obtained by passing through the projection point as the normal line of the object plane. The tangent vector refers to the tangent vector perpendicular to the normal vector, which is used to determine the direction of the interpolation template. See. Figure 2 , which is a schematic diagram of a high-precision GCIB method interpolation template provided in an embodiment of the present application. Figure 2 Medium g = (x g ,y g ) is a virtual unit, and P0 is its projection point on the object plane. The object plane normal through P0, the inward unit normal vector and the corresponding tangent vector are denoted as n and With P0 = (x0, y0) as the center, the long and short radii are r a ,r b Draw an ellipse, and the part falling inside the flow field is denoted as S P0, which is the interpolation template. It can be seen that most of the data information required for interpolation is concentrated near the normal direction.

[0088] According to the interpolation template S P0 Based on the flow field information in the virtual unit, a quadratic polynomial is constructed to approximate the value of the objective function q(x, y) at the virtual unit. In a specific implementation, the quadratic polynomial of the objective function q(x, y) can be expressed as follows:

[0089] (8)

[0090] is the quadratic polynomial of the objective function q(x,y), k and l are indices, a kl are polynomial coefficients, P0 = (x0, y0) is the center, △x and △y are the uniform mesh sizes in the x and y directions.

[0091] It should be noted that in order to make this polynomial accurately reflect the interpolation template S P0 The flow field information inside needs to ensure that the Vandermonde matrix A (A·a=b, a is the unknown vector in formula (8), b is S P0 The known flow field information values ​​in the matrix are M×6 columns of full rank, that is, the column vectors of the matrix are linearly independent. This usually requires that the flow field data points in the interpolation template are no less than 9 (because the two-variable quadratic polynomial has 9 unknown coefficients).

[0092] S106: Using a constrained least squares method to fit a quadratic polynomial of the objective function, a third-order approximation at the virtual unit is obtained.

[0093] Once the Vandermonde matrix and the known flow field information value vector are determined, the unknown coefficient vector can be solved by solving a set of linear equations. Specifically, in order to obtain the best approximation of the objective function q(x,y), the least squares fitting method can be used, that is, the optimal coefficient a is obtained by minimizing the following weighted cost function. Among them, the least squares fitting method is a mathematical optimization technique that finds the best function match for the data by minimizing the sum of squares of the error. In this step, the least squares method is used to fit the binary quadratic polynomial of the objective function, that is, to find a polynomial that minimizes the sum of squares of the error on the flow field data points within the interpolation template. The cost function (also called the loss function or objective function) is a function that measures the difference between the fitted polynomial and the actual flow field information. In some specific implementations, the cost function is shown in the following formula (9):

[0094] (9)

[0095] Where J(a) is the cost function, M is the number of training data sets, represents the true value q(x i ,y i ) and the model prediction value The square of the difference between is the distance |P g P i |Related weight function. The weight function usually assigns weights according to the distance of the data point from the interpolation center, so that the data points with closer distance have a greater impact on the fitting result.

[0096] At the same time, the fitting function is required Satisfy the following physical surface conditions (referring to the conditions that the fitting function must satisfy on the boundary), namely Among them, g b (x0,y0) is the given boundary value. α D =1,β N = 0 corresponds to Dirichlet-type boundary conditions, while α D =0,β N =1 corresponds to Neumann-type boundary conditions.

[0097] Due to the constraints of the object surface conditions, the Lagrange multiplier method is used to solve the problem. The Lagrange multiplier method is a mathematical method used to deal with constrained optimization problems. It introduces Lagrange multipliers to transform the constraints into part of the objective function, thereby transforming it into an unconstrained optimization problem. The formula of the Lagrange multiplier method is shown in the following formula (10):

[0098] (10)

[0099] in, is the augmented Lagrangian function, used to solve problems with constraints, J(a) is the cost function, is the Lagrange multiplier, is a constraint, is a matrix, a is a model parameter vector, g b is a constant vector. The least squares solution to the above problem satisfies .in, The KKT matrix is ​​an important matrix that appears when solving constrained optimization problems. It contains information such as the gradient of the objective function and the gradient of the constraints. In this step, the KKT matrix is ​​required to have certain properties (such as positive definiteness or full rank) to ensure that the constrained least squares problem has a unique solution. Therefore, by solving the constrained least squares problem, the fitting polynomial is obtained. The coefficients a determine the shape and size of the polynomial, allowing the virtual cell P to be estimated based on the flow field information within the interpolation template.g The third-order approximation of the objective function at It can be understood that this third-order approximation is obtained based on the flow field information and object surface conditions in the interpolation template, so it should be able to better reflect the actual flow field conditions at the virtual unit.

[0100] S107: When there is a large gradient or discontinuity in the flow field near the object surface, the WENO type extrapolation method is used, combined with the geometric characteristics of the object surface, to take equidistant interpolation points in the normal direction.

[0101] In fluid dynamics simulations, the flow field near the object surface is often very complex, and there may be large gradients or discontinuities (such as shock waves, boundary layer separation, etc.). These phenomena place high demands on the stability and accuracy of numerical methods. When high-order interpolation polynomials are used to derive the extrapolated values ​​of virtual cells, if there are large gradients or discontinuities in the flow field, the interpolation polynomials may produce false oscillations in these areas. False oscillations refer to unnecessary fluctuations in the numerical solution near the true solution, which affects the accuracy and stability of the numerical solution.

[0102] In order to overcome the problem of false oscillations caused by high-order interpolation polynomials in large gradient or discontinuous areas of the flow field, step S107 adopts the weighted essentially non-oscillatory (WENO) type extrapolation method. The WENO method is a high-order numerical method that combines the weighted average of multiple low-order interpolation polynomials to construct a high-order interpolation polynomial. By dynamically adjusting the weights, the WENO method can maintain high-order accuracy in smooth areas of the flow field while suppressing false oscillations in large gradient or discontinuous areas of the flow field.

[0103] When using the WENO type extrapolation method, step S107 also combines the geometric characteristics of the object surface. This means that when deriving the extrapolated value of the virtual unit, factors such as the shape, curvature and direction of the object surface are considered. This helps to more accurately capture the flow field characteristics near the object surface and improve the accuracy of the numerical solution. Specifically, it is necessary to combine the geometric characteristics of the object surface and take equidistant interpolation points in the normal direction of the object surface n = (cosθ, sinθ), that is, the following formula (11):

[0104] (11)

[0105] in, This approach helps maintain the stability and accuracy of the interpolation polynomial, especially when there are large gradients or discontinuities in the flow field.

[0106] S108: Perform nonlinear weighting based on the equidistant interpolation points to obtain a third-order approximation at the virtual unit, wherein a weight coefficient of the nonlinear weighting is determined according to a smoothing factor at the equidistant interpolation points.

[0107] After obtaining the equidistant interpolation points, these equidistant interpolation points are nonlinearly weighted to obtain the third-order approximation at the virtual unit. Among them, the weight coefficient of the nonlinear weighting is determined according to the smoothness factor at the equidistant interpolation point. The smoothness factor is an indicator to measure the smoothness of the flow field near the point. If the flow field is smooth near a certain point, the smoothness factor of the point will be smaller; if the flow field has a large gradient or discontinuity near a certain point, the smoothness factor of the point will be larger. In some specific implementations, the weight coefficient of the nonlinear weighting can be shown as follows (12):

[0108] (12)

[0109] in, is the smoothness factor of the two-point template. Formula (12) can be proved by using the mean inequality , ensuring the rationality and effectiveness of the weight coefficients of nonlinear weighting.

[0110] In some specific implementations, for Dirichlet-type boundary conditions, it is assumed that the known information is a first-order approximation of the virtual cell value. Then, the nonlinear weighting of the first-order approximation and the third-order approximation of the virtual cell value is defined by the following formula (13).

[0111] (13)

[0112] Therefore, q(P g ) is the nonlinear weighted result of the first-order approximation and the third-order approximation of the virtual unit value, is the weight value, P g = (x g ,y g ) is a virtual unit. The design of formula (13) enables the virtual unit value to achieve third-order accuracy when the flow field near the object surface is smooth. This is because in the smooth area, the weight coefficient tends to select those interpolation points that contribute more to the third-order approximation. However, when there are large gradients or discontinuities near the object surface, the weight coefficient is automatically adjusted to reduce the dependence on unstable interpolation points, thereby ensuring the stability of the format. Finally, by performing Taylor expansion on the virtual unit value in the normal direction near the object surface, it can be proved that when the flow field near the object surface is smooth, .

[0113] In summary, the embodiment of the present application provides a numerical simulation method based on a high-precision immersed boundary method, the method comprising: obtaining the left eigenvector and the right eigenvector of the flux at the target half-node; according to the left eigenvector, converting the flux value at the template node into a feature space variable; nonlinearly weighting the feature space variable at the template node by a WENO-type nonlinear interpolation method, obtaining the left state value and the right state value of the feature space variable at the half-node, wherein the weighting coefficient of the nonlinear weighting is determined by the Z-type nonlinear weight; according to the left state value, the right state value and the right eigenvector, the numerical approximation of the flux at the target half-node is obtained; for the nodes near the object surface outside the flow field, a constrained weighted least squares virtual unit immersed boundary method is proposed to handle the complex object surface boundary well. Therefore, the present application introduces a WENO-type nonlinear interpolation method to perform nonlinear weighting on the feature space variables at the template node, which can significantly reduce numerical oscillation while maintaining numerical stability, thereby improving the resolution and accuracy of the numerical approximation of the flux at the target half-node. In addition, the virtual cell immersed boundary method can solve the inviscid compressible flow problems around irregular geometries with high accuracy, efficiency and robustness.

[0114] See also Figure 3 , which is a schematic diagram of a numerical simulation device based on a high-precision immersed boundary method provided in an embodiment of the present application. The numerical simulation device 300 includes: a vector acquisition module 301, a variable acquisition module 302, a state acquisition module 303 and a numerical approximation module 304. Among them, the vector acquisition module 301 is used to obtain the left eigenvector and the right eigenvector of the flux at the target half node;

[0115] A variable acquisition module 302 is used to convert the positive flux value and the negative flux value at the template node into a feature space variable according to the left eigenvector;

[0116] The state acquisition module 303 is used to perform nonlinear weighting on the feature space variables at the template nodes by using a high-precision weighted essential non-oscillating WENO type nonlinear interpolation method to obtain the left state value and the right state value of the feature space variables at the half node, wherein the nonlinear weighting coefficient is determined by the Z-type nonlinear weight;

[0117] The numerical approximation module 304 is used to obtain the numerical approximation of the flux at the target half node by converting the feature space into the physical space according to the left state value, the right state value and the right eigenvector.

[0118] In some specific implementations, the weighting coefficient of the nonlinear weighting is determined by the Z-type nonlinear weight and the global smoothing factor, wherein the formula of the global smoothing factor is as follows:

[0119]

[0120] in, is the global smoothing factor, β n is the sub-smoothing factor of the nth sub-template.

[0121] In some specific implementations, the vector acquisition module 301 specifically includes: a first acquisition submodule and a second acquisition submodule;

[0122] The first acquisition submodule is used to obtain the Jacobian matrix of the flux, where the Jacobian matrix is ​​the matrix of all partial derivatives of the flux with respect to its independent variable, and the independent variable is U=(ρ,ρu,ρv,E) T , where ρ is the fluid density, u and v are the velocities of the fluid in the x and y directions respectively, and E is the total energy per unit volume;

[0123] The second acquisition submodule is used to perform eigendecomposition on the Jacobian matrix to obtain the left eigenvector and the right eigenvector of the flux at the target half-node, wherein the Jacobian matrix is ​​the product of a diagonal matrix, the left eigenvector and the right eigenvector of the flux at the target half-node.

[0124] In some specific implementations, the numerical simulation device 300 further includes: a polynomial construction module and a polynomial fitting module;

[0125] A polynomial construction module is used to construct a two-variable quadratic polynomial of the objective function at the virtual unit for the nodes near the object surface outside the flow field according to the flow field information in the interpolation template, wherein the interpolation template is determined based on the projection point of the virtual unit on the object surface, the object surface normal, the inward unit normal vector and the corresponding tangent vector;

[0126] The polynomial fitting module is used to fit the bivariate quadratic polynomial of the objective function by using the constrained least square method to obtain the third-order approximation at the virtual unit.

[0127] In some specific implementations, the numerical simulation device 300 further includes: an interpolation point acquisition module and an approximate value acquisition module;

[0128] The interpolation point acquisition module is used to obtain equidistant interpolation points in the normal direction by using the WENO type extrapolation method when there is a large gradient or discontinuity in the flow field near the object surface, combined with the geometric characteristics of the object surface;

[0129] The approximate value acquisition module is used to perform nonlinear weighting based on equidistant interpolation points to obtain a third-order approximate value at a virtual unit, wherein the weight coefficient of the nonlinear weighting is determined according to a smoothing factor at the equidistant interpolation points.

[0130] In summary, the embodiment of the present application provides a numerical simulation device based on the high-precision immersed boundary method. The present application introduces a WENO-type nonlinear interpolation method to perform nonlinear weighting on the feature space variables at the template nodes. While maintaining numerical stability, it can significantly reduce numerical oscillations, thereby improving the resolution and accuracy of the numerical approximation of the flux at the target half-node.

[0131] It should be noted that each embodiment in this specification is described in a progressive manner, and the same or similar parts between the embodiments can be referred to each other, and each embodiment focuses on the differences from other embodiments. In particular, for the device and system embodiments, since they are basically similar to the method embodiments, the description is relatively simple, and the relevant parts can be referred to the partial description of the method embodiments. The device and system embodiments described above are merely schematic, wherein the units described as separate components may or may not be physically separated, and the components indicated as units may or may not be physical units, that is, they may be located in one place, or they may be distributed on multiple network units. Some or all of the modules may be selected according to actual needs to achieve the purpose of the scheme of this embodiment. Ordinary technicians in this field can understand and implement it without paying creative labor.

[0132] The above is only a specific implementation of the present application, but the protection scope of the present application is not limited thereto. Any changes or substitutions that can be easily thought of by a person skilled in the art within the technical scope disclosed in the present application should be included in the protection scope of the present application. Therefore, the protection scope of the present application should be based on the protection scope of the claims.

Claims

1. A numerical simulation method based on a high-precision immersed boundary method, characterized in that: The method comprises: Get the Roe average of the state value at the target half node; According to the Roe average of the state value at the target half-node, the Jacobian matrix of the flux is obtained, and the Jacobian matrix is ​​the matrix of all partial derivatives of the flux with respect to its independent variable, and the independent variable is U=(ρ,ρu,ρv,E) T , where ρ is the fluid density, u and v are the velocities of the fluid in the x and y directions respectively, and E is the total energy per unit volume; Performing eigendecomposition on the Jacobian matrix to obtain a left eigenvector and a right eigenvector of the flux at a target half-node, wherein the Jacobian matrix is ​​a product of a diagonal matrix, a left eigenvector and a right eigenvector of the flux at a target half-node; According to the left eigenvector, the positive flux value and the negative flux value at the template node are converted into feature space variables; The feature space variables at the template nodes are nonlinearly weighted by a high-precision weighted essential non-oscillating WENO-type nonlinear interpolation method to obtain the left state value and the right state value of the feature space variables at the half-node, wherein the weighting coefficient of the nonlinear weighting is determined by the Z-type nonlinear weight and the global smooth factor, wherein the formula of the global smooth factor is as follows: ； in, is the global smoothing factor, β n is the sub-smoothing factor of the nth sub-template, and the global smoothing factor , , ; According to the left state value, the right state value and the right eigenvector, a numerical approximation of the flux at the target half-node is obtained by converting the eigenspace into the physical space.

2. The numerical simulation method based on the high-precision immersed boundary method according to claim 1 is characterized in that: The method further comprises: For nodes near the object surface outside the flow field, a two-variable quadratic polynomial of the objective function at the virtual unit is constructed according to the flow field information in the interpolation template, wherein the interpolation template is determined based on the projection point of the virtual unit on the object surface, the object surface normal, the inward unit normal vector and the corresponding tangent vector; The constrained least square method is used to fit the bivariate quadratic polynomial of the objective function to obtain a third-order approximation at the virtual unit.

3. The numerical simulation method based on the high-precision immersed boundary method according to claim 1 is characterized in that: The method further comprises: When there is a large gradient or discontinuity in the flow field near the object surface, the WENO type extrapolation method is used, and equidistant interpolation points are taken in the normal direction in combination with the geometric characteristics of the object surface; Nonlinear weighting is performed based on the equidistant interpolation points to obtain a third-order approximation at the virtual unit, wherein the weight coefficient of the nonlinear weighting is determined according to the smoothing factor at the equidistant interpolation points.

4. A numerical simulation device based on a high-precision immersed boundary method, characterized in that: The device comprises: a vector acquisition module, a variable acquisition module, a state acquisition module and a numerical approximation module; The vector acquisition module is used to obtain the Roe average of the state value at the target half-node; according to the Roe average of the state value at the target half-node, the Jacobian matrix of the flux is obtained, and the Jacobian matrix is ​​the matrix of all partial derivatives of the flux with respect to its independent variable, and the independent variable is U=(ρ,ρu,ρv,E) T , where ρ is the fluid density, u and v are the velocities of the fluid in the x and y directions respectively, and E is the total energy per unit volume; perform eigendecomposition on the Jacobian matrix to obtain the left eigenvector and the right eigenvector of the flux at the target half-node, where the Jacobian matrix is ​​the product of a diagonal matrix, the left eigenvector and the right eigenvector of the flux at the target half-node; The variable acquisition module is used to convert the positive flux value and the negative flux value at the template node into a feature space variable according to the left eigenvector; The state acquisition module is used to perform nonlinear weighting on the feature space variables at the template node by a high-precision weighted essential non-oscillating WENO type nonlinear interpolation method to obtain the left state value and the right state value of the feature space variable at the half node, wherein the weighting coefficient of the nonlinear weighting is determined by the Z-type nonlinear weight and the global smooth factor, wherein the formula of the global smooth factor is as follows: ； in, is the global smoothing factor, β n is the sub-smoothing factor of the nth sub-template, and the global smoothing factor , , ; The numerical approximation module is used to obtain the numerical approximation of the flux at the target half-node by converting the feature space into the physical space according to the left state value, the right state value and the right eigenvector.

5. The numerical simulation device based on the high-precision immersed boundary method according to claim 4, characterized in that: The device also includes: a polynomial construction module and a polynomial fitting module; The polynomial construction module is used to construct a two-variable quadratic polynomial of the objective function at the virtual unit for the nodes near the object surface outside the flow field according to the flow field information in the interpolation template, wherein the interpolation template is determined based on the projection point of the virtual unit on the object surface, the object surface normal, the inward unit normal vector and the corresponding tangent vector; The polynomial fitting module is used to fit the binary quadratic polynomial of the objective function by using the constrained least square method to obtain a third-order approximation at the virtual unit.

6. The numerical simulation device based on the high-precision immersed boundary method according to claim 4, characterized in that: The device further comprises: an interpolation point acquisition module and an approximate value acquisition module; The interpolation point acquisition module is used to adopt the WENO type extrapolation method and take equidistant interpolation points in the normal direction in combination with the geometric characteristics of the object surface when there is a large gradient or discontinuity in the flow field near the object surface; The approximate value acquisition module is used to perform nonlinear weighting based on the equidistant interpolation points to obtain a third-order approximate value at the virtual unit, wherein the weight coefficient of the nonlinear weighting is determined according to the smoothness factor at the equidistant interpolation points.

Images