Interface flow field information generation method suitable for ducted fan momentum source method

Abstract

The invention particularly relates to an interface flow field information generation method suitable for a ducted fan momentum source method, which comprises the following steps of: S1, recording a front interface of a movable blade as a front interface, the front interface is divided into a front propeller hub boundary layer area, a front duct boundary layer area, a wing boundary layer area, a low-energy airflow area, a duct lip boundary layer area, a duct lip acceleration area and a front uniform transition area; respectively generating speed distribution according to the divided areas, and assembling and combining to obtain speed distribution on a front interface; s2, marking the rear interface of the stationary blade as a rear interface, and dividing the rear interface into a rear propeller hub boundary layer region, a duct index distribution transition region and a stationary blade wake region; and then generating speed distribution according to the divided areas respectively, and assembling and merging to obtain speed distribution on the rear interface. The method can get rid of dependence on flow field calculation of multiple reference coordinate systems, can also obtain accurate momentum source distribution characteristics, and improves the calculation efficiency.

Description

Technical Field

[0001] The present invention belongs to the field of design methods for distributed ducted fan aircraft, and particularly relates to a method for generating interface flow field information applicable to the ducted fan momentum source method. Background Art

[0002] A ducted fan is a power device driven by electricity and arranged inside a duct. Compared with an open-type electric-driven propeller, it has the advantages of a compact structural layout, low noise level, high safety, and easy internal embedding. In recent years, the electric-driven distributed ducted fan-wing coupled layout has gradually become a key development direction for various research institutions. Conducting aerodynamic design research on such aircraft has important engineering practical significance and application value.

[0003] The characteristics of an electric-driven distributed ducted fan aircraft are that multiple ducted fans are arranged on the wing, and there is a relatively close aerodynamic coupling effect between the fan, the wing, and the duct lip. Through research, the flow field information on the front and rear interfaces has a relatively obvious impact on the aerodynamic performance of each component. The main influencing factors are reflected in two aspects, namely the flow rate and the velocity distribution characteristics. Therefore, accurately predicting the aerodynamic performance of an aircraft is a challenging task. Using the traditional multiple reference coordinate system method has a high calculation cost, and when using the traditional momentum source method, due to differences in the flow rate and velocity distribution characteristics, the prediction error of the drag characteristics is relatively large.

[0004] In the prior art, Chinese Patent CN117371344A proposes a flow field simulation method for a distributed ducted fan-wing coupled layout. In this method, the flow field velocity information at the front and rear interfaces of the ducted fan is obtained by calculating the multiple reference coordinate system flow field for calculating the momentum source term. Although this method has certain advantages in terms of efficiency and accuracy, the acquisition of the flow field velocity information depends on the multiple reference coordinate system flow field, and the multiple reference frame calculation is actually based on the RANS method to simulate the flow field. For the aerodynamic target evaluation method in aircraft design, this method still has the defect of low calculation efficiency. Summary of the Invention

[0005] In order to solve the problems existing in the prior art, the present invention provides a method for generating interface flow field information applicable to the ducted fan momentum source method. The front and rear interfaces are processed separately. First, the average velocity at the fan under different oncoming flow velocities and different thrusts is calculated through momentum theory, then the velocity distributions at the front and rear interfaces considering the influence of multiple components are generated, and finally the obtained velocity distributions are assembled and merged to obtain the velocity distributions on the front and rear interfaces. This method can get rid of the dependence on the calculation of the multiple reference coordinate system flow field, and can also obtain accurate momentum source distribution characteristics, improving the calculation efficiency.

[0006] The present invention is realized through the following technical solutions:

[0007] An interface flow field information generation method applicable to the ducted fan momentum source method, comprising the following steps:

[0008] Step S1: Denote the interface in front of the ducted fan rotor blades as the front interface, and divide the front interface into a front hub boundary layer region, a front duct boundary layer region, a wing boundary layer region, a low-energy airflow region, a duct lip boundary layer region, a duct lip acceleration region, and a front uniform transition region; then generate velocity distributions for each of the above-mentioned regions respectively, and then assemble and merge them to obtain the velocity distribution on the front interface;

[0009] Step S2: Denote the interface behind the ducted fan stator blades as the rear interface, and divide the rear interface into a rear hub boundary layer region, a duct exponential distribution transition region, and a stator wake region; then generate velocity distributions for each of the above-mentioned regions respectively, and then assemble and merge them to obtain the velocity distribution on the rear interface.

[0010] Furthermore, for the method of obtaining the front and rear interfaces, refer to the method of obtaining the interface in front of the rotor blades and the interface behind the stator blades in Chinese Patent CN117371344A.

[0011] Furthermore, in the step S1, the front hub boundary layer region refers to the annular region on the front interface from the hub wall surface to the outer boundary of the hub boundary layer; the front duct boundary layer region refers to the annular region on the front interface from the inner duct wall surface to the outer boundary of the duct boundary layer; the wing boundary layer region refers to the rectangular region from the upper surface of the wing to the outer boundary of the wing boundary layer; the low-energy airflow region refers to the region from the outer boundary of the wing boundary layer to the position of the average flow velocity of the interface; the duct lip boundary layer region refers to the region from the duct lip wall surface to the outer boundary of the duct lip boundary layer; the duct lip acceleration region refers to the region obtained by taking the intersection of the circumferential action range region and the radial action range region. The circumferential action range region is: respectively with as the central angle, and the range of each central angle is used as the circumferential action range. The radial action range region is: taking the intersection of each central angle ray and the outer duct ring as the radial center position, and the range of 1 / 10 to 1 / 5 of the duct radius around each radial center position is used as the radial action range; the front uniform transition region is the remaining region on the front interface after removing the above-mentioned partitions.

[0012] Furthermore, the step S1 includes the following sub-steps:

[0013] Step S11: Calculate the average velocity at the front interface under different oncoming flow velocities and different thrust conditions based on the momentum theory;

[0014] The momentum theory is a well-known method in the art. The specific process of calculating the average velocity at the front interface is as follows:

[0015] First, according to the momentum theory, the fan thrust T P is expressed as:

[0016]

[0017] In the formula, T P represents the fan thrust, ρ is the atmospheric density, A 1 is the flow channel area at the front interface, V 0 is the free stream velocity, V 2 is the duct exit velocity. From the above formula, the velocity at the duct exit under different free stream velocities and different thrust conditions can be obtained as:

[0018]

[0019] The average velocity V 1 at the front interface is calculated by mass conservation:

[0020]

[0021] In the formula, A 2 is the duct exit cross-sectional area;

[0022] Step S12: Based on the thickness formula of the laminar boundary layer around a flat plate in the downstream direction and the assumption of the velocity power distribution, calculate the thickness and two-dimensional velocity distribution of the front hub boundary layer region; based on the thickness formula of the turbulent boundary layer around a flat plate in the downstream direction and the assumption of the velocity power distribution, calculate the thickness and two-dimensional velocity distribution of the front duct boundary layer region; based on the thickness and two-dimensional velocity distribution of the front hub boundary layer region and the thickness and two-dimensional velocity distribution of the front duct boundary layer region, interpolate to calculate the two-dimensional velocity distribution of the front uniform transition region;

[0023] The thickness formula of the laminar boundary layer around a flat plate in the downstream direction and the assumption of the velocity power distribution, as well as the thickness formula of the turbulent boundary layer around a flat plate in the downstream direction and the assumption of the velocity power distribution, are all well-known technologies in the art. The specific process of this step is as follows:

[0024] ① Use the calculation formula for the thickness of the laminar boundary layer of a flat plate placed in the downstream direction to calculate the thickness of the front hub boundary layer region. The specific expression is:

[0025]

[0026] In the formula, δ f_hub is the front hub boundary layer thickness, x hub is the chordwise length of the front cone of the hub, and Re x_f_hub is the Reynolds number calculated based on the chordwise length of the front cone of the hub;

[0027] Based on the thickness δ f_hub, the internal velocity distribution in the front hub boundary layer region is calculated using the power function distribution, i.e.:

[0028]

[0029] In the formula, u f_hub is the velocity at different normal positions inside the front hub boundary layer region, and z f_hub represents the coordinates of different normal positions inside the front hub boundary layer. V f_hub is the velocity at the outer edge of the front hub boundary layer and is calculated using the following formula:

[0030] V f_hub = ζV 1

[0031] In the formula, ζ is the correction coefficient of the velocity at the outer edge of the front hub boundary layer, which characterizes the acceleration effect of the hub surface on the airflow and takes values between 1.05 and 1.1;

[0032] ② The thickness of the front duct boundary layer region is calculated using the formula for the thickness of the turbulent boundary layer around a flat plate placed along the flow direction. The specific expression is:

[0033]

[0034] In the formula, δ f_duct is the thickness of the front duct boundary layer region, x f_duct is the chord length of the duct, and Re x_f_duct is the Reynolds number calculated based on the chord length of the duct;

[0035] Based on the thickness δ f_duct of the front duct boundary layer region, the internal velocity distribution in the front duct boundary layer region is calculated using the power function distribution, i.e.:

[0036]

[0037] In the formula, u f_duct is the velocity distribution at different normal positions inside the front duct boundary layer region, and z f_duct represents the coordinates of different normal positions inside the front duct boundary layer. V f_duct is the velocity at the outer edge of the front duct boundary layer; it is calculated using the following formula:

[0038] V f_duct = ζ 2 V 1

[0039] In the formula, ζ 2 is the correction coefficient of the velocity at the outer edge of the front duct boundary layer, which characterizes the acceleration effect of the duct surface on the airflow and takes values between 1.2 and 1.3;

[0040] ③Based on the calculated thickness and velocity distribution of the front hub boundary layer region and the front duct boundary layer region, the two-dimensional velocity distribution u of the front uniform transition region is obtained by using the linear interpolation method. gd ;

[0041] Specifically, the outer edge coordinate values of the hub boundary layer and the front duct boundary layer are obtained through the thickness of the front hub boundary layer region and the front duct boundary layer region. A linear interpolation function is established between the outer edges of the front hub boundary layer region and the front duct boundary layer region to interpolate the velocity, so that the velocity smoothly transitions between the outer edges of the front hub boundary layer region and the front duct boundary layer region, and the two-dimensional velocity distribution u of the front uniform transition region is obtained. gd ;

[0042] Step S13: Assemble the two-dimensional velocity distribution of the front hub boundary layer region, the two-dimensional velocity distribution of the front uniform transition region, and the two-dimensional velocity distribution of the front duct boundary layer region obtained in step S12 according to the spatial position to obtain the two-dimensional radial velocity distribution on the front interface. Rotate the obtained two-dimensional radial velocity distribution one week around the fan axis in the circumferential direction to generate the "base" of the three-dimensional velocity distribution on the front interface.

[0043] Step S14: Use the Karman-Pohlhausen method to calculate the thickness and two-dimensional velocity distribution of the wing boundary layer region, use xfoil to calculate the thickness of the low-energy air flow region, and based on the thickness and two-dimensional velocity distribution of the wing boundary layer region and the thickness of the low-energy air flow region, use the cubic spline distribution to calculate the two-dimensional velocity distribution of the low-energy air flow region; Assemble the two-dimensional velocity distribution of the wing boundary layer region and the two-dimensional velocity distribution of the low-energy air flow region and perform spanwise extension to obtain the three-dimensional velocity distribution of the wing boundary layer region and the low-energy air flow region, and perform "bowing" processing on the obtained three-dimensional velocity distribution.

[0044] The Karman-Pohlhausen approximation method is a well-known method in the art, and xfoil is a well-known calculation program in the art. The specific process of this step is as follows:

[0045] ①Use the classical Karman-Pohlhausen approximation method to calculate the thickness and velocity distribution of the wing boundary layer region. This method establishes the relationship between the wall shear stress and the velocity distribution through the definition of the momentum thickness, combines the no-slip condition in the wall normal direction and the far-field free-stream velocity condition, and obtains the wing boundary layer thickness δ wing_bl and velocity distribution u wing_bl , specifically refer to Section 3, Chapter 5 of "Viscous Fluid Mechanics" (written by Zhang Zixiong and Dong Zengnan); Multiply the wing boundary layer velocity distribution obtained by Karman-Pohlhausen by the ratio of the average velocity V 1 at the front interface to the free-stream velocity V 0 to obtain the boundary layer velocity distribution considering the fan suction effect.

[0046] ② First, calculate the momentum loss thickness of the low-energy airflow region based on XFOIL, and then calculate the nominal boundary layer thickness of the low-energy airflow region based on the momentum loss thickness, and subtract the thickness of the wing boundary layer to obtain the thickness of the low-energy airflow region. Specifically:

[0047] The momentum loss thickness δ of the low-energy airflow region 2 The calculation expression is

[0048]

[0049] In the formula, δ is the nominal boundary layer thickness of the low-energy airflow region, Λ is a dimensionless quantity, and its expression is:

[0050]

[0051] In the formula, υ is the kinematic viscosity coefficient, is the velocity gradient along the wing surface;

[0052] From the above two formulas, we get:

[0053]

[0054] Solve the above formula to obtain the nominal boundary layer thickness δ of the low-energy airflow region, and subtract the thickness δ of the wing boundary layer from the calculated δ wing_bl That is, the thickness δ of the low-energy airflow region is obtained wing_lp ;

[0055] ③ Based on the obtained thickness of the wing boundary layer region, obtain the coordinates of the outer edge of the wing boundary layer. Based on the thickness δ of the low-energy airflow region wing_lp , obtain the coordinates of the outer edge of the low-energy airflow region, and make the outer edge velocity value of the low-energy airflow region equal to the average velocity V of the front interface calculated in step S11 1 , and then according to the coordinates of the outer edge of the wing boundary layer, the coordinates of the outer edge of the low-energy airflow region, and adding the tangency constraints of the velocity distribution of the low-energy airflow region with the velocity distribution of the wing boundary layer and the velocity distribution of the front uniform transition region, use cubic spline distribution to calculate and obtain the velocity distribution of the low-energy airflow region;

[0056] ④ Combine the obtained velocity distribution of the low-energy airflow region with the velocity distribution of the wing boundary layer; obtain the two-dimensional velocity distribution of the wing boundary layer and the low-energy airflow region, and extend the two-dimensional velocity distribution along the span direction to the outer edge of the duct to generate a three-dimensional rectangular velocity distribution;

[0057] ⑤ Perform "bowing" processing on the obtained three-dimensional velocity distribution. The specific process is as follows:

[0058] The calculation method of the z coordinate in the height direction of the bowing is:

[0059]

[0060] Among them, the coordinate system is defined as follows: the X-axis is along the free-stream direction, the Y-axis is to the right from the pilot's perspective, and the Z-direction is perpendicular to the XY plane and upward, satisfying the right-hand system; the calculation object in the formula is the velocity distribution scatter points at different spanwise positions in the three-dimensional rectangular velocity distribution, where z j_i_moi represents the corrected coordinate of the i-th point at the j-th spanwise position, z j_i_ori represents the initial coordinate of the i-th point at the j-th spanwise position, h j represents the z-direction height at the j-th spanwise position in the bow region, δ wing_lp is the thickness of the low-energy airflow region;

[0061] Step S15: Assemble the velocity distributions of the wing boundary layer region and the low-energy airflow region after the "bowing" process obtained in step S14 with the three-dimensional velocity distribution "base" on the front interface, and perform smoothing on the regions where the wing boundary layer region and the low-energy airflow region meet the "base" to obtain the velocity distribution of the front interface flow field considering the hub boundary layer, duct boundary layer, wing boundary layer, and low-energy airflow region;

[0062] Step S16: Based on the velocity distribution of the front interface flow field obtained in step S15, use the linear function scaling method to correct the velocity distributions of the duct lip boundary layer region and the duct lip acceleration region to obtain the final velocity distribution of the front interface flow field.

[0063] Furthermore, the assembly method for obtaining the two-dimensional radial velocity distribution on the front interface in step S13 is to merge the two-dimensional matrices containing coordinate information and velocity information, as shown in the following formula:

[0064]

[0065] In the formula, z f_hub is the position coordinate corresponding to the velocity distribution u f_hub of the front hub boundary layer region, z gd is the position coordinate corresponding to the velocity distribution u gd of the front uniform transition region, z f_duct is the position coordinate corresponding to the velocity distribution u f_duct of the front duct boundary layer region; z 2D represents the coordinate distribution after assembling the front hub boundary layer region, the front uniform transition region, and the front duct boundary layer region, and u 2D represents the velocity distribution after assembling the front hub boundary layer region, the front uniform transition region, and the front duct boundary layer region.

[0066] Furthermore, the specific process of assembling and smoothing the velocity distributions of the wing boundary layer region and the low-energy airflow region after the "bowing" process with the three-dimensional velocity distribution "base" on the front interface in step S15 is as follows:

[0067] Delete the region in the "base" of the three-dimensional velocity distribution where the height is lower than the top outer edge of the wing boundary layer region and the low-energy airflow region, and replace it with the velocity distribution of the wing boundary layer region and the low-energy airflow region after "bowing" treatment. For the step of the velocity distribution at the junction position, use the method of linear interpolation to smooth and correct the velocity step region, and obtain the velocity distribution of the front interface flow field considering the hub boundary layer, duct boundary layer, wing boundary layer and low-energy airflow region.

[0068] Furthermore, the specific process of step S16 is as follows:

[0069] ① Establish a linear scaling factor function. For the linear scaling function attenuation in the duct lip boundary layer region and the duct lip acceleration region, it mainly includes two parts: radial attenuation and circumferential attenuation. The formulas are as follows:

[0070] Radial scaling factor r factor :

[0071]

[0072] In the formula, r i represents the radius of the current grid point, r represents the duct radius, and σ r represents the parameter set to control the radial scaling range;

[0073] Circumferential scaling factor θ factor :

[0074]

[0075] In the formula, θ i -θ center represents the circumferential angle difference between the current grid point and the center of the duct lip boundary layer, and σ θ represents the parameter set to control the circumferential scaling range;

[0076] ② Multiply the velocity distribution of the front interface flow field obtained in step S15 by the corresponding scaling function to obtain the corrected velocity distribution of the front interface flow field for the duct lip boundary layer region and the duct lip acceleration region; specifically:

[0077] For the duct lip boundary layer region, which is a region where the velocity decreases, use the product of the circumferential and radial scaling factors as the total scaling factor sc factor_bl :

[0078] sc factor_bl = r factor ·θ factor

[0079] Based on the total scaling factor of the duct lip boundary layer region, use the following expression to perform the first scaling correction on the velocity on the front interface:

[0080] u isc1 = u i ·(1 - sc factor_bl )

[0081] Expanded as:

[0082]

[0083] In the formula, u i is the flow field velocity at each point on the front interface, and u isc1 is the flow field velocity at each point on the front interface after the first scaling correction of the velocity on the front interface;

[0084] For the ducted lip acceleration region, which is a region where the velocity increases, the following expression is used as the total scaling factor sc for the ducted lip acceleration region factor_ac :

[0085] sc factor_ac = 1 + (ms - 1)·r factor ·θ factor

[0086] In the formula, ms is the maximum value of velocity scaling. The curve at the maximum curvature of the ducted lip can be used as the leading edge. After supplementing it to a complete airfoil, the potential flow method is used to obtain the airflow velocity ratio at the relative position of the front interface. According to empirical values, it is taken as 1.2 - 1.4;

[0087] Based on the total scaling factor of the ducted lip acceleration region, the following expression is used for the second scaling correction of the velocity on the front interface:

[0088] u isc2 = u isc1 ·sc factor_ac

[0089] Expanded as:

[0090]

[0091] In the formula, u isc2 is the flow field velocity at each point on the front interface after the second scaling correction of the velocity on the front interface, that is, the final flow field velocity distribution of the front interface.

[0092] Furthermore, in step S2, the rear hub boundary layer region refers to the annular region from the hub wall surface to the outer boundary of the hub boundary layer on the rear interface; the stator wake region refers to obtaining the central contour of the stator wake region by projecting the shape of the stator trailing edge onto the rear interface, and then expanding the central contour of the stator wake region circumferentially The area of the angular range; the ducted exponential distribution transition area refers to the remaining area on the rear interface after removing the above-mentioned partition.

[0093] Furthermore, the step S2 includes the following sub-steps:

[0094] Step S21: Calculate the average velocity V at the rear interface under different oncoming flow velocities and different thrust conditions based on the momentum theory b1 ;

[0095] The momentum theory is a well-known method in the field. For the specific process of calculating the average velocity at the rear interface, refer to the calculation of the average velocity at the front interface.

[0096] Step S22: Calculate the thickness and two-dimensional velocity distribution of the rear hub boundary layer area based on the thickness formula of the turbulent boundary layer around a downstream-facing flat plate and the assumption of velocity power-law distribution, and calculate the two-dimensional velocity distribution of the ducted exponential distribution transition area based on the assumption of the velocity distribution law in a circular pipe;

[0097] The thickness formula of the turbulent boundary layer around a downstream-facing flat plate, the assumption of velocity power-law distribution, and the assumption of the velocity distribution law in a circular pipe are all well-known technologies in the field. The specific process of this step is as follows:

[0098] ① Calculate the thickness of the rear hub boundary layer area using the calculation formula for the thickness of the flat plate turbulent boundary layer:

[0099]

[0100] In the formula, l S is the length of the stator blade domain in the single-ducted fan-wing segment coupling configuration, and Re b_hub is the Reynolds number calculated with l S as the characteristic length. The calculation formula is:

[0101]

[0102] In the formula, μ is the dynamic viscosity coefficient, and V b1 is the average velocity at the rear interface;

[0103] Based on the thickness δ b_hub of the rear hub boundary layer area, use the power function distribution to calculate the internal velocity distribution of the rear hub boundary layer area, that is:

[0104]

[0105] In the formula, u b_hub is the velocity at different normal positions inside the rear hub boundary layer area, z b_hub represents the coordinates of different normal positions inside the rear hub boundary layer, and V b_maxis the maximum velocity in the two-dimensional radial velocity distribution on the rear interface, taking the velocity at the outer edge of the rear hub boundary layer, V b_max The flow velocity distribution law in a circular pipe is used for calculation, and the specific expression is as follows:

[0106]

[0107] In the formula, the constant n takes 7, V b1 is the average velocity at the rear interface;

[0108] From the coordinate point information z in the rear hub boundary layer region b_hub and the velocity information u b_hub the velocity information matrix of the rear hub boundary layer region is obtained:

[0109] V b_hub =[z b_hub , u b_hub

[0110] ② The velocity distribution in the transition zone of the duct index distribution is calculated by using the exponential formula of the flow velocity distribution in a circular pipe, that is:

[0111]

[0112] In the formula, u b_duct is the velocity at different normal positions in the transition zone of the duct index distribution, r - δ b_hub is the length from the outer ring to the outer edge of the rear hub boundary layer, z b_duct is the coordinate of different normal positions in the transition zone of the duct index distribution, and the constant m takes the value of 10;

[0113] From the coordinate point information z in the transition zone of the duct index distribution b_duct and the velocity information u b_duct the velocity information matrix of the transition zone of the duct index distribution is obtained:

[0114] V b_duct =[z b_duct , u b_duct

[0115] Step S23: Assemble the two-dimensional velocity distribution in the rear hub boundary layer region and the two-dimensional velocity distribution in the transition zone of the duct index distribution obtained in step S22 according to the spatial position to obtain the two-dimensional radial velocity distribution on the rear interface, and rotate the obtained two-dimensional radial velocity distribution one week along the circumferential direction around the fan axis to generate the three-dimensional velocity distribution "base" on the rear interface;

[0116] ​​Step S24: Generate the central contour of the stator wake region based on the shape of the trailing edge of the stator blade. Specifically, project the shape of the trailing edge of the stator blade onto the rear interface to obtain the centerline of the stator wake. On the rear interface, extend both ends of the centerline of the stator wake to the hub of the rear interface and the duct of the rear interface as the central contour of the stator wake region; based on the central contour of the stator wake region, obtain the circumferential angle θ of the central contour of the stator wake region. center ;

[0117] Step S25: Based on the three-dimensional velocity distribution "base" on the rear interface obtained in step S23 and the central contour of the stator wake region obtained in step S24, use the linear function scaling method to correct the velocity distribution in the transition region of the duct index distribution to obtain the final velocity distribution of the flow field on the rear interface.

[0118] Furthermore, the assembly method for obtaining the two-dimensional radial velocity distribution on the rear interface in step S23 is to merge two-dimensional matrices containing coordinate information and velocity information, as shown in the following formula:

[0119]

[0120] In the formula, z b_2D represents the coordinate distribution after assembling the rear hub boundary layer region and the transition region of the duct index distribution, and u b_2D represents the velocity distribution after assembling the rear hub boundary layer region and the transition region of the duct index distribution.

[0121] Furthermore, the specific process of step S25 is as follows:

[0122] The circumferential scaling factor θ of the wake wj_factor is:

[0123]

[0124] In the formula, θ i -θ center represents the circumferential angle difference between the current grid point and the central contour of the stator wake region, and σ θ represents the parameter for setting the control circumferential scaling range;

[0125] Based on the circumferential scaling factor θ of the wake wj_factor , use the following expression to scale and correct the velocity on the rear interface:

[0126] u isc_b = u i_b ·(1 - θ wj_factor )

[0127] In the formula, u i_b is the flow field velocity at each grid point on the "base" distribution of the rear interface, and u isc_bThe flow field velocity at each grid point of the rear interface after scaling and correcting the velocity at the rear interface is the final flow field velocity distribution of the rear interface.

[0128] Beneficial effects

[0129] The present invention provides an interface flow field information generation method applicable to the ducted fan momentum source method, which separately processes the front and rear interfaces. First, the average velocity at the fan under different oncoming flow velocities and different thrusts is calculated through momentum theory, and then the velocity distributions of the front and rear interfaces considering the influence of multiple components are generated. Finally, the obtained velocity distributions are assembled and merged to obtain the velocity distributions on the front and rear interfaces. This method can get rid of the dependence on the flow field calculation of the multiple reference coordinate systems, can also obtain relatively accurate momentum source distribution characteristics, and the calculation efficiency is improved. Description of the drawings

[0130] Figure 1 It is a flowchart of an embodiment of the present invention;

[0131] Figure 2 It is a single ducted fan-wing segment coupling configuration with a ducted fan in an embodiment of the present invention;

[0132] Figure 3 It is a single ducted fan-wing segment coupling configuration calculated using the momentum source in an embodiment of the present invention;

[0133] Figure 4 It is a schematic diagram of the partition of the front interface in an embodiment of the present invention;

[0134] Figure 5 It is a schematic diagram of the partition of the rear interface in an embodiment of the present invention;

[0135] Figure 6 It is a schematic diagram of the velocity distribution and partition of the central cross-section on the front interface in an embodiment of the present invention;

[0136] Figure 7 It is a schematic diagram of the velocity distribution and partition of the lower central cross-section on the front interface in an embodiment of the present invention;

[0137] Figure 8 It is a schematic diagram of the velocity distribution and partition of the central cross-section on the rear interface in an embodiment of the present invention;

[0138] Figure 9 It is the two-dimensional radial velocity distribution on the front interface in an embodiment of the present invention;

[0139] Figure 10 It is the "base" generated based on the two-dimensional radial velocity distribution on the front interface in an embodiment of the present invention;

[0140] Figure 11 It is the velocity distribution of the wing boundary layer and the low-energy air flow region in an embodiment of the present invention;

[0141] Figure 12 The "bowing" operation effect of the embodiment of the present invention;

[0142] Figure 13 The "bowing" region of the embodiment of the present invention is assembled with the "substrate";

[0143] Figure 14 Comparison of the front interface velocity distribution generated by the method of the present invention and the momentum source distribution calculated from the front interface velocity distribution extracted from the MRF flow field;

[0144] Among them, Figure 14 (a) is the momentum source distribution calculated based on the front interface velocity distribution generated by the method of the present invention, Figure 14 (b) is the momentum source distribution calculated based on the front interface velocity distribution extracted from the MRF flow field;

[0145] Figure 15 The two-dimensional radial velocity distribution of the rear interface of the embodiment of the present invention;

[0146] Figure 16 The "substrate" generated based on the two-dimensional radial velocity distribution of the rear interface of the embodiment of the present invention;

[0147] Figure 17 The wake center contour of the rear interface of the embodiment of the present invention;

[0148] Figure 18 Comparison of the rear interface velocity distribution generated by the method of the present invention and the momentum source distribution calculated from the rear interface velocity distribution extracted from the MRF flow field;

[0149] Among them, Figure 18 (a) is the momentum source distribution calculated based on the front interface velocity distribution generated by the method of the present invention, Figure 18 (b) is the momentum source distribution calculated based on the front interface velocity distribution extracted from the MRF flow field;

[0150] Figure 19 The lift coefficient, drag coefficient, and relative error of the drag coefficient with respect to the MRF method calculated by the method of the present invention at different oncoming flow velocities;

[0151] Figure 20 The lift coefficient, drag coefficient, and relative error of the drag coefficient with respect to the MRF method calculated by the method of the present invention at different fan rotational speeds;

[0152] In the figure: 1 - duct outer shell; 2 - rotor blade; 3 - hub; 4 - duct lip; 5 - exhaust duct; 6 - stator blade; 7 - wing; 8 - front interface; 9 - rear interface; 11 - duct lip boundary layer region; 12 - duct lip acceleration region; 13 - wing boundary layer region; 14 - low-energy airflow region; 15 - front hub boundary layer region; 16 - front uniform transition region; 17 - front duct boundary layer region; 21 - rear hub boundary layer; 22 - stator blade wake region; 23 - duct exponential distribution transition region. Detailed implementation manner

[0153] In order to make the technical problems, technical solutions and beneficial effects solved by the present invention clearer and more understandable, and enable those skilled in the art to better understand the present invention, the following further details and fully describe the present invention in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not used to limit the present invention.

[0154] This embodiment adopts an interface flow field information generation method suitable for the momentum source method of a ducted fan. Taking the wing section with a ducted fan at the trailing edge of the upper surface of the wing as an example, to verify the effectiveness and practicability of the method of the present invention, the wing section model with a ducted fan at the trailing edge of the upper surface of the wing is as Figure 2 shown, and the wing section model with front and rear interfaces at the trailing edge of the upper surface of the wing calculated by the momentum source method is as Figure 3 shown.

[0155] As Figure 1 shown, an interface flow field information generation method suitable for the momentum source method of a ducted fan in this embodiment includes the following steps:

[0156] Step S1: Denote the interface in front of the rotor blade of the ducted fan as the front interface, and divide the front interface into a front hub boundary layer region, a front duct boundary layer region, a wing boundary layer region, a low-energy airflow region, a duct lip boundary layer region, a duct lip acceleration region, and a front uniform transition region, as Figure 4 shown; then generate the velocity distributions in the above-mentioned divided regions respectively, and then assemble and merge them to obtain the velocity distribution on the front interface;

[0157] The front hub boundary layer region refers to the annular region from the hub wall surface to the outer boundary of the hub boundary layer on the front interface; the front duct boundary layer region refers to the annular region from the inner wall surface of the duct to the outer boundary of the duct boundary layer on the front interface; the wing boundary layer region refers to the rectangular region from the upper surface of the wing to the outer boundary of the wing boundary layer; the low-energy airflow region refers to the region from the outer boundary of the wing boundary layer to the average flow velocity position of the interface; the duct lip boundary layer region refers to the region from the duct lip wall surface to the outer boundary of the duct lip boundary layer; the duct lip acceleration region refers to the region obtained by taking the intersection of the circumferential action range region and the radial action range region. The circumferential action range region is: respectively with is the central angle, and the range is used as the circumferential acting range. The radial acting range area is: taking the intersection of each central angle ray and the outer ring of the duct as the radial center position, and the range of 1 / 10 to 1 / 5 of the duct radius around each radial center position is used as the radial acting range; the front uniform transition area is the remaining area on the front interface after removing the above-mentioned partitions;

[0158] Step S2: Denote the interface behind the stator blades of the ducted fan as the rear interface, and divide the rear interface into the rear hub boundary layer area, the duct exponential distribution transition area, and the stator blade wake area, as Figure 5 shown; then generate the velocity distributions in the above-mentioned divided areas respectively, and then assemble and merge them to obtain the velocity distribution on the rear interface.

[0159] The rear hub boundary layer area refers to the annular area on the rear interface from the hub wall surface to the outer boundary of the hub boundary layer; the stator blade wake area refers to obtaining the central contour of the stator blade wake area by projecting the shape of the trailing edge of the stator blade onto the rear interface, and then expanding the central contour of the stator blade wake area circumferentially by the angular range; the duct exponential distribution transition area refers to the remaining area on the rear interface after removing the above-mentioned partitions;

[0160] In this embodiment, the method for obtaining the front and rear interfaces refers to the method for obtaining the interface in front of the rotor blade and the interface behind the stator blade in Chinese Patent CN117371344A.

[0161] In this embodiment, step S1 includes the following sub-steps:

[0162] Step S11: Calculate the average velocity at the front interface under different oncoming flow velocities and different thrust conditions based on the momentum theory;

[0163] The momentum theory is a well-known method in the art. The specific process of calculating the average velocity at the front interface is as follows:

[0164] First, according to the momentum theory, the expression for the fan thrust T P is obtained:

[0165]

[0166] In the formula, T P represents the thrust of the fan, ρ is the atmospheric density, A 1 is the flow channel area at the front interface, V 0 is the free oncoming flow velocity, V 2 is the duct exit velocity. From the above formula, the velocity at the duct exit under different oncoming flow velocities and different thrust conditions can be obtained as:

[0167]

[0168] The average velocity V at the front interface is obtained by mass conservation calculation 1 :

[0169]

[0170] where A 2 is the cross-sectional area of the duct outlet;

[0171] In this embodiment, the duct diameter of the fan is 0.15 m, the hub diameter is 0.06 m, the fan thrust T P is 10.06 N, the atmospheric density is 1.11166 kg / m 3 , the flow area A at the front interface 1 is 0.014844 m 2 , the free stream velocity V 0 is 41.667 m / s, the duct outlet area A 2 is 0.0158 m 2 , therefore, the average velocity V at the front interface 1 is 57.85 m / s.

[0172] Step S12: Based on the thickness formula of the laminar boundary layer around a downstream-facing flat plate and the assumption of velocity power distribution, calculate the thickness and two-dimensional velocity distribution of the front hub boundary layer region; based on the thickness formula of the turbulent boundary layer around a downstream-facing flat plate and the assumption of velocity power distribution, calculate the thickness and two-dimensional velocity distribution of the front duct boundary layer region; interpolate and calculate the two-dimensional velocity distribution of the front uniform transition region based on the thickness and two-dimensional velocity distribution of the front hub boundary layer region and the thickness and two-dimensional velocity distribution of the front duct boundary layer region;

[0173] The thickness formula of the laminar boundary layer around a downstream-facing flat plate and the assumption of velocity power distribution, and the thickness formula of the turbulent boundary layer around a downstream-facing flat plate and the assumption of velocity power distribution are all well-known techniques in the art. The specific process of this step is as follows:

[0174] ① Use the calculation formula for the thickness of the laminar boundary layer of a flat plate placed downstream to calculate the thickness of the front hub boundary layer region. The specific expression is:

[0175]

[0176] where δ f_hub is the thickness of the front hub boundary layer, x hub is the chordwise length of the front cone of the hub, and Re x_f_hub is the Reynolds number calculated based on the chordwise length of the front cone of the hub;

[0177] In this embodiment, the chordwise length of the front cone of the hub is 0.04 m, and the Reynolds number calculated from the chordwise length of the front cone of the hub is 1.4484×105 , the thickness of the boundary layer of the front propeller hub is 5.2552×10 -4 m.

[0178] Based on the thickness δ of the boundary layer region of the front propeller hub f_hub , the internal velocity distribution in the boundary layer region of the front propeller hub is calculated using a power function distribution, that is:

[0179]

[0180] In the formula, u f_hub is the velocity at different normal positions inside the boundary layer region of the front propeller hub, z f_hub represents the coordinates of different normal positions inside the boundary layer of the front propeller hub, V f_hub is the velocity at the outer edge of the boundary layer of the front propeller hub, and is calculated using the following formula:

[0181] V f_hub = ζV 1

[0182] In the formula, ζ 1 is the correction coefficient of the velocity at the outer edge of the boundary layer of the front propeller hub, which characterizes the acceleration effect of the propeller hub surface on the air flow, and takes 1.08;

[0183] ② Calculate the thickness of the boundary layer region of the front duct using the calculation formula for the thickness of the turbulent boundary layer around a flat plate placed along the flow direction. The specific expression is:

[0184]

[0185] In the formula, δ f_duct is the thickness of the boundary layer region of the front duct, x f_duct is the chord length of the duct, and Re x_f_duct is the Reynolds number calculated based on the chord length of the duct;

[0186] In this embodiment, the chord length of the duct is 0.02 m, and the Reynolds number calculated from the chord length of the duct is 7.2418×10 4 , and the thickness of the boundary layer of the front propeller hub is 7.076×10 -4 m.

[0187] Based on the thickness δ of the boundary layer region of the front duct f_duct , the internal velocity distribution in the boundary layer region of the front duct is calculated using a power function distribution, that is:

[0188]

[0189] In the formula, u f_duct is the velocity distribution at different normal positions inside the boundary layer region of the front duct, z f_duct represents the coordinates of different normal positions inside the boundary layer of the front duct, V f_ductThe velocity at the outer edge of the front duct boundary layer is calculated using the following formula:

[0190] V f_duct = ζ 2 V 1

[0191] In the formula, ζ 2 is the correction coefficient of the velocity at the outer edge of the front duct boundary layer, which characterizes the acceleration effect of the duct surface on the air flow and takes the value of 1.24;

[0192] ③ Based on the calculated thickness and velocity distributions of the front hub boundary layer region and the front duct boundary layer region, the two-dimensional velocity distribution u of the front uniform transition region is obtained by using the linear interpolation method gd ;

[0193] Specifically, the outer edge coordinate values of the hub boundary layer and the front duct boundary layer are obtained from the thicknesses of the front hub boundary layer region and the front duct boundary layer region. A linear interpolation function is established between the outer edges of the front hub boundary layer region and the front duct boundary layer region to interpolate the velocity, so that the velocity smoothly transitions between the outer edges of the front hub boundary layer region and the front duct boundary layer region, and the two-dimensional velocity distribution u of the front uniform transition region is obtained gd ;

[0194] Step S13: Assemble the two-dimensional velocity distributions of the front hub boundary layer region, the front uniform transition region, and the front duct boundary layer region obtained in step S12 according to their spatial positions to obtain the two-dimensional radial velocity distribution on the front interface. Rotate the obtained two-dimensional radial velocity distribution one week around the fan axis to generate the three-dimensional velocity distribution "base" on the front interface;

[0195] The assembly method for obtaining the two-dimensional radial velocity distribution on the front interface is to merge the two-dimensional matrices containing coordinate information and velocity information, as shown in the following formula:

[0196]

[0197] In the formula, z gd is the position coordinate corresponding to the velocity distribution u gd of the front uniform transition region, z 2D represents the coordinate distribution after assembling the front hub boundary layer region, the front uniform transition region, and the front duct boundary layer region, and u 2D represents the velocity distribution after assembling the front hub boundary layer region, the front uniform transition region, and the front duct boundary layer region, as Figure 9 shown; the generated three-dimensional velocity distribution "base" on the front interface is as Figure 10 shown;

[0198] Step S14: Use the Karman-Pohlhausen method to calculate the thickness and two-dimensional velocity distribution of the wing boundary layer region, use xfoil to calculate the thickness of the low-energy airflow region, and based on the thickness and two-dimensional velocity distribution of the wing boundary layer region and the thickness of the low-energy airflow region, use cubic spline distribution to calculate the two-dimensional velocity distribution of the low-energy airflow region; Assemble the two-dimensional velocity distribution of the wing boundary layer region and the two-dimensional velocity distribution of the low-energy airflow region and perform spanwise expansion to obtain the three-dimensional velocity distribution of the wing boundary layer region and the low-energy airflow region, and perform "bowing" processing on the obtained three-dimensional velocity distribution;

[0199] The Karman-Pohlhausen approximation method is a well-known method in the art, and xfoil is a well-known calculation program in the art. The specific process of this step is as follows:

[0200] ① Use the classical Karman-Pohlhausen approximation method to calculate the thickness and velocity distribution of the wing boundary layer region. This method establishes the relationship between the wall shear stress and the velocity distribution through the definition of the momentum thickness, combines the no-slip condition in the wall normal direction and the far-field free-stream velocity condition, and obtains the wing boundary layer thickness δ wing_bl and velocity distribution u wing_bl by iteratively solving the boundary layer integral equation. For specific reference, see Section 3, Chapter 5 of "Viscous Fluid Mechanics" (written by Zhang Zixiong and Dong Zengnan); Multiply the wing boundary layer velocity distribution obtained by Karman-Pohlhausen by the ratio of the average velocity V 1 at the front interface to the free-stream velocity V 0 to obtain the boundary layer velocity distribution considering the fan suction effect; In this embodiment, the wing boundary layer thickness δ wing_bl is 7.84×10 -4 m;

[0201] ② First, use XFOIL to calculate the momentum loss thickness of the low-energy airflow region, and then calculate the nominal thickness of the boundary layer of the low-energy airflow region based on the momentum loss thickness and subtract the thickness of the wing boundary layer as the thickness of the low-energy airflow region; Specifically:

[0202] The momentum loss thickness δ 2 of the low-energy airflow region is calculated as

[0203]

[0204] where δ is the nominal thickness of the boundary layer of the low-energy airflow region, Λ is a dimensionless quantity, and its expression is:

[0205]

[0206] where υ is the kinematic viscosity coefficient, is the velocity gradient along the wing surface;

[0207] From the above two equations, we get:

[0208]

[0209] Solving the above equation gives the nominal thickness δ of the boundary layer in the low-energy airflow region. Subtracting the thickness δ of the wing boundary layer from the calculated δ wing_bl yields the thickness δ of the low-energy airflow region wing_lp ; in this embodiment, the nominal thickness of the boundary layer is 0.0203 m, and the thickness δ of the low-energy airflow region wing_lp is 0.019516 m;

[0210] ③ Based on the obtained thickness of the wing boundary layer region, obtain the coordinates of the outer edge of the wing boundary layer. Based on the thickness δ of the low-energy airflow region wing_lp , obtain the coordinates of the outer edge of the low-energy airflow region. Let the outer edge velocity value of the low-energy airflow region be equal to the average velocity V of the front interface calculated in step S11 1 , and then, according to the coordinates of the outer edge of the wing boundary layer and the coordinates of the outer edge of the low-energy airflow region, together with the tangency constraints of the velocity distribution in the low-energy airflow region with the velocity distribution in the wing boundary layer and the velocity distribution in the front uniform transition region, use cubic spline distribution to calculate and obtain the velocity distribution in the low-energy airflow region;

[0211] ④ Combine the obtained velocity distribution in the low-energy airflow region with the velocity distribution in the wing boundary layer; obtain the two-dimensional velocity distribution of the wing boundary layer and the low-energy airflow region, as Figure 11 shown, and extend the two-dimensional velocity distribution along the span direction to the outer edge of the duct. In this embodiment, the length of the extension along the span direction is 0.10262 m to generate a three-dimensional rectangular velocity distribution;

[0212] ⑤ Perform "bowing" processing on the obtained three-dimensional velocity distribution. The specific process is as follows:

[0213] The calculation method of the z coordinate in the height direction of the bowing is as follows:

[0214]

[0215] where the coordinate system is defined as: the X axis is along the free-stream direction, the Y axis is to the right from the pilot's perspective, and the Z direction is perpendicular to the XY plane upward, satisfying the right-hand system; the calculation object in the formula is the velocity distribution scatter points at different span positions in the three-dimensional rectangular velocity distribution. Among them, z j_i_moi represents the corrected coordinate of the i-th point at the j-th span position, z j_i_ori represents the initial coordinate of the i-th point at the j-th span position, h j represents the z-direction height at the j-th span position in the bow region, δ wing_lp is the thickness of the low-energy airflow region, as Figure 12 shown;

[0216] Step S15: Assemble the velocity distributions of the wing boundary layer region and the low-energy airflow region after the "bowing" process obtained in Step S14 with the three-dimensional velocity distribution "base" obtained in Step S13, and perform smoothing on the regions where the wing boundary layer region and the low-energy airflow region meet the "base" to obtain the velocity distribution of the front interface flow field considering the hub boundary layer, duct boundary layer, wing boundary layer, and low-energy airflow region;

[0217] The specific process is as follows:

[0218] Delete the regions in the three-dimensional velocity distribution "base" where the height is lower than the top outer edge of the wing boundary layer region and the low-energy airflow region, and replace them with the velocity distributions of the wing boundary layer region and the low-energy airflow region after the "bowing" process. For the velocity step at the junction position, use linear interpolation to smooth and correct the velocity step region to obtain the velocity distribution of the front interface flow field considering the hub boundary layer, duct boundary layer, wing boundary layer, and low-energy airflow region, as Figure 13 shown.

[0219] Step S16: Based on the velocity distribution of the front interface flow field obtained in Step S15, use the linear function scaling method to correct the velocity distributions of the duct lip boundary layer region and the duct lip acceleration region to obtain the final velocity distribution of the front interface flow field;

[0220] The specific process is as follows:

[0221] ① Establish a linear scaling factor function. For the linear scaling function attenuation of the duct lip boundary layer region and the duct lip acceleration region, it mainly includes two parts: radial attenuation and circumferential attenuation. The formulas are as follows:

[0222] Radial scaling factor r factor :

[0223]

[0224] In the formula, r i represents the radius of the current grid point, r represents the duct radius, and σ r represents the parameter set to control the radial scaling range;

[0225] Circumferential scaling factor θ factor :

[0226]

[0227] In the formula, θ i -θ center represents the circumferential angle difference between the current grid point and the center of the duct lip boundary layer, and σ θ represents the parameter set to control the circumferential scaling range;

[0228] ② Multiply the flow field velocity distribution at the front interface obtained in step S15 by the corresponding scaling function to obtain the corrected flow field velocity distribution at the front interface for the boundary layer region of the duct lip and the acceleration region of the duct lip. Specifically:

[0229] For the boundary layer region of the duct lip, which is a region where the velocity decreases, the product of the circumferential and radial scaling factors is used as the total scaling factor sc for the boundary layer region of the duct lip factor_bl :

[0230] sc factor_bl = r factor · θ factor

[0231] Based on the total scaling factor for the boundary layer region of the duct lip, the following expression is used to perform the first scaling correction on the velocity at the front interface:

[0232] u isc1 = u i · (1 - sc factor_bl )

[0233] Expanded as:

[0234]

[0235] In the formula, u i is the flow field velocity at each point on the front interface, and u isc1 is the flow field velocity at each point on the front interface after the first scaling correction of the velocity at the front interface;

[0236] For the acceleration region of the duct lip, which is a region where the velocity increases, the following expression is used as the total scaling factor sc for the acceleration region of the duct lip factor_ac :

[0237] sc factor_ac = 1 + (ms - 1) · r factor · θ factor

[0238] In the formula, ms is the maximum value of velocity scaling. The curve at the maximum curvature of the duct lip can be used as the leading edge. After supplementing it to a complete airfoil, the potential flow method is used to obtain the air flow velocity ratio at the relative position of the front interface. According to experience, it is taken as 1.3;

[0239] Based on the total scaling factor for the acceleration region of the duct lip, the following expression is used to perform the second scaling correction on the velocity at the front interface:

[0240] u isc2 = u isc1 · sc factor_ac

[0241] Expanded as:

[0242]

[0243] wherein, u isc2 is the flow field velocity at each point on the front interface after the second scaling correction of the velocity on the front interface, that is, the final flow field velocity distribution on the front interface; in this embodiment, the comparison between the front interface momentum source distribution generated based on the front interface velocity distribution obtained by the method of the present invention and the result of MRF calculation is as Figure 14 shown.

[0244] In this embodiment, the step S2 includes the following sub-steps:

[0245] Step S21: Calculate the average velocity V at the rear interface under different oncoming flow velocities and different thrust conditions based on the momentum theory b1 ;

[0246] The momentum theory is a well-known method in the art. The specific process of calculating the average velocity at the rear interface refers to the calculation of the average velocity at the front interface; in this embodiment, since the flow channel area of the front interface and the flow channel area of the rear interface are the same, therefore, V b1 = V 1 ;

[0247] Step S22: Calculate the thickness and two-dimensional velocity distribution of the rear hub boundary layer region based on the thickness formula of the turbulent boundary layer around the downstream flat plate and the assumption of the velocity power distribution, and calculate the two-dimensional velocity distribution of the duct index distribution transition region based on the assumption of the flow velocity distribution law in the circular pipe;

[0248] The thickness formula of the turbulent boundary layer around the downstream flat plate, the assumption of the velocity power distribution, and the assumption of the flow velocity distribution law in the circular pipe are all well-known technologies in the art. The specific process of this step is as follows:

[0249] ① Calculate the thickness of the rear hub boundary layer region using the calculation formula of the flat plate turbulent boundary layer thickness:

[0250]

[0251] wherein, l S is the length of the stator blade domain in the single-duct fan-wing segment coupling configuration, and Re b_hub is the Reynolds number calculated with l S as the characteristic length. The calculation formula is:

[0252]

[0253] wherein, μ is the dynamic viscosity coefficient, and V b1 is the average velocity at the rear interface; in this embodiment, δ b_hub = 0.0023

[0254] Based on the thickness δ of the rear hub boundary layer region b_hub , the internal velocity distribution in the rear hub boundary layer region is calculated using a power function distribution, that is:

[0255]

[0256] In the formula, u b_hub is the velocity at different normal positions inside the rear hub boundary layer region, z b_hub represents the coordinates of different normal positions inside the rear hub boundary layer, V b_max is the maximum velocity in the two-dimensional radial velocity distribution at the rear interface, taking the velocity at the outer edge of the rear hub boundary layer, V b_max is calculated using the flow velocity distribution law in a circular pipe, and the specific expression is as follows:

[0257]

[0258] In the formula, the constant n takes 7, V b1 is the average velocity at the rear interface. In this embodiment, V b_max is 66.133 m / s;

[0259] From the coordinate point information z b_hub and velocity information u b_hub of the rear hub boundary layer region, the velocity information matrix of the rear hub boundary layer region is obtained:

[0260] V b_hub =[z b_hub , u b_hub

[0261] ② The velocity distribution in the transition zone of the duct index distribution is calculated using the exponential formula of the flow velocity distribution in a circular pipe, that is:

[0262]

[0263] In the formula, u b_duct is the velocity at different normal positions in the transition zone of the duct index distribution, r - δ b_hub is the length from the outer ring to the outer edge of the rear hub boundary layer, z b_duct is the coordinate of different normal positions in the transition zone of the duct index distribution, and the constant m takes the value of 10;

[0264] From the coordinate point information z b_duct and velocity information u b_duct of the transition zone of the duct index distribution, the velocity information matrix of the transition zone of the duct index distribution is obtained:

[0265] V b_duct =[z b_duct , u b_duct ​​

[0266] Step S23: Assemble the two-dimensional velocity distribution in the rear hub boundary layer region and the two-dimensional velocity distribution in the transition region of the duct index distribution obtained in Step S22 according to the spatial position to obtain the two-dimensional radial velocity distribution on the rear interface. Rotate the obtained two-dimensional radial velocity distribution one week along the circumferential direction around the fan axis to generate the "base" of the three-dimensional velocity distribution on the rear interface;

[0267] The assembly method for obtaining the two-dimensional radial velocity distribution on the rear interface is to merge the two-dimensional matrices containing coordinate information and velocity information, as shown in the following formula:

[0268]

[0269] In the formula, z b_2D represents the coordinate distribution after assembling the rear hub boundary layer region and the transition region of the duct index distribution, and u b_2D represents the velocity distribution after assembling the rear hub boundary layer region and the transition region of the duct index distribution, as Figure 15 shown; the "base" of the three-dimensional velocity distribution generated on the rear interface is as Figure 16 shown;

[0270] Step S24: Generate the central contour of the stator wake region based on the shape of the stator trailing edge. Specifically, obtain the stator wake centerline by projecting the shape of the stator trailing edge onto the rear interface, and extend both ends of the stator wake centerline to the hub of the rear interface and the duct of the rear interface on the rear interface as the central contour of the stator wake region, as Figure 17 shown; based on the central contour of the stator wake region, obtain the circumferential angle θ center of the central contour of the stator wake region;

[0271] Step S25: Based on the "base" of the three-dimensional velocity distribution on the rear interface obtained in Step S23 and the central contour of the stator wake region obtained in Step S24, use the linear function scaling method to correct the velocity distribution in the transition region of the duct index distribution to obtain the final velocity distribution of the flow field on the rear interface;

[0272] The specific process is as follows:

[0273] The wake circumferential scaling factor θ wj_factor is:

[0274]

[0275] In the formula, θ i -θ center represents the circumferential angle difference between the current grid point and the central contour of the stator wake region, and σ θ represents the parameter for setting the control circumferential scaling range;

[0276] Based on the wake circumferential scaling factor θ wj_factor , the velocity on the rear interface is scaled and corrected using the following expression:

[0277] u isc_b = u i_b ·(1 - θ wj_factor )

[0278] In the formula, u i_b is the flow field velocity at each grid point on the "base" distribution of the rear interface, and u isc_b is the flow field velocity at each grid point on the rear interface after scaling and correcting the velocity on the rear interface, that is, the final flow field velocity distribution of the rear interface; in this embodiment, the comparison between the rear interface momentum source distribution generated based on the rear interface velocity distribution obtained by the method of the present invention and the result of MRF calculation is as Figure 18 shown.

[0279] Based on Figure 3 the bladeless grid configuration of the single-ducted fan-wing segment fusion configuration shown, the momentum source generation method developed in this patent is verified. The flight conditions and model geometric parameters are shown in Table 1;

[0280] Table 1 Calculation conditions and geometric parameters of the single-ducted fan-wing segment fusion configuration

[0281]

[0282]

[0283] Based on the Latin hypercube method, 25 sample points are randomly selected in the design space from 35 m / s to 60 m / s, and this velocity range covers the velocity from climb to cruise. First, the MRF method is used to calculate the flow field of the single-ducted fan-wing segment fusion configuration at a 4° angle of attack under different oncoming flow velocities, and the aerodynamic performance such as thrust, lift and drag coefficients is obtained. Then, according to the oncoming flow velocity and thrust of each sample point, the front and rear interface momentum source distributions are calculated using the front and rear interface velocity distributions obtained by the method of the present invention, and the lift and drag coefficients are calculated using a flow field simulation method for the distributed ducted fan-wing coupling layout in Chinese Patent CN117371344A, and compared with the results of the MRF method.

[0284] As Figure 19As shown, within the research scope, at different oncoming flow velocities, the average relative error of the lift coefficient is 2.61%, the error range for 25 sample points is from 1.13% to 3.72%, and the standard deviation is 0.77%; the average relative error of the drag coefficient is 2.44%, the error range is from 0.75% to 3.94%, and the standard deviation is 0.87%; the average relative error of the pitching moment coefficient is 1.01%, the error range is within 2.88%, and the standard deviation is 0.7%. Data analysis shows that the method of calculating the momentum source distribution at the front and rear interfaces using the velocity distributions at the front and rear interfaces obtained by the method of the present invention still exhibits high calculation accuracy under different oncoming flow velocity conditions;

[0285] Based on the Latin hypercube method, 25 sample points are randomly selected within the design space of 9000 RPM ≤ ω ≤ 12000 RPM. The MRF method is used to calculate the flow field of a single-ducted fan-wing segment fusion configuration at a 4° angle of attack under different fan speeds to obtain aerodynamic performances such as thrust and lift-drag coefficients. Then, based on the oncoming flow velocity and thrust at each sample point, the momentum source distribution at the front and rear interfaces is calculated using the velocity distributions at the front and rear interfaces obtained by the method of the present invention, and the lift-drag coefficients are calculated using a flow field simulation method for a distributed ducted fan-wing coupling layout in Chinese Patent CN117371344A, and compared with the results of the MRF method.

[0286] As Figure 20 shown, within the research scope, at different fan speeds, the average relative error of the lift coefficient is 2.93%, the error range for 25 sample points is from 1.71% to 4.06%, and the standard deviation is 0.63%; the average relative error of the drag coefficient is 1.9%, the error range is from 1.6% to 2.2%, and the standard deviation is 1.15%; the average relative error of the pitching moment coefficient is 0.72%, the error range is within 1.52%, and the standard deviation is 0.41%. Data analysis shows that the method of calculating the momentum source distribution at the front and rear interfaces using the velocity distributions at the front and rear interfaces obtained by the method of the present invention still exhibits high calculation accuracy under different fan speed conditions.

[0287] In summary, in order to solve the problem of the decrease in calculation efficiency caused by the new momentum source method relying on MRF flow field information, an interface flow field information generation method applicable to the ducted fan momentum source method proposed by the present invention can improve the prediction efficiency of the aerodynamic performance of a distributed ducted fan-wing fusion configuration. After verification under different oncoming flow velocities and different fan speed conditions, the method of the present invention improves the prediction efficiency of the aerodynamic performance while ensuring high calculation accuracy.

[0288] Although the embodiments of the present invention have been shown and described above, it can be understood that the above embodiments are exemplary and should not be construed as limiting the present invention. Those of ordinary skill in the art can make changes, modifications, substitutions, and variations to the above embodiments within the scope of the present invention without departing from the principles and spirit of the present invention.

Claims

1. A method for generating interface flow field information applicable to a ducted fan momentum source method, characterized in that: The following steps are involved: Step S1: the interface in front of the ducted fan rotor blade is recorded as the front interface, and the front interface is divided into a front hub boundary layer area, a front duct boundary layer area, a wing boundary layer area, a low-energy airflow area, a duct lip boundary layer area, a duct lip acceleration area, and a front uniform transition area; then, velocity distributions are generated respectively according to the above-divided areas, and then assembled and merged to obtain a velocity distribution on the front interface; Step S2: The rear interface of the ducted fan stator blade is recorded as the rear interface, and the rear interface is divided into a rear hub boundary layer area, a ducted index distribution transition area and a stator blade wake area; then, velocity distributions are generated respectively according to the above-divided areas, and then assembled and merged to obtain the velocity distribution on the rear interface.

2. The method for generating interface flow field information applicable to the ducted fan momentum source method according to claim 1, characterized in that: In step S1, the front hub boundary layer area refers to the annular area from the hub wall surface on the front interface to the outer boundary of the hub boundary layer; the front duct boundary layer area refers to the annular area from the duct inner wall surface on the front interface to the outer boundary of the duct boundary layer; the wing boundary layer area refers to the rectangular area from the upper surface of the wing to the outer boundary of the wing boundary layer; the low-energy airflow area refers to the area from the outer boundary of the wing boundary layer to the average flow velocity position of the interface; the duct lip boundary layer area refers to the area from the duct lip wall surface to the outer boundary of the duct lip boundary layer; the duct lip acceleration area refers to the area obtained by the intersection of the circumferential action range area and the radial action range area, and the circumferential action range area is: respectively is the center angle, and each center angle The range is taken as the circumferential action range, and the radial action range is: the intersection of each central angle ray and the outer ring of the duct is taken as the radial center position, and the 1 / 10 to 1 / 5 duct radius range around each radial center position is taken as the radial action range; the front uniform transition zone is the remaining area on the front interface after removing the above-mentioned partition.

3. The method for generating interface flow field information applicable to the ducted fan momentum source method according to claim 1, characterized in that: The step S1 comprises the following sub-steps: Step S11: calculating the average velocity at the front interface under different incoming flow velocities and different thrust conditions based on momentum theory; Step S12: Based on the thickness formula of the laminar boundary layer around the downstream flat plate and the assumption of power distribution of velocity, the thickness and two-dimensional velocity distribution of the front hub boundary layer are calculated; based on the thickness formula of the turbulent boundary layer around the downstream flat plate and the assumption of power distribution of velocity, the thickness and two-dimensional velocity distribution of the front duct boundary layer are calculated; based on the thickness and two-dimensional velocity distribution of the front hub boundary layer, the thickness and two-dimensional velocity distribution of the front duct boundary layer are interpolated to obtain the two-dimensional velocity distribution of the front uniform transition zone; Step S13: assembling the two-dimensional velocity distribution of the front hub boundary layer area, the two-dimensional velocity distribution of the front uniform transition area, and the two-dimensional velocity of the front duct boundary layer area obtained in step S12 according to the spatial position to obtain the two-dimensional radial velocity distribution on the front interface, and rotating the obtained two-dimensional radial velocity distribution around the fan axis in the circumferential direction for one circle to generate the three-dimensional velocity distribution "base" on the front interface; Step S14: using the Karman-Bohausen method to calculate the thickness and two-dimensional velocity distribution of the wing boundary layer region, using xfoil to calculate the thickness of the low-energy airflow region, and using cubic spline distribution to calculate the two-dimensional velocity distribution of the low-energy airflow region based on the thickness and two-dimensional velocity distribution of the wing boundary layer region and the thickness of the low-energy airflow region; assembling the two-dimensional velocity distribution of the wing boundary layer region and the two-dimensional velocity distribution of the low-energy airflow region and performing spanwise expansion to obtain the three-dimensional velocity distribution of the wing boundary layer region and the low-energy airflow region, and performing "bow-shaped" processing on the obtained three-dimensional velocity distribution; The specific process of "bow-forming" the obtained three-dimensional velocity distribution is as follows: The calculation method of the z coordinate in the direction of the bow height is: The coordinate system is defined as follows: the X axis is along the free flow direction, the Y axis is to the right from the pilot's perspective, and the Z direction is perpendicular to the XY plane and upward, satisfying the right-hand system; the calculation object is the velocity distribution scattered points at different spanwise positions in the three-dimensional rectangular velocity distribution, where z j_i_moi represents the corrected coordinates of the i-th point at the j-th spanwise position, z j_i_ori represents the initial coordinates of the i-th point at the j-th spanwise position, h j represents the z-height of the jth spanwise position in the arch region, δ wing_lp is the thickness of the low-energy airflow zone; Step S15: Assemble the velocity distributions of the wing boundary layer region and the low-energy airflow region after the "bow-forming" treatment obtained in step S14 with the three-dimensional velocity distribution "base" on the front interface obtained in step S13, and perform smoothing treatment on the area where the wing boundary layer region and the low-energy airflow region intersect with the "base" to obtain the velocity distribution of the front interface flow field that takes into account the hub boundary layer, the duct boundary layer, the wing boundary layer and the low-energy airflow region; Step S16: Based on the velocity distribution of the front interface flow field obtained in step S15, a linear function scaling method is used to correct the velocity distribution of the duct lip boundary layer area and the duct lip acceleration area to obtain the final velocity distribution of the front interface flow field.

4. The method for generating interface flow field information applicable to the ducted fan momentum source method according to claim 3, characterized in that: The assembly method for obtaining the two-dimensional radial velocity distribution on the front interface in step S13 is to merge the two-dimensional matrices containing coordinate information and velocity information, as shown in the following formula: In the formula, z f_hub is the velocity distribution u in the boundary layer area with the front hub f_hub The corresponding position coordinate, z gd is the velocity distribution u in the uniform transition zone before gd The corresponding position coordinate, z f_duct is the velocity distribution u in the boundary layer area of ​​the front duct f_duct The corresponding position coordinates; z 2D represents the coordinate distribution after assembling the front hub boundary layer area, the front uniform transition area, and the front duct boundary layer area, u 2D It represents the velocity distribution after assembling the front hub boundary layer area, the front uniform transition area, and the front duct boundary layer area.

5. The method for generating interface flow field information applicable to the ducted fan momentum source method according to claim 3, characterized in that: The specific process of assembling and smoothing the velocity distribution of the wing boundary layer region and the low-energy airflow region after the "bow-shaped" treatment with the three-dimensional velocity distribution "base" on the front interface in step S15 is as follows: The area in the "base" of the three-dimensional velocity distribution that is lower than the top outer edge of the wing boundary layer area and the low-energy airflow area is deleted, and the velocity distribution of the wing boundary layer area and the low-energy airflow area after "bow" processing is used to replace it. For the step of velocity distribution at the junction position, the linear interpolation method is used to smooth the velocity step area, and the velocity distribution of the front interface flow field that takes into account the hub boundary layer, duct boundary layer, wing boundary layer and low-energy airflow area is obtained.

6. The method for generating interface flow field information applicable to the ducted fan momentum source method according to claim 3, characterized in that: The specific process of step S16 is as follows: ① Establish a linear scaling factor function. The linear scaling function attenuation in the boundary layer area and acceleration area of ​​the duct lip mainly includes radial attenuation and circumferential attenuation. The formulas are: Radial scaling factor r factor : In the formula, r i represents the radius of the current grid point, r represents the duct radius, σ r Indicates the parameters set to control the radial zoom range; Circumferential scaling factor θ factor : In the formula, θ i -θ center Represents the circumferential angle difference between the current grid point and the boundary layer center of the duct lip, σ θ Indicates the parameters for setting the control circumferential scaling range; ② Multiply the front interface flow field velocity distribution obtained in step S15 by the corresponding scaling function to obtain the front interface flow field velocity distribution after correction of the duct lip boundary layer area and the duct lip acceleration area; specifically: For the duct lip boundary layer area, which is the area where the velocity decreases, the product of the circumferential and radial scaling factors is used as the total scaling factor sc of the duct lip boundary layer area. factor_bl : sc factor_bl =r factor ·θ factor Based on the total scaling factor of the boundary layer area at the duct lip, the following expression is used to perform the first scaling correction on the velocity on the front interface: u isc1 =u i ·(1-sc factor_bl ) Expands to: In the formula, u i is the flow field velocity at each point on the front interface, u isc1 The velocity of the flow field at each point on the front interface after the first scaling correction of the velocity on the front interface; For the duct lip acceleration zone, which is the area where the speed increases, the following expression is used as the total scaling factor sc of the duct lip acceleration zone factor_ac : sc factor_ac =1+(ms-1)·r factor ·θ factor In the formula, ms is the maximum value of velocity scaling. The curve at the maximum curvature of the duct lip can be used as the leading edge. After being supplemented into a complete airfoil, the potential flow method is used to obtain the airflow velocity ratio at the relative position of the front interface. The empirical value is 1.2 to 1.

4. Based on the total scaling factor of the acceleration zone of the duct lip, the velocity on the front interface is corrected for the second time using the following expression: u isc2 =u isc1 ·sc factor_ac Expands to: In the formula, u isc2 The flow field velocity at each point on the front interface after the second scaling correction of the velocity on the front interface is the final flow field velocity distribution on the front interface.

7. The method for generating interface flow field information applicable to the ducted fan momentum source method according to claim 1, characterized in that: In step S2, the rear hub boundary layer area refers to the annular area from the hub wall surface to the outer boundary of the hub boundary layer on the rear interface; the stator blade wake area refers to the central contour of the stator blade wake area obtained by projecting the shape of the stator blade trailing edge onto the rear interface, and then expanding the central contour of the stator blade wake area along the circumferential direction. The area within the angle range; the duct index distribution transition zone refers to the remaining area on the rear interface after removing the above-mentioned partitions.

8. The method for generating interface flow field information applicable to the ducted fan momentum source method according to claim 1, characterized in that: The step S2 comprises the following sub-steps: Step S21: calculating the average velocity at the rear interface under different incoming flow velocities and different thrust conditions based on momentum theory; Step S22: Based on the thickness formula of the turbulent boundary layer around the downstream flat plate and the assumption of velocity power distribution, the thickness and two-dimensional velocity distribution of the rear hub boundary layer area are calculated, and based on the assumption of the velocity distribution law in the circular tube, the two-dimensional velocity distribution of the duct index distribution transition zone is calculated; Step S23: assembling the two-dimensional velocity distribution of the rear hub boundary layer region and the two-dimensional velocity distribution of the duct index distribution transition region obtained in step S22 according to the spatial position to obtain the two-dimensional radial velocity distribution on the rear interface, and rotating the obtained two-dimensional radial velocity distribution around the fan axis in the circumferential direction to generate a three-dimensional velocity distribution "base" on the rear interface; Step S24: Based on the shape of the stator blade trailing edge, generate the center contour of the stator blade wake area. Specifically, the shape of the stator blade trailing edge is projected onto the rear interface to obtain the center line of the stator blade wake. On the rear interface, the two ends of the stator blade wake center line are extended to the hub of the rear interface and the duct of the rear interface as the center contour of the stator blade wake area. Based on the center contour of the stator blade wake area, obtain the circumferential angle θ of the center contour of the stator blade wake area center ; Step S25: Based on the three-dimensional velocity distribution "base" on the rear interface obtained in step S23 and the central contour of the stator wake area obtained in step S24, the velocity distribution in the transition zone of the duct index distribution is corrected by a linear function scaling method to obtain the final velocity distribution of the rear interface flow field.

9. The method for generating interface flow field information applicable to the ducted fan momentum source method according to claim 8, characterized in that: The assembly method of obtaining the two-dimensional radial velocity distribution on the rear interface in step S23 is to merge the two-dimensional matrices containing coordinate information and velocity information, as shown in the following formula: In the formula, z b_2D represents the coordinate distribution after assembling the rear hub boundary layer area and the ductance index distribution transition area, u b_2D It represents the velocity distribution after assembling the rear hub boundary layer area and the duct index distribution transition area.

10. The method for generating interface flow field information applicable to the ducted fan momentum source method according to claim 8, characterized in that: The specific process of step S25 is: Wake circumferential scaling factor θ wj_factor for: In the formula, θ i -θ center represents the circumferential angle difference between the current grid point and the central contour of the stator wake area, σ θ Indicates the parameters for setting the control circumferential scaling range; Based on the wake circumferential scaling factor θ wj_factor , the following expression is used to scale and correct the velocity on the rear interface: in isc_b =in i_b ·(1-θ wj_factor ) In the formula, u i_b is the flow field velocity at each grid point on the "base" distribution of the rear interface, u isc_b The flow field velocity at each grid point on the rear interface after scaling and correcting the velocity on the rear interface is the final flow field velocity distribution on the rear interface.

Images