Wing boundary simulation method, terminal device and storage medium based on airfoil curved surface normal
By simulating the wing boundary based on the normal direction of the airfoil curved surface, and using the Simpson rule and Newton iteration method to calculate the arc length and Gaussian integration points, the problem of numerical integration accuracy of the airfoil curved boundary in unstructured grids is solved, thereby improving the accuracy of the simulation results and the computational efficiency.
Patent Information
- Application Number
- CN202510767049.8
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- Filing Date
- 2025-06-10
- Publication Date
- 2025-09-23
- Estimated Expiration
- 2045-06-10
AI Technical Summary
In the existing unstructured grid finite volume numerical simulation method, it is difficult to accurately reflect the flight state of the wing, especially the numerical integration accuracy at the curved boundary of the airfoil is insufficient, which affects the accuracy of the calculation results.
A wing boundary simulation method based on the normal direction of the airfoil curved surface is adopted. By obtaining the wing curved boundary information, the arc length is calculated using the Simpson rule and the Simpson 3/8 rule. The Gaussian integration point is determined in combination with the Newton iteration method, and the normal direction of the curved boundary is calculated to achieve high-precision numerical integration.
The accuracy of numerical integration calculation and the precision of simulation results are improved, the calculation process is simplified, the calculation complexity is reduced, and the accuracy of wing flow field calculation is improved.
Smart Images

Figure CN120277929B_ABST
Abstract
Description
Technical Field
[0001] The present invention relates to aircraft wing simulation technology, in particular to a wing boundary simulation method based on the normal direction of an airfoil curved surface, a terminal device and a storage medium. Background Art
[0002] The airfoil, commonly known as the wing profile or blade cross-section, is one of the core factors affecting the overall performance of an aircraft. The airfoil is a fundamental element in the design of aerodynamic surfaces such as aircraft wings, tail planes, helicopter rotors, propellers, and wind turbine blades. It directly influences the aerodynamic performance of an aircraft, including key aerodynamic characteristics such as lift, drag, and pitching moment. The aerodynamic characteristics of the airfoil also directly affect the maneuverability and stability of an aircraft. A well-designed airfoil can provide sufficient control torque, ensuring stability and maneuverability under various flight conditions, thereby improving flight safety.
[0003] In addition, airfoils are closely related to energy conservation, emission reduction, and consumption reduction of aircraft. Optimizing airfoil design is one of the effective measures to reduce drag, energy consumption, and fuel consumption. Airfoil design determines the lift-to-drag ratio of an aircraft [Li Yunpeng, Han Yongzhi. Research Progress on Design of Lift-Enhancing Devices Based on Laminar Flow Wings [J]. Progress in Aeronautical Engineering, 2021, 12(4): 1-11.], that is, the ratio of lift to drag. The greater the lift-to-drag ratio, the higher the flight efficiency, which can provide a longer range under limited fuel conditions while reducing environmental pollution (Allison, E., Kroo, I., Sturdza, P., Suzuki, Y., & Martins-Rivas, H. (2010). Aircraft conceptual design with natural laminar flow. InProceedings of the 27th Congress of the International Council of theAeronautical Sciences (pp. 428–436). Nice, France: ICAS.). Excellent airfoil design can achieve a high lift-to-drag ratio, effectively reduce fuel consumption, lower carbon emissions, increase range and endurance, which is crucial to improving the flight efficiency and economy of aircraft.
[0004] In order to economically and efficiently obtain excellent airfoils with strong stability and good drag reduction effects, computational fluid dynamics methods are often used in engineering applications to achieve the design, simulation, and optimization of airfoils. Computational fluid dynamics is a numerical simulation method that discretizes the physical space into a computational space composed of grids, and then uses time and space numerical formats to simulate the complex flow problems of the wing during flight. The process of discretizing the physical space is called grid generation. Commonly used grids in engineering applications include structured grids, unstructured grids, and Cartesian grids. Structured grids have high generation efficiency in simple geometric shapes, but grid generation becomes increasingly difficult in complex situations. Cartesian grids are simple to generate and can provide high computational accuracy within a regular computational domain, but complex geometric shapes are prone to mismatch between the grid and the boundary, requiring special processing methods such as the immersed boundary method, which increases the complexity of implementation. Unstructured grids have strong automatic generation capabilities and can flexibly adapt to arbitrarily complex geometric shapes. They have been widely used in practical engineering problems in recent decades. Numerical schemes for unstructured grids include finite volume (FV), spectral volume (SV), discontinuous Galerkin (DG), PNPM, and correction procedure via reconstruction (CPR). Finite volume schemes are relatively simple to derive and implement, and can be easily extended to high-order accuracy. In recent years, the introduction of the compact finite volume (CFV) scheme and mature shock wave capture techniques have continuously increased the influence of finite volume schemes.
[0005] In order for unstructured grid finite volume numerical simulation methods to truly reflect the flight state of an airfoil, two aspects require special attention. First, the mesh edges of unstructured grids generated in engineering applications are all line segments, while airfoils generally have curved boundaries. Replacing curved boundaries with line segments results in a slight difference between the numerical shape of the wing and the actual shape, which may lead to inaccurate numerical simulations. Therefore, it is necessary to express the airfoil with a high-order curved boundary. Second, to improve computational accuracy and solution efficiency, high-precision numerical formats are often used for solutions in engineering applications. However, high-precision formats are often used for calculations within the calculation region. Numerical formats on the curved boundaries of airfoils often have difficulty obtaining high-order accuracy. The values at the airfoil boundaries directly affect the calculation results of the entire calculation region. In other words, numerical simulations of airfoils with curved boundaries cannot achieve high-order accuracy, which will inevitably affect the overall accuracy.
[0006] Finite volume schemes discretize conservation law equations in integral form. Therefore, achieving high-order accuracy relies on matching the scheme's accuracy with high-order numerical integration, such as surface integrals of fluxes and volume integrals of source terms. High-order numerical integration is relatively straightforward for straight-edge elements. Accurate results can be obtained by following the integration points and weights described in the literature (C. Ollivier-Gooch, A. Nejat, K.Michalak, Obtaining and Verifying High-Order Unstructured Finite Volume Solutions to the Euler Equations, AIAA Journal 47 (2009) 2105–2120. https: / / doi.org / 10.2514 / 1.40585.). However, for elements with curved edges, particularly those with curved airfoils, the numerical integration process presents numerous challenges. Two common approaches exist. The first is the isoparametric element method, commonly used in the finite element field. This method establishes a mapping between physical space and isoparametric space to achieve equivalent computations on isoparametric elements. However, this method requires defining the distribution of parameter points, solving and storing the Lagrange interpolation coefficients and Jacobians, so the process is somewhat cumbersome, especially in high-order isoparametric elements, which will further increase the computational complexity. In 2006, Krivodonova and Berger proposed a new curved element processing method for the DG format. It is based on straight-edge element calculations and uses the normal on the surface at the integration point, which significantly reduces the complexity of curved element processing, but the overall accuracy cannot be guaranteed. Li and Nishikawa used this method in their work and achieved high-order accuracy in the DG / FV and EB3 (Edge-BasedThird-Oirder) formats, but the accuracy of the calculation results is not as good as the isoparametric method, and the accuracy can be further improved. Summary of the Invention
[0007] The technical problem to be solved by the present invention is to provide a wing boundary simulation method, terminal device and storage medium based on the normal of the airfoil curved surface, so that the unstructured grid finite volume numerical simulation method can truly reflect the flight state of the wing in view of the shortcomings of the existing technology.
[0008] To solve the above technical problems, the technical solution adopted by the present invention is: a wing boundary simulation method based on the normal direction of the airfoil curved surface, comprising the following steps:
[0009] S1. Obtain wing curve boundary information;
[0010] S2. Determine the integral degree of the curved unit based on the wing curved boundary information. If the integral degree is 3, use Simpson's rule to calculate the arc length L; if the integral degree is 4, use Simpson's 3 / 8 rule to calculate the arc length L;
[0011] S3. Use Newton iteration method to calculate Gaussian integration points of wing curvature boundary;
[0012] S4, calculate the normal direction of the Gaussian integration point at the wing curvature boundary;
[0013] S5. Substitute the normal of the Gaussian integral point into the straight boundary discretization method formula to calculate the numerical integral:
[0014] ;
[0015] in, It represents the average value of wing speed or temperature in the unit cell. represents the wing unit volume, represents the number of unit faces, represents the number of unit area points, is the surface integral coefficient, represents the convection and viscous flux at the area integral point k, represents the number of unit volume points, represents the volume fraction coefficient, Represents the heat source at volume point q.
[0016] During the simulation of aircraft wing boundaries, inaccurate numerical integration can lead to inaccurate calculations of the wing's flow field, which in turn causes erroneous simulation results. This method obtains wing curved boundary information, determines the integration degree of the curved unit, determines the arc length based on the different integration degrees, and then determines the Gaussian integration points based on the arc length. Finally, it obtains the numerical integral, greatly improving the accuracy of the numerical integration calculation and reducing the computational complexity.
[0017] The wing boundary curve calculated using Simpson's rule is in the interval The expression of the arc length L on is:
[0018] ;
[0019] in, , n is The number of subintervals, 、 、 They are the airfoil curve function f(x) in 、 、 The first derivative at , 、 、 They represent the horizontal coordinates of point a, the i-th integration point, and point b respectively.
[0020] The wing boundary curve calculated using Simpson's rule is in the interval The expression of the arc length L on is:
[0021] ;
[0022] in, , n is The number of subintervals.
[0023] In step S3, the specific implementation process of using the Newton iteration method to solve the Gaussian integral point of the airfoil curved boundary includes:
[0024] 1) Set two parameters: ; , and are the positions of the two integration points, 、 is the target arc length corresponding to the two integration point positions, and are the arc lengths from the starting point to the two integration points respectively;
[0025] 2) When and When the error is greater than the tolerance, the positions of the two integration points are updated using the following formula: and ;
[0026] 3) Repeat step 2) until and If the error is less than or equal to the error tolerance, the updated integration point position is output.
[0027] In step S4, the normal calculation formula of the Gaussian integral point at the wing curve boundary is:
[0028] ;
[0029] in, 、 are the first derivatives of the wing boundary curve at the two integration points respectively.
[0030] As an inventive concept, the present invention also provides a terminal device, including a memory, a processor, and a computer program stored in the memory; the processor executes the computer program to implement the steps of the above method.
[0031] As an inventive concept, the present invention also provides a computer-readable storage medium having a computer program / instruction stored thereon; the computer program / instruction implements the steps of the above method when executed by a processor.
[0032] As an inventive concept, the present invention also provides a computer program product, comprising a computer program / instruction; when the computer program / instruction is executed by a processor, the steps of the above method are implemented.
[0033] Compared with the prior art, the present invention has the following beneficial effects:
[0034] 1. The method of the present invention does not require the definition of the distribution of parameter points, and the calculation process is simple;
[0035] 2. After obtaining the wing curved boundary information, the present invention judges the integration degree of the curved unit, determines the arc length according to different integration degrees, and then determines the Gaussian integration point according to the arc length, and finally obtains the numerical integral, which greatly improves the accuracy of the numerical integration calculation, and indirectly improves the calculation accuracy of the simulation results. BRIEF DESCRIPTION OF THE DRAWINGS
[0036] Figure 1 It is a quadrilateral mesh NACA0012 airfoil;
[0037] Figure 2 Schematic diagram of spatial discretization of high-order finite volume format;
[0038] Figure 3 Schematic diagram of the calculation method for straight-edge elements based on curved surface normals; (a) curved boundary and curved normal, (b) straight boundary and curved normal;
[0039] Figure 4 This is a flow chart of a method according to an embodiment of the present invention;
[0040] Figure 5 The distribution of three sets of grids near the airfoil surface; (a) the sparsest, (b) the medium, and (c) the densest;
[0041] Figure 6 is the negative pressure coefficient (-C) of the airfoil surface obtained by different methods at three different attack angles p ); (a) 3 degrees angle of attack, (b) 4 degrees angle of attack, (c) 5 degrees angle of attack;
[0042] Figure 7 is the entropy error distribution of the wall unit; (a) is the sparsest, (b) is the densest;
[0043] Figure 8 The meshes for the flow around a viscous NACA0012 airfoil are: (a) the sparsest (14200 meshes), (b) the medium (15620 meshes), and (c) the densest (17040 meshes).
[0044] Figure 9 Three sets of grid local positions reconstruct the meshes near the airfoil curve and wall; (a) sparsest position 1, (b) medium position 1, (c) densest position 1, (d) sparsest position 2, (e) medium position 2, (f) densest position 2;
[0045] Figure 10 The pressure coefficient C obtained by different methods p ; (a) sparsest, (b) medium, (c) densest;
[0046] Figure 11 The integral C of the resistance coefficient obtained by different methods D ; (a) sparsest, (b) medium, (c) densest;
[0047] in, Figure 6 、 Figure 7 、 Figure 9 、 Figure 10 The horizontal axis is the horizontal coordinate of the wing boundary point, Figure 9 The vertical axis is the vertical coordinate of the wing boundary point. DETAILED DESCRIPTION
[0048] To make the objectives, technical solutions, and advantages of the embodiments of the present invention more clear, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the accompanying drawings 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 making creative efforts shall fall within the scope of protection of the present invention. Example
[0049] The airfoil curved boundary expression requires the use of a curve function to represent the boundary of the airfoil's geometric model, which involves the mapping from the geometric model to the computational model.
[0050] Mapping the geometric model to the computational model requires three steps: generating a mesh, building a coordinate system, and calculating coordinates. Figure 1 Taking the NACA0012 airfoil shown in the figure as an example, the three steps are explained.
[0051] 1. Generate a mesh
[0052] like Figure 1 As shown in the figure, firstly, the geometric model boundary of the NACA0012 airfoil is filled with points for discretization of the geometric boundary. These points are called grid points, and adjacent grid points are connected to form grid edges. After the grid edges on the geometric boundary are formed, the mesh generation software can be used to generate the mesh of the external calculation area of the NACA0012 airfoil, as shown in the figure. Figure 1As shown, a quadrilateral mesh is generated, and a triangular mesh can also be generated. Although the mesh edges of unstructured meshes are composed of line segments, the airfoil's curved boundary can be accurately expressed using the high-order curved boundary representation method (Gao H, Wang Z, Liu Y. A study of curved boundary representations for 2D high-order Euler solvers [J]. Journal of Scientific Computing. 2010, 44: 323–336.). This provides the basis for constructing a high-order precision format for airfoil de-elementing and performing high-order numerical integration.
[0053] 2. Constructing the coordinate system
[0054] You need to construct the XOY plane. Set the leading edge of the NACA0012 airfoil to point O. Set the x-axis along the chord line of the NACA0012 airfoil, the line connecting the leading and trailing edges, with the positive x-axis direction. Set the y-axis perpendicular to the x-axis through point O, with the positive y-axis direction pointing upward.
[0055] 3. Calculate coordinates
[0056] The construction of the numerical format is related to the grid point coordinates. For each grid point in the calculation model, the grid coordinates of the grid point in the XOY plane are calculated according to the grid scale. .
[0057] After completing the above three steps, the high-order precision numerical format of the airfoil curved element can be constructed.
[0058] The following introduces the discretization method of airfoil curved unit and straight edge unit.
[0059] (1) High-order accuracy unstructured finite volume scheme
[0060] The high-order unstructured finite volume scheme requires solving the integral form of the Navier-Stokes (NS) equations. The integral conservation law equations are:
[0061]
[0062] Where u is the conserved variable, s is the source term, and F c With F v are the convection and viscous fluxes respectively. The above formula can be converted to:
[0063]
[0064] in, is the conserved variable in the control volume The average value on and are the convection and viscous fluxes along the normal direction outside the unit surface, respectively. According to the discretization scheme of the high-order finite volume format, on the control volume element i, the equation The semi-discrete form of can be expressed as:
[0065]
[0066] Combine Figure 1 It can be seen that represents the number of element faces of the control volume element i, and represent the number of surface integral and volume integral points respectively, and and are the surface integral and volume integral coefficients. Depending on the accuracy requirements of the format, the number, location, and integral weight coefficients of the surface integral and volume integral points will vary.
[0067] In addition, the formula in represents the numerical flux at the integration point k, which is given by the convection and stickiness It consists of two parts:
[0068]
[0069] The numerical convection flux can be calculated using the Roe format:
[0070]
[0071] in, is the Jacobian based on the Roe averaging method, and are the left and right values of the variables at the integration point. In contrast, the viscous numerical flux calculation needs to rely on the variable gradient at the interface. , can be calculated using the α-damping format:
[0072]
[0073] Wherein, α is a constant, and the calculation in the embodiment of the present invention is , and are the length and unit vector of the line connecting the reference points of element i and its adjacent element j:
[0074]
[0075] (2) Straight-edge unit discretization method
[0076] When integrating curved airfoil elements, isoparametric elements or straight-edge elements are often used. The isoparametric element method offers a rigorous mathematical derivation and accurately derives the surface integral of the flux and volume integral of the source term over the curved element, providing the foundation for achieving high-order accuracy. However, the overall process is somewhat cumbersome. Based on this current situation, Krivodonova and Berger proposed a simplified scheme for the DG format that avoids the construction of isoparametric elements and the calculation and storage of Jacobians.
[0077] like Figure 3 As shown, for airfoil elements with curved edges, they can be directly calculated as straight-edge elements. That is, the curved boundary is treated as a straight edge on the curved boundary, and the corresponding boundary length is directly solved according to formula (7). Then, according to formula (3), a high-order finite volume format is obtained on the curved boundary of the airfoil. This method is simple to implement and does not require much special treatment of the curved elements. According to the results obtained in the literature (L. Krivodonova, M. Berger, High-order accurate implementation of solid wall boundary conditions in curved geometries, Journal of ComputationalPhysics 211 (2006) 492–512. https: / / doi.org / 10.1016 / j.jcp.2005.05.029.), using this method in the DG format, high-order accuracy results are obtained for inviscid problems.
[0078] When using the straight-edge element discretization method to calculate the numerical flux at the Gaussian integration point on the curved edge, the velocity component still needs to be calculated using the outer normal of the curved boundary at the fitting location. However, based on the high-order curved boundary expression of the airfoil, the straight-edge element discretization method cannot clearly select the outer normal of the object surface of the curved boundary of the airfoil. Krivodonova did not reconstruct the high-order expression of the curved object surface. Instead, he constructed an approximate circular arc based entirely on the mesh information of the straight edge using two element vertices located on the curved boundary. He then calculated the average value of the arc radius based on the adjacent vertices on the left or right side, and thus obtained the outer normal of the Gaussian integration point on the straight boundary at the projection point on the approximate circular arc curved boundary.
[0079] In his work, Li applied the curved unit processing concept of Krivodonova and Berger to the DG / FV format. However, Li differed from this approach by using cubic spline interpolation to construct a curved boundary representation, thereby overcoming the problem of directly approximating the true curved boundary with circular arcs. He also achieved high-order precision results for the curved boundary problem. However, analysis revealed three key issues with this concept. First, Li's work did not address how to determine the calculated position of the normal vector outside the curved boundary based on the reconstructed standard cubic spline interpolation curve. Second, using the standard cubic spline interpolation curve is not an ideal choice. Even after improving it using a piecewise approach, the reconstructed curve still exhibits significant geometric errors compared to the improved polynomial curve. Finally, this curved boundary processing scheme has only been verified in the high-order DG and DG / FV formats. Whether it remains effective within the framework of the high-order FV format requires further verification.
[0080] The following introduces the straight-edge element discretization method based on the normal direction of the airfoil curved surface.
[0081] In response to the three key issues mentioned above, the example of the present invention is based on the airfoil boundary represented by the high-order curved boundary expression method, improves the straight-edge unit discretization method directly used on the airfoil curved unit from three aspects, and proposes a straight-edge unit discretization method based on the normal of the airfoil curved surface.
[0082] With respect to the curved boundary unit external normal, Gaussian integration point and arc length required for solving the airfoil curved unit integral in formula (5), formula (3) and formula (7), the straight edge unit discretization method based on the curved surface normal of the airfoil proposed in the embodiment of the present invention has made three optimizations compared with the direct straight edge unit discretization of the curved boundary of the airfoil: first, an airfoil curved boundary arc length calculation algorithm is proposed, which no longer directly calculates the length of the straight edge, but accurately obtains the arc length of the curved boundary of the airfoil; second, a new Gaussian integration point calculation algorithm is proposed, and the Gaussian integration point is no longer selected on the straight edge, but on the curved boundary of the airfoil; third, a method for calculating the unit external normal on the curved surface is proposed, which no longer calculates the surface normal on the straight edge, but obtains the exact surface normal of the Gaussian integration point on the curved boundary.
[0083] (3) Calculation method of arc length of airfoil curved boundary
[0084] The position of the Gaussian integration point is determined by the arc length. This is easier for straight-edge elements. Assuming that under two-dimensional conditions, the two endpoints of the surface to be integrated are , If the number of integration points is 2, two integration points and The location is:
[0085] (8)
[0086] If the number of integration points is 3, we have:
[0087] (9)
[0088] However, it is not easy to determine the location of the integration point according to the arc length of the unit curve on the curved boundary. For example, in the interval The arc length L on can be obtained according to the following analytical formula:
[0089] (10)
[0090] However, when the curve expression is relatively complex, it is impossible to directly obtain the exact integral of the arc length, and a numerical method is needed. The most common calculation method is the trapezoidal method, which is to use the curve Divide into n subintervals, the width of each subinterval is:
[0091] (11)
[0092] At this point, according to the trapezoidal method, we have:
[0093] (12)
[0094] However, the integral degree of the trapezoidal method is only 1, which will introduce low-order truncation errors when the grid points on the curved boundary are sparsely distributed.
[0095] The present invention proposes a method for calculating the curved side length of an airfoil curved element based on the Simpson and Simpson 3 / 8 rules. The Simpson and Simpson 3 / 8 rules offer higher integration degrees than the trapezoidal method, at 3 and 4, respectively, resulting in higher accuracy. This makes them more compatible with third- and fourth-order finite volume schemes.
[0096] Simpson's rule approximates the integral region by a quadratic polynomial, but requires The interval is divided into an even number of uniform subintervals:
[0097] (13)
[0098] At this time, the arc length L can be approximated as:
[0099] (14)
[0100] On this basis, the more accurate Simpson 3 / 8 rule requires dividing the interval into multiples of 3:
[0101] (15)
[0102] The arc length L is:
[0103] (16)
[0104] The following introduces the calculation method of Gaussian integration points on the curved boundary of the airfoil.
[0105] After calculating the arc length of the curved side of the airfoil curved unit using the Simpson and Simpson 3 / 8 rules, the Gaussian integral point on the surface can be found based on the arc length, that is, the calculation position of the outer normal component on the surface can be determined.
[0106] Based on the Simpson and Simpson 3 / 8 rules, the arc length of the curved edge of the airfoil curved element can be calculated using the Newton iteration method to find the integration points. Taking two Gaussian integration points as an example, assuming that the starting and end points of the curved edge are , the position of the integration point to be found is and The target arc length corresponding to the two integration point positions is:
[0107] (17)
[0108] Let's take two parameters:
[0109] (18)
[0110] in and The arc lengths from the starting point to the current position can be calculated using the Simpson and Simpson 3 / 8 rules. The error tolerance set in the embodiment of the present invention is ,when and When it is greater than the error tolerance, Newton iteration is performed:
[0111] (19)
[0112] Until two integration point positions that meet the error tolerance are found.
[0113] The following introduces the calculation method of the unit external normal of the airfoil curved boundary.
[0114] On the basis of accurately obtaining the curved boundary integral point, the unit external normal vector at the corresponding position can be directly obtained according to the reconstructed interpolation polynomial. Figure 2 The relative positions of the integration points are shown. The tangent vectors at the two integration points on the surface are,
[0115] (20)
[0116] The unit normals pointing outward are
[0117] (twenty one)
[0118] The above three methods form the core of the high-order discretization method for straight-edge elements based on the direction of the curved object surface. From a procedural perspective, this method is relatively simple to implement. Although it requires the use of the Simpson and Simpson 3 / 8 rules with high integration degrees to determine the position of the calculated external normal vector on the curved object surface, it avoids the construction of high-order isoparametric elements and does not require separate storage of Jacobians.
[0119] The flow chart of the high-order discretization method of straight-edge elements based on the direction of the airfoil curved surface is as follows: Figure 4 shown.
[0120] 1. First, read the high-order precision airfoil curve boundary information;
[0121] 2. According to the requirements of the high-order format, determine the integration degree of the airfoil curved unit, and change the integration degree from 3 to 3 and from 4 to 4;
[0122] 3. Use Simpson's rule to calculate arc length, turn 5;
[0123] 4. Use Simpson's 3 / 8 rule to calculate arc length, turn 5;
[0124] 5. Use Newton iteration method to solve the Gaussian integral points of the airfoil curved boundary;
[0125] 6. Calculate the normal direction of the Gaussian integration point at the curved boundary of the airfoil;
[0126] 7. Substitute the normal of the Gaussian integral point into the straight boundary discretization method formula to calculate the numerical integral;
[0127] 8. End.
[0128] The high-order curved boundary airfoil expression method is used to test the calculation results of the curved unit discretization method based on the quintic polynomial curve.
[0129] The following is an example of flow around an inviscid NACA0012 airfoil.
[0130] Basic information of the example: This example is the flow around NACA0012 airfoil, and the incoming flow Mach number is , using four sets of quadrilateral grids from sparse to dense, with the number of elements being 8000, 9400, 12720 and 14300 respectively. The four grids are named Coa, Med, Fin, Vfin Figure 5 The figure first shows the distribution of the first two sets of grids and the last set of grids near the airfoil surface.
[0131] Accuracy test and correctness test:
[0132] Figure 6 The figure shows the distribution of airfoil surface pressure coefficients obtained using different methods in third-order and fourth-order accuracy formats. The third-order format - cubic isoparametric element and the third-order format - quintic isoparametric element are respectively third-order cubic and quintic isoparametric element discretization methods, the fourth-order format - cubic isoparametric element and the fourth-order format - quintic isoparametric element are respectively fourth-order cubic and quintic isoparametric element discretization methods according to the embodiments of the present invention, and the third-order format - simplified curved edge element and the fourth-order format - simplified curved edge element are respectively third-order and fourth-order high-order discretization methods based on the straight edge element of the curved object surface direction of the airfoil according to the embodiments of the present invention.
[0133] It can be seen that Figure 6 The negative pressure coefficient is shown in , because when analyzing lift, we usually want to see the negative pressure contribution on the upper surface of the airfoil. When , the negative pressure contribution of the upper surface is positive. The positive pressure of the lower surface is negative, which can show the contribution distribution of lift more clearly.
[0134] Judging from the results obtained, under the three angles of attack conditions, the results obtained by different curved unit discretization methods in the third-order and fourth-order formats are very close. As the angle of attack continues to increase, the negative pressure provided by the upper surface of the airfoil and the positive pressure provided by the lower surface gradually increase, providing higher lift for the airfoil.
[0135] Figure 7 The figure shows that on the sparsest and densest grids, the angle of attack is The entropy error distribution of the wall element.
[0136] It can be seen that the entropy error distribution of the wall element obtained under different curved element discretization schemes is almost the same. As the format accuracy is improved and the grid is encrypted, the entropy error will further decrease. Table 1 further shows At the same angle of attack, different methods on all wall elements Error statistics.
[0137] Table 1 Statistics of the L2 errors of different methods on all wall elements and the calculation accuracy between the last two sets of dense grids
[0138]
[0139] The correctness of the improved polynomial curved boundary expression combined with the straight-edge element discretization method based on the curved boundary normal is demonstrated by the flow around the inviscid NACA0012 airfoil at three different attack angles.
[0140] The following is an example of flow around a viscous NACA0012 airfoil.
[0141] Basic situation of the example: This part is a viscous NACA0012 flow problem. The given incoming flow Mach number is , the Reynolds number is .like Figure 8 As shown in the figure, this part uses three sets of units, which are successively encrypted near the wing boundary. The number of grids is 14200, 15620 and 17040 respectively. We name them as the sparsest, medium and densest respectively. The height of the first layer grid is 5×10 -4 .
[0142] Because the first layer of mesh height is relatively small in high-Reynolds-number viscous problems, before presenting the specific calculation results of the different methods, we first examine the relative positional relationship between the wall curve reconstructed based on the quintic polynomial and the mesh near the wall in some areas of the airfoil on the three sets of meshes.
[0143] It can be seen that the improved quintic polynomial curve can accurately express the airfoil wall. When the first-layer grid height is small, the polynomial curve does not pass through the first-layer grid.
[0144] In the post-processing of numerical experiments, even if the high-order discretization method of straight-edge elements based on the direction of the curved surface of the airfoil can ensure the flow results of corresponding accuracy, the calculation of the pressure and drag coefficients involves the actual integration point positions and the integration along the boundary surface. At this time, the integration points and the numerical integration process are both along the straight line boundary, so the calculation results may be significantly different from those of the high-order discretization method of straight-edge elements based on the direction of the curved surface of the airfoil.
[0145] Therefore, the post-processing process of the high-order discretization method of the straight-edge unit based on the direction of the airfoil curved surface is improved, and the process of finding the calculation position of the normal vector outside the curved boundary is based on the Simpson 3 / 8 rule to ensure that the calculation position of the pressure coefficient is located at C G2 The integration of the drag coefficient is performed along the curved boundary, and the arc length of the curved boundary calculated based on the Simpson 3 / 8 rule is used. The calculation results obtained based on this method are named simplified curved edge element (improved).
[0146] Figure 9 The pressure coefficient C is shown in D The curve of change with the number of iteration steps, where the cubic isoparametric unit and the quintic isoparametric unit are respectively the cubic and quintic optimized isoparametric unit curved unit discretization methods of the embodiment of the present invention, the simplified curved edge unit is the high-order discretization method of the straight edge unit based on the direction of the curved object surface of the airfoil in the embodiment of the present invention, the simplified curved edge unit (improved) is the result of improving the post-processing of the high-order discretization method of the straight edge unit based on the direction of the curved object surface of the airfoil in the embodiment of the present invention, and the straight edge unit refers to the straight edge unit discretization method.
[0147] Figure 10The pressure coefficient distribution of different methods on three sets of grids is first demonstrated. Among them, the cubic isoparametric unit and the quintic isoparametric unit are the cubic and quintic optimized isoparametric unit curved unit discretization methods respectively. The simplified curved edge unit is a high-order discretization method of the straight edge unit based on the curved surface direction of the airfoil in an embodiment of the present invention. The simplified curved edge unit (improved) is the result of improving the post-processing of the high-order discretization method of the straight edge unit based on the curved surface direction of the airfoil in an embodiment of the present invention. The straight edge unit refers to the straight edge unit discretization method.
[0148] Judging from the results obtained, selecting the calculation point of the pressure coefficient on the curved boundary at this time does not bring obvious improvement. The results are still close to the discretization method based entirely on straight-edge elements, and there is a clear deviation between the discretization method and the high-order discretization method of straight-edge elements based on the direction of the curved surface of the airfoil.
[0149] Unlike the pressure coefficient, the drag coefficient integral obtained by the improved post-processing step using a high-order discretization method for straight-edge elements in the direction of the airfoil's curved surface is significantly improved, and is closer to the results of the optimized isoparametric curved element discretization method, especially on denser grids. We analyzed that the possible cause of this phenomenon is that the pressure coefficient is only the result at each x position. Therefore, for high Reynolds number airfoil flows, the shape and area difference between straight-edge elements and curved-edge elements is very small due to the higher aspect ratio of the wall grid itself, and the difference in results before and after the improved calculation position cannot be directly reflected. The drag coefficient output is the integral result along the entire airfoil surface, which leads to cumulative errors.
[0150] Through the problem of flow around a viscous NACA0012 airfoil, it is further explained that performing the integration process along the straight edge and the curved edge will lead to large deviations, and that the use of a high-order discretization method of straight-edge elements based on the direction of the curved surface of the airfoil can achieve high accuracy in airfoil design and is more economical than the isoparametric element method. Example
[0151] Embodiment 2 of the present invention provides a terminal device corresponding to the above-mentioned embodiment 1. The terminal device can be a processing device for a client, such as a mobile phone, a laptop computer, a tablet computer, a desktop computer, etc., to execute the method of the above-mentioned embodiment.
[0152] The terminal device of this embodiment includes a memory, a processor, and a computer program stored in the memory; the processor executes the computer program in the memory to implement the steps of the method in the above-mentioned embodiment 1.
[0153] In some implementations, the memory may be a high-speed random access memory (RAM), and may also include a non-volatile memory, such as at least one disk storage.
[0154] In other implementations, the processor may be a central processing unit (CPU), a digital signal processor (DSP), or other general-purpose processors, which are not limited herein. Example
[0155] Embodiment 3 of the present invention provides a computer-readable storage medium corresponding to the above-mentioned embodiment 1, on which a computer program / instruction is stored. When the computer program / instruction is executed by a processor, the steps of the method of the above-mentioned embodiment 1 are implemented.
[0156] Computer readable storage media can be tangible devices that hold and store instructions used by instruction execution devices. Computer readable storage media can be, for example, but not limited to, electronic storage devices, magnetic storage devices, optical storage devices, electromagnetic storage devices, semiconductor storage devices, or any combination thereof.
[0157] Those skilled in the art will appreciate that the embodiments of the present application may be provided as methods, systems, or computer program products. Therefore, the present application may take the form of a fully hardware embodiment, a fully software embodiment, or an embodiment combining software and hardware. Furthermore, the present application may take the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code. The solutions in the embodiments of the present application may be implemented in various computer languages, such as the object-oriented programming language Java and the interpreted scripting language JavaScript.
[0158] The present application is described with reference to the flowcharts and / or block diagrams of the methods, devices (systems), and computer program products according to the embodiments of the present application. It should be understood that each process and / or block in the flowchart and / or block diagram, as well as the combination of processes and / or blocks in the flowchart and / or block diagram, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing device to produce a machine, so that the instructions executed by the processor of the computer or other programmable data processing device generate instructions for implementing the processes in the flowchart and / or block diagram. Figure 1 a process or multiple processes and / or boxes Figure 1 A device that provides the functions specified in a block or multiple blocks.
[0159] These computer program instructions can also be loaded onto a computer or other programmable data processing device so that a series of operational steps are executed on the computer or other programmable device to produce a computer-implemented process, thereby providing the instructions executed on the computer or other programmable device for implementing the process. Figure 1 a process or multiple processes and / or boxes Figure 1A step that specifies a function in one or more boxes.
[0160] Although the preferred embodiments of the present application have been described, those skilled in the art may make additional changes and modifications to these embodiments once they have learned the basic creative concept. Therefore, the appended claims are intended to be interpreted as including the preferred embodiments and all changes and modifications that fall within the scope of the present application.
[0161] Obviously, those skilled in the art may make various changes and modifications to this application without departing from the spirit and scope of this application. Thus, if these modifications and variations of this application fall within the scope of the claims of this application and their equivalents, this application is intended to include these modifications and variations.
Claims
1. A wing boundary simulation method based on the normal direction of the airfoil curved surface, characterized in that: The following steps are involved: S1. Obtain wing curve boundary information; S2. Determine the integral degree of the curved unit according to the wing curved boundary information. If the integral degree is 3, use Simpson's rule to calculate the arc length L; if the integral degree is 4, use Simpson's 3 / 8 rule to calculate the arc length L; S3. Using the arc length L, calculate the Gaussian integral point of the wing curvature boundary using the Newton iteration method; S4, calculate the normal direction of the Gaussian integration point at the wing curvature boundary; S5. Substitute the normal of the Gaussian integral point into the straight boundary discretization method formula to calculate the numerical integral: ; in, It represents the average value of wing speed or temperature in the unit cell. represents the wing unit volume, represents the number of unit faces, represents the number of unit area points, is the surface integral coefficient, represents the convection and viscous flux at the area integral point k, represents the number of unit volume points, represents the volume fraction coefficient, represents the heat source at the volume point q; In step S3, the specific implementation process of using the Newton iteration method to solve the Gaussian integral point of the airfoil curved boundary includes: 1) Set two parameters: ; , and are the positions of the two integration points, 、 is the target arc length corresponding to the two integration point positions, and are the arc lengths from the starting point to the two integration points respectively; 2) When and When the error is greater than the tolerance, the positions of the two integration points are updated using the following formula: and ; 、 are the first derivatives of the wing boundary curve at the two integration points respectively; 3) Repeat step 2) until and If the error is less than or equal to the error tolerance, the updated integration point position is output.
2. The wing boundary simulation method based on the normal direction of the airfoil curved surface according to claim 1, characterized in that: The wing boundary curve calculated using Simpson's rule is in the interval The expression of the arc length L on is: ; in, , n is The number of subintervals, 、 、 They are the airfoil curve function f(x) in 、 、 The first derivative at , 、 、 They represent the horizontal coordinates of point a, the i-th integration point, and point b respectively.
3. The wing boundary simulation method based on the normal direction of the airfoil curved surface according to claim 2, characterized in that: The wing boundary curve calculated using Simpson's 3 / 8 rule is in the interval The expression of the arc length L on is: ; in, , n is The number of subintervals.
4. The wing boundary simulation method based on the normal direction of the airfoil curved surface according to claim 1, characterized in that: In step S4, the normal calculation formula of the Gaussian integration point at the wing curve boundary is: 。 5. A terminal device comprising a memory, a processor, and a computer program stored in the memory; characterized in that: The processor executes the computer program to implement the steps of the method according to any one of claims 1 to 4.
6. A computer-readable storage medium having a computer program / instruction stored thereon; characterized in that: When the computer program / instructions are executed by a processor, the steps of the method according to any one of claims 1 to 4 are implemented.
7. A computer program product comprising a computer program / instructions; characterized in that When the computer program / instruction is executed by a processor, the steps of the method according to any one of claims 1 to 4 are implemented.
Citation Information
Patent Citations
NURBS curve bi-directional adaptive interpolation algorithm based on S-curve acceleration and deceleration algorithm
CN107817764A
A numerical simulation method for three-dimensional non-stick low-speed streaming based on curved surface boundary conditions
CN109726433A