Numerical Investigation on Aerodynamic Characteristics of Damaged Infinite Wings With Variation in Penetration Angle

: Aerodynamic performance of a full span NACA 64 1 -412 airfoil with a circular-shaped damage at various attack directions has been numerically investigated in this study. To assess the aerodynamic effects of different penetration angles in which threats such as projectiles can pass through the wings, attack directions of 30°, 60°, –30° and –60° relative to the normal axis of the chord line has been studied and compared with attack direction of 0°. To validate with published studies about damaged wing, the 200 mm chord airfoil was simulated with the damage hole diameter of 20% chord at the midspan and midchord location in Reynolds number of 500,000. Quantitative and qualitative results of this numerical study had a good agreement with published experimental data due to appropriate structured mesh and turbulence modelling. In addition to lift, drag and pitching moment coefficient, surface pressure distribution around the damage hole has been studied. Results show that, if the penetration angle becomes more negative, aerodynamics performance of the wing will be further decreased; therefore, attack directions of threat mechanisms such as “ahead and above” or “below from the rear” have severe negative impact than other directions on aerodynamic performance of the damaged infinite wing. improved aerodynamic characteristics of the damaged airfoil. Experimental and numerical analysis on NACA 64 1 -412 airfoil with triangular and star-shaped damage in Reynolds number of 180,000 was published by Etemadi et al . (2012); the results indicated that for each damage type, by increasing incidence, flow disturbance was increased with the damage jet effect. As a result, it had negative impact on aerodynamic characteristics. Furthermore, pressure coefficient magnitude in lower and upper surface was dependent on damage shape and jet velocity variation at the edge of damages. Render and Pickhaver (2012) experimentally investigated the influence of orientation on the aerodynamics of battle-damaged wings. Their wind tunnel tests were carried out on a NASA LS (1)-0417MOD airfoil with a circular hole simulating gunfire damage. It was mentioned that the negative obliquity orientations, where upper surface damage hole was closer to the leading edge and the lower hole to the trailing edge, had more unfavorable impacts on aerodynamic coefficients comparing to positive obliquity orientations. Adding skew damage where the upper and lower holes were offset in a spanwise orientation had a little effect on aerodynamic coefficients; however, it had an adverse effect when combined with greatest negative obliquity. with negative attack directions were contrariwise. Data from the present investigation in Reynolds number of 5e+05 showed good consistency with the results of Irwin and Render (2000) at the same Reynolds number. ω problem of Numerical point in the lower surface. In fact, usage of the self-supplying air-jet vortex generators on the airfoil delays the flow separation at higher incidences. It seems that +60° and +70° obliquity angles are applied as self-supplying air-jet vortex generators in higher incidences, undermined the damage jet flow, returned downstream flow to undisturbed surface flow velocities, and delayed surface flow separation.


DAMAGE MODELLING
The NACA 64 1 -412 airfoil was selected for the study in order to validate the numerical results with Irwin's published experimental works in the first place. The airfoil model had the chord of 200 mm and the span of 450 mm. This size selection was done to match with dimensions of the utilized wind tunnel test section and reduce the number of future experimental works.
The Reynolds number set for all simulations was 500,000; indeed, the flow velocity for the simulation was 35 m/s. Irwin et al. (1995) showed that the greatest effects for simulated gunfire damage occurred for damage located at quarter-and half-chord. Hence, the hole axis in the current study was at the wing midspan and half-chord position. Render et al. (2007) studied the influence of sharp edges as a star-shaped hole to simulate battle damage in a more realistic way and showed that, by assessing damage flow characteristics and changes in aerodynamic coefficients, using a circular hole is a reasonable representative to simulate gunfire damage if the circular hole diameter is close to the maximum width of the simulated damage. Thus, the circular-shaped damage has been opted for gunfire damage modelling in the present study and key assumptions were made that all simulated gunfire damage cases were modeled by a circular shape with 20% chord diameter through both upper and lower wing surfaces.
To allow a systematic study about the variation in penetration angle, it was considered that the damaged wings with different penetration angles have been simulated to compare with the one with an attack direction of 0° relative to the normal axis of the chord line. These models would have entry and exit circular holes at different chord-wise locations, i.e., with an attack direction of 30°, 60°, -30° and -60° relative to the normal axis of the chord line in Fig. 1. In models with positive attack directions, upper surface holes were closer to the airfoil's trailing edge and lower surface holes were closer to the airfoil's leading edge. The models with negative attack directions were contrariwise. Data from the present investigation in Reynolds number of 5e+05 showed good consistency with the results of Irwin and Render (2000) at the same Reynolds number.

NUMERICAL METHOD
To perform the numerical study, a commercial fluid dynamics program, Fluent TM , was used that utilizes a pressure based finite volume method within incompressible flows.

GRID GENERATION
As a Fluent's TM preprocessor, GAMBIT TM was used for geometry and grid generation. To have a precise modelling of viscous effects close to the solid boundaries, flow field has been divided by a multi-block structured C-type grid in Fig. 2. Structured grid improves the accuracy in capturing the flow parameters within the domain (Steinthorsson and Modiano 1995), more precision in modelling boundary layer regions and have a better estimation of lift and drag coefficients. Moreover, better control can be achieved by structured grid in some critical parts such as the leading edge, trailing edge, boundary viscous sublayer, strong pressure gradients, and vortical wake in the damage hole and wing upper surface. Therefore, fine mesh was used in some important zones including leading edge, trailing edge, damage hole and wake regions. Figure 3(a) depicts a sample view of computational grid near the damaged wing body in different attack directions while Fig. 3(b) shows the grid inside the circular damage hole for attack directions of +30°, -30°, +60° and -60°. Table 1 shows the number of mesh cells applied for each case in this study. It can be seen that the significant changes exist in total hexahedral cell numbers from 2,256,680 in undamaged airfoil to 3,995,720 in obliquity damaged airfoil with attack directions of +60°.   ,873,600 Obliquity damaged airfoil with the angles of -60°Th ree-dimensional grid independence study for the straight through damaged airfoil with attack direction of 0° is shown in Table 2. It indicates that aerodynamic coefficients are more sensitive to cell numbers in y direction than other directions; thus, more precision has been applied to this direction. To evaluate the grid quality, the value of dimensionless wall distance, so-called Y + , was estimated for the first cell height near the solid wall. Y + is nondimensional distance of first node from the wall, which depends on friction velocity and kinematic viscosity in the turbulent boundary layer (Schlichting and Gersten 2000). Clearly, by increasing the incidence angle and having a large separation region on the upper surface, maximum turbulent viscosity near the wall is increased, Therefore, the value of Y + on the most part of the upper surface is lower than 1, which meets the criteria of selected turbulence model that will be discussed in the following section. Y + contour is shown in Fig. 4 for incidence angle of 10.5° in the damaged airfoil with attack directions of -60°.

BOUNDARY AND INITIAL CONDITIONS
As shown in Fig. 5, three types of boundaries were assumed in this numerical study: 1) constant velocity at the inlet surfaces of the flow field (blue area in the domain), 2) constant static pressure at the outlet surfaces (red area in the domain), and 3) solid walls at both the upper and lower surfaces of the airfoil as well as the lateral surfaces of the span (white area in the domain). Computational flow domain dimensions, as shown in Fig. 2, were designated great enough to alleviate turbulences before reaching the domain boundaries. To be more accurate, in upstream direction, the domain was extended about 10 times the chord length from the wing leading edge while in downstream direction the extension was about 21 times the chord length from the wing trailing edge. The flow velocity is 35 m/s. Operating pressure and temperature of the flow are 101,325 Pa and 290 K, respectively. Conditions for the numerical investigation is defined in a way to be the same as the conditions in the Irwin's experiments (Irwin and Render 2000).
x z y Figure 5. Three types of boundaries with red, blue and white colors around the infinite wing.

NUMERICAL SOLUTION CHARACTERISTICS
Given that the flow is stationary, single phase, and turbulent without any heat transfer, Reynolds averaged Navier-Stokes (RANS) equations were considered as the governing equations of the numerical analysis. To model the turbulent viscosity term, the Menter's shear stress transport (SST) k-ω turbulence model was utilized in the present study. The SST k-ω turbulence model is a two-equation eddy-viscosity model; the two variables calculated are usually interpreted as "k" that is the turbulent kinetic energy and "ω" that is the rate of dissipation of eddies (Wilcox 2006). This method is applied because of its remarkable performance in estimating adverse pressure gradients and flow separating at the vicinity of the wall boundaries. The model combines the k-ω turbulence model and k-ε turbulence model such that the k-ω is used in the inner region of the boundary layer and switches to the k-ε in the free shear flow. The use of a k-ω formulation in the inner parts of the boundary layer makes the model directly usable all the way down to the wall through the viscous sublayer. According to Menter (1994), the SST model switches to a k-ε set of equation in the free-stream (far enough from the wall boundaries) to avoid the common k-ω problem of intense sensitivity to the inlet free-stream turbulence properties.
As shown by Kader (1981), enhancement wall functions were applied for the near wall areas because of having sufficed fine mesh to resolve the laminar sub-layer.
Pressure equation was discretized by the second order method and the momentum equation, turbulent kinetic energy equation and the turbulent dissipation rate equation were discretized by the second-order upwind method. Pressure and velocity were coupled by a SIMPLEC method. Absolute criterion for continuity, velocity, kinetic energy and dissipation rate is 1e-05.

NUMERICAL VALIDATION
Prior to study the damage cases, the validation of the current undamaged model data was investigated by comparing them with published data from Irwin and Render (2000) with the same profile and Reynolds number. Figures 6(a), (b) and (c) show the comparison of lift, drag and pitching moment coefficients for the present undamaged model and Irwin's data at 500,000 Reynolds number. Lift curve slope is 0.0833 for Irwin experimental result and 0.0866 for current numerical study, an error of 3.96%. A rapid drag increase are seen in high incidence angles for both results and the differences are small in pitching moment results. The last 1.5% of chord length around the trailing edge in Irwin's wing profile might be modelled inaccurately due to manufacture limitations, and it led to much more differences with current numerical results in high incidence angles.
Figures 6(d), (e) and (f) illustrates the comparison of lift, drag and pitching moment coefficients for the straight through damaged model in the current study with Irwin's data (1999). There are acceptable trend and value compatibility between experimental and numerical data in lower incidences while in higher incidences as the through damage jet becomes stronger, the numerical data slightly diverged from experimental data. The relative error between the numerical data and Irwin's experimental results is up to 10%. Lift coefficient error in high incidence angles can be originated from ignoring the laminar separation bubble in numerical solution and drag coefficient error is a consequence of wake influences in high incidence ranges.

STRAIGHT THROUGH DAMAGE
To quantify the damage effects, aerodynamic coefficient discrepancies caused by damage presence are indicated by subtracting the undamaged wing aerodynamic coefficient from the same coefficient of the damaged wing. Meanwhile, the reference area of the models in calculating the aerodynamic coefficients was defined by multiplying the wing chord to span dimensions. The aerodynamic coefficient discrepancies are calculated by Eqs. 1-3: The presence of flow passing through the wing results in decreasing the pressure difference between the lower and upper surfaces of the airfoil that can lower the aerodynamic performance. Figure 7 exhibits the lift coefficient discrepancy in current numerical study. The zero-lift coefficient incidence is near -3° (Fig. 6a) showing that the flow passing through the damage is negligible at this incident angle. For incidences lower than -3°, the through flow direction is reversed from the upper to lower surface; therefore, higher lift coefficient increments were observed. In positive lift incidences range, the through flow direction is from the lower surface to the upper surface. In such cases, by increasing the incidence, more lift loss is occurred. As shown by Irwin and Render (2000), in the lower positive angle of attacks, the small part of flow passed through damage produces a weak jet, which leads to a smaller affected area behind the damage hole than in the higher positive angle of attacks. Figure 8(a) depicts the weak jet structure in the upper surface flow for straight through damage at 0°. Surface flow-structures of damage are similar to flow observations of Irwin and Render (2000). At the front edge of the damage hole where the surface cross flow and damage jet intersect, an adverse pressure gradient causes a separation region with two lines shown by forward "1" and secondary "2" separation lines. These lines are the boundaries of a vortex that forms a greater horse-shoe vortex "3" around the damage hole determining the extent of the affected area behind the damage hole. At the rear edge of the damage hole, there are two contra-rotating vortices "4", moving upward around the damage edge strengthening the jet due to increasing in the incidence. However, these two contra-rotating vortices are symmetrical around the hole center line in all incidence ranges.
As the incidence is increased above 0°, the pressure difference between the upper and lower surface is increased, thus, a higher flow mass passes through the damage. Figure 8(b) shows this type of flow structure called "strong jet" in 8.5° incidence. In comparison with the weak jet pattern, the forward "1" and secondary "2" separation lines in the strong jet are located farther to the damage front edge. This can be explained as the ability of strong jet to push the upstream flow further comparing to the weak jet. This chord wise movement of separation lines continues with a larger span wise horse-shoe vortex "3" and a greater reverse flow region "4" that comes from a strong jet with more penetration in the freestream. Two small contra-rotating vortices "5" at the rear edge of damage are located far from each other relative to their previous locations in the weak jet. Moreover, there are two large contra-rotating vortices "6" near the trailing edge; the occurrence is due to the fluid entrainment into the damage wake from the surrounding upper surface flow in combination with lower wing surface flow at the trailing edge. Figure 9 shows the contours of velocity magnitude from the lateral view for both weak jet and strong jet structures at mid span station. As discussed above, due to the pressure difference between upper and lower surfaces, flow passes through damage and normally penetrates into upper surface flow behind the damage hole. However, normal penetration of the weak jet in 0° incidence shows a smaller vertical height than higher positive angles of attacks (Irwin 1999). In 8.5° incidence, the strong jet emerging from the damage is intersected with upstream flow and fully detached from the upper surface making a larger separation region between them. It starts immediately behind the rear edge of damage and extending up to trailing edge. In this region, there is significant region of three-dimensional reverse flows mainly due to drawn flow from the lower surface at the trailing edge. 6.66e+01 6.33e+01 5.99e+01 5.66e+01 5.33e+01 4.99e+01 4.66e+01 4.33e+01 4.00e+01 3.66e+01 3.33e+01 3.00e+01 2.66e+01 2.33e+01 2.00e+01 1.66e+01 1.33e+01 9.99e+00 6.66e+00 3.33e+00 0.00e+00 6.65e+01 6.32e+01 5.98e+01 5.65e+01 5.32e+01 4.99e+01 4.65e+01 4.32e+01 3.99e+01 3.66e+01 3.32e+01 2.99e+01 2.66e+01 2.33e+01 1.99e+01 1.66e+01 1.33e+01 9.97e+00 6.65e+00 3.32e+00 0.00e+00 3 3 1 1 2 2 5 6 4 4 weak jet structure in 0° incidence (a) (b) strong jet structure in 8.5° incidence Putting a magnifier inside the damage hole at the midspan station in Fig. 10 shows that the through flow mainly exits from the rear portion of the damage hole. The exit flow cross section is less than half of the damage hole cross section even in higher incidences. Due to the pressure difference between surfaces, flow streamlines are drawn into the damage hole in the lower surface, collide with rear wall of the damage hole make the highest-pressure region. At this stagnation point, the flow divides into two parts: part of the flow turns around the bottom edge of the rear wall and continues its way to the trailing edge in the lower surface while the other, being the main part of the flow, passes along the damage rear wall with lower static pressure up to the hole exit. Above the damage hole, there is a high static pressure area due to collision of upstream flow and the damage jet flow which leads to another flow division itself.
After the collision, part of the flow continues to join damage wakes, but the other part changes its direction into the damage hole and makes two groups of contra-rotating vortices in the hole. The lower group of vortices rotates counter-clockwise while the upper group rotates clockwise. Figure 11 illustrates how these two groups of vortices are formed inside the damage hole and the red-lines indicate the cores of them. It clearly shows how the streamlines exit from the hole. At this point, one group of vortices on top have the same velocity direction to its adjacent group of vortices at the bottom with a tiny high velocity area between them as indicated with light blue color in Fig. 9 (b). As the angle of attack increases up to 8.5° incidence ( Fig.10 b), the pressure difference between surfaces increases, the jet passing through the damage hole becomes stronger and its cross section becomes wider inside the damage hole. Consequently, this strong jet pushes the two groups of vortices toward the front wall as well as draws them upward by itself.  As shown in Fig. 7, between -4.5° and 0° incidence, approximately, constant lift loss is seen because the flow disruption on the upper surface of the model is at the minimum level in the weak jet pattern. The sharp inclination in the rate of dC l is started from 0° onwards due to the transition of the flow from weak jet to strong jet. In the strong jet region, the more incidence increases, the more lift loss occurs.
Correspondingly, Fig. 12 indicates that the presence of damage leads to drag coefficient increments in all incidence ranges; after a slight slope in weak jet pattern between -4.5° to 0°, above 0°, the rate of dC d changes increases sharply with incidence while at 10.5°, the difference between damaged and undamaged data was at the minimum level in all of the strong jet region.
As shown in Fig. 13, at higher incidence with damage presence, the nose down pitching moment increment occurs. In other words, if the incidence has bigger positive values, the rate of change of dC m will have higher negative values, especially above 2.5° incidence where strong jet region exists. This negative trend of pitching moment increments can be attributed to the difference in the surface pressure coefficient differential at the forward and aft ward of the wing. Figure 14 exhibits the pressure coefficient versus chord locations in x coordinates for damaged and undamaged airfoils at 8.5° incidence in different span-wise locations. As the flow over undamaged model is two-dimensional, the chord-wise pressure distribution is only illustrated for the midspan station (center-line). For damaged case, the chord-wise pressure distributions are shown for five span-wise stations; the damage center-line, 0.5R, 1.5R, 2.5R and 5R (Fig.15); where R is the radius of the circle shaped damage.  For undamaged wing, the pressure difference between the upper and lower surface at the damage location can define the extent of through flow from the lower to the upper surface with the presence of damage (Fig.15). For damaged wing, it generally shows a reduction in absolute value of pressure coefficients for both upper and lower surfaces, except at some regions at the rear part of the upper surface. Figure 15 also depicts that the presence of damage does not change the stagnation point position at different span-wise stations. In the front part of the damaged wing, the pressure coefficient reduction in upper surface suction peak has some discrepancies in different span-wise stations.
The center line and 0.5R damaged profile are almost close to each other with a small back and forth, confirming the occurrence of through flow in the damage hole. However, it should be considered that, compared with 0.5R, the upper surface flow at 1.5R station meets the through flow from the damage hole with some delay and, on the other hand, the rear edge of the damage jet exits sooner at 1.5R span-wise station. The 1.5 to 5R stations respectively shows the C p profile as getting closer to undamaged case, although, none of them are not fully recovered towards the undamaged C p profile. In the lower surface prior to the forward damage edge, the pressure coefficient is decreased since the flow is accelerated and suctioned into the damage hole. At the rear edge, damage hole wall blocks the through flow, which initially leads to a sudden increase of the pressure coefficient, and, after the stagnation point, a small amount of flow opens its way to the lower surface with acceleration and the rest of it accelerates and overflows from the upper surface toward the damage wake.
Downstream of the damage, almost none of the profiles fully recover the undamaged values up to trailing edge; however, 2.5R and 5R span-wise stations are closer to undamaged profile than the others. At 8.5° incidence, strong jet characteristics lead to the upper surface flow separation behind the damage that makes an increase in the pressure coefficient negative values. After that, the C p negative values are also greater than the undamaged profile up to trailing edge. Behind the damage in the lower surface, all the span-wise stations depict a slight rise in positive C p values toward the trailing edge. In a small region at the trailing edge location, fluid from the lower surface is drawn to the upper surface reverse flow with acceleration and a decrease in C p values might happen.
For the damaged wing, the differential pressure coefficient between the upper and lower surface at the rearward of the wing is more negative than the same calculation at the forward of the wing. It results into positive lift at the rearward and negative lift at the forward of the wing, which is responsible for nose down pitching moment in the damaged wing.
All the aforementioned quantitative and qualitative results are consistent with experimental results of Irwin and Render (2000), such as damage effect on aerodynamic coefficients, pressure distribution and upper surface oil flow visualization for the same NACA 64 1 -412 airfoil with circular shape quarter chord damage and the same Reynolds number.

THE INFLUENCE OF OBLIQUITY
All obliquity cases in low incidences show the typical characteristics of weak jet, which is a small scale through flow in damage holes due to the less differential pressure between surfaces in low incidences. To acquire a better understanding of variation in penetration angle, strong jet characteristics for each case at 8.5° incidence are shown for the upper surface flow as a velocity path line in Fig. 16. It shows that obliquity angle has a significant influence on the quality of flow passing the damage hole and the surface flow structure is affected by the damage. This collection of obliquity angles depicts that as the obliquity angle becomes more negative, the damage jet strength and its wake extent would tend to rise. According to Fig. 15, the discrepancy between the locations of damage hole centers at the undamaged upper and lower wing surfaces can be related to the pressure coefficient.
As described in the previous section about strong jet characteristics for straight through damage in 8.5° incidence, for obliquity angles, when the upper surface flow meets the damage jet, its velocity is decreased, and its pressure is increased. The behavior leads to an adverse pressure gradient forming the primary and secondary separation lines ahead of the damage hole. The jet flow is pushed upstream and generates surface flow vorticity behind the primary separation line. This vorticity starts to form a great horse-shoe vortex around the damage that eventuates to the two major contrarotating vortices as a result of the entrainment of free stream and the lower surface flow into the damage wake. The horse-shoe vortex is like an isosceles triangle with its head at the forward of damage edge and the equilateral sides at the boundaries between damage wake and free stream, out of the affected area of the damage.  Figure 16. The path line colored by velocity magnitude for different obliquity angles at 8.5° incidence.
As shown in Fig. 16, the horse-shoe vortex is larger in negative obliquity angles than positive ones at the same obliquity magnitude. It explains that the through damage with the negative obliquity angle has a wide extent influence on the upper surface flow structure, located behind the damage hole with greater affected zone. Some of these influences are: bringing the primary and secondary separation line closer to forward damage edge, creating greater contrarotating vortices at the trailing edge, making small contrarotating vortices farther from each other around the rear damage edge, and promoting stronger reverse flow from the lower surface around the trailing edge that reduces the flow velocity in the damage wake at the rear area of the wing.

±60° OBLIQUITY ANGLE
To compare the positive and negative obliquity angles, Figs. 17 and 18 show the contours of velocity magnitude and static pressure with damage hole path lines from lateral view at 8.5° incidence for +60° and -60° obliquity angles. Due to the pressure difference between the upper and lower surface, the flow passes through the damage from lower to the upper surface. After stagnation point at the leading edge, there are clear high velocity zones on the upper surfaces of the airfoils. When upper surface flow meets the damage through flow in front of the damage edge, the upcoming flow velocity decreases and separation lines occur (Fig.16 d and e). Moreover, flow in lower surface ahead of the damage is drawn into the damage hole with a velocity raise around the damage edge, especially in +60°. In -60° obliquity angle, the direction of through flow turns into the leading edge and has a strong penetration into the upstream flow making a great separation region with a large-scale low speed reverse flow underneath. Therefore, the affected area behind the damage jet is vertically and horizontally wider in -60° obliquity angle than +60° obliquity angle. In the lower surface of +60° behind the damage, there is only a small low speed zone due to collision of streamlines to the rear edge; however, in -60° obliquity angle, the stagnation point is inside the damage at the rear wall dividing the flow into two parts; one part emerges from the upper surface as a strong jet, the other part exits from the rear edge in the lower surface while its velocity increases toward the trailing edge and after that turns into the reverse flow in the upper surface. Figure 18 (a) shows that the static pressure for +60° obliquity angle increases in two visible zones; one in the upper surface at the separation point in front of the damage and the other for the stagnation point at the damage rear edge in the lower surface. As it is seen, the jet flow passes through the damage toward downstream direction without any significant changes in velocity except near the front wall inside the damage at the lower part where jet velocity decreases after the flow acceleration around the forward edge and trapped in the form of a small vortex same as a recirculation region in the backward step flow.
For -60° obliquity angle, the stagnation point at the rear wall lower part of the damage determines the pressure rise at the intersection of streamlines with the damage wall; after that, the through flow changes its direction along with the rear wall toward the exit hole. When the damaged jet penetrates the freestream, it collides with the upcoming flow results in increased pressure and divides the flow into two parts: a small scale part returns into the damage hole made two groups of contra-rotating vortices in the damage hole; the other great part changes the direction again and continues its path toward the trailing edge with a wide extent separation region beneath. The upper group of vortices inside the damage rotates clockwise being originated from the upper surface streamlines and the lower group of vortices rotates counter-clockwise due to originating from lower surface streamlines. Having the same direction for streamlines between these two contra-rotating vortices results in a velocity increase as shown in Figure 17 (b) (light-blue color zone).
All span-wise pressure distribution results for obliquity angles show that centerline and 0.5R stations are close to each other for all obliquity cases as well as 2.5R and 5R stations for positive obliquity angles. Hence, for fully coverage of the chord-wise pressure distribution in different obliquity angles, centerline, 1.5R and 2.5R stations are chosen to illustrate. Figure 19 (a), (b) and Fig. 20 show the surface pressure distribution for ± 60° at centerline, 1.5R and 2.5R span-wise stations. For a better comparison, undamaged and straight through damage wings data are also added to charts at each span-wise station. 1.5R Figure 19. Chord-wise pressure distribution at two different span-wise stations for ±60°, 0° obliquity angles and undamaged wings for 8.5° incidence.  Figure 20. Chord-wise pressure distribution at 2.5R station for ± 60 °, 0° obliquity angles and undamaged wings for 8.5° incidence. Figure 19 (a) shows that all three damage cases in the centerline station of the upper surface experience a reduction in their negative pressure peak in front of the damage edge such that pressure coefficient decline in -60° obliquity angle has the steepest slope in comparison to straight through damage (0° obliquity) and + 60° obliquity angle. It is also noted that the front damage edge of obliquity -60° case in the upper surface is the nearest to the leading edge (x/c = 0.27) than obliquity 0° (x/c = 0.4) and obliquity +60° (x/c = 0.61) cases.
Pressure coefficient variations for the front wall of the damage in obliquity -60° case are defined in the chart as an area with ups and downs between x = 0.053 and x = 0.092. C p values of the front wall inside the damage hole are being affected by two groups of contrarotating vortices, which leads to a reduction in the C p values adjacent to them and an increase in the C p values in the area between them close to the front wall (Fig. 18 b). In the centerline station of lower surface, due to drawn flow around the damage edge, C p positive values decrease up to the damage edge; there are slight discrepancies between obliquity cases in the pressure decreasing process in this region except +60° obliquity case between x = 0.066 and x = 0.122. The sudden negative pressure peak relates to its front wall inside the damage; in this region in +60°, due to the high pressure difference between the center of lower and upper surface holes in +60°, flow is suctioned into the damage and its C p values suddenly rise up to negative values, which made a recirculation region adjacent to the front wall (Fig.18 a) similar to back-step flow in Westphal et al. (1984). After this separation, when flow approaching to the reattachment point, its negative pressure value gradually decreases.
Downstream of the damage holes in the upper surface, the onset of the strong jet from the rear edge of each obliquity case varies with the location of the obliquity cases rear edge, consequently the location of maximum pressure downstream of the damage cases are in consistent with the location of the exit jet in any case. Whatever the jet exit is closer to the trailing edge in any case, the negative pressure peak is nearer to the trailing edge as well. Inside the +60° obliquity case, the rear wall is between x = 0.11 and x = 0.138. Due to flow collision with the damage rear edge in the lower surface, C p values increase to positive values. After this collision area, the flow divides into two parts: in the lower surface, the flow accelerates toward to the trailing edge and static pressure decreases up to the trailing edge, for the other part in the upper surface, flow accelerates toward the jet exit adjacent to the rear wall inside the damage hole and suction pressure decreases up to the rear edge at the upper surface and after that flow emerges and separates and C p value increases.
Downstream of the damage in the lower surface, the C p profile depicts that the order of pressure reduction versus undamaged profile is +60° obliquity, 0° obliquity and the last one -60° obliquity case, respectively. For -60° obliquity case, the rear wall is between the x = 0.092 in the upper surface and x = 0.127 in the lower surface. Oncoming flow from the lower surface collides with the rear wall inside the damage at x = 0.12 that increases C p to positive values. In the following, a part of this-high pressure flow turns around the rear edge and emerges into the lower surface as separated and its velocity firstly decreases due to separation and after that increases toward the trailing edge. The other part of the flow divides from the collision point, goes toward the exit hole adjacent to the rear wall, simultaneously the C p values decrease from positive values and increase around the exit edge to the negative values.
At 1.5R distance from the centerline span-wise station, the C p values reduction trend also continued due to damage effects. -60° obliquity angle shows a greater reduction in negative values in front of the damage. The suction peak in the upper surface for each case is almost parallel to their rear part of damage hole where the strong jet exits. As it is seen in downstream of the damage cases, C p has almost the constant values due to its presence in the separation region. There is a slight difference between damage cases in the lower surface, but the location of the minimum positive C p values for each case are around their rear part of damage holes.
At the 2.5R span-wise station away from damage location, +60° obliquity angle profile depicts the most recovery of C p profile to the undamaged case in the whole chord length. It means that flow at this station in +60° returns to undisturbed surface flow C p values. Moreover, +0° obliquity angle profile has the second most recovery of C p profiles among the cases and -60° obliquity angle is the profile that has the biggest difference from the undamaged case.
Generally, Figs. 19 and 20 illustrate that, in all different span-wise stations, the -60° obliquity angle has more significant effect on C p profile and, subsequently on wing performance versus +0° and +60° obliquity angles. In all obliquity cases, by moving away from the damage location in span-wise direction, the extent of similarity to the undamaged C p profile will be increased. Not only has the +60° obliquity angle case the closest C p profile to the undamaged wing in all different span-wise stations, but also its trend in different angle of attacks is inconsistent comparing to other discussed cases.
According to Fig. 21 (a), (b) and (c), which indicate lift, drag and pitching moment, coefficients changes versus angle of attack and Irwin's experimental results (1999), and also, as discussed above, in velocity and static pressure contours, two distinct region can be observed: weak jet and strong jet. The lower incidences between -4.5° to 0° are the weak jet region for almost all cases with a mild slope in the change of aerodynamic coefficients. The absolute value of dC l from -4.5° onward is gradually increased with two different slopes in weak and strong jet region. Figure 21 shows that the strong jet area started after 2.5° incidence for all damage cases except for -60° obliquity case, because its aerodynamic coefficient changes seem to be occurred from 0° onward. The greater magnitudes of the coefficient discrepancies and the faster rate of lift loss, drag increase and negative pitching moment after 2.5° reflect the strong jet characteristics discussed in the qualitative results. The two main parameters of incidence angle and damage holes location define the ratio of freestream to through flow, which influences the type of jet present. In the cases of jet in cross flow (Andreopoulos and Rodi 1984), there is no interdependency between freestream and the jet velocity; however, these two values could not be set separately from the other in the damaged wing cases.
Quantitative results show that +60° obliquity case has the least negative effect on aerodynamic performance as shown in Figs. 17 (a) and 18 (a) in terms of qualitative results. Moreover, there is a slight sensible effect in pitching moment coefficient increments (Fig 21c).
For better assessment of the +60° trend, another positive obliquity angle slightly greater than +60°, such as +70° obliquity angle, is considered. Figure 22 shows the upper surface flow structure of +70° obliquity case. As can be seen, at 4.5° incidence, through flow emerges from the hole as a strong jet with separated flow downstream of the damage. At 6.5° incidence, since the wake size is obviously smaller, the exit jet is defined as a weak jet with attached wake downstream of the damage. The flow downstream of the damage almost remains as a weak jet in greater incidences such as 8.5° incidence. The cause of jet strength reduction in higher incidences could be similar to the performance of mechanism of self-supplying air-jet vortex generators in Krzysiak (2008). They remain inactive at low angles of attack and only become activated at the higher angles of attack. In higher angles of attack, the pressure difference between upper and lower surfaces increases, especially near the stagnation point in the lower surface. In fact, usage of the self-supplying air-jet vortex generators on the airfoil delays the flow separation at higher incidences. It seems that +60° and +70° obliquity angles are applied as self-supplying air-jet vortex generators in higher incidences, undermined the damage jet flow, returned downstream flow to undisturbed surface flow velocities, and delayed surface flow separation.  Figure 22. path line colored by velocity magnitude for +70° obliquity angle at 4.5°, 6.5° and 8.5° incidences.

±30° OBLIQUITY ANGLE
As mentioned above, to assess qualitatively the flow characteristics in positive and negative obliquity angles, it was decided to compare the aerodynamic performance of +30° and -30° obliquity angles at 8.5° incidence where damage jets were of the strong type. Figures 23 and 24 indicate the contours of velocity magnitude and static pressure with damage hole path lines from lateral view. As discussed, the upcoming flow on the upper surface crossed with damage through flow and its velocity decreased; flow in the lower surface ahead of the damage edge suctioned into the hole due to the pressure difference between surfaces; the jet flow emerged only from the rear parts of the damage hole and in front parts of the hole, two groups of contra-rotating vortices have been made with a high velocity area between them. Figure 24 (a) and (b) show that there are two regions with high pressure in +30° and -30° obliquity angles; firstly, due to the collision of through flow inside the damage hole with the bottom part of the rear wall; secondly, due to the jet flow collision and upcoming flow above the damage hole. However, in -30° obliquity angle, these high-pressure areas are wider and stronger than +30° obliquity angle. In -30° obliquity angle, at the lower surface of the rear edge, high-pressure flow turns around the edge, becomes separated and then accelerates toward the trailing edge. At the upper surface of the rear edge, flow emerges from the damage edge and penetrates to a broad separation area with a large amount of reverse flow beneath it. To put the upper and lateral view of +30° and -30° obliquity angles into a close perspective, it can be seen that vertical height and span-wise width of -30° obliquity wake is greater than +30° obliquity wake.   Fig. 26 show the surface pressure distribution in 8.5° incidence for undamaged, straight through damage and ±30° obliquity angles, at centerline, 1.5R and 2.5R span-wise stations. The previous discussion in the surface pressure distribution of straight through damage and ±60° could be applicable to ±30° obliquity angles. In Fig. 25 (a), at the centerline station, the order of C p values reduction ahead of the damage hole for different cases versus undamaged case is +30°, 0° and -30° obliquity angles, respectively.
As expected, according to previous results for ±60°, -30° obliquity angle has the most surface pressure decrease in front of the damage edge; its front wall inside the damage hole is between x = 0.071 in the upper surface and x = 0.084 in the lower surface; in between its C p values experience constant magnitude part due to neighboring with a small recirculation region at the top. In the following, due to neighboring with two groups of vortices adjacent to the front wall, C p values get initially increased and its velocity decreased near the top vortex and then the velocity gets increased and C p values decreased near the bottom vortex.
In the lower surface ahead of the damage hole, there is a slight difference in pressure reduction rate between different cases. The front wall of the +30° damage hole is between x = 0.075 in the lower surface and x = 0.089 in the upper surface that C p values gradually increase next to the front wall adjacent to the lower groups of vortices where the flow velocity decreases and after that C p values is fixed due to the small recirculation region near the forward edge at the top.   Figure 26. Chord-wise pressure distribution at 2.5R stations for ±30°, 0° obliquity angles and undamaged wings for 8.5° incidence.
Downstream of the damage, the rear wall of the +30° damage hole is between x = 0.117 in the lower surface and x = 0.127 in the upper surface that in bottom part due to the collision of flow to the rear wall, the positive increase in C p values is observed and, after that, due to flow acceleration toward the trailing edge, C p values become decreased. Inside the damage hole at the rear wall in +30°, the blockage effect of the passing flow results in decreasing C p values from the large positive value at the bottom part near the stagnation point to the large negative value at the top part near the jet exit, hence through flow emerges from the rear edge with the negative pressure reduction and makes the separated wake behind the damage.
The rear wall of the -30° damage hole is between x = 0.111 in the upper surface and x = 0.122 in the lower surface, the maximum C p values at this distance belongs to the flow collision point to the rear wall and after this point flow is divided into two parts; one part is suctioned around the rear edge at the bottom part of the hole then firstly separates in the remainder of the lower surface and secondly accelerates toward the trailing edge, the other part of the flow continues its way as a through flow adjacent to the rear wall and its pressure firstly decreases from the maximum positive value at the stagnation point (Fig. 24 b) and secondly increases to the negative values near the exit edge on the upper surface then emerged flow adds to the separated wake downstream of the hole, effect of separation on the C p profile appears as a negative pressure peak in the wake region. Figure 25 (b) depicts different C p profiles cases at 1.5R span-wise station. As evident on the upper surface, more negative obliquity leads to more pressure peak reduction in front of the damage hole and the more negative pressure value in the damage wake. For each obliquity case, the minimum magnitude of C p values on the upper surface aligns with the latter 30% of the damage hole diameter length where the jet exits. C p profiles for all the cases in Fig. 26 shows that the span-wise extent of damage flow effects still exists at the 2.5R span-wise station. +30° obliquity angle has also the close proximity to undamaged C p profile and -30° obliquity angle has the most different C p profile versus undamaged case. Different cases in the lower surface at 1.5R and 2.5R span-wise stations are quite similar.
As discussed above, -60° and -30° obliquity angles decreases the aerodynamic performance versus its equal positive obliquity angles. For better assessment of the results, it would be rational to compare the +60° and +30° obliquity angles and also -30° and -60° obliquity angles. Quantitative results depict that the closer the obliquity is to 0° in positive angles and the obliquity farther from 0° in negative angles, the more aerodynamic performance reduction has occurred. The upper view of all obliquity angles in figure 16 confirm that the disruptive influence of negative obliquity angles are much more than positive ones.
To compare the +60° and +30° obliquity angles, it was decided to survey chord-wise turbulence intensity slices over these two obliquity cases at 8.5° incidence and measure the damping distance of their turbulence intensities. Figure 27 illustrates the chord-wise turbulence intensity slices over these two obliquity cases at 8.5° incidence. Despite the fact that jet flow starting point in +60° obliquity angle is almost 6% of chord length closer to trailing edge versus the +30° obliquity angle, but turbulence intensity damping point for +60° obliquity angle is 3.75 times of the chord length behind the wing while this distance is at a further point for +30° obliquity angle such as 4.3 times of the chord length. This is also another reason that +30 ° obliquity angle led to stronger jet, which had more negative impact on aerodynamic performance.
To compare the -60° and -30° obliquity angles, as an example, Fig. 28 shows the chord-wise pressure distribution of -30° and -60° obliquity angles at 2.5R span-wise station for 8.5° incidence. It is also clear that -60° obliquity angle is the most different C p profile from the undamaged case upstream and downstream of the damage hole location. After that -30° and 0° obliquity angles, respectively, gain the second and third place versus undamaged C p profile. . Chord-wise pressure distribution at 2.5R station for -60°, -30° , 0° obliquity angles and undamaged wings for 8.5° incidence.

CONCLUSION
Numerical investigation on NACA 64 1 -412 airfoil has been done for the first time in order to assess the effects of variation in the projectile penetration angle. This airfoil has been chosen in order to validate the numerical methodology with published experimental data about its undamaged case and straight through damage one. The numerical results of this study estimated the aerodynamic coefficients in an acceptable agreement with previous experimental results due to appropriate structured mesh and turbulence modelling. Damage presence generally reduced the infinite wing aerodynamic performance; increasing the angle of attack for all damaged cases resulted in decreasing lift coefficient, increasing drag coefficient and more negative pitching moment. The amount of flow passing through the wing that was alongside of the vortices, which is formed inside the damage hole, played the crucial role in results.
The current numerical study illustrated the interactions between the through flow and vortical pattern inside the damage hole and downstream of the damage. It has been numerically shown in NACA 64 1 -412 airfoil that, as the obliquity angle became more negative, aerodynamic performance decreased. Conversely, whatever the obliquity angle was farther from the vertical straight line with positive angle, the damage effects on aerodynamic performance decreased as the chord-wise turbulence intensity slices also confirmed that.
According to current verified numerical results on NACA 64 1 -412 airfoil, if a threat mechanism attacks the midspan of current wings from the "ahead and above" or "below from the rear" directions, it may have more negative impact on the aerodynamic performance than other directions. However, further investigation about different span-wise locations in three-dimensional models of current wings is necessary, especially at wing root and wing tip locations.
Different span-wise stations of pressure distribution also confirmed the results described above; damage presence led to a reduction in the chord-wise pressure difference between lower and upper surfaces. However, this negative effect gradually decreased with getting away from damage in span-wise direction. In more negative obliquity angles, damage adverse effects continued farther away along the span which indicates the greater width of the affected area by negative obliquity angles in comparison with the smaller width for positive obliquity angles.

FUNDING
There are no funders to report for this submission.