Numerical Investigation of Oriﬁce Nearﬁeld Flow Development in Oleo-Pneumatic Shock Absorbers

: The ﬂow ﬁeld development through a simpliﬁed shock absorber oriﬁce geometry is investigated using a single phase Large Eddy Simulation. Hydraulic oil is used as the working ﬂuid with a constant inlet velocity and an open top boundary to allow the study to focus on the free shear layer and the ﬂow development in the vicinity of the main oriﬁce. The ﬂow ﬁeld is validated using standard mixing layer dynamics. The impact of the oriﬁce shape is discussed with regards to the initial free shear layer growth, boundary layer development and the potential appearance of cavitation bubbles. Observations are made regarding the presence of ﬂow ﬁeld disturbances upstream of and through the oriﬁce, thereby, leading to a notable turbulence intensity level in those regions. A.A.S.A.-S.; writing— review and editing, B.G., D.V., P.T., A.F.A. and M.S.; visualization, A.A.S.A.-S., B.G. and D.V.; supervi-sion, P.T., A.F.A. and M.S.; project administration, M.S.; funding acquisition, M.S. All authors have the


Introduction
Aircraft landing gear shock absorbers are used to protect the main aircraft structure from the loads and forces during the landing impact and other ground manoeuvres. Most large commercial airliners use oleo-pneumatic type shock absorbers, which have been the most popular type on jet aircraft for decades. These usually include hydraulic oil being forced through a main orifice creating hydrodynamic resistance to the motion and a pressure differential across the orifice in addition to some gas (often nitrogen), that is either in direct contact with the hydraulic oil or separated in its own chamber to provide a restoring force to an equilibrium position under static loading conditions on the ground.
The hydraulic resistance of fluid moving through the orifice of the shock absorber has been identified by previous studies as the main energy absorption mechanism during the operation of oleo-pneumatic shock absorbers [1,2]. This hydraulic resistance is directly associated with the free shear layer developing downstream of the main orifice in the shock absorber, which contributes to the impact absorption by dissipating the bulk flow energy through the turbulence cascade, as viscous heat from the smallest scales of turbulent eddies. Hence, the orifice design and its impact on the downstream near-field flow development is of particular importance when considering the internal dynamics of shock absorbers and their impact on performance.
Previous investigations into oleo-pneumatic shock absorber performance mostly focused on the use of dynamic systems models to simulate the system [1,3], while others concentrated on the estimate of a particular parameter, such as the discharge coefficient of the main orifice [4]. The task of predicting the performance of an oleo-pneumatic shock absorber usually requires simplifying assumptions due to the complexity of the field and the multiphysics nature of the problem. Furthermore, interest in the overall integral performance often leads to the adoption of lower fidelity methods and, hence, reduces the details available on the internal shock absorber dynamics.
The present study aims to investigate a specific feature of the underlying flow physics inside a typical shock absorber by conducting a single phase simulation through representative oleo-pneumatic shock absorber geometry, which was closely based on the geometry used by Milwitzky and Cook [1] in their drop down experiments. The flow development through the orifice is discussed using a range of flow variables to help identify some key elements of the design and their effects on the flow field.

Numerical Methods
OpenFoam (v2006) [5,6] was used to conduct time-dependent scale resolving simulations. OpenFoam uses an unstructured finite volume approach to simulate fluid flow, which facilitates the handling of complex geometries that are expected in industrial type applications, as in the case of the shock absorbers investigated here. Field properties are stored at the centroid of each cell and then calculated at each cell surface "face" by interpolating the values at the centroids of the two cells sharing that face. Conservation laws are then applied based on theses face values, in accordance with the chosen discretization scheme.
The PIMPLE algorithm [7] is selected in the present simulations, as it provides good stability and convergence control options. PIMPLE utilises an outer loop based on the PISO algorithm [8] to capture the field development in time, while an inner loop based on the SIMPLE algorithm [9] solves each field instance as a steady problem. The interpolation method used in the present study is a second order linear interpolation, which is a standard choice in many finite volume solvers. The convective term is discretized using a second order limited linear scheme, which provides a balance between accuracy and stability through limiting the nonlinear term. While the diffusive term is discretized using a second order linear scheme and temporal integration is conducted using a second order Crank-Nicolson method.
After the discretization process, the resulting linear system for the pressure field is solved for using a Preconditioned Conjugate Gradient (PCG) method, while the velocity field system is found using a Preconditioned Bi-Conjugate Gradient method (PBiCG). An adjustable time step is used to insure that the Courant-Friedrich-Lewy (CFL) number remains below 1.0 throughout the simulation, to aid in achieving a balance between efficiency and stability during these simulations. The WALE subgrid scale model [10] is used to model the eddy viscosity of the subgrid scales in the Large Eddy Simulation, (LES) framework. The flow field is allowed to develop by allowing multiple through-flows across the entire domain, then time averaging is conducted until statistical convergence is reached for all turbulent properties of interest.

Methodology and Computational Setup
The baseline orifice geometry, illustrated in Figure 1, is essentially similar to that of a convergent nozzle, with the throat being at the orifice exit. The rounded outer (upstream) edge of the orifice helps avoid any boundary layer separation on the orifice's inner surface, while the accelerating flow through the favourable pressure gradient through the orifice allows relatively strong flow turning in the flow direction so that the flow can be tangential to the axial (streamwise) direction at the orifice exit. An axisymmetric three dimensional simulation is conducted on a simplified geometry of the shock absorber that focuses on the orifice tip region and the developing shear layer in its vicinity.
Grids of varying resolution were generated and tested to establish grid independence and optimise the setup where possible. The resolution in the shear layer region was tested using a preliminary simulation on a grid with 185,000 nodes. An analysis was conducted to check the required resolution in the main simulation, including looking at the estimated mixing layer thickness near the orifice tip to determine the grid resolution required there. The boundary conditions are illustrated on Figure 1, while further information about the case is provided in Table 1.  The inlet boundary condition is taken to be a constant vertical flow through the shock absorber V = 1.0 m/s as a representative figure of an approximately time-averaged value over the stroke period from the Milwitzky and Cook [1] drop test. The stroke shape and value, corresponding to the drop test, can be seen in Figure 2, and these were obtained using a dynamic system model simulation of the drop test following the analysis in Milwitzky and Cook [1]. The initial delay in the response of the shock absorber, the flat region at time = 0.0 s, is due to the force required to deflect the tire and other internal effects (e.g., the gas pressure force), which leads to the shock absorber behaving as a rigid body during the early stages of impact; however, this is quickly overcome by the impact resulting in the subsequent stroke profile.
The flow in the present simulation is assumed to be incompressible, which is justified due to the low Mach number prevalent throughout the computational domain, and, since the focus of the study is on the free shear layer under a constant inlet velocity, we ignore the influence of the transient compression pressure waves associated with the movement of the piston in a real shock absorber. This assumption helps to increase the efficiency of the simulation by allowing the numerical stability to be determined by the hydrodynamic time-scale rather than the more restrictive acoustic time-scale.  Once the free shear layer becomes fully turbulent, it is expected to follow a nearly linear spreading rate, which is a standard result that has been reported in many studies [12]. There are two main methods of measuring the spreading rate of the mixing layer. The first tracks the width of the shear layer defined by the expansion of the lines through specific velocity ratios across the shear layer as detailed in Equations (1) to (3). The other method of calculating the spreading rate focuses on the growth of the shear layer's momentum thickness as defined in Equations (4) and (5).
The edge of the free shear layer for each profile was found locally based on the vorticity magnitude and streamwise velocity value, and then the high and low velocity values across the profile were defined, and their difference was used for normalisation of the velocity values across that specific profile. The self-similarity distance variable η is defined according to Equation (6) where the capital vertical coordinate Y has its origin at the orifice exit.
It should be noted that, in the present simulation, the streamwise direction is along the vertical y-direction with the corresponding velocity component and turbulent (Reynolds) stress component, v v . For the purposes of analysing the near-field free shear layer development, a two-dimensional (2D) section (x-y plane at z = 0) is extracted from the three-dimensional domain, where the mixing layer is quasi two-dimensional within a distance equivalent to a few orifice diameters downstream of the orifice exit.
In this region, the mixing layer is expected to exhibit standard behaviour that can be used to validate the simulation against similar mixing layer characteristics reported in experimental investigations, which acts as a substitute to the lack of internal shock absorber data available in the literature for validation. Noting that the Cylindrical and Cartesian coordinate systems are equivalent on the extracted x-y plane, the cross-flow Reynolds stress component is u u , and the spanwise stress component is w w relative to the mixing layer development in the quasi two-dimensional region.

Results and Discussion
The results from the three-dimensional and scale-resolved turbulence simulation of the shock absorber's shear layer are generally consistent with standard results for the mixing layer and round jet cases found in the literature. Transition occurs relatively quickly due to the strong primary instability associated with the discontinuity in the velocity profile at the orifice tip, thus, leading to an exponentially growing instability. Within about one orifice diameter downstream of the tip, the shear layer develops secondary and tertiary instabilities, and the larger coherent structures start breaking down to smaller eddy scales as the turbulent cascade progressively forms the inertial sub-range and the dissipation range at the smallest scales where energy is dissipated as viscous heating.
The instantaneous velocity contour in Figure 3a reflects these developments as the initially laminar sides of the shear layer become visibly turbulent with intermittent passage of patches of turbulent flow at the edges. Similarly, the static pressure contours, Figure 3c, show the periodic passage of pressure fluctuations associated with the instabilities developing in the shear layers. The time-averaged velocity contours in Figure 3b show the (expected) more smoothly defined edge of the free shear layers and jet as it spreads travelling downstream. The transition to turbulence is also evident from the turbulence shear stress contours plotted in Figure 3d, where the strong rise in the shear stress downstream of the orifice tip is clearly seen.  Turbulence can also be visualised using the Q-criterion, which is found by plotting the iso-surfaces of the second invariant of the velocity strain rate tensor as shown in Figure 4. The gradual development of increasingly complex flow structures in the shear layer is clearly demonstrated as initially laminar instabilities, such as ring vortices in Figure 4a, develop into more complex structures and breakdown to smaller scales until the full turbulent cascade is established in Figure 4c.
Further insight can be attained by more closely examining the development of the mixing layer downstream of the orifice tip. It is worth noting that the mixing layer in this case is not a standard one, due to the potential constraint imposed by the shock absorber casing, which has a blockage effect, similar to other geometrical details that might be included, such as an orifice support tube or bearing structures. This is in contrast to the unconstrained development region into the far field in most investigations dedicated to the study of mixing layer dynamics, where any physical boundaries are typically more than an order of magnitude away. Nevertheless, comparing the developing axisymmetric mixing layer with standard results could provide useful validation to establish confidence in the simulations. This comparison can also help reach a deeper understanding of the shock absorber flow, by identifying the differences between a standard mixing layer and in a shock absorber.
An initial important observation is the tendency of mixing layers towards an approximately self-similar state moving in the streamwise direction [13]. This state is reached once the free shear layer passes the initial transitional stage near the tip and enters the turbulent regime around one orifice diameter downstream of the exit, where the profiles of turbulent flow variables taken across the shear layer tend to collapse onto a single standard profile once they are normalised appropriately.  The velocity profiles are plotted in Figure 5a, in which the tendency towards a self similar profile is clearly manifested. It is seen in Figure 5a that the first profile plotted at Y/D ori f ice = 0.75 exhibits small deviation from the self similar profile, which is not surprising given that the turbulence flow is still developing at that location, as seen in the visualisation in Figures 3 and 4. The profiles downstream become increasingly more conforming to the self-similar profile.  Other turbulent properties are also expected to exhibit self-similarity. Profiles of the cross flow Reynolds stress component u u are also plotted in Figure 5b, where it is again demonstrated that the normalised profiles tend to collapse towards a single curve when moving downstream.
The stress profiles require a longer streamwise distance to reach approximate selfsimilarity when compared with the velocity profiles. This observation is consistent with the findings of previous mixing layer investigations. However, the less common observation, in the context of a standard mixing layer, is that the Reynolds stress profiles in Figure 5b seem to have a higher value on the negative η side, corresponding to the flow on the high velocity (jet) side of the free shear layer, which will be examined further later in the discussion.
The general pattern of growth is shown in Figure 6a for profiles reaching up to Y/D ori f ice ≈ 9; however, it is important to remember that the jet's potential core only usually extends to roughly five diameters downstream of the exit, and after that, there is strong interaction at the mixing layer's inner edge as the flow tends towards a new self-similar state corresponding to the round jet, which is only achieved much further downstream beyond the boundary of the current simulation domain.
The spreading rate value, SR δ = 0.0878, found here for the mixing layer region in Figure 6b falls within the experimental range of observations, 0.06 ≤ SR δ ≤ 0.11, which was reported for planar mixing layers in Dimotakis [14] and Pope [12]. While the mixing layer here is axisymmetric and constrained by the casing to some degree, these values could be regarded as indicative of the general range expected. The mixing layer's momentum thickness results presented in Figure 7 follow a similar pattern to the δ-based spreading rate. The profiles nearest to the orifice exit plane in Figure 7b are used to calculate the spreading rate, while those in the jet interaction region, more than five orifice diameters downstream, tend to be less reliable as the pattern in Figure 7a shows. The value found here, SR θ = 0.03718, is in good agreement with the spreading rate correlation with velocity ratio across the mixing layer that is available in the literature [15], which, for a single stream mixing layer, has a value of dθ/dx = 0.036. Having verified that the simulated mixing layer exhibits features that are generally consistent with the expected physical behaviour reported in laboratory measurements, it is possible to consider some of the finer details of the specific flow inside this shock absorber case. The contours of the vorticity magnitude and pressure gradient magnitude are shown in Figure 8 for an instantaneous developed flow snapshot.
The pressure gradient magnitude in Figure 8b reaches high intensity values in two regions; the first is in the initial stages of mixing layer development, under the influence of the strong primary instability; the second is on the inner surface of the convergent orifice, where the flow is turned and accelerated as it passes through the orifice leading to a low pressure region near the surface. It is this low pressure region that is of concern, due to its susceptibility to the inception of cavitation bubbles under some loading conditions. These cavitation bubbles usually burst violently when they are convected to higher pressure regions, potentially leading to high-stress transient loading on nearby surfaces. Hence, it can be beneficial to limit the orifice surface area that is subjected to a significantly low pressure value.
A key advantage of the orifice design here is its ability to gradually turn the flow on the upstream boundary of the orifice while avoiding boundary layer separation on the inner surface. The presence of a separation bubble inside the orifice could result in sudden changes in the effective flow area available and the mass flow through it, leading to undesirable variations in the performance under certain loading conditions. This convergent nozzle shape effectively aligns the flow through the orifice in the streamwise direction so that the flow at exit is approximately tangential to the orifice tip. These flow characteristics are confirmed by the instantaneous and mean velocity contours in Figure 3.
The vorticity magnitude contours plot in Figure 8a is consistent with the observations made earlier regarding the flow development in the mixing layer, where, indeed, the high vorticity magnitude is concentrated in the high turbulence intensity zones inside the mixing layer, particularly in the transitional region near the orifice tip. However, it is worth noting that even the potential core of the jet through the orifice has significant vorticity intensity levels and fluctuations within it, indeed there is some notable vorticity fluctuations even upstream of the orifice. Evidently, the orifice design in a shock absorber does not necessarily prioritise settling the flow within the orifice jet's potential core or having a laminar (quiet) flow at the orifice exit. It could, in fact, enhance mixing in the free shear layers and through the central orifice jet to allow more efficient energy dissipation and, hence, improved shock absorber performance in that respect. The flow development upstream of the orifice is the focus of Figure 9, where the contours and streamlines of an instantaneous snapshot of the flow are used to visualise the interactions and flow conditions upstream of the orifice. The spanwise velocity component contours in Figure 9a clearly demonstrate the relatively high intensity of its value, even upstream of the orifice. Given that the inlet velocity is simply set to a constant 1 m/s in the vertical y direction, by symmetry, the spanwise velocity component would be expected to maintain a negligible value until some three-dimensional flow mechanism generates local variations in the flow direction, similar to the ones observed in the axisymmetric jet development region downstream of the orifice, for example.
Streamline plots in Figure 9b,c indicate that the unsteadiness is associated with the flow convergence on approaching the orifice. The streamlines in Figure 9c are coloured by the pressure field, showing the strong pressure gradient across the orifice driving the flow through the orifice to the upper chamber, as expected. It is evident in these streamlines plots how the convergent orifice shape helps gradually turn the streamlines adjacent to the walls from a horizontal direction, the x direction, on the lower wall of the orifice, to the vertical direction, the y direction; however, the flow nearer the central axis of symmetry impinges on the radial flow converging in that region, causing a more violent flow turning in the streamwise direction and leading to the observed flow unsteadiness.
The Reynolds stress profiles, normalised in accordance with self-similarity, are presented in Figure 10, where the relatively high turbulence intensity in the potential core, in comparison with a standard mixing layer case, is observed in both stress components. The disparity between the high velocity and low velocity sides of the mixing layer is greatest in the stress component corresponding to the streamwise direction, v v , while the other components, u u and w w , have a lower level inside the potential core. All the Reynolds stress profiles show a trend of higher maximum amplitude values in the profiles near the orifice exit, a feature associated with the transition process and the strong instabilities present in that region. The profiles settle to a lower more consistent maximum amplitude value for the normalised profiles further downstream in the mixing layer.   The mixing layer development is aided by the convergent nozzle aligning the flow in the streamwise direction at exit. Different orifice design options could potentially reduce the area over which a high-intensity pressure gradient exists and its corresponding low static pressure region, which poses a cavitation risk near the surface. However, that will often involve increasing the streamlines' radius of curvature near the orifice exit and loosing the gradual flow turning/alignment to have the mixing layer flow exit in the streamwise direction, thus, affecting the mixing layer development, particularly in the transitional near field of the orifice, which could, in turn, have a significant impact on the shock absorber performance.

Conclusions
The flow development in a simplified shock absorber setup was investigated using a Large Eddy Simulation. The top boundary was taken to be an outlet, rather than having a closed domain, and a constant inlet velocity, approximately representative of an averaged drop test value, was used. The mixing layer properties, in particular self-similarity and the spreading rate, were found to be consistent with the available literature.
However, significant levels of unsteadiness and turbulence intensity were observed in the orifice jet's potential core and even in the upstream vicinity of the orifice, potentially due to the converging cross flow in that region. This relatively high unsteadiness is not only relevant to the understanding of the flow field around the orifice but can also be significant if a reduced domain LES is conducted with the computational domain starting at the orifice-in which case, representative inlet boundary conditions have to be specified, including a turbulence intensity value.
Further work will be conducted to identify the impact of the orifice shape and loading conditions on the extent and severity of the unsteadiness upstream of the orifice and inside the potential core, in addition to the tendency to generate cavitation regions near the orifice and the overall impact on the shock absorber performance.
The flow in an actual shock absorber is likely to be significantly more complex due to the highly unsteady nature of the loading and the influence of other physical factors, such as heat transfer and multiphase interaction between the hydraulic oil and the gas, which is expected to have a more pronounced impact on the later stages of the stroke, depending on the exact shock absorber design and loading characteristics.
The relevant multiphase phenomena include aeration and foaming that can influence the working fluid properties and, hence, the shock absorber performance as well as cavitation that could impact the service life of internal components. Future work will increase the complexity and scope of the simulations to allow a more accurate prediction of the shock absorber internal dynamics and overall performance.

Data Availability Statement:
The underlying data can be accessed at https://doi.org/10.17862 /cranfield.rd, for the data appearing in the figures. Raw data can be made available on request.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: 2D