A Numerical Simulation Method and Device for Contra-Rotating Open Rotor Considering Nacelle Coupling

Abstract

This application belongs to the field of aero-engine simulation technology, and particularly relates to a numerical simulation method and device for a contra-rotating open rotor considering nacelle coupling. The method includes: calculating the induced velocity vector generated by the nacelle at the lift line positions of two rows of propfan blades; calculating the total induced velocity of each propfan based on the induced velocity vector of each propfan, the self-induced velocity of each propfan, and the induced velocity between the two propfans, and then calculating the lift coefficient and drag coefficient based on the angle of attack; then determining the circulation distribution, and controlling the iterative cycle according to the convergence of the circulation distribution. During the cycle, the self-induced velocity of each propfan and the induced velocity between the two propfans are updated based on the circulation distribution. This application can calculate the characteristics of the contra-rotating propfan in the full operating condition range, with high calculation accuracy and fast calculation speed.

Description

Technical Field

[0001] This application belongs to the field of aeroengine simulation technology, and particularly relates to a numerical simulation method and device for a contra-rotating open rotor considering nacelle coupling. Background Art

[0002] As a new generation of high-efficiency and energy-saving green power, the open rotor engine has always been a research hotspot in recent years. During the research and development process of the open rotor engine, the overall performance design model is the basis for carrying out the overall engine matching research and component design. Therefore, in the preliminary design stage, establishing an open rotor engine overall performance design model with high accuracy, small computational amount, high integration, and stable numerical calculation helps to accelerate the design process and is an urgent need in the future aircraft power plant design field.

[0003] Currently, the overall simulation method of the open rotor engine is mainly based on the component-level zero-dimensional model. Although the zero-dimensional method based on the characteristic map can achieve the rapid calculation of the overall performance of the open rotor engine, this method has certain limitations: ① It is extremely dependent on the propfan characteristics, and there are few publicly published test characteristic maps, and the tested characteristic maps are also relatively sparse, making it difficult to ensure the accuracy of the model calculation. ② The semi-empirical model between the two rows of propfans cannot fully reflect the complex interference effects between the contra-rotating propfans, thus bringing certain errors. ③ For different propfan design parameters and different working conditions, it is usually necessary to scale and interpolate the characteristic map, thus introducing additional errors, especially at off-design point conditions, the errors are very large. ④ The zero-dimensional model cannot reflect the influence of design parameters such as blade geometry, number of blades, diameter, rotational speed, and pitch angle on the aerodynamic performance, and can only be used as a module for evaluating thrust and cannot be used for the design of contra-rotating propfans. In summary, the calculation accuracy of the contra-rotating propfan model based on the single-prop characteristic map is limited.

[0004] In order to improve the model accuracy, some scholars have used CFD technology to achieve three-dimensional calculations of contra-rotating propfans, but three-dimensional simulations require a large amount of computing resources and computing time, and it takes more than one hour to calculate a single point. Therefore, three-dimensional simulations are mainly used for mechanism analysis and are not suitable for the rapid design requirements in the preliminary design stage of open rotor engines. Summary of the Invention

[0005] In order to solve the above problems, this application provides a numerical simulation method and device for a contra-rotating open rotor considering nacelle coupling, taking into account both computational efficiency and computational accuracy.

[0006] The first aspect of this application provides a numerical simulation method for a contra-rotating open rotor considering nacelle coupling, mainly including:

[0007] Step S1: Divide each blade of the contra-rotating propfan into multiple blade elements along the span direction, construct a three-dimensional model of each blade according to the input airfoil parameters, and determine the lift line position coordinates of each propfan.

[0008] Step S2: Calculate the induced velocity vector generated by the nacelle on the blades of the two rows of contra-rotating fans at the lift line positions based on the air flow velocity at the inlet of the contra-rotating fans, the grid point coordinates of the nacelle geometry, and the position coordinates of all lift lines of the two rows of fans.

[0009] Step S3: Calculate the total induced velocity of each fan based on the air flow velocity at the inlet of the contra-rotating fans, the rotational speeds of the two rows of fans respectively, the induced velocity vectors of the nacelle on each fan, the self-induced velocity of each fan, and the induced velocity between the two fans, where the initial values of the self-induced velocity of each fan and the induced velocity between the two fans are set to 0.

[0010] Step S4: Perform coordinate transformation on the total induced velocity of the two fans to determine the resultant velocity, chordwise component velocity, and normal component velocity of each blade element.

[0011] Step S5: Determine the angle of attack distribution of each blade element of the two rows of contra-rotating fans based on the chordwise component velocity and normal component velocity of each blade element.

[0012] Step S6: Determine the lift coefficient and drag coefficient corresponding to the angle of attack of each blade element in the two-dimensional airfoil lift-drag characteristic database.

[0013] Step S7: Determine the circulation distribution of the blade elements based on the lift coefficient and the resultant velocity of the blade elements.

[0014] Step S8: Judge whether the circulation of each blade element converges. After the circulation of each blade element converges in the iterative calculation, calculate the performance of the contra-rotating fans according to the circulation distribution of the blade elements; otherwise, recalculate the new circulation distribution of the blade elements according to the relaxation factor.

[0015] Step S9: Determine the induced velocity generated by each vortex filament segment of each wake vortex of each blade at the control point of each blade element according to the new circulation of the blade elements, so as to calculate the induced velocity of each wake vortex of each fan in the axial direction of the fan.

[0016] Step S10: Determine the geometric coordinates of each vortex filament segment according to the induced velocity of the wake vortex filaments in the axial direction of the fan, the air flow velocity at the inlet of the contra-rotating fans, and the rotational speed of the rotating shaft.

[0017] Step S11: Calculate the self-induced geometric coefficient and mutual-induced geometric coefficient of each fan of the contra-rotating fans based on the lift line position coordinates of the two fans and the geometric coordinates of each vortex filament segment of the two fans.

[0018] Step S12: Determine the self-induced velocity of each fan and the induced velocity between the two fans based on the self-induced geometric coefficient, the mutual-induced geometric coefficient, and the new circulation distribution of the blade elements on the two rows of fans, and return to Step S3, and loop to execute Step S3 - Step S12 until the circulation of the blade elements converges in the iterative calculation.

[0019] Preferably, in step S3, the total induced velocity of each propeller fan is calculated by the following formula:

[0020] ;

[0021] Wherein, is the total induced velocity of the i-th element of the front row propeller fan, is the total induced velocity of the i-th element of the rear row propeller fan, is the air flow velocity at the inlet of the contra-rotating propeller fan, is the rotational speed of the front row propeller fan shaft, is the rotational speed of the rear row propeller fan shaft, is the distance from the i-th element of the front row propeller fan to the front row propeller fan shaft, is the distance from the i-th element of the rear row propeller fan to the rear row propeller fan shaft, is the induced velocity vector of the nacelle on the i-th element of the front row propeller fan, is the induced velocity vector of the nacelle on the i-th element of the rear row propeller fan, is the self-induced velocity of the front row propeller fan at its i-th element, is the induced velocity of the rear row propeller fan at the i-th element of the front row propeller fan, is the self-induced velocity of the rear row propeller fan at its i-th element, is the induced velocity of the front row propeller fan at the i-th element of the rear row propeller fan.

[0022] Preferably, in step S8, when the norm of the difference between the circulation of each element calculated in the current iteration and the circulation of each element calculated in the previous iteration is less than or equal to 10 -6 , it is determined that the circulation of the element calculated by the iterative calculation converges, otherwise the new circulation of the element is calculated by the following formula :

[0023] ;

[0024] Wherein, is the circulation of the element calculated in the previous iteration, is the circulation of the element calculated in the current iteration, is the relaxation factor.

[0025] Preferably, in step S8, the performance of the contra-rotating propeller fan is calculated by the following steps:

[0026] Step S81, calculate the lift and drag of each element according to the circulation of the element, the resultant velocity, the lift coefficient and the drag coefficient;

[0027] Step S82, calculate the component forces of each element in each direction according to the lift and drag of each element;

[0028] Step S83: Obtain the thrust and torque of the front row of propfans and the rear row of propfans through integration.

[0029] Preferably, in Step S10, the geometric coordinates of each vortex filament segment in the cylindrical coordinate system are determined by the following formula , , :

[0030] ;

[0031] wherein, is the radius position of the j-th segment of the k-th wake vortex filament, is the radius position of the i-th blade element, is the azimuth angle of the wake vortex filament, is the number of turns allowed for the development of the wake vortex filament, is the number of segments into which the wake vortex filament is divided, is the number of each segment into which the wake vortex filament is divided, is the rotational speed of the propfan shaft, is the induced velocity of the wake vortex filament in the axial direction of the propfan, is the parameter for controlling the grid level.

[0032] The second aspect of the present application provides a contra-rotating open rotor numerical simulation device considering nacelle coupling, mainly including:

[0033] Lift line determination module, configured to divide each blade of the contra-rotating propfan into multiple blade elements along the span direction, construct a three-dimensional model of each blade according to the input airfoil parameters, and determine the position coordinates of the lift line of each propfan;

[0034] Nacelle induced velocity vector calculation module, configured to calculate the induced velocity vector generated by the nacelle on the blades of the two rows of propfans at the lift line position according to the air flow velocity at the inlet of the contra-rotating propfan, the grid point coordinates of the nacelle geometry, and the position coordinates of all lift lines of the two rows of propfans;

[0035] Total induced velocity calculation module, configured to calculate the total induced velocity of each propfan according to the air flow velocity at the inlet of the contra-rotating propfan, the rotational speeds of the two rows of propfans respectively, the induced velocity vector of the nacelle on each propfan, the self-induced velocity of each propfan, and the induced velocity between the two propfans, wherein the initial values of the self-induced velocity of each propfan and the induced velocity between the two propfans are set to 0;

[0036] Blade element velocity determination module, configured to perform coordinate transformation on the total induced velocity of the two propfans to determine the resultant velocity, chordwise component velocity, and normal component velocity of each blade element;

[0037] The blade element angle of attack determination module is used to determine the angle of attack distribution of each blade element of the two rows of propfans based on the chordwise component velocity and the normal component velocity of each blade element;

[0038] The lift coefficient and drag coefficient determination module is used to determine the lift coefficient and drag coefficient corresponding to the angle of attack of each blade element in the two-dimensional airfoil lift-drag characteristic database;

[0039] The blade element circulation determination module is used to determine the blade element circulation distribution according to the lift coefficient and the resultant velocity of the blade element;

[0040] The loop control module is used to judge whether the circulation of each blade element converges. After the circulation of each blade element in the iterative calculation converges, the performance of the contra-rotating propfan is calculated according to the blade element circulation distribution. Otherwise, a new blade element circulation distribution is recalculated according to the relaxation factor;

[0041] The propfan axial direction induced velocity determination module is used to determine the induced velocity generated by each vortex filament segment of each wake vortex filament of each blade at the control point of each blade element according to the new blade element circulation, so as to calculate the induced velocity of each wake vortex filament of each propfan in the propfan axial direction;

[0042] The propfan wake vortex filament geometric coordinate determination module is used to determine the geometric coordinates of each vortex filament segment according to the induced velocity of the wake vortex filament in the propfan axial direction, the air flow velocity at the inlet of the contra-rotating propfan, and the rotational speed of the rotating shaft;

[0043] The geometric coefficient determination module is used to calculate the self-induced geometric coefficient and the mutual-induced geometric coefficient of each propfan of the contra-rotating propfan based on the lift line position coordinates of the two propfans and the geometric coordinates of each vortex filament segment of the two propfans;

[0044] The induced velocity update module is used to determine the self-induced velocity of each propfan and the induced velocity between the two propfans based on the self-induced geometric coefficient, the mutual-induced geometric coefficient, and the new blade element circulation distribution on the two rows of propfans.

[0045] Preferably, in the total induced velocity calculation module, the total induced velocity of each propfan is calculated by the following formula:

[0046] ;

[0047] Wherein, is the total induced velocity of the i-th blade element of the front row propfan, is the total induced velocity of the i-th blade element of the rear row propfan, is the air flow velocity at the inlet of the contra-rotating propfan, is the rotational speed of the rotating shaft of the front row propfan, is the rotational speed of the rotating shaft of the rear row propfan, is the distance from the i-th blade element of the front row propfan to the rotating shaft of the front row propfan, is the distance from the i-th blade element of the rear row propeller fan to the axis of the rear row propeller fan, is the induced velocity vector of the nacelle on the i-th blade element of the front row propeller fan, is the induced velocity vector of the nacelle on the i-th blade element of the rear row propeller fan, is the self-induced velocity of the front row propeller fan at its i-th blade element, is the induced velocity of the rear row propeller fan on the i-th blade element of the front row propeller fan, is the self-induced velocity of the rear row propeller fan at its i-th blade element, is the induced velocity of the front row propeller fan on the i-th blade element of the rear row propeller fan.

[0048] Preferably, in the circulation control module, the norm of the difference between the circulation of each blade element calculated in the current iteration and the circulation of each blade element calculated in the previous iteration is less than or equal to 10 -6 When it is determined that the circulation of the blade element in the iterative calculation converges, otherwise calculate the new circulation of the blade element through the following formula :

[0049] ;

[0050] wherein, is the circulation of the blade element calculated in the previous iteration, is the circulation of the blade element calculated in the current iteration, is the relaxation factor.

[0051] Preferably, the circulation control module includes:

[0052] A lift and drag calculation unit for calculating the lift and drag of each blade element according to the blade element circulation, the resultant velocity, the lift coefficient and the drag coefficient;

[0053] A component force calculation unit for calculating the component forces of each blade element in each direction according to the lift and drag of each blade element;

[0054] A thrust and torque calculation unit for obtaining the thrust and torque of the front row propeller fan and the rear row propeller fan through integration.

[0055] Preferably, in the blade element wake vortex filament geometric coordinate determination module, the geometric coordinates of each vortex filament segment in the cylindrical coordinate system are determined through the following formula 、 、 :

[0056] ;

[0057] wherein, is the radius position of the j-th segment of the k-th wake vortex filament, is the radius position of the i-th blade element, is the azimuth angle of the wake vortex filament, is the number of loops allowing the wake vortex filament to develop, is the number of segments into which the wake vortex filament is divided, is the number of each segment into which the wake vortex filament is divided, is the rotational speed of the shaft of the propeller fan, is the induced velocity of the wake vortex filament in the axial direction of the propeller fan, is a parameter for controlling the grid level.

[0058] This application improves the simulation accuracy of the contra-rotating propeller fan, increases the calculation speed, and can quickly obtain the characteristics of the contra-rotating propeller fan in the full operating condition range. BRIEF DESCRIPTION OF THE DRAWINGS

[0059] Figure 1 is a flowchart of a preferred embodiment of the numerical simulation method for a contra-rotating open rotor considering nacelle coupling in this application.

[0060] Figure 2 This application Figure 1 is a schematic diagram of the blade lift line positions of the front and rear row propeller fans of the embodiment shown in this application.

[0061] Figure 3 This application Figure 1 is a schematic diagram of the wake vortex filament structure and coordinate positions of the contra-rotating propeller fan of the embodiment shown in this application.

[0062] Figure 4 This application Figure 1 is a schematic diagram of the velocity triangle of the contra-rotating propeller fan of the embodiment shown in this application. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0063] To make the objectives, technical solutions, and advantages of the implementation of this application clearer, the technical solutions in the embodiments of this application will be described in more detail below with reference to the accompanying drawings in the embodiments of this application. In the drawings, the same or similar reference numerals denote the same or similar elements or elements with the same or similar functions from beginning to end. The described embodiments are some, but not all, of the embodiments of this application. The embodiments described below by referring to the drawings are exemplary and are intended to explain this application and should not be construed as limiting this application. All other embodiments obtained by those of ordinary skill in the art based on the embodiments in this application without creative efforts shall fall within the scope of protection of this application. The embodiments of this application will be described in detail below with reference to the drawings.

[0064] In the first aspect of this application, a numerical simulation method for a contra-rotating open rotor considering nacelle coupling is provided. Referring to Figure 1 , it mainly includes:

[0065] Step S1: Divide each blade of the contra-rotating propeller fan into multiple blade elements along the spanwise direction, construct a three-dimensional model of each blade according to the input airfoil parameters, and determine the position coordinates of the lift line of each propeller fan.

[0066] In this application, a three-dimensional model of the contra-rotating propeller fan is first constructed in Step S1. The curve distribution of the airfoil parameters of the contra-rotating propeller fan along the radial direction is used as the input. These airfoil parameters include: the designed lift coefficient , the maximum relative thickness ratio TOC, the chord length c, the twist angle θ, and the distribution curves of the sweep angle Λ along the radial direction. Each propeller fan blade is divided into N blade elements. The airfoil of each blade element is determined according to the above airfoil parameters, and the airfoil coordinates of the calculated blade element are obtained. First, stack the airfoil coordinates of each blade element by the centroid, then translate the centroid position according to the sweep angle Λ, then rotate around its centroid according to the twist angle θ, and finally install according to the given installation angle β to obtain the final propeller fan blade coordinates. Take the 1 / 4 chord length position of each blade element as the lift line position of the blade, as Figure 2 shown. Figure 2 The lift lines of the front and rear rows of propeller fans shown are displayed in the global xyz coordinate system. In this application, the global xyz coordinate system, the local scn coordinate system, and the cylindrical coordinate system are realized, and the conversion relationships of each coordinate system are clarified. Among them, the global xyz coordinate system is for the entire contra-rotating propeller fan, xy is the rotation plane of the propeller fan, and z is the rotation axis. The local scn coordinate system is for each blade element, s is the spanwise direction of the blade / blade element, c is the chord line direction of the blade element, and n is the normal direction.

[0067] Given the axial spacing between the two rows of propeller fans, translate the blade coordinates of the front row of propeller fans along the z-axis to obtain the blade coordinates of the rear row of propeller fans in all coordinate systems when the phase angle φ = 0°. If the phase angle φ between the front and rear rows of propeller fans is not equal to 0°, then the blade coordinates of the rear row of propeller fans translated along the z-axis can be further rotated and then coordinate-transformed according to the phase angle φ, and the blade coordinates of the rear row of propeller fans in all coordinate systems can also be obtained quickly.

[0068] Step S2: Calculate the induced velocity vector generated by the nacelle on the blades of the two rows of propeller fans at the lift line positions according to the air flow velocity at the inlet of the contra-rotating propeller fan, the grid point coordinates of the nacelle geometry, and the position coordinates of all lift lines of the two rows of propeller fans.

[0069] The air flow velocity in this step can be quickly calculated according to the atmospheric conditions (altitude, flight Mach number, temperature, etc.) at the given design point. For the nacelle, the induced effect of the two rows of propeller fans on the nacelle is not considered. Input the air flow velocity, the grid point coordinates of the nacelle geometry, and the position coordinates of all lift lines of the two rows of propeller fans into the Panair code developed by Boeing based on the panel method, and the induced velocity vectors generated by the nacelle on the blades of the two rows of propeller fans at the lift line positions can be obtained respectively and 。

[0070] Step S3: Calculate the total induced velocity of each paddle fan according to the inlet air velocity of the contra-rotating paddle fan, the rotational speeds of the two rows of paddle fans respectively, the induced velocity vectors of the nacelle on each paddle fan, the self-induced velocity of each paddle fan, and the induced velocity between the two paddle fans. Among them, the initial values of the self-induced velocity of each paddle fan and the induced velocity between the two paddle fans are set to 0.

[0071] It can be understood that in step S3, in order to calculate the total induced velocity, the present application not only considers the induced velocity of the nacelle on the two rows of paddle fans, but also considers the self-induced velocity of each paddle fan and the induced velocity between the two paddle fans, thereby improving the calculation accuracy.

[0072] In some alternative embodiments, in step S3, referring to Figure 4 the inlet velocity triangles of the two rows of paddle fans shown, calculate the total induced velocity of each paddle fan through the following formula:

[0073] ;

[0074] Among them, is the total induced velocity of the i-th blade element of the front row of paddle fans, is the total induced velocity of the i-th blade element of the rear row of paddle fans, is the air velocity at the inlet of the contra-rotating paddle fan, is the rotational speed of the shaft of the front row of paddle fans, is the rotational speed of the shaft of the rear row of paddle fans, is the distance from the i-th blade element of the front row of paddle fans to the shaft of the front row of paddle fans, is the distance from the i-th blade element of the rear row of paddle fans to the shaft of the rear row of paddle fans, is the induced velocity vector of the nacelle on the i-th blade element of the front row of paddle fans, is the induced velocity vector of the nacelle on the i-th blade element of the rear row of paddle fans, is the self-induced velocity of the front row of paddle fans at its i-th blade element, is the induced velocity of the rear row of paddle fans at the i-th blade element of the front row of paddle fans, is the self-induced velocity of the rear row of paddle fans at its i-th blade element, is the induced velocity of the front row of paddle fans at the i-th blade element of the rear row of paddle fans.

[0075] Step S4: Perform coordinate transformation on the total induced velocities of the two paddle fans to determine the resultant velocity, chordwise component velocity, and normal component velocity of each blade element.

[0076] Since the velocities at each blade element in the global xyz coordinate system are calculated in step S3, based on the coordinate systems constructed in step S1 and the conversion relationships between the coordinate systems, it is easy to obtain the resultant velocity of the cn airfoil plane of each blade element in the local scn coordinate system in step S4. and the chordwise component velocity and the normal component velocity . Here, the subscript i has the same meaning as in step S3 and is the blade element number.

[0077] Step S5: Determine the angle of attack distribution of each blade element of the two rows of paddle fans based on the chordwise component velocity and the normal component velocity of each blade element.

[0078] In this step, according to the formula , calculate the angle of attack distribution of the airflows of the two rows of paddle fans respectively, that is, calculate the angle of attack of the i-th blade element of the front row of paddle fans and the angle of attack of the i-th blade element of the front row of paddle fans.

[0079] Step S6: Determine the lift coefficient and drag coefficient corresponding to the angle of attack of each blade element in the two-dimensional airfoil lift-drag characteristic database.

[0080] This step makes use of the two-dimensional airfoil lift-drag characteristic database. Referring to Figure 1 , from the design lift coefficient , the maximum relative thickness ratio TOC, the angle of attack calculated in step S5, the Reynolds number Re, and the Mach number Ma, through the two-dimensional airfoil lift-drag characteristic database using kriging surrogate, obtain the lift coefficient and the drag coefficient on the i-th blade element. Here, the lift coefficients and drag coefficients of the NACA16 series and NACA65 series airfoil databases are specifically used, and the surrogate model is:

[0081] .

[0082] where shape is the airfoil type.

[0083] Step S7: Determine the blade element circulation distribution according to the lift coefficient and the resultant velocity of the blade element.

[0084] In this step, according to the circulation equation , calculate the circulation of each blade element.

[0085] where is the circulation of the i-th blade element, is the chord length of the i-th blade element, is the resultant velocity of the cn airfoil plane of the i-th blade element, is the lift coefficient of the i-th blade element.

[0086] Step S8: Determine whether the circulation of each blade element converges (it is defaulted that it does not converge in the first iteration). After the circulation of each blade element in the iterative calculation converges, calculate the performance of the contra-rotating propeller fan according to the distribution of the blade element circulation; otherwise, recalculate the distribution of the new blade element circulation according to the relaxation factor.

[0087] This step determines whether to perform simulation iteration through the convergence coefficient. For example, in some alternative embodiments, in step S8, when the norm of the difference between the circulation of each blade element calculated in the current iteration number and the circulation of each blade element calculated in the previous iteration number is less than or equal to 10 -6 , it is determined that the circulation of the blade element in the iterative calculation converges; otherwise, calculate the new blade element circulation through the following formula :

[0088] ;

[0089] wherein, is the circulation of the blade element calculated in the previous iteration number, is the circulation of the blade element calculated in the current iteration number, is the relaxation factor. In the first iteration, the induced velocity is set to 0 in step S3, and the circulation is obtained through steps S3 to S7, and it is defaulted that the circulation does not converge, and it is used as the for the second iteration to calculate. In the second iteration, the induced velocity in step S3 needs to be calculated according to and the wake vortex filaments.

[0090] Step S9: Determine the induced velocity generated by each vortex filament segment of each wake vortex filament of each blade at the control point of each blade element according to the new blade element circulation, so as to calculate the induced velocity of each wake vortex filament of each propeller fan in the axial direction of the propeller fan.

[0091] In this step, first, according to the Biot - Savart law, the induced velocity generated by the j - th vortex filament segment of the k - th wake vortex filament of the m - th blade at the control point of the i - th blade element is:

[0092] ;

[0093] wherein, the regularization factor is calculated according to the Lamb - Oseen vortex core model. is the vector from the end point A of the vortex filament segment to the control point P, is the vector from the end point B of the vortex filament segment to the control point P, is the vector from the end point A to the end point B of the vortex filament segment, is the distance from the vortex filament AB to the control point P, is the Lamb - Oseen constant, with a value of 1.25463, is the radius of the vortex core.

[0094] After that, based on the induced velocity generated by the j-th segment of the k-th wake vortex filament of the m-th blade at the control point of the i-th blade element Calculate the induced velocity of each wake vortex filament of each propeller fan in the axial direction of the propeller fan .

[0095] Specifically, by arranging the calculation formula of the induced velocity , it can be written as , is the geometric coefficient calculated according to the coordinates of the j-th segment of the k-th wake vortex filament of the m-th blade and the coordinates of the control point of the i-th blade element. Integrating all blades and all vortex filaments, the circulation on the k-th wake vortex filament is , then the induced velocity of all wake vortex filaments on the i-th blade element is , where is the circulation on the blade element when i = k, is the geometric coefficient of the k-th wake vortex filament on the i-th blade element. Finally, perform coordinate transformation, according to Obtain the induced velocity of each wake vortex filament of each propeller fan in the axial direction of the propeller fan .

[0096] In this step, the axial direction of the propeller fan is also the z direction in the global xyz coordinate system.

[0097] Step S10: Determine the geometric coordinates of each vortex filament segment according to the induced velocity of the wake vortex filament in the axial direction of the propeller fan, the air flow velocity at the inlet of the contra-rotating propeller fan, and the rotational speed of the rotating shaft.

[0098] In some alternative embodiments, in step S10, the geometric coordinates of each vortex filament segment in the cylindrical coordinate system are determined by the following formula , , :

[0099] ;

[0100] Among them, is the radius position of the j-th segment of the k-th wake vortex filament, is the radius position of the i-th blade element. Since the wake vortex filament is generated from the boundary of the blade element and the radius of each vortex filament is constant, the radii of all vortex filament segments of the k-th wake vortex filament are the same, and the radius position of the k-th wake vortex filament is determined by the average value of the radius positions of its adjacent two blade elements. is the azimuth angle of the wake vortex filament, is the number of turns allowed for the wake vortex filament to develop, is the number of segments into which the wake vortex filament is divided, is the number assigned to each segment into which the wake vortex filament is divided, is the rotational speed of the shaft of the propfan, is the induced velocity of the wake vortex filament in the axial direction of the propfan, is the parameter for controlling the grid level.

[0101] In this step, according to the specified wake model, the coordinates of all the wake vortex filament segments generated by all the propfan blades after discretization of the front and rear rows of propfans are calculated respectively, as Figure 3 shown.

[0102] Step S11: Calculate the self-induced geometric coefficient and the mutual-induced geometric coefficient of each propfan of the contra-rotating propfan based on the lift line position coordinates of the two propfans and the geometric coordinates of each vortex filament segment of the two propfans.

[0103] In this step, the isolated wake model described above is used for wake superposition for both rows of propfans, and the wake vortex filament structure changes with the relative phase angle φ of the two rows of propfans to consider the periodic influence. The mutual influence of the wake models between the two rows of propfans is considered by adding the axial components of the self-induced velocity and the mutual-induced velocity of the two rows of propfans (i.e., the induced velocity mentioned above).

[0104] Specifically, the self-induced geometric coefficient generated by the k-th vortex filament of all the blades of the front propfan on the i-th blade element of the front propfan is calculated from the lift line coordinates of the front propfan and the coordinates of the wake vortex filaments of the front propfan , the self-induced geometric coefficient generated by the k-th vortex filament of all the blades of the rear propfan on the i-th blade element of the rear propfan is calculated from the lift line coordinates of the rear propfan and the coordinates of the wake vortex filaments of the rear propfan , the induced geometric coefficient of the rear propfan on the front propfan generated by the k-th vortex filament of all the blades of the rear propfan on the i-th blade element of the front propfan is calculated from the lift line coordinates of the front propfan and the wake coordinates of the rear propfan , and the induced geometric coefficient of the front propfan on the rear propfan generated by the k-th vortex filament of all the blades of the front propfan on the i-th blade element of the rear propfan is calculated from the lift line coordinates of the rear propfan and the wake coordinates of the front propfan .

[0105] Step S12: Determine the self-induced velocity of each propfan and the induced velocity between the two propfans based on the self-induced geometric coefficient, the mutual-induced geometric coefficient, and the new blade element circulation distribution on the two rows of propfans, return to step S3, and loop through steps S3 - S12 until the blade element circulation calculated by iteration converges.

[0106] In this step, the following formula is used to calculate each induced velocity:

[0107] ;

[0108] where, as mentioned above, here is the circulation on the k-th blade element of the front row propeller fan, is the circulation on the k-th blade element of the rear row propeller fan. As mentioned above, is the circulation of the i-th blade element. Since the circulation on the k-th wake vortex filament is equal to , that is, the k-th wake vortex filament is determined by the difference in circulation between the k - 1-th and k-th blade elements. The circulations on different wake vortex filaments are different. The induced velocity on the i-th blade element is . In special cases, when k = 1 or k = N + 1, . Therefore, after expanding and arranging this formula, it becomes , where is the circulation on the blade element when i = k, which is used to calculate the induced velocity generated by the k-th wake vortex filament.

[0109] After obtaining the above-mentioned induced velocities, return to step S3, and iterate to update the wake geometry coordinates, induced velocity distribution, total induced velocity distribution, angle of attack distribution, lift / drag coefficient, etc. in turn, and determine whether to terminate the iteration by whether the norm of the difference between the new blade element circulation and the previous step circulation is less than 10 -6 .

[0110] In some alternative embodiments, in step S8, the performance of the contra-rotating propeller fan is calculated through the following steps:

[0111] Step S81: Calculate the lift and drag of each blade element according to the blade element circulation, resultant velocity, lift coefficient, and drag coefficient;

[0112] Step S82: Calculate the component forces of each blade element in each direction according to the lift and drag of each blade element;

[0113] Step S83: Obtain the thrust and torque of the front row propeller fan and the rear row propeller fan through integration.

[0114] In this embodiment, according to the circulation distribution, velocity distribution, and lift coefficient distribution, calculate the lift distribution and drag distribution of the propeller fan blade, and then calculate the component forces in the local coordinate system:

[0115] ;

[0116] ;

[0117] In the above formula, is the angle of attack of the i-th blade element, is the lift of the i-th blade element, is the drag of the i-th blade element, is the resultant velocity of the i-th blade element on the cn airfoil plane, is the chord length of the i-th blade element, is the length of the lift line segment of the i-th blade element, is the air density. are the component forces in the local scn coordinate system respectively.

[0118] For different propfans, different angles of attack are adopted. As mentioned above, the angle of attack of the ith blade element of the front row of propfans and the angle of attack of the ith blade element of the front row of propfans ; the lift and the drag are similar.

[0119] Convert the component forces in the local coordinate system to the component forces in the cylindrical coordinate system , and , and finally perform integration. Calculate the thrust and torque of the front and rear row propfans respectively according to the following formulas, and finally obtain the total power of the contra-rotating propfans , the total thrust and the efficiency .

[0120] .

[0121] In the formula, T is the thrust, Q is the torque, is the number of blades, is the distance from the position of the ith blade element to the axis of the propfan.

[0122] It should be noted that in the spatial period of the contra-rotating propfans of this application, P relative positions of the propfans are divided, then P wake geometric structures of the contra-rotating propfans are obtained, corresponding to P relative phase angles. Change the relative phase angle φ in turn, rotate the blades and wake geometric coordinates of the rear row of propfans, repeat the steps of this application P times, the thrust and torque of the contra-rotating propfans at different phases can be obtained, and finally the average value is calculated after integrating the thrust and torque obtained at all phases, and the thrust and torque of the contra-rotating propfans considering the periodic influence can be obtained.

[0123] This application can calculate the induction effect of the nacelle on the two rows of propfans and obtain the nacelle installation performance of the contra-rotating propfans.

[0124] The lift line model of the contra-rotating propfans developed in this application can calculate the characteristics of the contra-rotating propfans in the full operating range, can reflect the influence of design parameters such as blade geometry, number of blades, diameter, rotational speed and pitch angle on the aerodynamic performance, can directly obtain the real characteristics, thus avoiding the error caused by the characteristic scaling characteristics, with high calculation accuracy, fast calculation speed, good convergence and stability. Compared with the wind tunnel test data of the F7-A7 contra-rotating propfans, the error is only about 0.5%, and it only takes about 20 s to calculate a working point with an ordinary computer.

[0125] The developed lifting-line model of the contra-rotating propeller fan in this application can be used as an independent characteristic diagram calculation program. By specifying the geometric shape of the propeller fan blade, the characteristic diagrams under the full operating conditions for different combinations of the pitch angles of the two rows of blades can be obtained, solving the problem of the scarcity of the characteristic diagrams of the contra-rotating propeller fan. It can also be used for the design of the contra-rotating propeller fan, parameter sensitivity analysis, and the study of the matching characteristics of the two rows of propeller fans, etc.

[0126] The second aspect of this application provides a numerical simulation device for a contra-rotating open rotor considering nacelle coupling corresponding to the above method, mainly including:

[0127] A lifting-line determination module, which is used to divide each blade of the contra-rotating propeller fan into multiple blade elements along the span direction, construct a three-dimensional model of each blade according to the input blade profile parameters, and determine the position coordinates of the lifting line of each propeller fan;

[0128] A nacelle induced velocity vector calculation module, which is used to calculate the induced velocity vector generated by the nacelle on the blades of the two rows of propeller fans at the lifting line position according to the air flow velocity at the inlet of the contra-rotating propeller fan, the grid point coordinates of the nacelle geometry, and the position coordinates of all the lifting lines of the two rows of propeller fans;

[0129] A total induced velocity calculation module, which is used to calculate the total induced velocity of each propeller fan according to the air flow velocity at the inlet of the contra-rotating propeller fan, the rotational speeds of the two rows of propeller fans respectively, the induced velocity vector of the nacelle on each propeller fan, the self-induced velocity of each propeller fan, and the induced velocity between the two propeller fans. Among them, the initial values of the self-induced velocity of each propeller fan and the induced velocity between the two propeller fans are set to 0;

[0130] A blade element velocity determination module, which is used to perform coordinate transformation on the total induced velocity of the two propeller fans to determine the resultant velocity, chordwise component velocity, and normal component velocity of each blade element;

[0131] A blade element angle of attack determination module, which is used to determine the angle of attack distribution of each blade element of the two rows of propeller fans based on the chordwise component velocity and normal component velocity of each blade element;

[0132] A lift coefficient and drag coefficient determination module, which is used to determine the lift coefficient and drag coefficient corresponding to the angle of attack of each blade element in the two-dimensional airfoil lift-drag characteristic database;

[0133] A blade element circulation determination module, which is used to determine the blade element circulation distribution according to the lift coefficient and the resultant velocity of the blade element;

[0134] A circulation control module, which is used to judge whether the circulation of each blade element converges. When the circulation of each blade element converges after iterative calculation, calculate the performance of the contra-rotating propeller fan according to the blade element circulation distribution. Otherwise, recalculate the new blade element circulation distribution according to the relaxation factor;

[0135] The axial induced velocity determination module of the contra-rotating propeller fan is used to calculate the induced velocity generated by each filament segment of each wake vortex filament of each blade at the control point of each blade element according to the new blade element circulation, so as to calculate the induced velocity of each wake vortex filament of each contra-rotating propeller fan in the axial direction of the contra-rotating propeller fan;

[0136] The geometric coordinate determination module of the wake vortex filament of the contra-rotating propeller fan is used to determine the geometric coordinates of each filament segment according to the induced velocity of the wake vortex filament in the axial direction of the contra-rotating propeller fan, the air flow velocity at the inlet of the contra-rotating propeller fan, and the rotational speed of the rotating shaft;

[0137] The geometric coefficient determination module is used to calculate the self-induced geometric coefficient and the mutual-induced geometric coefficient of each propeller fan of the contra-rotating propeller fan based on the lift line position coordinates of the two propeller fans and the geometric coordinates of each filament segment of the two propeller fans;

[0138] The induced velocity update module is used to determine the self-induced velocity of each propeller fan and the induced velocity between the two propeller fans based on the self-induced geometric coefficient, the mutual-induced geometric coefficient, and the new blade element circulation distribution on the two rows of propeller fans.

[0139] In some alternative embodiments, in the total induced velocity calculation module, the total induced velocity of each propeller fan is calculated by the following formula:

[0140] ;

[0141] wherein, is the total induced velocity of the i-th blade element of the front row propeller fan, is the total induced velocity of the i-th blade element of the rear row propeller fan, is the air flow velocity at the inlet of the contra-rotating propeller fan, is the rotational speed of the rotating shaft of the front row propeller fan, is the rotational speed of the rotating shaft of the rear row propeller fan, is the distance from the i-th blade element of the front row propeller fan to the rotating shaft of the front row propeller fan, is the distance from the i-th blade element of the rear row propeller fan to the rotating shaft of the rear row propeller fan, is the induced velocity vector of the nacelle on the i-th blade element of the front row propeller fan, is the induced velocity vector of the nacelle on the i-th blade element of the rear row propeller fan, is the self-induced velocity of the front row propeller fan at its i-th blade element, is the induced velocity of the rear row propeller fan at the i-th blade element of the front row propeller fan, is the self-induced velocity of the rear row propeller fan at its i-th blade element, is the induced velocity of the front row propeller fan at the i-th blade element of the rear row propeller fan.

[0142] In some alternative embodiments, in the loop control module, the norm of the difference between the circulation of each blade element calculated for the current iteration number and the circulation of each blade element calculated for the previous iteration number is less than or equal to 10 -6 When this occurs, it is determined that the circulation of the blade element calculated by the iterative calculation converges; otherwise, the new circulation of the blade element is calculated using the following formula :

[0143] ;

[0144] where is the circulation of the blade element calculated for the previous iteration number, is the circulation of the blade element calculated for the current iteration number, is the relaxation factor

[0145] In some alternative embodiments, the loop control module includes:

[0146] a lift and drag calculation unit configured to calculate the lift and drag of each blade element based on the circulation of the blade element, the resultant velocity, the lift coefficient, and the drag coefficient;

[0147] a component force calculation unit configured to calculate the component forces of each blade element in each direction based on the lift and drag of each blade element;

[0148] a thrust and torque calculation unit configured to obtain the thrust and torque of the front row of paddle fans and the rear row of paddle fans through integration

[0149] In some alternative embodiments, in the blade tip trailing vortex filament geometric coordinate determination module, the geometric coordinates of each vortex filament segment in the cylindrical coordinate system are determined using the following formula , , :

[0150] ;

[0151] where is the radius position of the j-th segment of the k-th trailing vortex filament, is the radius position of the i-th blade element, is the azimuth angle of the trailing vortex filament, is the number of turns allowed for the development of the trailing vortex filament, is the number of segments into which the trailing vortex filament is divided, is the number of each segment into which the trailing vortex filament is divided, is the rotational speed of the axis of the paddle fan, is the induced velocity of the trailing vortex filament in the axial direction of the paddle fan, is the parameter for controlling the grid level

[0152] As described above, it is only the specific implementation manner of the present application. However, the protection scope of the present application is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the technical scope disclosed in the present application should be covered within the protection scope of the present application. Therefore, the protection scope of the present application shall be subject to the protection scope of the claims described above.

Claims

1. A numerical simulation method for a contra-rotating open rotor considering nacelle coupling, characterized in that Comprising: Step S1: Divide each blade of the contra-rotating propeller fan into multiple blade elements along the spanwise direction, construct a three-dimensional model of each blade according to the input airfoil parameters, and determine the position coordinates of the lift lines of each propeller fan; Step S2: Calculate the induced velocity vector generated by the nacelle on the blades of the two rows of propeller fans at the lift line positions according to the air flow velocity at the inlet of the contra-rotating propeller fan, the grid point coordinates of the nacelle geometry, and the position coordinates of all lift lines of the two rows of propeller fans; Step S3: Calculate the total induced velocity of each propeller fan according to the air flow velocity at the inlet of the contra-rotating propeller fan, the rotational speeds of the two rows of propeller fans respectively, the induced velocity vector of the nacelle on each propeller fan, the self-induced velocity of each propeller fan, and the induced velocity between the two propeller fans, wherein the initial values of the self-induced velocity of each propeller fan and the induced velocity between the two propeller fans are set to 0; Step S4: Perform coordinate transformation on the total induced velocity of the two propeller fans to determine the resultant velocity, chordwise component velocity, and normal component velocity of each blade element; Step S5: Determine the angle of attack distribution of each blade element of the two rows of propeller fans based on the chordwise component velocity and normal component velocity of each blade element; Step S6: Determine the lift coefficient and drag coefficient corresponding to the angle of attack of each blade element in the two-dimensional airfoil lift-drag characteristic database; Step S7: Determine the circulation distribution of the blade elements according to the lift coefficient and the resultant velocity of the blade elements; Step S8: Judge whether the circulation of each blade element converges. After the circulation of each blade element in the iterative calculation converges, calculate the performance of the contra-rotating propeller fan according to the circulation distribution of the blade elements. Otherwise, recalculate the new circulation distribution of the blade elements according to the relaxation factor; Step S9: Determine the induced velocity generated by each vortex filament segment of each wake vortex of each blade at the control point of each blade element according to the new circulation of the blade elements, so as to calculate the induced velocity of each wake vortex of each propeller fan in the axial direction of the propeller fan; Step S10: Determine the geometric coordinates of each vortex filament segment according to the induced velocity of the wake vortex filament in the axial direction of the propeller fan, the air flow velocity at the inlet of the contra-rotating propeller fan, and the rotational speed of the rotating shaft; Step S11: Calculate the self-induced geometric coefficient and mutual-induced geometric coefficient of each propeller fan of the contra-rotating propeller fan based on the position coordinates of the lift lines of the two propeller fans and the geometric coordinates of each vortex filament segment of the two propeller fans; Step S12: Determine the self-induced velocity of each propeller fan and the induced velocity between the two propeller fans based on the self-induced geometric coefficient, mutual-induced geometric coefficient, and the new circulation distribution of the blade elements of the two rows of propeller fans, and return to Step S3, and loop to execute Step S3 - Step S12 until the circulation of the blade elements in the iterative calculation converges.

2. The numerical simulation method for a contra-rotating open rotor considering nacelle coupling as described in claim 1, characterized in that, In Step S3, the total induced velocity of each propeller fan is calculated by the following formula: ； Among them, is the total induced velocity of the ith blade element of the front row propeller fan, is the total induced velocity of the ith blade element of the rear row propeller fan, is the air flow velocity at the inlet of the contra-rotating propeller fan, is the rotational speed of the rotating shaft of the front row propeller fan, is the rotational speed of the rotating shaft of the rear row propeller fan, is the distance of the ith blade element of the front row propeller fan from the rotating shaft of the front row propeller fan, is the distance of the ith blade element of the rear row propeller fan from the rotating shaft of the rear row propeller fan, is the induced velocity vector of the nacelle on the ith blade element of the front row propeller fan, is the induced velocity vector of the nacelle on the ith blade element of the rear row propeller fan, is the self-induced velocity of the front row propeller fan at its ith blade element, is the induced velocity of the rear row propeller fan at the ith blade element of the front row propeller fan, is the self-induced velocity of the rear row propeller fan at its ith blade element, is the induced velocity of the front row propeller fan at the ith blade element of the rear row propeller fan.

3. The contra-rotating open rotor numerical simulation method considering nacelle coupling according to claim 1, characterized in that In step S8, when the norm of the difference between the circulation of each blade element calculated in the current iteration and the circulation of each blade element calculated in the previous iteration is less than or equal to 10 -6 , it is determined that the circulation of the blade element calculated by the iterative calculation converges; otherwise, the new circulation of the blade element is calculated by the following formula : ； Among them, is the circulation of the blade element calculated for the previous iteration number, is the circulation of the blade element calculated for the current iteration number, is the relaxation factor.

4. The numerical simulation method of a contra-rotating open rotor considering nacelle coupling according to claim 1, characterized in that, In Step S8, the performance of the contra-rotating propeller fan is calculated by the following steps: Step S81: Calculate the lift and drag of each blade element according to the circulation of the blade element, resultant velocity, lift coefficient, and drag coefficient; Step S82: Calculate the component forces of each blade element in each direction according to the lift and drag of each blade element; Step S83: Obtain the thrust and torque of the front row propeller fan and the rear row propeller fan through integration.

5. The numerical simulation method for a contra-rotating open rotor considering nacelle coupling as claimed in claim 1, characterized in that In step S10, the geometric coordinates of each vortex filament segment in the cylindrical coordinate system are determined by the following formula , , : ； Among them, is the radius position of the j-th segment of the k-th wake vortex filament, is the radius position of the i-th blade element, is the azimuth angle of the wake vortex filament, is the number of turns allowed for the development of the wake vortex filament, is the number of segments into which the wake vortex filament is divided, is the numbering of each segment into which the wake vortex filament is divided, is the rotational speed of the shaft of the propeller fan, is the induced velocity of the wake vortex filament in the axial direction of the propeller fan, is a parameter for controlling the grid level.

6. A contra-rotating open rotor numerical simulation device considering nacelle coupling, characterized in that Comprising: A lift line determination module, configured to divide each blade of the contra-rotating propeller fan into multiple blade elements along the spanwise direction, construct a three-dimensional model of each blade according to the input airfoil parameters, and determine the position coordinates of the lift lines of each propeller fan; The nacelle induced velocity vector calculation module is used to calculate the induced velocity vector generated by the nacelle on the blades of the two rows of contra-rotating propfans at the lift line positions according to the air flow velocity at the inlet of the contra-rotating propfans, the grid point coordinates of the nacelle geometry, and the position coordinates of all lift lines of the two rows of propfans; The total induced velocity calculation module is used to calculate the total induced velocity of each propfan according to the air flow velocity at the inlet of the contra-rotating propfans, the rotational speeds of the two rows of propfans respectively, the induced velocity vector of the nacelle on each propfan, the self-induced velocity of each propfan, and the induced velocity between the two propfans. Among them, the initial values of the self-induced velocity of each propfan and the induced velocity between the two propfans are set to 0; The blade element velocity determination module is used to perform coordinate transformation on the total induced velocity of the two propfans to determine the resultant velocity, chordwise component velocity, and normal component velocity of each blade element; The blade element angle of attack determination module is used to determine the angle of attack distribution of each blade element of the two rows of propfans based on the chordwise component velocity and normal component velocity of each blade element; The lift coefficient and drag coefficient determination module is used to determine the lift coefficient and drag coefficient corresponding to the angle of attack of each blade element in the two-dimensional airfoil lift-drag characteristic database; The blade element circulation determination module is used to determine the blade element circulation distribution according to the lift coefficient and the resultant velocity of the blade element; The circulation control module is used to judge whether the circulation of each blade element converges. When the circulation of each blade element in the iterative calculation converges, calculate the performance of the contra-rotating propfans according to the blade element circulation distribution, otherwise recalculate the new blade element circulation distribution according to the relaxation factor; The propfan axial direction induced velocity determination module is used to determine the induced velocity generated by each vortex filament segment of each wake vortex filament of each blade at the control point of each blade element according to the new blade element circulation, so as to calculate the induced velocity of each wake vortex filament of each propfan in the axial direction of the propfan; The propfan wake vortex filament geometric coordinate determination module is used to determine the geometric coordinates of each vortex filament segment according to the induced velocity of the wake vortex filament in the axial direction of the propfan, the air flow velocity at the inlet of the contra-rotating propfans, and the rotational speed of the rotating shaft; The geometric coefficient determination module is used to calculate the self-induced geometric coefficient and mutual-induced geometric coefficient of each propfan of the contra-rotating propfans based on the lift line position coordinates of the two propfans and the geometric coordinates of each vortex filament segment of the two propfans; The induced velocity update module is used to determine the self-induced velocity of each propfan and the induced velocity between the two propfans based on the self-induced geometric coefficient, mutual-induced geometric coefficient, and the new blade element circulation distribution on the two rows of propfans.

7. The contra-rotating open rotor numerical simulation device considering nacelle coupling according to claim 6, characterized in that In the total induced velocity calculation module, the total induced velocity of each propfan is calculated by the following formula: ； Among them, is the total induced velocity of the i-th element of the front row propeller fan, is the total induced velocity of the i-th element of the rear row propeller fan, is the air flow velocity at the inlet of the contra-rotating propeller fan, is the rotational speed of the rotating shaft of the front row propeller fan, is the rotational speed of the rotating shaft of the rear row propeller fan, is the distance from the i-th element of the front row propeller fan to the rotating shaft of the front row propeller fan, is the distance from the i-th element of the rear row propeller fan to the rotating shaft of the rear row propeller fan, is the induced velocity vector of the nacelle on the i-th element of the front row propeller fan, is the induced velocity vector of the nacelle on the i-th element of the rear row propeller fan, is the self-induced velocity of the front row propeller fan at its i-th element, is the induced velocity of the rear row propeller fan at the i-th element of the front row propeller fan, is the self-induced velocity of the rear row propeller fan at its i-th element, is the induced velocity of the front row propeller fan at the i-th element of the rear row propeller fan.

8. The counter-rotating open rotor numerical simulation device considering nacelle coupling according to claim 6, characterized in that, In the loop control module, the norm of the difference between the circulation of each blade element calculated for the current iteration number and the circulation of each blade element calculated for the previous iteration number is less than or equal to 10 -6 When it is determined that the circulation of the blade element calculated by the iterative calculation converges, otherwise calculate the new circulation of the blade element through the following formula : ； Among them, the circulation of the blade element calculated for the previous iteration number, the circulation of the blade element calculated for the current iteration number, is the relaxation factor.

9. The contra-rotating open rotor numerical simulation device considering nacelle coupling according to claim 6, characterized in that, The circulation control module includes: The lift and drag calculation unit is used to calculate the lift and drag of each blade element according to the blade element circulation, resultant velocity, lift coefficient, and drag coefficient; The component force calculation unit is used to calculate the component forces of each blade element in each direction according to the lift and drag of each blade element; The thrust and torque calculation unit is used to obtain the thrust and torque of the front row propfan and the rear row propfan through integration.

10. The contra-rotating open rotor numerical simulation device considering nacelle coupling according to claim 6, characterized in that, In the propeller fan wake vortex filament geometric coordinate determination module, the geometric coordinates of each vortex filament segment in the cylindrical coordinate system are determined by the following formula , , : ； Among them, is the radius position of the j-th segment of the k-th wake vortex filament, is the radius position of the i-th blade element, is the azimuth angle of the wake vortex filament, is the number of turns allowing the development of the wake vortex filament, is the number of segments into which the wake vortex filament is divided, is the numbering of each segment into which the wake vortex filament is divided, is the rotational speed of the shaft of the propeller fan, is the induced velocity of the wake vortex filament in the axial direction of the propeller fan, is the parameter for controlling the grid level.

Images