# Influence of Roughness on the Asymmetric Flow Field of a Slender Body

Authors: Kai Wei; Shaosong Chen; Yihang Xu; Dongdong Tang
Journal: International Journal of Aerospace Engineering
 (2022)

---

## Abstract

Low speed and high angle of attack are problems that must be faced in vertical launching missiles. A natural asymmetric vortex phenomenon occurs at a low speed and high angle of attack of a slender body. In this paper, the Detached Eddy Simulation (DES) method is adopted to simulate the asymmetric flow of a slender body at a high angle of attack. The influence of roughness is analysed from the flow field and pressure distribution. And the DES method is verified by comparing it with the wind tunnel test results. The flow fields of the four models with the surface roughness of 0, 0.8, 5, and 100 are compared. The downstream flow still appears asymmetric when the surface is smooth, indicating that the stability of the flow field is not enough to maintain the symmetry of the flow field. The simulation results show that the adverse pressure gradient increases in the region where the boundary layer separates with the roughness greatly increasing, but the structure of the flow field at the head is slightly different. On the whole, the surface pressure and lateral force per unit length (Cfzpm) transform alternately along the axial direction, and the period of alternating reverses increases with the increase of roughness. Finally, the pressure tends to be in equilibrium, and Cfzpm approaches zero. It should be noted that the distribution of Cfzpm is slightly different on the head, indicating that the asymmetry of pressure on the cylinder section is an important factor controlling the magnitude and the direction of lateral force. The influence of roughness on the flow around a slender body is acquired in this paper, and it has reference significance to the roughness problem of the actual missile.


---

## Body

## 1. Introduction

Low speed and high angle of attack are problems that vertical launching missiles must face. At low speed and high angle of attack, a completely symmetrical slender body generates an asymmetric vortex flow at a high angle of attack, and the lateral force generated by the asymmetric vortex leads to lateral instability. However, the reason for the asymmetry flow around the slender body at a high angle of attack is still unresolved. Siclari and Marconi [1] confirmed that an asymmetric flow exists in reality, rather than being caused by numerical simulation errors. The causes of an asymmetric flow field with a high angle of attack around a slender body include hydrodynamic instability [2] and asymmetry of flow separation [3]. Jiménez-Varona et al. [4] analysed the vorticity distribution of a slender body from the aspect of hydrodynamic stability by comparing numerical simulations and experiments and concluded that the entire flow field was unstable. In terms of the asymmetric study of flow separation, Obeid et al. [5] captured that a slight asymmetry of the separation point would lead to the inconsistency of shear flow on both sides, which in turn would cause the asymmetry of the separation vortex. Furthermore, Taligoski et al. [6] explained that the boundary layer affects the asymmetry of shear flow and leads to the asymmetric phenomenon of the leeward vortex; they also investigated the sensitivity of the boundary layer and overall lateral aerodynamic force to the rolling attitude when there were small machining defects on the head.The transition and separation of the boundary layer affect the asymmetry of the shear flow, whereas roughness is an important factor affecting the transition and separation of the boundary layer. Jeong et al. [7] conducted experiments on the influence of roughness on the boundary layer in a low-speed wind tunnel and captured that roughness would promote the transition, and the rough surface at the back of the transition region would enhance the turbulent mixing of a low-momentum flow near the wall and a high-momentum flow outside the boundary layer, resulting in momentum loss. Noureddine et al. [8] simulate the cone-cylinder flow with a rough surface using the transition Shear Stress Transfer (transition SST) turbulence model. It was concluded that the existence of surface roughness would promote the transition of the boundary layer and reduce the lateral force. In terms of numerical simulation, Bimbato et al. [9] developed a roughness model in two dimensions by combining the Lagrange vortex method with a large eddy simulation (LES) to analyse the influence of roughness.On a smaller scale, Karn et al. [10] used electron microscopy to confirm the asymmetry of the flow induced by small machining defects on the head. Moskovitz et al. [11] also proposed that the influence of roughness on the sharp head controls the asymmetry. A significant number of experiments [12–14] and numerical simulations [15–18] have confirmed that the initial asymmetric disturbance located at the upstream tip of the nose triggers strong asymmetry in the downstream flow field. In addition, surface defects and roughness cause secondary vortices in the shear layer [19]. Overall, roughness has a significant effect on the flow patterns around a cylinder, which was confirmed by Buresti [20]. The asymmetrical pressure distribution of the cone head was also experimentally and numerically determined by Meng et al. [21], Jia et al. [22], and Zheng et al. [23]. In the control of asymmetric flow, head spoiler [24], nose cone bluntness [25], and plasma generator [26] all play a role in reducing the asymmetric flow at high angles of attack.Current research on the effect of roughness on an asymmetric flow around slender bodies mostly focuses on the causes of asymmetric flow; there is a lack of analysis of the effect of roughness on the overall flow field structure. In this study, the effect of roughness on the flow field structure on the leeward side of a slender body with a high angle of attack was analysed from the flow field structure of the entire surface and near-wall surface. To study the effect of roughness on the asymmetric flow field near a slender body at a high angle of attack, a numerical simulation method based on the Detached Eddy Simulation (DES) method was used to simulate the flow near the slender body at a high angle of attack, and the simulation results were compared with those of Particle Image Velocimetry (PIV) experiments. Numerical simulations were used to analyse the flow field structure in the leeward zone, the separation characteristics of the two sides of the head, and the surface pressure distribution in the leeward zone of the entire missile under four types of roughness.


## 2. Model and Experiment

### 2.1. Model and PIV Experimental Conditions

As shown in Figure1, the experimental model lathed with hard alloy was a body of revolution with a large slenderness ratio. The head of the model was blunt with a diameter of d=2 mm. The model’s diameter is D=43.2 mm, and the length is L=531.2 mm. The definition of roughness is shown in Figure 2. Taking the average height of the surface of the sampling length as the reference line, the area enclosed between the surface curve below the reference line and the reference line is negative Sblack, and the area enclosed between the surface curve above the reference line and the reference line is positive Sblack. The sum of the positive and negative areas in the sampling area is zero. And there is a straight line at Ra above the reference line. The area Sred enclosed by the Ra line and the reference line in the sampling area is equal to the sum of the absolute values of the positive and negative areas Sblack in the sampling area. That is, ∑x=1nSblackx=0and∑x=1nSblackx=Sred. Ra is defined as roughness. The surface roughness of the model was Ra=0.8 μm, which was compared with that of the machined roughness contrast block. The surface roughness of the model was Ra=0.8 μm. The PIV experiment was conducted under the following experimental conditions: room temperature, 25°C; humidity, 56%; pressure, 100293 Pa; wind speed v=3.4 m/s, and model angle of attack α=55°. The selection of the angle of attack takes into account the conclusion of literature [27]: the maximum lateral force is obtained at the angle of attack of 55°, which is convenient for further analysis. The experimental Reynolds number was calculated by Re=ρvD/μ and Re=9240. The right-hand rule was used to establish a Cartesian coordinate system, with the X-axis pointing from the tail of the body along the axis of the body and the Y-axis perpendicular to the X-axis vertically upward.Figure 1
Schematic of experimental model and partially enlarged drawing of head.Figure 2
The definition ofRa.


### 2.2. PIV Experimental Case

The PIV experiment was conducted in the HG-1 low-speed wind tunnel of NJUST. The wind tunnel is an open-jet closed-circuit low-speed wind tunnel with a turbulence of less than 1%. The layout scheme of the experimental system is shown in Figure3.Figure 3
Schematic of PIV experiment.The PIV experiment uses a 532 nm solid pulse laser at 10 Hz as the light source, and the acquisition frequency is also 10 Hz. In the experiment, the imaging particles were scattered in the wind tunnel. After the particles were evenly distributed in the experiment section, the CCD camera and the laser were exposed twice at an interval of 200μs to obtain two frames of particle images with 2048∗2048 pixels. Two frames of images are divided into multiple query areas in the image processing software, and cross-correlation statistical analysis is used to obtain the size and direction of particle displacement in the query area. The time interval of two frames of images is the set interval time of the pulse, and the particle velocity vector can be obtained. The velocity vector of the query area can be calculated by the statistical average of all particles in the query area, and the velocity vector field of the whole solution area can be calculated by the above determination and statistics of all the query areas. In the experiment, an impulse disturbance was applied to shift the asymmetric flow of the head eddy to a constant direction. After the flow field stabilized, the cross-correlation algorithm was used to obtain the two-dimensional velocity distribution diagram of the four cross-section positions within 5 s, and the vorticity was further calculated. The magnitudes of the vorticity contours of the four sections obtained from the PIV experiment were compared with those from the numerical simulation results. The magnitude of vorticity indicates the strength of the vortex at this location. The positive direction of the X-axis is the vorticity ωx (X-vorticity), and its value is obtained as follows:
(1)ωx=∂Uz∂y−∂Uy∂z.


### 2.3. Force Measurement

For quantitative comparison, force measurement experiments in the low-speed wind tunnel were carried out and used the TP-18B strain balance measurement system. The measuring bridge of the balance was a double matrix type. The performance indexes of balance calibration are shown in Table1.Table 1
The performance indexes of balance calibration.
Normal forcePitching momentAxial forceRoll momentLateral forceYaw momentDesign load (N, N·m)603.0201.2301.5Absolute error (N, N·m)0.060.0030.030.00240.030.0015Limiting error (%)0.30.30.450.450.30.3Comprehensive precision (%)0.030.030.060.100.080.07Under the condition of 3.4 m/s inflow, the aerodynamic force of the experimental model is on the same order of magnitude as the absolute error of the balance. Therefore, it is considered to compare the force measurement results under 15.3 m/s flow with DES simulation results under the same conditions to verify the credibility of the numerical algorithm. Supplementary experiments were carried out in the HG-1 wind tunnel, and the wind speed was 15.3 m/s. Wind speed was measured using the complete set of hot-wire anemometer systems of Dantec Company. The calibration error of hot-wire anemometer was no more than 2%, and the measured turbulence of inlet flow was 1.4%. Ambient relative humidity was 59% when the experiment was carried out. The experimental model is fixed on the base of the balance with an adjustable angle of sideslip and angle of attack. The base’s sideslip angle adjustment accuracy is 0.025°, and the angle of attack adjustment accuracy is 0.1°. The experiment was carried out under the conditions of a sideslip angle of 55° and angle of attack of 0°. Before the experiment, the attitude of zero angle of attack with zero lift force and zero lateral force was corrected. During the experiment, 10000 samples were taken within 10 s, and the sampling results were averaged arithmetically to obtain the experimental results. The experimental results were compared with the numerical simulation values after the coordinate system conversion.The average value of aerodynamic force is obtained after the calculation of balance measurement value, coordinate system conversion, and arithmetic average, as shown in Table2.Table 2
Average aerodynamic force.
Axial forceNormal forceLateral forcePitching momentRoll momentYaw momentEXP result (N, N·m)-0.02122.50300.4584-0.66830.0000-0.1047


### 2.4. Numerical Simulation

A computational fluid dynamic (CFD) method was used to simulate the flow of slender bodies at high angles of attack to study the influence of different roughness values on the flow field structure. The flow around models of four kinds of roughnessRa=0, 0.8, 5, and 100 μm was calculated under the condition of the angle of attack α=55°, as shown in Table 3. The results of roughness Ra=0.8 μm were compared with the experimental results to verify the feasibility of the numerical simulation method.Table 3
Models of 4 kinds of roughness.
Model 1Model 2Model 3Model 4RoughnessRa=0 μmRa=0.8 μmRa=5 μmRa=100 μm


### 2.5. Numerical Simulation of Turbulence Model

The DES method based onk−ωSST was adopted for numerical simulation [28, 29]. The DES method is a turbulent solution for 3D unsteady and large separated flows. The Reynolds-averaged Navier-Stokes (RANS) model was used in the near-wall region, where the grid resolution was not suitable for LES. The LES was used in the region away from the wall where the grid resolution was not suitable for RANS, thus ensuring the calculation accuracy, significantly reducing the grid requirements, and saving computing resources.In the Navier-Stokes equation (N-S equation), the velocityui can be decomposed into the average quantity ui¯ and fluctuation quantity ui′ as follows:
(2)ui=ui¯+ui′.The continuity and momentum equations are expressed as follows:(3)∂ρ∂t+∂ρuj∂xj=0,∂ρui∂t+∂ρuiuj∂xj=−∂p∂xi+∂∂xjμ+μt∂ui∂xj+∂uj∂xi−23∂ui∂xjδij.The initial SST model uses thek−ω model in the near-wall area, along with the k−ε model in the boundary layer edge and free shear layer, which are transited by the mixing function. The two-equation turbulence model of k−ω SST can be expressed as follows:
(4)∂ρk∂t+∂ρuik∂xi=Pk~−β∗ρkω+∂∂xiμ+σkμt∂k∂xik,∂ρω∂t+∂ρuiω∂xi=αρS2−βρω2+∂∂xiμ+σωμt∂k∂xj+21−F1ρσω21ω∂k∂xi∂ω∂xi,where F1 denotes the following first mixing function:
(5)F1=tanhminmaxkβ∗ωy,500vy2ω,4ρσω2kCDkωy24,where
(6)CDkω=max2ρσω21ω∂k∂xi∂ω∂xi,10−10,where y denotes the distance from the nearest wall. The function F1 represents the k−ε model at the wall and k−ω is the model in the boundary layer. Turbulence vortex viscosity is obtained as follows:
(7)vt=a1kmaxa1ω,SF2,whereS denotes an invariant of the strain rate and F2 is the second mixing function obtained as follows:
(8)F2=tanhmax2kβ∗ωy,500vy2ω2.A production limiter exists to prevent the accumulation of turbulence in the stagnation zone:(9)Pk~=minPk,10⋅β∗ρkω,Pk=μt∂ui∂xj∂ui∂xj+∂uj∂xi,α=α1F+α21−F.All constants were calculated by mixing the constants of thek−ε and k−ω models. The constants were α1=5/9, α2=0.44, β∗=0.09, β2=0.0828, σk1=0.85, σk2=1, σω1=0.5, and σω2=0.856.The DES method based on SST was used to modify the turbulent dissipation term in thek equation as follows:
(10)ρε=β∗ρkω⟶β∗ρkω⋅FDES.Correction term:(11)FDES=maxLtCDESΔ,1,where ε denotes the turbulence dissipation rate, the maximum local grid size is Δ=maxΔx,Δy,Δz in the Cartesian coordinate system, the turbulence characteristic length size is Lt=k/β∗ω, and constant CDES=0.61.


### 2.6. Numerical Method Verification

To verify the correctness of the numerical simulation algorithm, the DES method was used to simulate the smooth cylinder in literature [30], and the results were compared with the wind tunnel experiment results in the literature [30]. In the reference, the diameter of the cylinder is 57 mm, the incoming flow velocity is 3 m/s, and the Reynolds number is Ra=1.17×104. The experimental results obtained in the reference are the average value of the pressure coefficient distribution on one side of the cylinder; the pressure coefficient was defined as CP=P−P∞/1/2ρv2. A three-dimensional cylindrical model with the same diameter was used for numerical simulation. The incoming flow velocity was set at 3 m/s, the incoming flow pressure was set at 101325 Pa, the Reynolds number was calculated at 11706, and the time step of unsteady flow was set as 2E−5 s. The numerical simulation results were compared with the experimental values in literature by taking the mean pressure distribution within 0.17 s, that is, the mean pressure distribution within three cycles. The results are shown in Figure 4, which is in good agreement with the experimental measurement results on the whole cylinder. Therefore, the numerical method is reliable.Figure 4
Comparison diagram of numerical simulation results and experimental values of smooth cylinder mean pressure coefficient distribution at 3 m/s.Furthermore, the flow around a cylinder with a roughness of 500μm in the literature [30] was simulated by the DES method to verify the roughness case. The calculation was carried out under standard atmospheric conditions of air pressure 101325 kPa and air temperature 288.15 K. And the flow with incoming velocity 15.2 m/s and Reynolds number 59326 is simulated. After averaging the pressure distribution within 0.24 s, the pressure distribution of the rough cylinder is shown in Figure 5. The simulation results are in good agreement with the experimental values.Figure 5
Comparison diagram of numerical simulation results and experimental values of rough cylinder mean pressure coefficient distribution at 15 m/s.


### 2.7. Solution Method and Grid Verification

Because the study involved a low-speed incompressible unsteady flow, a pressure-based unsteady solver was selected, and the coupling algorithm was used for the pressure-velocity coupling solution. A second-order implicit scheme was used for the time discretization method. For the spatial discretization method, the least-squares cell-based scheme was used for gradient discretization, the second-order scheme was used for pressure discretization, and the second-order upwind scheme was used for density, momentum, turbulent kinetic energy, turbulent dissipation rate, and energy discretization. In addition, higher-order term relaxation was used to enhance the convergence performance of the calculation. The calculation domain was set according to a previous study [27] and slightly expanded. As shown in Figures 6 and 7, the computational domain is a cylinder with a diameter of 1600 mm. The head vertex is approximately 23D away from the forward field. The height of the first layer of the grid near the wall was approximately 0.000463D. Wall y plus y+=1. In the experiment, the model is fixed on the datum platform through the support rod and bolt. Because the support structure is inevitable, the measurement will be affected by the support structure. Therefore, in the numerical simulation, the influence of the support structure is considered, and the support rod at the back of the model is included in the calculation model, which is reflected in the schematic diagram of the calculation domain (Figures 6 and 7). The surface mesh which is evenly distributed along the circumference is shown in Figure 8. The coordinate system used in the numerical simulation was different from the experimental coordinate system. The Cartesian coordinate system used in the numerical simulation was established according to the right-hand rule. The X-axis points from the head to the tail along the missile axis, the Y-axis is vertical, and the X-axis is vertical and upward.Figure 6
Diagram of boundary conditions.Figure 7
Grid diagram.Figure 8
A view of the surface mesh and a partial enlargement.Considering the randomness of the initial asymmetric flow direction, the full-multigrid initialization method is used in the numerical simulation to make the initial conditions produce the same asymmetric flow direction as the experiment. In this way, the initial disturbance of the incoming flow turbulence is not enough to determine the offset direction of the asymmetric flow field. Due to the existence of a bistable state, the calculation will converge to a near steady state, that is, one side of the initial offset. If the completely symmetric initialization method is adopted, the initial disturbance of incoming turbulence will play a decisive role in the direction of asymmetric back vortex migration of the flow field and thus will converge randomly to one of the two steady states.Grid independence verification is required for numerical simulation. Grid files of Grid 1, Grid 2, and Grid 3 with different grid quantities are selected for grid independence verification, and the uniform surface roughness isRa=0.8 μm. The numerical simulation results of grid number and time-average aerodynamic coefficient of the three grids are shown in Table 4. The elevation direction of the pitching moment coefficient is defined as positive; the reference length is 0.0432 mm, and the reference area is 0.001464 m2. As can be seen from Table 4, there is little difference between the simulation results of Grid 2 and Grid 3, while there is a great difference between the simulation results of Grid 1 and the other two sets of grids. It can be seen that the continuous increase of the grid quantity has little influence on the simulation results.Table 4
Calculation results of aerodynamic coefficients for three kinds of grids (no bottom resistance is included).
Grid quantityAxial force coefficientNormal force coefficientPitching moment coefficientLateral force coefficientYaw moment coefficientRoll moment coefficientGrid 13596092-0.110213.0865-95.51991.1024-16.35720.0008Grid 25059776-0.204714.8344-98.3868-3.518712.96500.0046Grid 36436364-0.200814.8987-99.0282-3.556512.77930.0021Furthermore, the average vorticity contours of the typical cross-section of asymmetric flow were selected and compared with the PIV experimental results. As shown in Figures9–12, it can be seen that the simulated results of vorticity distribution of Grid 2 and Grid 3 are close to the experimental results of PIV, showing the same asymmetry as the experimental results, while the simulated results of Grid 1 differ greatly from the experimental values, showing no obvious asymmetry. To save computing resources, the grid file with 5059776 grid quantity was selected for the next calculation.Figure 9
X-vorticity contour comparison at section 1: (a) PIV experiment, (b) result of Grid 1, (c) result of Grid 2, and (d) result of Grid 3.
(a)(b)(c)(d)Figure 10
X-vorticity contour comparison at section 2: (a) PIV experiment, (b) result of Grid 1, (c) result of Grid 2, and (d) result of Grid 3.
(a)(b)(c)(d)Figure 11
X-vorticity contour comparison at section 3 (a) PIV experiment, (b) result of Grid 1, (c) result of Grid 2, and (d) result of Grid 3.
(a)(b)(c)(d)Figure 12
X-vorticity contour comparison at section 4: (a) PIV experiment, (b) result of Grid 1, (c) result of Grid 2, and (d) result of Grid 3.
(a)(b)(c)(d)To verify the unsteady performance, delta time verification is required. Grid 2 is used for delta time verification, and three delta times,1E−4 s, 2E−5 s, and 5E−6 s are selected for delta time verification. The results are shown in Table 5. It can be seen that when the delta time is smaller than 2E−5 s, the aerodynamic coefficient does not change much. Therefore, the delta time 2E−5 s is selected for calculation.Table 5
Calculation results of aerodynamic coefficients for three kinds of delta time (no bottom resistance is included).
Time stepAxial force coefficientNormal force coefficientPitching moment coefficientLateral force coefficientYaw moment coefficientRoll moment coefficient1E−4 s-0.219414.6317-96.2298-3.702314.48900.00662E−5 s-0.204714.8344-98.3868-3.518712.96500.00465E−6 s-0.205314.8577-98.4543-3.517312.94010.0044In order to verify the unsteady characteristics of the flow around the slender body, a solid-state continuous laser was used to replace the pulse laser in the experiment, and a high-speed camera was used to replace the CCD camera to obtain the smoke flow diagram of the instantaneous flow. The frame rate of the high-speed camera was set as 1000 frames/s, and the aperture was f1.4. As shown in Figure13, the images obtained by the high-speed camera show that the structure of the flow field is steady within 2 s under experimental conditions. As shown in Figure 14, both the smoke flow diagram obtained in the experiment and the numerical simulation results show that there is no obvious difference in the flow field structure. Therefore, the time-average method is adopted in the subsequent analysis to compare the time-average flow field within 0.01 s.Figure 13
Flow structure is steady with time increase.Figure 14
There is no obvious difference in flow field structure.


### 2.8. Comparison of Aerodynamic Coefficients

The DES method was used to calculate the flow with velocity at 15.3 m/s, pressure at 102.56 kPa, and angle of attack at 55°. The surface roughness of the model was set to 0.8μm. The experimental results are compared with the DES simulation results after the coordinate system transformation, as shown in Table 6. The reference area is 1.4657E−3 m2, and the reference length is 0.0432 m.Table 6
The comparison between DES simulation results and experimental results after coordinate transformation.
Axial force coefficientNormal force coefficientLateral force coefficientPitching moment coefficientRoll moment coefficientYaw moment coefficientDES result-0.108411.30612.4612-70.19990.0011-13.0665EXP result-0.103712.27262.2477-75.85460.0042-11.8878Relative error4.57%7.87%9.50%7.45%73.78%9.91%As shown in Table6, except for the rolling moment coefficient, the errors of other force coefficients are all within 10%, while the rolling moment coefficient is quite close to 0. The calculation result of the DES method is reliable. The above content is added to the solution model and grid validation section.


## 2.1. Model and PIV Experimental Conditions

As shown in Figure1, the experimental model lathed with hard alloy was a body of revolution with a large slenderness ratio. The head of the model was blunt with a diameter of d=2 mm. The model’s diameter is D=43.2 mm, and the length is L=531.2 mm. The definition of roughness is shown in Figure 2. Taking the average height of the surface of the sampling length as the reference line, the area enclosed between the surface curve below the reference line and the reference line is negative Sblack, and the area enclosed between the surface curve above the reference line and the reference line is positive Sblack. The sum of the positive and negative areas in the sampling area is zero. And there is a straight line at Ra above the reference line. The area Sred enclosed by the Ra line and the reference line in the sampling area is equal to the sum of the absolute values of the positive and negative areas Sblack in the sampling area. That is, ∑x=1nSblackx=0and∑x=1nSblackx=Sred. Ra is defined as roughness. The surface roughness of the model was Ra=0.8 μm, which was compared with that of the machined roughness contrast block. The surface roughness of the model was Ra=0.8 μm. The PIV experiment was conducted under the following experimental conditions: room temperature, 25°C; humidity, 56%; pressure, 100293 Pa; wind speed v=3.4 m/s, and model angle of attack α=55°. The selection of the angle of attack takes into account the conclusion of literature [27]: the maximum lateral force is obtained at the angle of attack of 55°, which is convenient for further analysis. The experimental Reynolds number was calculated by Re=ρvD/μ and Re=9240. The right-hand rule was used to establish a Cartesian coordinate system, with the X-axis pointing from the tail of the body along the axis of the body and the Y-axis perpendicular to the X-axis vertically upward.Figure 1
Schematic of experimental model and partially enlarged drawing of head.Figure 2
The definition ofRa.


## 2.2. PIV Experimental Case

The PIV experiment was conducted in the HG-1 low-speed wind tunnel of NJUST. The wind tunnel is an open-jet closed-circuit low-speed wind tunnel with a turbulence of less than 1%. The layout scheme of the experimental system is shown in Figure3.Figure 3
Schematic of PIV experiment.The PIV experiment uses a 532 nm solid pulse laser at 10 Hz as the light source, and the acquisition frequency is also 10 Hz. In the experiment, the imaging particles were scattered in the wind tunnel. After the particles were evenly distributed in the experiment section, the CCD camera and the laser were exposed twice at an interval of 200μs to obtain two frames of particle images with 2048∗2048 pixels. Two frames of images are divided into multiple query areas in the image processing software, and cross-correlation statistical analysis is used to obtain the size and direction of particle displacement in the query area. The time interval of two frames of images is the set interval time of the pulse, and the particle velocity vector can be obtained. The velocity vector of the query area can be calculated by the statistical average of all particles in the query area, and the velocity vector field of the whole solution area can be calculated by the above determination and statistics of all the query areas. In the experiment, an impulse disturbance was applied to shift the asymmetric flow of the head eddy to a constant direction. After the flow field stabilized, the cross-correlation algorithm was used to obtain the two-dimensional velocity distribution diagram of the four cross-section positions within 5 s, and the vorticity was further calculated. The magnitudes of the vorticity contours of the four sections obtained from the PIV experiment were compared with those from the numerical simulation results. The magnitude of vorticity indicates the strength of the vortex at this location. The positive direction of the X-axis is the vorticity ωx (X-vorticity), and its value is obtained as follows:
(1)ωx=∂Uz∂y−∂Uy∂z.


## 2.3. Force Measurement

For quantitative comparison, force measurement experiments in the low-speed wind tunnel were carried out and used the TP-18B strain balance measurement system. The measuring bridge of the balance was a double matrix type. The performance indexes of balance calibration are shown in Table1.Table 1
The performance indexes of balance calibration.
Normal forcePitching momentAxial forceRoll momentLateral forceYaw momentDesign load (N, N·m)603.0201.2301.5Absolute error (N, N·m)0.060.0030.030.00240.030.0015Limiting error (%)0.30.30.450.450.30.3Comprehensive precision (%)0.030.030.060.100.080.07Under the condition of 3.4 m/s inflow, the aerodynamic force of the experimental model is on the same order of magnitude as the absolute error of the balance. Therefore, it is considered to compare the force measurement results under 15.3 m/s flow with DES simulation results under the same conditions to verify the credibility of the numerical algorithm. Supplementary experiments were carried out in the HG-1 wind tunnel, and the wind speed was 15.3 m/s. Wind speed was measured using the complete set of hot-wire anemometer systems of Dantec Company. The calibration error of hot-wire anemometer was no more than 2%, and the measured turbulence of inlet flow was 1.4%. Ambient relative humidity was 59% when the experiment was carried out. The experimental model is fixed on the base of the balance with an adjustable angle of sideslip and angle of attack. The base’s sideslip angle adjustment accuracy is 0.025°, and the angle of attack adjustment accuracy is 0.1°. The experiment was carried out under the conditions of a sideslip angle of 55° and angle of attack of 0°. Before the experiment, the attitude of zero angle of attack with zero lift force and zero lateral force was corrected. During the experiment, 10000 samples were taken within 10 s, and the sampling results were averaged arithmetically to obtain the experimental results. The experimental results were compared with the numerical simulation values after the coordinate system conversion.The average value of aerodynamic force is obtained after the calculation of balance measurement value, coordinate system conversion, and arithmetic average, as shown in Table2.Table 2
Average aerodynamic force.
Axial forceNormal forceLateral forcePitching momentRoll momentYaw momentEXP result (N, N·m)-0.02122.50300.4584-0.66830.0000-0.1047


## 2.4. Numerical Simulation

A computational fluid dynamic (CFD) method was used to simulate the flow of slender bodies at high angles of attack to study the influence of different roughness values on the flow field structure. The flow around models of four kinds of roughnessRa=0, 0.8, 5, and 100 μm was calculated under the condition of the angle of attack α=55°, as shown in Table 3. The results of roughness Ra=0.8 μm were compared with the experimental results to verify the feasibility of the numerical simulation method.Table 3
Models of 4 kinds of roughness.
Model 1Model 2Model 3Model 4RoughnessRa=0 μmRa=0.8 μmRa=5 μmRa=100 μm


## 2.5. Numerical Simulation of Turbulence Model

The DES method based onk−ωSST was adopted for numerical simulation [28, 29]. The DES method is a turbulent solution for 3D unsteady and large separated flows. The Reynolds-averaged Navier-Stokes (RANS) model was used in the near-wall region, where the grid resolution was not suitable for LES. The LES was used in the region away from the wall where the grid resolution was not suitable for RANS, thus ensuring the calculation accuracy, significantly reducing the grid requirements, and saving computing resources.In the Navier-Stokes equation (N-S equation), the velocityui can be decomposed into the average quantity ui¯ and fluctuation quantity ui′ as follows:
(2)ui=ui¯+ui′.The continuity and momentum equations are expressed as follows:(3)∂ρ∂t+∂ρuj∂xj=0,∂ρui∂t+∂ρuiuj∂xj=−∂p∂xi+∂∂xjμ+μt∂ui∂xj+∂uj∂xi−23∂ui∂xjδij.The initial SST model uses thek−ω model in the near-wall area, along with the k−ε model in the boundary layer edge and free shear layer, which are transited by the mixing function. The two-equation turbulence model of k−ω SST can be expressed as follows:
(4)∂ρk∂t+∂ρuik∂xi=Pk~−β∗ρkω+∂∂xiμ+σkμt∂k∂xik,∂ρω∂t+∂ρuiω∂xi=αρS2−βρω2+∂∂xiμ+σωμt∂k∂xj+21−F1ρσω21ω∂k∂xi∂ω∂xi,where F1 denotes the following first mixing function:
(5)F1=tanhminmaxkβ∗ωy,500vy2ω,4ρσω2kCDkωy24,where
(6)CDkω=max2ρσω21ω∂k∂xi∂ω∂xi,10−10,where y denotes the distance from the nearest wall. The function F1 represents the k−ε model at the wall and k−ω is the model in the boundary layer. Turbulence vortex viscosity is obtained as follows:
(7)vt=a1kmaxa1ω,SF2,whereS denotes an invariant of the strain rate and F2 is the second mixing function obtained as follows:
(8)F2=tanhmax2kβ∗ωy,500vy2ω2.A production limiter exists to prevent the accumulation of turbulence in the stagnation zone:(9)Pk~=minPk,10⋅β∗ρkω,Pk=μt∂ui∂xj∂ui∂xj+∂uj∂xi,α=α1F+α21−F.All constants were calculated by mixing the constants of thek−ε and k−ω models. The constants were α1=5/9, α2=0.44, β∗=0.09, β2=0.0828, σk1=0.85, σk2=1, σω1=0.5, and σω2=0.856.The DES method based on SST was used to modify the turbulent dissipation term in thek equation as follows:
(10)ρε=β∗ρkω⟶β∗ρkω⋅FDES.Correction term:(11)FDES=maxLtCDESΔ,1,where ε denotes the turbulence dissipation rate, the maximum local grid size is Δ=maxΔx,Δy,Δz in the Cartesian coordinate system, the turbulence characteristic length size is Lt=k/β∗ω, and constant CDES=0.61.


## 2.6. Numerical Method Verification

To verify the correctness of the numerical simulation algorithm, the DES method was used to simulate the smooth cylinder in literature [30], and the results were compared with the wind tunnel experiment results in the literature [30]. In the reference, the diameter of the cylinder is 57 mm, the incoming flow velocity is 3 m/s, and the Reynolds number is Ra=1.17×104. The experimental results obtained in the reference are the average value of the pressure coefficient distribution on one side of the cylinder; the pressure coefficient was defined as CP=P−P∞/1/2ρv2. A three-dimensional cylindrical model with the same diameter was used for numerical simulation. The incoming flow velocity was set at 3 m/s, the incoming flow pressure was set at 101325 Pa, the Reynolds number was calculated at 11706, and the time step of unsteady flow was set as 2E−5 s. The numerical simulation results were compared with the experimental values in literature by taking the mean pressure distribution within 0.17 s, that is, the mean pressure distribution within three cycles. The results are shown in Figure 4, which is in good agreement with the experimental measurement results on the whole cylinder. Therefore, the numerical method is reliable.Figure 4
Comparison diagram of numerical simulation results and experimental values of smooth cylinder mean pressure coefficient distribution at 3 m/s.Furthermore, the flow around a cylinder with a roughness of 500μm in the literature [30] was simulated by the DES method to verify the roughness case. The calculation was carried out under standard atmospheric conditions of air pressure 101325 kPa and air temperature 288.15 K. And the flow with incoming velocity 15.2 m/s and Reynolds number 59326 is simulated. After averaging the pressure distribution within 0.24 s, the pressure distribution of the rough cylinder is shown in Figure 5. The simulation results are in good agreement with the experimental values.Figure 5
Comparison diagram of numerical simulation results and experimental values of rough cylinder mean pressure coefficient distribution at 15 m/s.


## 2.7. Solution Method and Grid Verification

Because the study involved a low-speed incompressible unsteady flow, a pressure-based unsteady solver was selected, and the coupling algorithm was used for the pressure-velocity coupling solution. A second-order implicit scheme was used for the time discretization method. For the spatial discretization method, the least-squares cell-based scheme was used for gradient discretization, the second-order scheme was used for pressure discretization, and the second-order upwind scheme was used for density, momentum, turbulent kinetic energy, turbulent dissipation rate, and energy discretization. In addition, higher-order term relaxation was used to enhance the convergence performance of the calculation. The calculation domain was set according to a previous study [27] and slightly expanded. As shown in Figures 6 and 7, the computational domain is a cylinder with a diameter of 1600 mm. The head vertex is approximately 23D away from the forward field. The height of the first layer of the grid near the wall was approximately 0.000463D. Wall y plus y+=1. In the experiment, the model is fixed on the datum platform through the support rod and bolt. Because the support structure is inevitable, the measurement will be affected by the support structure. Therefore, in the numerical simulation, the influence of the support structure is considered, and the support rod at the back of the model is included in the calculation model, which is reflected in the schematic diagram of the calculation domain (Figures 6 and 7). The surface mesh which is evenly distributed along the circumference is shown in Figure 8. The coordinate system used in the numerical simulation was different from the experimental coordinate system. The Cartesian coordinate system used in the numerical simulation was established according to the right-hand rule. The X-axis points from the head to the tail along the missile axis, the Y-axis is vertical, and the X-axis is vertical and upward.Figure 6
Diagram of boundary conditions.Figure 7
Grid diagram.Figure 8
A view of the surface mesh and a partial enlargement.Considering the randomness of the initial asymmetric flow direction, the full-multigrid initialization method is used in the numerical simulation to make the initial conditions produce the same asymmetric flow direction as the experiment. In this way, the initial disturbance of the incoming flow turbulence is not enough to determine the offset direction of the asymmetric flow field. Due to the existence of a bistable state, the calculation will converge to a near steady state, that is, one side of the initial offset. If the completely symmetric initialization method is adopted, the initial disturbance of incoming turbulence will play a decisive role in the direction of asymmetric back vortex migration of the flow field and thus will converge randomly to one of the two steady states.Grid independence verification is required for numerical simulation. Grid files of Grid 1, Grid 2, and Grid 3 with different grid quantities are selected for grid independence verification, and the uniform surface roughness isRa=0.8 μm. The numerical simulation results of grid number and time-average aerodynamic coefficient of the three grids are shown in Table 4. The elevation direction of the pitching moment coefficient is defined as positive; the reference length is 0.0432 mm, and the reference area is 0.001464 m2. As can be seen from Table 4, there is little difference between the simulation results of Grid 2 and Grid 3, while there is a great difference between the simulation results of Grid 1 and the other two sets of grids. It can be seen that the continuous increase of the grid quantity has little influence on the simulation results.Table 4
Calculation results of aerodynamic coefficients for three kinds of grids (no bottom resistance is included).
Grid quantityAxial force coefficientNormal force coefficientPitching moment coefficientLateral force coefficientYaw moment coefficientRoll moment coefficientGrid 13596092-0.110213.0865-95.51991.1024-16.35720.0008Grid 25059776-0.204714.8344-98.3868-3.518712.96500.0046Grid 36436364-0.200814.8987-99.0282-3.556512.77930.0021Furthermore, the average vorticity contours of the typical cross-section of asymmetric flow were selected and compared with the PIV experimental results. As shown in Figures9–12, it can be seen that the simulated results of vorticity distribution of Grid 2 and Grid 3 are close to the experimental results of PIV, showing the same asymmetry as the experimental results, while the simulated results of Grid 1 differ greatly from the experimental values, showing no obvious asymmetry. To save computing resources, the grid file with 5059776 grid quantity was selected for the next calculation.Figure 9
X-vorticity contour comparison at section 1: (a) PIV experiment, (b) result of Grid 1, (c) result of Grid 2, and (d) result of Grid 3.
(a)(b)(c)(d)Figure 10
X-vorticity contour comparison at section 2: (a) PIV experiment, (b) result of Grid 1, (c) result of Grid 2, and (d) result of Grid 3.
(a)(b)(c)(d)Figure 11
X-vorticity contour comparison at section 3 (a) PIV experiment, (b) result of Grid 1, (c) result of Grid 2, and (d) result of Grid 3.
(a)(b)(c)(d)Figure 12
X-vorticity contour comparison at section 4: (a) PIV experiment, (b) result of Grid 1, (c) result of Grid 2, and (d) result of Grid 3.
(a)(b)(c)(d)To verify the unsteady performance, delta time verification is required. Grid 2 is used for delta time verification, and three delta times,1E−4 s, 2E−5 s, and 5E−6 s are selected for delta time verification. The results are shown in Table 5. It can be seen that when the delta time is smaller than 2E−5 s, the aerodynamic coefficient does not change much. Therefore, the delta time 2E−5 s is selected for calculation.Table 5
Calculation results of aerodynamic coefficients for three kinds of delta time (no bottom resistance is included).
Time stepAxial force coefficientNormal force coefficientPitching moment coefficientLateral force coefficientYaw moment coefficientRoll moment coefficient1E−4 s-0.219414.6317-96.2298-3.702314.48900.00662E−5 s-0.204714.8344-98.3868-3.518712.96500.00465E−6 s-0.205314.8577-98.4543-3.517312.94010.0044In order to verify the unsteady characteristics of the flow around the slender body, a solid-state continuous laser was used to replace the pulse laser in the experiment, and a high-speed camera was used to replace the CCD camera to obtain the smoke flow diagram of the instantaneous flow. The frame rate of the high-speed camera was set as 1000 frames/s, and the aperture was f1.4. As shown in Figure13, the images obtained by the high-speed camera show that the structure of the flow field is steady within 2 s under experimental conditions. As shown in Figure 14, both the smoke flow diagram obtained in the experiment and the numerical simulation results show that there is no obvious difference in the flow field structure. Therefore, the time-average method is adopted in the subsequent analysis to compare the time-average flow field within 0.01 s.Figure 13
Flow structure is steady with time increase.Figure 14
There is no obvious difference in flow field structure.


## 2.8. Comparison of Aerodynamic Coefficients

The DES method was used to calculate the flow with velocity at 15.3 m/s, pressure at 102.56 kPa, and angle of attack at 55°. The surface roughness of the model was set to 0.8μm. The experimental results are compared with the DES simulation results after the coordinate system transformation, as shown in Table 6. The reference area is 1.4657E−3 m2, and the reference length is 0.0432 m.Table 6
The comparison between DES simulation results and experimental results after coordinate transformation.
Axial force coefficientNormal force coefficientLateral force coefficientPitching moment coefficientRoll moment coefficientYaw moment coefficientDES result-0.108411.30612.4612-70.19990.0011-13.0665EXP result-0.103712.27262.2477-75.85460.0042-11.8878Relative error4.57%7.87%9.50%7.45%73.78%9.91%As shown in Table6, except for the rolling moment coefficient, the errors of other force coefficients are all within 10%, while the rolling moment coefficient is quite close to 0. The calculation result of the DES method is reliable. The above content is added to the solution model and grid validation section.


## 3. Results and Discussion

### 3.1. Comparison between Numerical Simulation and Experimental Results

TheX-vorticity contour of the surface with roughness Ra=0.8 μm obtained by numerical simulation was compared with the PIV experimental values in four sections to verify the accuracy of the numerical simulation method. The location of the four sections was shown in Table 7. The X-vorticity contour with surface roughness Ra=0.8 μm was compared with that of the PIV experiment shown in Figures 15–17. The simulation results of the DES method based on k−ωSST are consistent with the experimental results of PIV in the flow-field structure. The values of the numerical simulation results of the main vortex intensity on the right side of section 1 are larger than those of the experimental results, and the main vortex intensities in the other three sections are close to the experimental values. The vorticity is slightly different in the free shear layer near the wall, but the positions of the main and additional vortices are consistent with the experimental values.Table 7
The location of the four sections.
Section 1Section 2Section 3Section 4Locationsx/D=1x/D=2x/D=3x/D=4Figure 15
Comparison contour of experimental and numerical simulation results for surface roughnessRa=0.8 μm at section 1.Figure 16
Comparison contour of experimental and numerical simulation results for surface roughnessRa=0.8 μm at section 2.Figure 17
Comparison contour of experimental and numerical simulation results for surface roughnessRa=0.8 μm at section 3.A comparison of the calculation and experimental results shows that the DES method can simulate the body well in a high angle of attack separation flow. The results of the numerical simulation of the flow-field structure are consistent with those of the experimental results. Although the near-wall flow field has a specific error in theX-direction vortex intensity, overall, the DES method simulation results can reflect the rule of flow, and we believe that the DES method computing result is credible.


### 3.2. Analysis of Numerical Simulation Results of Different Roughness Values

#### 3.2.1. Vorticity and Flow Separation

To compare the effect of a higher roughness value on the flow field with a high angle of attack, the DES method was used to simulate the high angle of attack flow of the spiral body with four roughness typesRa=0, 0.8, 5, and 100 μm.Figures18–21 shows the X-vorticity contour in four sections. We observe that for the left main vortex, with an increase in roughness, the high vorticity area of the left main vortex decreases slightly, but the vortex influence area increases visibly, and the free shear layer is further away from the wall surface. Simultaneously, the additional vortex in the same direction as the main vortex on the left is reduced, which may be caused by the decrease in the velocity gradient with the increase in roughness. From the streamline, the influence of the right main vortex on the right wall decreases with an increase in the roughness. Furthermore, the downstream flow still appears asymmetric when the surface is smooth, indicating that the stability of the flow field is not enough to maintain the symmetry of the flow field.Figure 18
Contrast contour ofX-vorticity with different roughness values at section 1: (a) Model 1, (b) Model 2, (c) Model 3, and (d) Model 4.
(a)(b)(c)(d)Figure 19
Contrast contour ofX-vorticity with different roughness values at section 2: (a) Model 1, (b) Model 2, (c) Model 3, and (d) Model 4.
(a)(b)(c)(d)Figure 20
Contrast contour ofX-vorticity with different roughness values at section 3: (a) Model 1, (b) Model 2, (c) Model 3, and (d) Model 4.
(a)(b)(c)(d)Figure 21
Contrast contour ofX-vorticity with different roughness values at section 4: (a) Model 1, (b) Model 2, (c) Model 3, and (d) Model 4.
(a)(b)(c)(d)A three-dimensional vortex intensity after time-averageω⇀ can be defined as follows:
(12)ω⇀=∂Uz∂y−∂Uy∂zi⇀+∂Ux∂z−∂Uz∂xj⇀+∂Uy∂x−∂Ux∂yk⇀.To facilitate the display in the contour, the norm ofω⇀ is ω⇀ in the cross-section contour, and it is obtained as follows:
(13)ω⇀=∂Uz∂y−∂Uy∂zi⇀+∂Ux∂z−∂Uz∂xj⇀+∂Uy∂x−∂Ux∂yk⇀.From Figures22 and 23, it is observed from the vorticity contour of multiple cross-sections that the back vortex is initially symmetrical and that the free shear layer generated from the lateral edge of the body enhanced the head vortex during the backward development process. The difference between the left and right sides is that the free shear layer generated at the left-side edge mainly strengthens the vorticity magnitude directly, and the vortex sheet rolled up at the right-side edge mainly increases the vortex region on the right side by merging small additional vortices with the main vortex. In addition, because the free shear layer on the left is closer to the wall, the velocity gradient of the shear layer is larger than that on the right, and the additional vortex is stronger than that on the right.Figure 22
Comparison of vorticity and streamline on the left side with different values of roughness: (a) Model 1, (b) Model 2, (c) Model 3, and (d) Model 4.
(a)(b)(c)(d)Figure 23
Comparison of vorticity and streamline on the right side with different values of roughness: (a) Model 1, (b) Model 2, (c) Model 3, and (d) Model 4.
(a)(b)(c)(d)Figure24 shows the streamlines, viewed from the positive to the negative direction of the Y-axis. The streamlines start from the head and extend backward as spirals. The spiral radius becomes larger further downstream, and with the increase in roughness, the pitch of the streamlines in the negative direction of the Z-axis also increases. By comparing the streamlines in the negative direction of the z-axis in Figures 22 and 23, it can be observed that while the vortex sheet rolled up on the left side enhancing the strength of the main vortex, it also changes the position of the main vortex core, such that the streamlines originally located near the vortex core are thrown out of the vicinity of the vortex core. The streamline in the positive z-direction does not exhibit this change.Figure 24
Comparison of streamline with different roughness values: (a) Model 1, (b) Model 2, (c) Model 3, and (d) Model 4.
(a)(b)(c)(d)To compare the separation of both sides of the head with different roughness values, the relative velocity angle was used to identify the separation position of the appendage flow. As shown in Figure25, the boundary of the leeward region in the figure represents the sudden reversal of the attached flow, implying that the attached flow is separated at this point, and several separation lines in the leeward region are caused by the main and additional vortices. Figure 26 shows the separation lines on both sides extracted from Figure 25.Figure 25
Comparison of surface relative velocity angle with different roughness values: (a) Model 1, (b) Model 2, (c) Model 3, and (d) Model 4.Figure 26
Comparison of separation lines with different roughness values on both sides.By comparing the separation lines on both sides of the same model, it is observed that the separation lines on both sides of the head are close, but the separation line on the left is clearly higher than that on the right in the back, implying that the separation on the left side of the low vortex occurs later than that on the right side of the high vortex, and this is consistent with the experimental phenomenon. Comparing different models of the separation line, we observe that with the increase in surface roughness, the separation line on the left side of the wave becomes smooth. In addition, the left separation lines of the four roughness models are similar. In the front segment of the head, the model separation line with roughnessRa=100 μm was the highest while the model separation line with roughness Ra=100 μm was the highest at the back segment close to the cylinder.To explain the aforementioned changes of different roughness separation points, the velocity and pressure contour ofx/D=1 and x/D=3 were selected for comparison, as shown in Figures 27–29. Figures 27 and 28 show the velocity contour of the left separation area at x/D=1 and x/D=3, respectively, and Figure 29 shows the pressure contour of the left separation area at x/D=3. It is observed from Figure 27 that there is no significant difference in the velocity contour of the separation area, and the low-velocity layer at the bottom of the boundary layer becomes thinner with an increase in the surface roughness. In Figure 28, it is clear that the velocity contour of the separation zone is different. A large increase in roughness leads to a higher velocity in the separation zone, and the shear layer exhibits a slight outward deviation when the roughness is relatively small. The pressure distribution shown in Figure 29 can explain the changes in the separation lines with different roughness values. At the head, the separation point is mainly affected by the turbulence of the boundary layer, and the increase in roughness causes the separation point to move backward. At the rear of the pointed arch section, the separation point was mainly affected by the adverse pressure gradient. As the roughness increased, the adverse pressure gradient increased, leading to premature separation.Figure 27
Comparison of closed-wall velocity with different roughness values at cross-section 1: (a) Model 1, (b) Model 2, (c) Model 3, and (d) Model 4.
(a)(b)(c)(d)Figure 28
Comparison of closed-wall velocity with different roughness values at cross-section 3: (a) Model 1, (b) Model 2, (c) Model 3, and (d) Model 4.
(a)(b)(c)(d)Figure 29
Comparison of closed-wall pressure with different roughness values at cross-section 3: (a) Model 1, (b) Model 2, (c) Model 3, and (d) Model 4.
(a)(b)(c)(d)


#### 3.2.2. Surface Pressure

The influence of roughness on the flow field in the head is limited owing to the presence of strong vortices in the head. In the rear region of the body, the influence of roughness is significantly enhanced with an increase in the length of the body and the separation, fragmentation, and dissipation of vortices.Figure30 shows the surface pressure distribution contour in the leeward area of the entire missile, Figure 31 shows the surface pressure distribution contour in the leeward area of the head region, and Figure 32 shows the surface pressure distribution contour in the leeward area of the rear part of the missile body. By comparison, the pressure distribution on the head does not change significantly with the roughness, and the pressure on the positive Z-axis side is higher than that on the negative Z-axis side. From the head of the cylinder section, the high- and low-pressure areas on both sides began to appear alternately along the axial direction. It is clearly observed from Figure 30 that the period of pressure alternation increases with the increase in roughness, and finally, the pressure on both sides tends to be consistent. This phenomenon leads to different effects of different roughness values on the lateral force on the body with finite length.Figure 30
Comparison of leeward pressure with different roughness values: (a) Model 1, (b) Model 2, (c) Model 3, and (d) Model 4.
(a)(b)(c)(d)Figure 31
Comparison of leeward pressure at the head of body with different roughness values: (a) Model 1, (b) Model 2, (c) Model 3, and (d) Model 4.
(a)(b)(c)(d)Figure 32
Comparison of leeward pressure at the middle and body tail with different roughness values: (a) Model 1, (b) Model 2, (c) Model 3, and (d) Model 4.
(a)(b)(c)(d)As shown in Figure30, Model 1 with roughness Ra=0 μm pressure alternation changes at the leeward area five times, and Model 4 with roughness Ra=100 μm pressure alternation changes only 2 times at the leeward area; the pressure distribution on both sides of the body in the leeward leads to lateral force changes, as shown in Table 8. The stress distribution of the cylindrical section controls the size of the full play to the lateral force and direction.Table 8
Lateral force coefficient for 4 kinds of models.
ModelRa=0 μmRa=0.8 μmRa=5 μmRa=100 μmLateral force coefficient-2.4900-2.3673-3.50251.3310To study the lateral force distribution, the lateral force coefficient per unit length is defined as follows:(14)Cfzpm=Fzpm1/2ρv2D,where Fzpm is the lateral force per unit length.The lateral force coefficient per unit lengthCfzpm corresponds to the axial alternating change of pressure. The lateral force coefficient per unit length Cfzpm of the four models shows the same axial alternating reverse phenomenon as the pressure, as shown in Figure 33. The four models in the head section (x/D<5) have relatively little difference in Cfzpm, and the main difference is in the body part (x/D>5). The effect of roughness is relatively small because the main vortex of the head is obviously asymmetric and increasing. Shown in Figures 34–40 is the pressure coefficient distribution of the four models with different roughness on different sections. The intersection point between the angle of attack plane and the lower surface of the body is defined as 0°, and the point of rotation around the counterclockwise direction is 360°. From Figures 34–36, the distribution of the circumferential pressure coefficient from 1D to 3D is given, corresponding to the downward phase of the curve in Figure 33. It can be seen that the asymmetry of the distribution of the circumferential pressure coefficient increases significantly from the nose to the axial 3D position. Cfzpm of the four models all reached a minimum value near the 3D position, and the distribution trend of the circumferential pressure coefficient of different roughness was consistent; although there was a gap, the gap was not significant. The position of axial distance from nose 3D to 6.5D corresponds to the curve rising stage in Figure 33, as shown in Figures 36–39. At this stage, the asymmetry of the distribution of the circumferential pressure coefficient decreases, and the asymmetry is not obvious at position 5D. The distribution law of the circumferential pressure coefficient at 6.5D is obviously contrary to that before 5D. In addition, from the 4D position, the circumferential pressure coefficients of the model with Ra=100 μm begin to show obvious differences. After 5D from the nose, Cfzpm of the four models showed obvious differences. The roughness Ra=0 μm, 0.8 μm, and 5 μm all experienced multiple positive and negative Cfzpm changes, and the positive Cfzpm interval length was relatively short, and the peak value was low. Cfzpm of the Ra=0 μm model has a higher peak value at the second peak, which means the vortex intensity attenuation is slow. The Cfzpm peak values of Ra=0.8 μm and Ra=5 μm are both lower and close to 0 after the second peak value. Model 4 with roughness Ra=100 μm will produce a large unit lateral force coefficient (Cfzpm) in the z-square direction, which gradually decreases to near 0 at the tail. As shown in Figure 40, at position 8D away from the nose, the model with roughness Ra=100 μm still has obvious asymmetry. According to the integration of Cfzpm in the x-direction of the four models, the lateral force coefficients of the four models were obtained. Obviously, the integral value of Ra=100 μm was the largest, and the integral value of Ra=5 μm was the smallest, which was consistent with the results in Figure 33.Figure 33
Axial variation of lateral force coefficient per unit length.Figure 34
Pressure coefficient distribution of different models at1D.Figure 35
Pressure coefficient distribution of different models at2D.Figure 36
Pressure coefficient distribution of different models at3D.Figure 37
Pressure coefficient distribution of different models at4D.Figure 38
Pressure coefficient distribution of different models at5D.Figure 39
Pressure coefficient distribution of different models at6.5D.Figure 40
Pressure coefficient distribution of different models at8D.The circumferential pressure coefficient distribution of different sections of the same model is compared, as shown in Figure41. In Model 1 with Ra=5 μm, the curve asymmetry increases significantly from 1D to 3D, and the pressure coefficient on the left side is significantly lower than that on the right side. The pressure coefficient difference on both sides of the 3D section is the largest, and the pressure coefficient on the left side is higher than that on the right side in the 6D section. At 9D and 10D sections, the pressure coefficients on both sides reverse again and finally reach equilibrium at the 12D section. The change in Model 2 with Ra=0.8 μm and Model 3 with Ra=5 μm was consistent with Model 1 with Ra=0 μm, as shown in Figures 42 and 43, but the pressure distribution of Model 2 and Model 3 tended to be symmetric faster. Compared with Model 3, the pressure distribution asymmetry at the 6D section is more obvious in Model 2. The pressure coefficient distribution of Model 4 with Ra=100 μm is shown in Figure 44. After, the pressure coefficient distribution of the section from 3D to 6D changes from low left and high right to high left and low right, and there is no lower pressure coefficient on the left, which is consistent with the trend in Figure 33. In addition, the pressure coefficient on the left side of both 6D section and 9D section shows local drastic changes. According to the analysis of vorticity contour in Figure 45, in these two sections, the secondary vortex near the wall is the cause of local drastic changes in the surface pressure coefficient.Figure 41
The pressure coefficient of different sections varies with the roll angle at Model 1 withRa=0 μm.Figure 42
The pressure coefficient of different sections varies with the roll angle at Model 2 withRa=0.8 μm.Figure 43
The pressure coefficient of different sections varies with the roll angle at Model 3 withRa=5 μm.Figure 44
The pressure coefficient of different sections varies with the roll angle at Model 4 withRa=100 μm.Figure 45
Vorticity at sectionx=6D and x=9D in Model 4.


## 3.1. Comparison between Numerical Simulation and Experimental Results

TheX-vorticity contour of the surface with roughness Ra=0.8 μm obtained by numerical simulation was compared with the PIV experimental values in four sections to verify the accuracy of the numerical simulation method. The location of the four sections was shown in Table 7. The X-vorticity contour with surface roughness Ra=0.8 μm was compared with that of the PIV experiment shown in Figures 15–17. The simulation results of the DES method based on k−ωSST are consistent with the experimental results of PIV in the flow-field structure. The values of the numerical simulation results of the main vortex intensity on the right side of section 1 are larger than those of the experimental results, and the main vortex intensities in the other three sections are close to the experimental values. The vorticity is slightly different in the free shear layer near the wall, but the positions of the main and additional vortices are consistent with the experimental values.Table 7
The location of the four sections.
Section 1Section 2Section 3Section 4Locationsx/D=1x/D=2x/D=3x/D=4Figure 15
Comparison contour of experimental and numerical simulation results for surface roughnessRa=0.8 μm at section 1.Figure 16
Comparison contour of experimental and numerical simulation results for surface roughnessRa=0.8 μm at section 2.Figure 17
Comparison contour of experimental and numerical simulation results for surface roughnessRa=0.8 μm at section 3.A comparison of the calculation and experimental results shows that the DES method can simulate the body well in a high angle of attack separation flow. The results of the numerical simulation of the flow-field structure are consistent with those of the experimental results. Although the near-wall flow field has a specific error in theX-direction vortex intensity, overall, the DES method simulation results can reflect the rule of flow, and we believe that the DES method computing result is credible.


## 3.2. Analysis of Numerical Simulation Results of Different Roughness Values

### 3.2.1. Vorticity and Flow Separation

To compare the effect of a higher roughness value on the flow field with a high angle of attack, the DES method was used to simulate the high angle of attack flow of the spiral body with four roughness typesRa=0, 0.8, 5, and 100 μm.Figures18–21 shows the X-vorticity contour in four sections. We observe that for the left main vortex, with an increase in roughness, the high vorticity area of the left main vortex decreases slightly, but the vortex influence area increases visibly, and the free shear layer is further away from the wall surface. Simultaneously, the additional vortex in the same direction as the main vortex on the left is reduced, which may be caused by the decrease in the velocity gradient with the increase in roughness. From the streamline, the influence of the right main vortex on the right wall decreases with an increase in the roughness. Furthermore, the downstream flow still appears asymmetric when the surface is smooth, indicating that the stability of the flow field is not enough to maintain the symmetry of the flow field.Figure 18
Contrast contour ofX-vorticity with different roughness values at section 1: (a) Model 1, (b) Model 2, (c) Model 3, and (d) Model 4.
(a)(b)(c)(d)Figure 19
Contrast contour ofX-vorticity with different roughness values at section 2: (a) Model 1, (b) Model 2, (c) Model 3, and (d) Model 4.
(a)(b)(c)(d)Figure 20
Contrast contour ofX-vorticity with different roughness values at section 3: (a) Model 1, (b) Model 2, (c) Model 3, and (d) Model 4.
(a)(b)(c)(d)Figure 21
Contrast contour ofX-vorticity with different roughness values at section 4: (a) Model 1, (b) Model 2, (c) Model 3, and (d) Model 4.
(a)(b)(c)(d)A three-dimensional vortex intensity after time-averageω⇀ can be defined as follows:
(12)ω⇀=∂Uz∂y−∂Uy∂zi⇀+∂Ux∂z−∂Uz∂xj⇀+∂Uy∂x−∂Ux∂yk⇀.To facilitate the display in the contour, the norm ofω⇀ is ω⇀ in the cross-section contour, and it is obtained as follows:
(13)ω⇀=∂Uz∂y−∂Uy∂zi⇀+∂Ux∂z−∂Uz∂xj⇀+∂Uy∂x−∂Ux∂yk⇀.From Figures22 and 23, it is observed from the vorticity contour of multiple cross-sections that the back vortex is initially symmetrical and that the free shear layer generated from the lateral edge of the body enhanced the head vortex during the backward development process. The difference between the left and right sides is that the free shear layer generated at the left-side edge mainly strengthens the vorticity magnitude directly, and the vortex sheet rolled up at the right-side edge mainly increases the vortex region on the right side by merging small additional vortices with the main vortex. In addition, because the free shear layer on the left is closer to the wall, the velocity gradient of the shear layer is larger than that on the right, and the additional vortex is stronger than that on the right.Figure 22
Comparison of vorticity and streamline on the left side with different values of roughness: (a) Model 1, (b) Model 2, (c) Model 3, and (d) Model 4.
(a)(b)(c)(d)Figure 23
Comparison of vorticity and streamline on the right side with different values of roughness: (a) Model 1, (b) Model 2, (c) Model 3, and (d) Model 4.
(a)(b)(c)(d)Figure24 shows the streamlines, viewed from the positive to the negative direction of the Y-axis. The streamlines start from the head and extend backward as spirals. The spiral radius becomes larger further downstream, and with the increase in roughness, the pitch of the streamlines in the negative direction of the Z-axis also increases. By comparing the streamlines in the negative direction of the z-axis in Figures 22 and 23, it can be observed that while the vortex sheet rolled up on the left side enhancing the strength of the main vortex, it also changes the position of the main vortex core, such that the streamlines originally located near the vortex core are thrown out of the vicinity of the vortex core. The streamline in the positive z-direction does not exhibit this change.Figure 24
Comparison of streamline with different roughness values: (a) Model 1, (b) Model 2, (c) Model 3, and (d) Model 4.
(a)(b)(c)(d)To compare the separation of both sides of the head with different roughness values, the relative velocity angle was used to identify the separation position of the appendage flow. As shown in Figure25, the boundary of the leeward region in the figure represents the sudden reversal of the attached flow, implying that the attached flow is separated at this point, and several separation lines in the leeward region are caused by the main and additional vortices. Figure 26 shows the separation lines on both sides extracted from Figure 25.Figure 25
Comparison of surface relative velocity angle with different roughness values: (a) Model 1, (b) Model 2, (c) Model 3, and (d) Model 4.Figure 26
Comparison of separation lines with different roughness values on both sides.By comparing the separation lines on both sides of the same model, it is observed that the separation lines on both sides of the head are close, but the separation line on the left is clearly higher than that on the right in the back, implying that the separation on the left side of the low vortex occurs later than that on the right side of the high vortex, and this is consistent with the experimental phenomenon. Comparing different models of the separation line, we observe that with the increase in surface roughness, the separation line on the left side of the wave becomes smooth. In addition, the left separation lines of the four roughness models are similar. In the front segment of the head, the model separation line with roughnessRa=100 μm was the highest while the model separation line with roughness Ra=100 μm was the highest at the back segment close to the cylinder.To explain the aforementioned changes of different roughness separation points, the velocity and pressure contour ofx/D=1 and x/D=3 were selected for comparison, as shown in Figures 27–29. Figures 27 and 28 show the velocity contour of the left separation area at x/D=1 and x/D=3, respectively, and Figure 29 shows the pressure contour of the left separation area at x/D=3. It is observed from Figure 27 that there is no significant difference in the velocity contour of the separation area, and the low-velocity layer at the bottom of the boundary layer becomes thinner with an increase in the surface roughness. In Figure 28, it is clear that the velocity contour of the separation zone is different. A large increase in roughness leads to a higher velocity in the separation zone, and the shear layer exhibits a slight outward deviation when the roughness is relatively small. The pressure distribution shown in Figure 29 can explain the changes in the separation lines with different roughness values. At the head, the separation point is mainly affected by the turbulence of the boundary layer, and the increase in roughness causes the separation point to move backward. At the rear of the pointed arch section, the separation point was mainly affected by the adverse pressure gradient. As the roughness increased, the adverse pressure gradient increased, leading to premature separation.Figure 27
Comparison of closed-wall velocity with different roughness values at cross-section 1: (a) Model 1, (b) Model 2, (c) Model 3, and (d) Model 4.
(a)(b)(c)(d)Figure 28
Comparison of closed-wall velocity with different roughness values at cross-section 3: (a) Model 1, (b) Model 2, (c) Model 3, and (d) Model 4.
(a)(b)(c)(d)Figure 29
Comparison of closed-wall pressure with different roughness values at cross-section 3: (a) Model 1, (b) Model 2, (c) Model 3, and (d) Model 4.
(a)(b)(c)(d)


### 3.2.2. Surface Pressure

The influence of roughness on the flow field in the head is limited owing to the presence of strong vortices in the head. In the rear region of the body, the influence of roughness is significantly enhanced with an increase in the length of the body and the separation, fragmentation, and dissipation of vortices.Figure30 shows the surface pressure distribution contour in the leeward area of the entire missile, Figure 31 shows the surface pressure distribution contour in the leeward area of the head region, and Figure 32 shows the surface pressure distribution contour in the leeward area of the rear part of the missile body. By comparison, the pressure distribution on the head does not change significantly with the roughness, and the pressure on the positive Z-axis side is higher than that on the negative Z-axis side. From the head of the cylinder section, the high- and low-pressure areas on both sides began to appear alternately along the axial direction. It is clearly observed from Figure 30 that the period of pressure alternation increases with the increase in roughness, and finally, the pressure on both sides tends to be consistent. This phenomenon leads to different effects of different roughness values on the lateral force on the body with finite length.Figure 30
Comparison of leeward pressure with different roughness values: (a) Model 1, (b) Model 2, (c) Model 3, and (d) Model 4.
(a)(b)(c)(d)Figure 31
Comparison of leeward pressure at the head of body with different roughness values: (a) Model 1, (b) Model 2, (c) Model 3, and (d) Model 4.
(a)(b)(c)(d)Figure 32
Comparison of leeward pressure at the middle and body tail with different roughness values: (a) Model 1, (b) Model 2, (c) Model 3, and (d) Model 4.
(a)(b)(c)(d)As shown in Figure30, Model 1 with roughness Ra=0 μm pressure alternation changes at the leeward area five times, and Model 4 with roughness Ra=100 μm pressure alternation changes only 2 times at the leeward area; the pressure distribution on both sides of the body in the leeward leads to lateral force changes, as shown in Table 8. The stress distribution of the cylindrical section controls the size of the full play to the lateral force and direction.Table 8
Lateral force coefficient for 4 kinds of models.
ModelRa=0 μmRa=0.8 μmRa=5 μmRa=100 μmLateral force coefficient-2.4900-2.3673-3.50251.3310To study the lateral force distribution, the lateral force coefficient per unit length is defined as follows:(14)Cfzpm=Fzpm1/2ρv2D,where Fzpm is the lateral force per unit length.The lateral force coefficient per unit lengthCfzpm corresponds to the axial alternating change of pressure. The lateral force coefficient per unit length Cfzpm of the four models shows the same axial alternating reverse phenomenon as the pressure, as shown in Figure 33. The four models in the head section (x/D<5) have relatively little difference in Cfzpm, and the main difference is in the body part (x/D>5). The effect of roughness is relatively small because the main vortex of the head is obviously asymmetric and increasing. Shown in Figures 34–40 is the pressure coefficient distribution of the four models with different roughness on different sections. The intersection point between the angle of attack plane and the lower surface of the body is defined as 0°, and the point of rotation around the counterclockwise direction is 360°. From Figures 34–36, the distribution of the circumferential pressure coefficient from 1D to 3D is given, corresponding to the downward phase of the curve in Figure 33. It can be seen that the asymmetry of the distribution of the circumferential pressure coefficient increases significantly from the nose to the axial 3D position. Cfzpm of the four models all reached a minimum value near the 3D position, and the distribution trend of the circumferential pressure coefficient of different roughness was consistent; although there was a gap, the gap was not significant. The position of axial distance from nose 3D to 6.5D corresponds to the curve rising stage in Figure 33, as shown in Figures 36–39. At this stage, the asymmetry of the distribution of the circumferential pressure coefficient decreases, and the asymmetry is not obvious at position 5D. The distribution law of the circumferential pressure coefficient at 6.5D is obviously contrary to that before 5D. In addition, from the 4D position, the circumferential pressure coefficients of the model with Ra=100 μm begin to show obvious differences. After 5D from the nose, Cfzpm of the four models showed obvious differences. The roughness Ra=0 μm, 0.8 μm, and 5 μm all experienced multiple positive and negative Cfzpm changes, and the positive Cfzpm interval length was relatively short, and the peak value was low. Cfzpm of the Ra=0 μm model has a higher peak value at the second peak, which means the vortex intensity attenuation is slow. The Cfzpm peak values of Ra=0.8 μm and Ra=5 μm are both lower and close to 0 after the second peak value. Model 4 with roughness Ra=100 μm will produce a large unit lateral force coefficient (Cfzpm) in the z-square direction, which gradually decreases to near 0 at the tail. As shown in Figure 40, at position 8D away from the nose, the model with roughness Ra=100 μm still has obvious asymmetry. According to the integration of Cfzpm in the x-direction of the four models, the lateral force coefficients of the four models were obtained. Obviously, the integral value of Ra=100 μm was the largest, and the integral value of Ra=5 μm was the smallest, which was consistent with the results in Figure 33.Figure 33
Axial variation of lateral force coefficient per unit length.Figure 34
Pressure coefficient distribution of different models at1D.Figure 35
Pressure coefficient distribution of different models at2D.Figure 36
Pressure coefficient distribution of different models at3D.Figure 37
Pressure coefficient distribution of different models at4D.Figure 38
Pressure coefficient distribution of different models at5D.Figure 39
Pressure coefficient distribution of different models at6.5D.Figure 40
Pressure coefficient distribution of different models at8D.The circumferential pressure coefficient distribution of different sections of the same model is compared, as shown in Figure41. In Model 1 with Ra=5 μm, the curve asymmetry increases significantly from 1D to 3D, and the pressure coefficient on the left side is significantly lower than that on the right side. The pressure coefficient difference on both sides of the 3D section is the largest, and the pressure coefficient on the left side is higher than that on the right side in the 6D section. At 9D and 10D sections, the pressure coefficients on both sides reverse again and finally reach equilibrium at the 12D section. The change in Model 2 with Ra=0.8 μm and Model 3 with Ra=5 μm was consistent with Model 1 with Ra=0 μm, as shown in Figures 42 and 43, but the pressure distribution of Model 2 and Model 3 tended to be symmetric faster. Compared with Model 3, the pressure distribution asymmetry at the 6D section is more obvious in Model 2. The pressure coefficient distribution of Model 4 with Ra=100 μm is shown in Figure 44. After, the pressure coefficient distribution of the section from 3D to 6D changes from low left and high right to high left and low right, and there is no lower pressure coefficient on the left, which is consistent with the trend in Figure 33. In addition, the pressure coefficient on the left side of both 6D section and 9D section shows local drastic changes. According to the analysis of vorticity contour in Figure 45, in these two sections, the secondary vortex near the wall is the cause of local drastic changes in the surface pressure coefficient.Figure 41
The pressure coefficient of different sections varies with the roll angle at Model 1 withRa=0 μm.Figure 42
The pressure coefficient of different sections varies with the roll angle at Model 2 withRa=0.8 μm.Figure 43
The pressure coefficient of different sections varies with the roll angle at Model 3 withRa=5 μm.Figure 44
The pressure coefficient of different sections varies with the roll angle at Model 4 withRa=100 μm.Figure 45
Vorticity at sectionx=6D and x=9D in Model 4.


## 3.2.1. Vorticity and Flow Separation

To compare the effect of a higher roughness value on the flow field with a high angle of attack, the DES method was used to simulate the high angle of attack flow of the spiral body with four roughness typesRa=0, 0.8, 5, and 100 μm.Figures18–21 shows the X-vorticity contour in four sections. We observe that for the left main vortex, with an increase in roughness, the high vorticity area of the left main vortex decreases slightly, but the vortex influence area increases visibly, and the free shear layer is further away from the wall surface. Simultaneously, the additional vortex in the same direction as the main vortex on the left is reduced, which may be caused by the decrease in the velocity gradient with the increase in roughness. From the streamline, the influence of the right main vortex on the right wall decreases with an increase in the roughness. Furthermore, the downstream flow still appears asymmetric when the surface is smooth, indicating that the stability of the flow field is not enough to maintain the symmetry of the flow field.Figure 18
Contrast contour ofX-vorticity with different roughness values at section 1: (a) Model 1, (b) Model 2, (c) Model 3, and (d) Model 4.
(a)(b)(c)(d)Figure 19
Contrast contour ofX-vorticity with different roughness values at section 2: (a) Model 1, (b) Model 2, (c) Model 3, and (d) Model 4.
(a)(b)(c)(d)Figure 20
Contrast contour ofX-vorticity with different roughness values at section 3: (a) Model 1, (b) Model 2, (c) Model 3, and (d) Model 4.
(a)(b)(c)(d)Figure 21
Contrast contour ofX-vorticity with different roughness values at section 4: (a) Model 1, (b) Model 2, (c) Model 3, and (d) Model 4.
(a)(b)(c)(d)A three-dimensional vortex intensity after time-averageω⇀ can be defined as follows:
(12)ω⇀=∂Uz∂y−∂Uy∂zi⇀+∂Ux∂z−∂Uz∂xj⇀+∂Uy∂x−∂Ux∂yk⇀.To facilitate the display in the contour, the norm ofω⇀ is ω⇀ in the cross-section contour, and it is obtained as follows:
(13)ω⇀=∂Uz∂y−∂Uy∂zi⇀+∂Ux∂z−∂Uz∂xj⇀+∂Uy∂x−∂Ux∂yk⇀.From Figures22 and 23, it is observed from the vorticity contour of multiple cross-sections that the back vortex is initially symmetrical and that the free shear layer generated from the lateral edge of the body enhanced the head vortex during the backward development process. The difference between the left and right sides is that the free shear layer generated at the left-side edge mainly strengthens the vorticity magnitude directly, and the vortex sheet rolled up at the right-side edge mainly increases the vortex region on the right side by merging small additional vortices with the main vortex. In addition, because the free shear layer on the left is closer to the wall, the velocity gradient of the shear layer is larger than that on the right, and the additional vortex is stronger than that on the right.Figure 22
Comparison of vorticity and streamline on the left side with different values of roughness: (a) Model 1, (b) Model 2, (c) Model 3, and (d) Model 4.
(a)(b)(c)(d)Figure 23
Comparison of vorticity and streamline on the right side with different values of roughness: (a) Model 1, (b) Model 2, (c) Model 3, and (d) Model 4.
(a)(b)(c)(d)Figure24 shows the streamlines, viewed from the positive to the negative direction of the Y-axis. The streamlines start from the head and extend backward as spirals. The spiral radius becomes larger further downstream, and with the increase in roughness, the pitch of the streamlines in the negative direction of the Z-axis also increases. By comparing the streamlines in the negative direction of the z-axis in Figures 22 and 23, it can be observed that while the vortex sheet rolled up on the left side enhancing the strength of the main vortex, it also changes the position of the main vortex core, such that the streamlines originally located near the vortex core are thrown out of the vicinity of the vortex core. The streamline in the positive z-direction does not exhibit this change.Figure 24
Comparison of streamline with different roughness values: (a) Model 1, (b) Model 2, (c) Model 3, and (d) Model 4.
(a)(b)(c)(d)To compare the separation of both sides of the head with different roughness values, the relative velocity angle was used to identify the separation position of the appendage flow. As shown in Figure25, the boundary of the leeward region in the figure represents the sudden reversal of the attached flow, implying that the attached flow is separated at this point, and several separation lines in the leeward region are caused by the main and additional vortices. Figure 26 shows the separation lines on both sides extracted from Figure 25.Figure 25
Comparison of surface relative velocity angle with different roughness values: (a) Model 1, (b) Model 2, (c) Model 3, and (d) Model 4.Figure 26
Comparison of separation lines with different roughness values on both sides.By comparing the separation lines on both sides of the same model, it is observed that the separation lines on both sides of the head are close, but the separation line on the left is clearly higher than that on the right in the back, implying that the separation on the left side of the low vortex occurs later than that on the right side of the high vortex, and this is consistent with the experimental phenomenon. Comparing different models of the separation line, we observe that with the increase in surface roughness, the separation line on the left side of the wave becomes smooth. In addition, the left separation lines of the four roughness models are similar. In the front segment of the head, the model separation line with roughnessRa=100 μm was the highest while the model separation line with roughness Ra=100 μm was the highest at the back segment close to the cylinder.To explain the aforementioned changes of different roughness separation points, the velocity and pressure contour ofx/D=1 and x/D=3 were selected for comparison, as shown in Figures 27–29. Figures 27 and 28 show the velocity contour of the left separation area at x/D=1 and x/D=3, respectively, and Figure 29 shows the pressure contour of the left separation area at x/D=3. It is observed from Figure 27 that there is no significant difference in the velocity contour of the separation area, and the low-velocity layer at the bottom of the boundary layer becomes thinner with an increase in the surface roughness. In Figure 28, it is clear that the velocity contour of the separation zone is different. A large increase in roughness leads to a higher velocity in the separation zone, and the shear layer exhibits a slight outward deviation when the roughness is relatively small. The pressure distribution shown in Figure 29 can explain the changes in the separation lines with different roughness values. At the head, the separation point is mainly affected by the turbulence of the boundary layer, and the increase in roughness causes the separation point to move backward. At the rear of the pointed arch section, the separation point was mainly affected by the adverse pressure gradient. As the roughness increased, the adverse pressure gradient increased, leading to premature separation.Figure 27
Comparison of closed-wall velocity with different roughness values at cross-section 1: (a) Model 1, (b) Model 2, (c) Model 3, and (d) Model 4.
(a)(b)(c)(d)Figure 28
Comparison of closed-wall velocity with different roughness values at cross-section 3: (a) Model 1, (b) Model 2, (c) Model 3, and (d) Model 4.
(a)(b)(c)(d)Figure 29
Comparison of closed-wall pressure with different roughness values at cross-section 3: (a) Model 1, (b) Model 2, (c) Model 3, and (d) Model 4.
(a)(b)(c)(d)


## 3.2.2. Surface Pressure

The influence of roughness on the flow field in the head is limited owing to the presence of strong vortices in the head. In the rear region of the body, the influence of roughness is significantly enhanced with an increase in the length of the body and the separation, fragmentation, and dissipation of vortices.Figure30 shows the surface pressure distribution contour in the leeward area of the entire missile, Figure 31 shows the surface pressure distribution contour in the leeward area of the head region, and Figure 32 shows the surface pressure distribution contour in the leeward area of the rear part of the missile body. By comparison, the pressure distribution on the head does not change significantly with the roughness, and the pressure on the positive Z-axis side is higher than that on the negative Z-axis side. From the head of the cylinder section, the high- and low-pressure areas on both sides began to appear alternately along the axial direction. It is clearly observed from Figure 30 that the period of pressure alternation increases with the increase in roughness, and finally, the pressure on both sides tends to be consistent. This phenomenon leads to different effects of different roughness values on the lateral force on the body with finite length.Figure 30
Comparison of leeward pressure with different roughness values: (a) Model 1, (b) Model 2, (c) Model 3, and (d) Model 4.
(a)(b)(c)(d)Figure 31
Comparison of leeward pressure at the head of body with different roughness values: (a) Model 1, (b) Model 2, (c) Model 3, and (d) Model 4.
(a)(b)(c)(d)Figure 32
Comparison of leeward pressure at the middle and body tail with different roughness values: (a) Model 1, (b) Model 2, (c) Model 3, and (d) Model 4.
(a)(b)(c)(d)As shown in Figure30, Model 1 with roughness Ra=0 μm pressure alternation changes at the leeward area five times, and Model 4 with roughness Ra=100 μm pressure alternation changes only 2 times at the leeward area; the pressure distribution on both sides of the body in the leeward leads to lateral force changes, as shown in Table 8. The stress distribution of the cylindrical section controls the size of the full play to the lateral force and direction.Table 8
Lateral force coefficient for 4 kinds of models.
ModelRa=0 μmRa=0.8 μmRa=5 μmRa=100 μmLateral force coefficient-2.4900-2.3673-3.50251.3310To study the lateral force distribution, the lateral force coefficient per unit length is defined as follows:(14)Cfzpm=Fzpm1/2ρv2D,where Fzpm is the lateral force per unit length.The lateral force coefficient per unit lengthCfzpm corresponds to the axial alternating change of pressure. The lateral force coefficient per unit length Cfzpm of the four models shows the same axial alternating reverse phenomenon as the pressure, as shown in Figure 33. The four models in the head section (x/D<5) have relatively little difference in Cfzpm, and the main difference is in the body part (x/D>5). The effect of roughness is relatively small because the main vortex of the head is obviously asymmetric and increasing. Shown in Figures 34–40 is the pressure coefficient distribution of the four models with different roughness on different sections. The intersection point between the angle of attack plane and the lower surface of the body is defined as 0°, and the point of rotation around the counterclockwise direction is 360°. From Figures 34–36, the distribution of the circumferential pressure coefficient from 1D to 3D is given, corresponding to the downward phase of the curve in Figure 33. It can be seen that the asymmetry of the distribution of the circumferential pressure coefficient increases significantly from the nose to the axial 3D position. Cfzpm of the four models all reached a minimum value near the 3D position, and the distribution trend of the circumferential pressure coefficient of different roughness was consistent; although there was a gap, the gap was not significant. The position of axial distance from nose 3D to 6.5D corresponds to the curve rising stage in Figure 33, as shown in Figures 36–39. At this stage, the asymmetry of the distribution of the circumferential pressure coefficient decreases, and the asymmetry is not obvious at position 5D. The distribution law of the circumferential pressure coefficient at 6.5D is obviously contrary to that before 5D. In addition, from the 4D position, the circumferential pressure coefficients of the model with Ra=100 μm begin to show obvious differences. After 5D from the nose, Cfzpm of the four models showed obvious differences. The roughness Ra=0 μm, 0.8 μm, and 5 μm all experienced multiple positive and negative Cfzpm changes, and the positive Cfzpm interval length was relatively short, and the peak value was low. Cfzpm of the Ra=0 μm model has a higher peak value at the second peak, which means the vortex intensity attenuation is slow. The Cfzpm peak values of Ra=0.8 μm and Ra=5 μm are both lower and close to 0 after the second peak value. Model 4 with roughness Ra=100 μm will produce a large unit lateral force coefficient (Cfzpm) in the z-square direction, which gradually decreases to near 0 at the tail. As shown in Figure 40, at position 8D away from the nose, the model with roughness Ra=100 μm still has obvious asymmetry. According to the integration of Cfzpm in the x-direction of the four models, the lateral force coefficients of the four models were obtained. Obviously, the integral value of Ra=100 μm was the largest, and the integral value of Ra=5 μm was the smallest, which was consistent with the results in Figure 33.Figure 33
Axial variation of lateral force coefficient per unit length.Figure 34
Pressure coefficient distribution of different models at1D.Figure 35
Pressure coefficient distribution of different models at2D.Figure 36
Pressure coefficient distribution of different models at3D.Figure 37
Pressure coefficient distribution of different models at4D.Figure 38
Pressure coefficient distribution of different models at5D.Figure 39
Pressure coefficient distribution of different models at6.5D.Figure 40
Pressure coefficient distribution of different models at8D.The circumferential pressure coefficient distribution of different sections of the same model is compared, as shown in Figure41. In Model 1 with Ra=5 μm, the curve asymmetry increases significantly from 1D to 3D, and the pressure coefficient on the left side is significantly lower than that on the right side. The pressure coefficient difference on both sides of the 3D section is the largest, and the pressure coefficient on the left side is higher than that on the right side in the 6D section. At 9D and 10D sections, the pressure coefficients on both sides reverse again and finally reach equilibrium at the 12D section. The change in Model 2 with Ra=0.8 μm and Model 3 with Ra=5 μm was consistent with Model 1 with Ra=0 μm, as shown in Figures 42 and 43, but the pressure distribution of Model 2 and Model 3 tended to be symmetric faster. Compared with Model 3, the pressure distribution asymmetry at the 6D section is more obvious in Model 2. The pressure coefficient distribution of Model 4 with Ra=100 μm is shown in Figure 44. After, the pressure coefficient distribution of the section from 3D to 6D changes from low left and high right to high left and low right, and there is no lower pressure coefficient on the left, which is consistent with the trend in Figure 33. In addition, the pressure coefficient on the left side of both 6D section and 9D section shows local drastic changes. According to the analysis of vorticity contour in Figure 45, in these two sections, the secondary vortex near the wall is the cause of local drastic changes in the surface pressure coefficient.Figure 41
The pressure coefficient of different sections varies with the roll angle at Model 1 withRa=0 μm.Figure 42
The pressure coefficient of different sections varies with the roll angle at Model 2 withRa=0.8 μm.Figure 43
The pressure coefficient of different sections varies with the roll angle at Model 3 withRa=5 μm.Figure 44
The pressure coefficient of different sections varies with the roll angle at Model 4 withRa=100 μm.Figure 45
Vorticity at sectionx=6D and x=9D in Model 4.


## 4. Conclusion

PIV experiment, continuous smoke flow experiment, force measurement experiment, and DES numerical simulation are carried out in this paper. The numerical simulation results of the DES method are compared with those of three experiments, and the flow field structure and calculation results are verified. It shows that the DES method can accurately simulate asymmetric flow with a high angle of attack at low speed. Based on the analysis of the experimental and numerical simulation results, the following conclusions were drawn in low Reynolds number and steady asymmetric flow with a high angle of attack:(1)
The downstream flow still appears asymmetry when the surface is smooth, indicating that the stability of the flow field is not enough to maintain the symmetry of the flow field(2)
There are one or more smaller size vortices due to the large velocity gradient at the side of the lower vortex. The small additional separation vortex produces and then merged with the main vortex after separation at the side of the higher vortex(3)
The adverse pressure gradient increases in the region where the boundary layer separates with the roughness greatly increasing and further affecting the separation point(4)
The head is the main part to generate lateral force, and the influence of roughness change on the head is slight. The influence of roughness is mainly to increase the fluctuation period of pressure along the axial direction. And the change in the pressure fluctuation period in the cylindrical section is the reason for the change in lateral force(5)
The secondary vortices near the wall will lead to drastic changes in local pressure coefficientsThe presence of roughness is closer to that of a real missile, and the results of this paper are useful for practical engineering applications.



---

*Source: 6590152-2022-07-27.xml*