Technology of 3D Simulation of High-Speed Damping Processes in the Hydraulic Brake Device

This chapter describes a three-dimensional simulation technology for physical processes in concentric hydraulic brakes with a throttling-groove partly filled hydraulic cylinder. The technology is based on the numerical solution of a system of Navier – Stokes equations. Free surface tracking is provided by the volume of fluid (VOF) method. Recoiling parts are simulated by means of moving transformable grids. Numerical solution of the equations is based on the finite-volume discretization on an unstructured grid. Our technology enables simulations of the whole working cycle of the hydraulic brake. Results of hydraulic brake simulations in the counter-recoil regime are reported. The results of the simulations are compared with experimental data obtained on JSC “ KBP ” test benches. The calculated and the experimental sets of data are compared based on the piston velocity as a function of distance. The performance of the hydraulic brake is studied as a function of the fluid mass and firing elevation of the gun.


Introduction
This chapter describes our three-dimensional simulation technology for physical processes in hydraulic brakes with a throttling-groove partly filled hydraulic cylinder. The technology is based on the numerical solution of a system of Navier-Stokes equations with an additional transport equation to track the working fluid/ free volume interface by the volume of fluid (VOF) method [1,2]. As a solver, we use an iterative algorithm, PISO [1,3], combined with a SLAE solver based on the algebraic multigrid method [4]. To improve the accuracy of solution, we use interface capturing schemes, HRIC [1,4], near the interface between the phases. Moving parts are simulated by means of moving deforming meshes [5], the motion of which is represented in the initial equations with the Lagrange-Eulerian approximation [6]. Our technology enables modeling of the whole working cycle of the hydraulic brake.
To demonstrate the performance of our technology, we are considering simulation of a hydraulic brake with a throttling-groove partly filled hydraulic cylinder. The simulation outputs are compared with experimental data obtained on the test benches at JSC "KBP." The calculated and the experimental sets of data are compared based on the piston velocity as a function of displacement. The performance of the hydraulic brake is studied as a function of the fluid mass and hydraulic brake angle.

Problem definition
The hydraulic brake serves to absorb the energy of moving parts. Its structure consists of a case with a grooved inner surface, a ring piston mounted directly on the outer surface of the barrel and having holes matching the guide rods with stranded counter-recoil springs uniformly distributed around the barrel in the hydraulic brake chamber, and front and rear packings with circular-cross-section rubber O-rings to seal the hydraulic brake chamber filled with hydraulic fluid (Figure 1).
The hydraulic brake starts to react as soon as there is no more free space in the recoil chamber (space between the piston and the rear packing); the fluid is sprinkled through the gaps from the recoil chamber to the counter-recoil chamber (space between the piston and the front packing) at a growing rate and stops the moving parts.
Recoil of the moving parts is countered by the energy stored by the counterrecoil springs while braking. In the counter-recoil phase, the fluid from the counterrecoil chamber is sprinkled through the gaps into the recoil chamber providing the required counter-recoil braking law. Both for recoil and for counter recoil, the areas of fluid flow in the hydraulic brake are generated by the gaps between the piston and the grooved case surface defining the law of motion, and between the piston holes and the guide rods.
Numerical simulation of the hydraulic brake includes the following physical processes: turbulent flow of fluid encountering sudden contractions, expansions, and axisymmetric gaps; multi-phase flow with different equations of state (it is reasonable to use constant-density approximation for fluid and perfect gas law for air); and flow with moving solids.
The scope of this work included the development of a mathematical model and numerical method to simulate these physical processes.

Mathematical model and numerical scheme
Simulation of hydraulic brake operation involves multi-phase flow simulation with two phases. Each phase can have its individual equation of state. Let us assume that the flow is isothermal within one operation cycle and that the phases have the same field of velocity. These assumptions allow us to construct our numerical model using the VOF method [1,7]. Considering the assumptions above, let us describe the flow of matter by a system of equations, comprised of the mass and momentum equations and the volume fraction transfer equation: where t is the time, u i = {u 1 , u 2 , u 3 } = {u, v, w} is the velocity, x i is the spatial vector component, τ ij is the viscous stress tensor, g i is the gravity acceleration vector, ξ is the phase index, α ξ is the volume fraction of the ξ-th phase, and ρ is the resulting density defined as an average over all the phases: where N is the number of phases. The volume fraction transfer equation is written in terms of the transfer of the quantity α ξ : Summing Eq. (3) over all the phases and considering the equality ∑ ξ α ξ ¼ 1, we obtain the following mass equation with respect to the velocity divergence: The momentum equation is written in a half-divergent form. This stabilizes the iterative procedure of numerical solution of this equation [3,4]: An important feature of this problem is that we need to consider the motion of moving parts, which travel to a distance exceeding the half length of the hydraulic brake during computation. To incorporate their motion, we decided to use a topology-preserving mesh deformation method. To prevent unnecessary cell deformations, the mesh nodes located at fixed elements (case and guide rods) follow the moving piston preserving their geometry. The mesh deformation was implemented by the IDW method [5].
To allow for the mesh motion in Eq. (1), we have rewritten the equations of volume fraction transfer and momentum of the phases in a moving frame of coordinate in accordance with the known law [6]: where d * φ dt is the substantial derivative φ with respect to the moving frame of reference and v i is the mesh displacement velocity vector. Using Eq. (3), the volume fraction transfer equation can be written as: Here d * α ξ dt denotes the substantial derivative on the moving mesh. The momentum equation is also defined with respect to the moving frame of reference considering Eq. (3): The mass equation is defined with respect to the velocity in the moving frame of reference: This form of the equations is easy to implement within the framework of finitevolume discretization [8][9][10][11]. The chosen way of taking into account the mesh motion is optimal, because it does not require topology reconstruction, though conservative with respect to the major quantities.
The experimental data on the hydraulic brake operation make it possible to estimate the Reynolds number, which equals Re = 10 5 -10 6 in the cylindrical gap during recoil of the moving parts. The most reasonable way to include turbulent flow components with such a Reynolds number is to use the RANS approach based on the solution of the Reynolds-averaged system Eq. (2) [8,9]. The averaged system of equations is closed using the SST turbulence model, which has proven its efficiency in real-life problem simulations [11][12][13]. In the SST model, the (k-ε) model is formulated in terms of (k-ω) and is focused on resolving small-scale turbulence in the outer flow region. A (k-ω) model intended to describe large-scale turbulence is used in the boundary layer. The combination of these models together is carried out with the help of a function that ensures the proximity of the total model to the (k-ε) model far from the solid walls and to the (k-ω) model in the near-wall flow area.
The resulting system of equations is solved by numerical integration on the finite-volume mesh. The relationship of pressure and velocity while solving Eq. (1) is determined by the SIMPLE algorithm [2]. The SIMPLE algorithm allows to find the pressure field for closure of the continuity equation. Let us write the equation of conservation of momentum at time discretization according to the Euler scheme: where n is the solution from the previous iteration, and j is the solution from the previous time step. To solve this equation, the pressure and velocity are supposed to have the form: Here 0 ≤ α p ≤ 1 is the parameter of relaxation. Substitution of the first expression in Eq. (11) into Eq. (10) yields a preliminary estimate of the velocity value in the next step from the equation: The molecular and turbulent components of the tensor of tangential stress in Eq. (12) are also calculated using u * i . In the second stage, the full speed value at the (n + 1)th iteration is calculated by adjustment using the correction to pressure: The pressure correction is found from Eq. (13) using the continuity condition for u nþ1 i . Taking the derivative of both sides of Eq. (13), we obtain the Poisson equation for pressure: This iterative procedure allows to obtain fields of velocity and pressure satisfying the system of Eq. (1).
Discretization of the equations is performed by a finite-volume technology. Let us show it on the example of the equation of transfer of a scalar quantity φ: The time discretization of Eq. (15) can be carried out using one of the known schemes, for example, the Euler scheme, which is used in the work to solve nonstationary problems using the RANS approach: Here j is the number of the time step.
Let us integrate Eq. (16) over volume and proceed to integration over the area for the convective and diffusion terms: For approximation on a finite-volume grid, the convective term can be written in the form: þ where F k is the volume flow through the face k. The value of φ k on the face k is determined by the applied convective scheme.
The discrete analog of the diffusion term can be written in the following form: þ where n k is the normal to face k.
The value of ∂φ ∂x i k on the face k is found by linear interpolation on the values of the gradient in the cells, which are determined by one of the known methods, for example, the Gauss method: Using this discretization, Eq. (15) is replaced by a system of linear algebraic equations written for each calculation cell: Time discretization of the equations is performed by a three-layer second-order scheme [10]. Discretization of the convective terms in the equation of motion, equation of transfer of turbulent parameters, is performed by the upwind scheme LUD [10], and in the volume fraction transfer equation, by the HRIC scheme [1], which prevents excessive numerical diffusion of the phase interface. The force of gravity was included using a bulk force fitting algorithm [14].
These methods and models have been implemented in the software package LOGOS. LOGOS is a 3D multi-physics code for convective heat and mass transfer, aerodynamic and hydrodynamic simulations on parallel computers [11,[13][14][15][16][17][18]. LOGOS has been successfully verified and demonstrated with sufficiently high efficiency on a number of various hydrodynamic tests, including propagation of gravity waves on a free surface (tsunami) [2,14,18] and industry-specific simulations [11,15]. Speedup of computations on highly parallel computers is provided by an original implementation of the algebraic multigrid method [11,19].

Numerical simulation of hydraulic brake operation
Let us consider the counter-recoil phase in the operation of the hydraulic brake, when the moving parts together with the piston (see Figure 1b) are driven by the springs from their rightmost position to the leftmost one resisting the hydrodynamic force acting on the fluid side.
To solve the problem numerically, an unstructured hexahedral-dominant mesh was constructed using the pre-and post-processing tools of the LOGOS. The mesh is refined near the piston and in the gaps between the piston and the case, where the flow of fluid is particularly strong. Figure 2 shows mesh fragments.
When constructing the mesh, we placed the piston in the middle between its rightmost initial and leftmost end positions of counter recoil. This minimizes mesh deformation during computation. The problem was simulated in the unsteady setup; the time step was 5 Â 10 À4 t 0 , where t 0 is the total time of piston counterrecoil motion. The forces acting on the piston, its velocity and displacement were calculated at the end of each time step after finding the solution of the Navier-Stokes equation.
During computation, the piston together with the moving parts is accelerated by the springs and reaches its maximum velocity during the second half of its travel. The piston then enters the grooved region of the case with a decreasing gap, as a result of which the hydrodynamic resistance grows and slows down the moving parts. The computation ends when the piston is at d m = 0.01 m from its end position, which is attributed to the maximum admissible mesh deformation in the counter-recoil chamber.
When the piston is moving, in addition to the force exerted by the springs and hydrodynamic resistance, the moving parts are exposed to friction arising in the packings, between the guides and the piston, and inside the stranded springs during their compression or tension. The resultant of the forces exerted by the springs and friction can be expanded into two components. The first does not depend on the velocity of the moving parts and can be easily measured experimentally by pulling the moving parts at a slow rate. The second depends on the velocity of the moving parts and is mostly related to friction in the turns of the stranded springs. Empirical estimates show that the second component of the resultant force can be expressed empirically as where u is the velocity, M is the mass of the moving parts, η is an empirical coefficient, which can vary in the range of η = 0-5 s À1 for the springs we use.
Computations were conducted with η = 0, 3, 4, and 5. Figure 3 shows the field of volume fraction of the fluid at different times for η = 0 s À1 .
At t = 0.2 t 0 , where t 0 is the total time of counter recoil, a wave of fluid rises upstream of the piston. Pressure in the counter-recoil chamber grows. Air is actively bled over through the upper part of the gap. A flow of fluid emerges in the bottom part of the gap. Its rate is 10-15 times slower than the rate of the air flow ( Figure 4). As the piston continues to move, the fluid is sprinkled through the gaps at a growing rate. At t = 0.95 t 0 , most of the fluid is already pressed out of the counterrecoil chamber; it forms a nearly homogeneous gas-liquid mixture in the opposite chamber. Pressure in the counter-recoil chamber reaches P = 150 P 0 , where P 0 = 1 atm is the initial pressure.
Fluid velocity in the gap grows as high as V = 50-70 m/s ( Figure 5). After the fluid has flown over through the gap, the bulk gas content in the mixture increases steeply as a result of a pressure drop. This leads to a further growth of mixture velocity to V = 90 m/s ( Figure 5). This effect is less pronounced in the bottom part of the gap because of the lower gas content in the mixture sprinkled through the gap. Figure 6 shows a plot of piston velocity as a function of displacement for different values of η.
On the whole, our results match the experimental data both qualitatively and quantitatively. All the simulated cases consistently predict the state of piston acceleration. In the high-velocity phase of piston motion, the parameter η starts to play a significant role. The best agreement with experiment is provided by η = 3 s À1 . This does not mean that the value above is close to the actual one, but is indicative of the fact that our numerical model, input data, and mesh are capable of describing the process of piston motion with required accuracy, if we use this value of η in our computational model.   The counter-recoil phase occurs at zero angle (the hydraulic brake is placed horizontally). Let us consider the dependence of piston motion on the accuracy of horizontal positioning of the hydraulic brake. For this purpose, let us carry out calculations with angles α = 1°and α = À1°. The calculation results are shown in Figure 7.
The plots demonstrate that the piston velocity grows with increasing angle. This is quite predictable, because the amount of fluid in the counter-recoil chamber decreases as the angle increases. One should also note that the gravity force acting on the moving parts in the direction opposite to that of their motion at α > 0 or in the direction of their motion at α < 0 has no measurable effect on the result.
The angle of the hydraulic brake also controls the static pressure in the counterrecoil chamber (Figure 8). The first pressure maximum is attributed to a decrease in the volume fraction of air in the counter-recoil chamber as the piston moves, and the subsequent drop, to a growth in the gap size between the piston and the grooved case region. As the angle increases, the initial amount of fluid in the counter-recoil chamber decreases, and the maximum shifts toward the range of greater piston   displacements, which allows the piston to speed up to a higher velocity. This increase in the piston velocity leads to a higher pressure maximum in the region of strong slowdown (the second pressure maximum in Figure 6). Thus, an increase in the amount of fluid in the counter-recoil chamber leads to a higher and earlier first pressure maximum and a lower second pressure maximum, which is quite consistent with process physics.
The results of piston motion analysis as a function of initial fluid level in the hydraulic brake (h = h 0 , h = h 0 + 0.01 m, h = h 0 -0.01 m) are shown in Figure 9.
The fluid level has a considerable effect on the piston velocity. A 1-cm increase in the level results in a 10% drop in the average piston velocity, which is explained  by a decrease in the force of hydrodynamic resistance acting on the moving parts. Note that the slopes of the velocity plot are indicative of the same trend: as the fluid level in the counter-recoil chamber rises, the first pressure maximum increases and the second one decreases.
An important characteristic of the counter-recoil piston motion is piston velocity V extr at 159 mm to the counter-recoil end (ΔZ = 153 mm), where the shell extractor is located, the performance of which depends on the velocity of the moving parts. Another important characteristic is piston velocity close to the end position V end , which will govern the force exerted by the piston on the packing in the counterrecoil chamber. The values of these characteristics are presented in Table 1.
The table demonstrates that a change in the fluid level within 1 cm leads to a 15% change in the velocity V extr . A change in the angle to 1 degree results in a 5% change in the velocity. Note that the level and angle variations have a minor effect on the velocity V end , the oscillations of which do not exceed 2-3%.
The results demonstrate that the governing quantity, which controls the performance of the hydraulic brake, is the initial fluid level in the counter-recoil chamber. An increase in the fluid level results in a lower average piston velocity, smaller V extr and V end , shift and rise in the first pressure maximum, and drop in the second pressure maximum. The process pattern, however, does not depend on the reason  of the rise in the fluid level in the counter-recoil chamber: both general increase in the initial amount of fluid and negative angle of the hydraulic brake result in the same trends.

Conclusion
In this chapter, we have presented a three-dimensional numerical modeling technology for physical processes occurring in hydraulic brake devices. The technology is based on the numerical solution of a system of Reynolds-averaged Navier-Stokes equations. To track the free surface, we use the volume of fluid (VOF) method. Moving parts are simulated by means of moving deforming meshes. Numerical solution of the equations is based on the finite-volume discretization, which is used to model the counter-recoil phase in hydraulic brake operation. This technology enables simulations on three-dimensional unstructured meshes. The technology is implemented based on the program package LOGOS, which provides high parallel efficiency of simulations.
The technology has been used to model the counter-recoil phase in hydraulic brake operation. Results of studying the effect of the parameter related to friction, hydraulic brake angle, and initial fluid level are reported. Our numerical experiments have demonstrated that the initial fluid level in the counter-recoil chamber is the governing quantity in hydraulic brake operation. An increase in the fluid level in the counter-recoil chamber as a result of pouring more fluid or placing the hydraulic brake at a negative angle results in the same trend in the change in pressure and its maximums.