’HES'S illllUHIHHIIIHI\HMIWHIHHNI‘\l‘HIWIIHHW 3 1293 017074 This is to certify that the thesis entitled AN INVESTIGATION OF FLUID MOTION INSIDE SHOCK ABSORBERS BY NUMERICAL SIMULATION presented by Ganesha Ekanayake has been accepted towards fulfillment of the requirements for Master of Science degree in Mechanical, Engineering Major professor Date April 24, 1998 0-7639 MS U i: an Affirmative Action/Equal Opportunity Institution . ___——_ __.i _— .. __, . _ _ . _~__——._—_... LIBRARY M'cmgan State Unlverslty PLACE IN RETURN BOX to remove this checkout from your record. To AVOID FINE-3 return on or before date due. DATE DUE DATE DUE DATE DUE 1198 WWW“ AN INVESTIGATION OF FLUID MOTION INSIDE SHOCK ABSORBERS BY NUMERICAL SIMULATION By Ganesha Ekanayake A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of MASTER OF SCIENCE Department of Mechanical Engineering 1998 ABSTRACT AN INVESTIGATION OF FLUID MOTION INSIDE SHOCK ABSORBERS BY NUMERICAL SIMULATION By Ganesha Ekanayake Fluid behavior inside a shock absorber under various conditions is investigated using commercially available CFD code, STAR-CD. Incompressible laminar and turbulent models are used with two forcing functions (constant speed and sinusoidal). Moving boundaries are specified based on the mass conservation. Newtonian and a combination of Newtonian and Non-Newtonian (Power law and modified Bingham plastic) models are considered. The effects of the geometry and the forcing conditions on flow velocity and pressure fields are investigated. The results of the flow field visualization are presented in terms of velocity vector and pressure contour plots. Damping force on the piston is calculated and presented as functions of piston displacement and piston speed. The results agree quantitatively with that of experiments for Newtonian, incompressible, laminar and turbulent flows. ACKNOWLEDGMENT I would like to thank my advisor Prof. Mei Zhuang for guidance, along with my committee, Profs. Craig W. Somerton, and Andre Bénard. Gerald Schmidt and John Rogers at Adapco provided me with the assistance to overcome the numerous problems encountered in STAR-CD. Chassis and Vehicle System Department, Research and Development Center of General Motors Corporation, for their generous financial support, without which this study would not be a success. I would like to make a special thanks to Alex A. Alexandridis, Chandra Namaduri, Bob Foister, Bill Hamilton and Dave Rule of Research and Development Center, General Motors Corporation, for all the assistance and advice, and for giving me the opportunity to work on this study. iii TABLE OF CONTENTS LIST OF TABLES LIST OF FIGURES CHAPTER 1 INTRODUCTION CHAPTER 2 PROBLEM DESCRIPTION CHAPTER 3 NUMERICAL TECHNIQUES AND FORMULATIONS CHAPTER 4 RESULTS AND DISCUSSION CHAPTER 5 SUMMARY AND CONCLUSIONS APPENDD( A USER SUBROUTINE FOR BOUNDARY DEFINITION APPENDIX B USER SUBROUTINE FOR MODIFIED BINGHAM PLASTIC MODEL APPENDD( C CHANGE GRID COMMANDS FOR MESH MOTION APPENDIX D USER SUBROUTINE FOR MESH MOTION BIBLIOGRAPHY iv vi vii 22 38 81 83 84 86 87 91 Table 1 Table 2 Table 3 Table 4 Table 5 Table 6 LIST OF TABLES Material properties and flow conditions used for various cases of simulation Default values for the standard k—e turbulence model Top and Outer boundary speeds for various constant speed piston forcing functions Parameters used in the power law mode] Parameters used in the modified Bingham plastic model Time steps used with various constant speed piston forcing functions 11 12 14 17 19 62 Figure 1 Figure 2 Figure 3a Figure 3b Figure 3c Figure 3d Figure 4 Figure 5 Figure 6 Figure 7 Figure 8 Figure 9 LIST OF FIGURES Sectional view of the controllable shock absorber considered for the study View of the piston and gap and an explanation of geometric parameters Mesh for the Entrance/Exit region of r/g = O Mesh for the Entrance/Exit region of r/g = l Mesh for the Entrance/Exit region of r/g = 5 Mesh for the Entrance/Exit region of r/g = 10 Top view of the shock absorber cylinder with the boundary locations Behavior of stress-strain rate variation for Non-Newtonian fluid models The effective viscosity vs. second invariant of the rate of deformation tensor for the Non-Newtonian models considered in the study Cell with centered node P and neighboring cell with centered Node N used in the finite volume method Node labeling convention for flux discretization for the finite volume method Arrangement of variables and notations for P180 algorithm Figure 10a Velocity field at the entrance to the gap for r/g = O at 900 Figure 10b Velocity field at the entrance to the gap for r/g = 1 at 900 vi 10 I6 19 27 31 34 39 39 Figure 10c Velocity field at the entrance to the gap for r/g = 5 at 900 Figure 10d Velocity field at the entrance to the gap for r/g = 10 at 900 Figure 11a Pressure contours at the entrance to the gap for r/g = 0 at 900 Figure 11b Pressure contours at the entrance to the gap for r/g = 1 at 900 Figure12 Pressure contours in the gap for r/g = 5 at 900 Figure 13a Velocity field at the exit from the gap for r/g = O at 900 Figure 13b Velocity field at the exit from the gap for r/g = 1 at 900 Figure 13c Velocity field at the exit from the gap for r/g = 5 at 900 Figure 13d Velocity field at the exit from the gap for r/g = 10 at 900 Figure 14a Pressure contours at the exit from the gap for r/g = 0 at 900 Figure 14b Pressure contours at the exit from the gap for r/g = 1 at 900 Figure 15a Velocity field at the exit from the gap for r/g = 0 at 1800 Figure 15b Velocity field at the exit from the gap for r/g = 1 at 1800 Figure 15c Velocity field at the exit from the gap for r/g = 5 at 1800 Figure 15d Velocity field at the exit from the gap for r/g = 10 at 1800 Figure 16 Pressure contours at the exit from the gap for r/g = 0 at 1800 Figure 17 Velocity field at the entrance to the gap at 1800 Figure18a Velocity field near the gap at 1780 just before the piston changes its direction of motion Figure18b Velocity field near the gap at 1820 just after the piston changes its direction of motion Figure 19a Pressure contours near the gap at 1780 just before the piston changes its direction of motion Figure 19b Pressure contours near the gap at 1820 just after the piston changes its direction of motion vii 4O 4O 41 41 42 43 43 45 46 47 47 48 48 49 49 50 51 51 52 Figure20a Figure20b Figure20c Figure20d Figure21 Velocity field at the entrance to the gap for r/g = 0 at 2700 Velocity field at the entrance to the gap for r/g = l at 2700 Velocity field at the entrance to the gap for r/g = 5 at 2700 Velocity field at the entrance to the gap for r/g = 10 at 2700 Velocity and pressure fields at the entrance to the gap for r/g = 0 at 270° Figure 22a Velocity field at the exit from the gap for r/g = 0 at 2700 Figure 22b Velocity field at the exit from the gap for r/g = 10 at 270° Figure 23a Velocity field at the exit from the gap for r/g = 0 at 3560 Figure 23b Velocity field at the exit from the gap for r/g = 1 at 3560 Figure 23c Velocity field at the exit from the gap for r/g = 5 at 3560 Figure 23d Velocity field at the exit from the gap for r/g = 10 at 3560 Figure 24 Figure 25 Figure 26 Figure 27 Figure 28 Figure 29 Figure 30 Figure 31 Pressure contours at the exit from the gap for r/g = O at 3560 The damping force vs. time for different constant speed forcing functions with the mesh motion defined using the Change Grid (CG) commands The damping force vs. time for 1.2 m/s constant speed forcing function with a time step of 1.47E-5 s. The number of P180 correctors vs. time for 1.2 m/s constant speed forcing function with the time step of 1.47E-5 s. The damping forces vs. piston displacement for a time step of 2.22E-4s. with the mesh motion defined by the CG commands The damping forces vs. piston displacement for a time step of 4.44E-4s. with the mesh motion defined by the CG commands. The damping forces vs. piston displacement for a time step of 8.87E-43. with the mesh motion defined by the CG commands. The damping forces vs. piston displacement for a time step of viii 53 53 54 54 55 56 56 57 58 58 59 59 61 63 65 66 66 Figure 32 Figure 33 Figure 34 Figure 35 Figure 36 Figure 37 Figure 38 Figure 39 Figure 40 Figure 41 Figure 42 Figure 43 Figure 44 Figure 45 1.13E-3s. with the mesh motion defined by the CG commands. The damping forces vs. piston speed for a time step of 1.77E-4s. with the mesh motion defined by the user-subroutine NEWXYZ. The damping forces vs. piston speed for a time step of 1.77E-53. with the mesh motion defined by the user-subroutine NEWXYZ. Comparison of experimental and numerical results of damping force vs. piston speed for Case 1 of Table 1. Comparison of experimental and numerical results of damping force vs. piston speed for Cases 2 and 3 of Table 1 under a laminar assumption. Comparison of experimental and numerical results of damping force vs. piston speed for Cases 2 and 3 of Table 1 under a turbulent assumption. Damping forces vs. piston speed for Cases 1 and 4 of Table 1 Non-dimensionalized damping force vs. non—dimensionalized piston speed for Cases 1 and 4 of Table 1 Comparison of the damping forces vs. piston speed for different cases of r/g. The damping force vs. piston speed for the power law model of Case 3 of Table 4. The damping force vs. piston speed for the power law model of Case 2 of Table 4. The damping force vs. piston speed for the power law model of Case 1 of Table 4 for a time step of 3.0 E—45. The damping force vs. piston speed for the power law model of Case 1 of Table 4 for a time step of 5.33 E75. The damping force vs. piston speed for the modified Bingham model of Case 2 of Table 5 for a time step of 1.0 E-4s. The damping force vs. piston speed for the modified Bingham model of Case 1 of Table 4 for a time step of 5.33 E75. 67 69 69 7O 71 72 73 73 74 76 76 77 77 79 8O Chapter 1 INTRODUCTION Computational Fluid Dynamic (CFD) is gaining popularity in fluid flow analysis and in this study CFD is used to analyze the flow inside a shock absorber. The study was initiated when the Chassis and Vehicle System department of General Motors Research and Development Center decided to do a feasibility study to evaluate the possibilities of using numerical simulation for the fluid motion inside a shock absorber under different forcing conditions. Fluid motion characteristics inside a shock absorber would ultimately dictate shock absorber behavior. If a design engineer had a design tool to study the fluid motion, then he or she would be able to correlate the characteristics of the shock absorber to the fluid motion inside and hence would be in a better position to modify the shock absorber characteristics. Apart from that, the simulation would allow the optimization of different geometric parameters, and usage of various fluid properties etc., while minimizing the need to obtain experimental data, thus reducing the cost and development time. The study was initiated when the writer was a summer intern at the Chassis and Vehicle System department. It was later extended as a research grant to the CFD group at Michigan State University. The analysis was carried out using the commercially available 1 2 thermo-fluid analysis code STAR-CD, on SUN UltraSparc workstations operating on SunOS 5.5. 1.1 Controllable shock absorbers In an ordinary shock absorber the characteristics are controlled by the valve behavior and by the fluid properties. The shock absorber characteristics for this design are controlled only by changing the fluid properties. When the shock absorber is in operation, the liquid passes through a very narrow annulus (0.735mm) of about 24mm in length and behaves as a Newtonian fluid. The viscosity of the liquid in this annular region can be changed to that of a Non-Newtonian by applying a magnetic field while elsewhere the fluid still behaves as a Newtonian fluid. The characteristics of the.shock absorber are therefore controlled by the magnetic field. The electroviscous effect, which indicates the reversible change in apparent viscosity of a fluid subjected to a magnetic field, was first reported by A.W. Duff [1896]. A patent was obtained in 1947, for the Electro-Rheological (ER) effect, which is considered to be a part of electroviscous effect by W.M.Winslow [1949]. Since that time, a number of attempts have been made to use this effect for clutches, brakes, hydraulic valves, and active damper systems. The main attractions of ER fluid devices are their wide variation of apparent viscosity and their fast response time. Under normal conditions ER fluid behaves as a Newtonian fluid, and an application of electric field increases the resistance to the flow, i.e. increase of viscosity, and the fluid behavior resembles that of Bingham plastic model 3 Stevens et a1. [1984] made a feasibility study of a practical ER damper suitable for vibration control. Morishita, et a1. [1990] have done extensive work on ER fluid controllable shock absorber devices for automobiles. Wylie, et al. [1989] have worked on numerical modeling of a controllable shock absorber. The ultimate objective of this study is to build a CFD tool to model the fluid behavior inside a controllable shock absorber. Though numerous studies have been done on controllable shock absorbers, most of the work concentrates on the controlling aspect. The emphasis of this study, however, is on inflow fluid motion and its effect on damping characteristics. Chapter 2 PROBLEM DESCRIPTION The objective of this study, as mentioned earlier, is to build a CFD model to visualize the fluid motion inside the controllable shock absorber and to predict the damping forces as a function of piston displacement and piston speed under various forcing conditions. A cross section of the shock absorber that was considered for this study is shown in Figure 1. The experimental work that supplemented the numerical results was done at the General Motors Research and Development center. Figure 1 Sectional view of the controllable shock absorber considered for the study 5 This study focuses on how the following factors, entrance and exit geometry, fluid properties, flow physics, Newtonian and Non-Newtonian fluids, affect the damping force behavior. 2.1 Mesh generation The mesh was generated using the Pre-Post processor of STAR-CD, PROSTAR. Since the geometry was simple, a patch was created to represent the cross sectional area of the shock absorber using vertices and splines. It was smoothed to minimize the non- orthogonality of cell. Due to the axis-symmetric nature of the problem, it was decided to use three layers, of one degree each in the azimuth direction. Therefore a three-degree wedge of the shock absorber was modeled. This patch was extruded around the centerline. A higher cell density was used in the gap and the entrance and exit regions. This enabled a higher accuracy and a better resolution in these regions. In order to maintain the continuum, cells of different cell densities had to be coupled using the arbitrary coupling. Because of the arbitrary coupling rather than the integral couplinga fraction of a cell may be coupled to another cell creating partial boundaries, and therefore the partial boundary facility had to be activated. Cells consisted of triangular prisms at the centerline and hexahedrons elsewhere. The mesh motion was defined in the z direction, as indicated in Figure 1. To study the effects of the entrance and exit geometry, various radii of curvature at the entrance and exit were considered. Magnified view of the piston and the gap is shown in Figure 2. The radius r was non—dimensionalized with the gap thickness g. A gap thickness of 0.735 mm was used in all the cases. A ratio of r/g = 6.8 was used as the 6 standard case. Four other cases of different r/g ratios for a given geometry were 0,1,5 and 10. The gap length 1 was kept constant at 23.82 mm. The same r/g ratio was used for both the entrance and the exit with r/g = 0 being the square edge entrance or exit. 1/ inner outer A‘ij/ Figure 2 View of the piston and gap and an explanation of geometric parameters The mesh at the entrance/exit for r/g = 0, 1, 5, and 10 are given in Figures 3a-d respectively. Figure 3a Mesh for the Entrance/Exit region of r/g = 0 Figure 3b Mesh for the Entrance/Exit region of r/g = 1 Figure 3c Mesh for the Entrance/Exit region of r/g = 5 Figure 3d Mesh for the Entrance/Exit region of r/g = 10 9 To study the effects of other parameters on the damping characteristics the standard geometry of r/g = 6.8 was used. In the cases of Non-Newtonian calculations, the mesh at the entrance and exit was refined in order to specify a gradual transition from Newtonian fluid to Non-Newtonian fluid behavior. 2.2 Boundary conditions Symmetry plane boundary conditions were imposed on the two surfaces in the azimuth direction. The normal component of the velocity and the normal gradient of all other variables were set to zero. The boundary locations are shown in Figure 4. The boundary conditions, for the outer wall was assigned through the user subroutine. It was defined as a no-slip, smooth, moving wall. The components of the velocity in r and 0 directions were set to zero. i.e. u = v = 0. The component in the z direction, depended on the forcing function. For the cases of constant speed forcing, it was simply assigned the same magnitude as the piston forcing speed. For the sinusoidal forcing function, the time derivative of the motion function was used and these are discussed in Section 2.8. For all the other walls, no-slip, smooth and stationary boundary conditions were imposed. Additionally, it was assumed all these walls behave as adiabatic. For the turbulent calculations, wall functions, which would be discussed later, were used. 2.3 Initial conditions All the velocity components were set to 1.0E—5 m/s for the initial conditions. This was necessary in order to ease the start up process and to eliminate possible discontinuities due to zero velocity. Experimental data showed that the Nitrogen chamber 10 was pressurized to a pressure of 2.241 MPa, when the piston was furthest from the chamber. It was assumed, the same pressure would be felt everywhere in the liquid, at stationary conditions, and hence the initial pressure was set 2.241 MPa. l outer wall plane of entry Figure 4 Top view of the shock absorber cylinder with the boundary locations 2.4 Material properties For the study it was assumed that the temperature variations were negligible and had no effect on the flow field. Hence, temperature solver was turned off. The density was considered as a constant and so was the laminar viscosity for Newtonian fluids. Various viscosity and density combinations with appropriate flow conditions were used in the simulation and are tabulated in Table l. Laminar and turbulent flow conditions were distinguished using the Reynolds number based on the gap thickness. The critical Reynolds number is 2100, if the flow is assumed to be in an annular region [Bird, 1960]. The Reynolds number was 2600 for Case 2 and 2400 for Case 3. The Case 4 listed in Table l is hypothetical. The density and the viscosity values are twice those of Case 1, but the Reynolds number remains the same. The purpose was to perform a dimensionless analysis for the damping characteristics. Table 1 Material properties and flow conditions used for various cases of simulations Case Density (kg/m3) Viscosity (Pas) Flow Condition 1 818 24.0 E -3 Laminar 2 2276 31.97 E -3 Transition 3 2964 45.9 E -3 Transition 4 1636 48.0 E -3 Laminar 2.5 Turbulent flow For the Turbulent flow calculations with Newtonian fluids, standard k-e model was used with a characteristic length of 0.735 mm which is the gap width g. All the other parameters were assigned the default values, shown in Table 2. 12 Table 2 Default values for the standard k-e turbulence model Cu 0.09 o, 1.0 o, 1.22 c,l 1.44 C,2 1.92 Cr, 1.44 C,4 -0.33 k 0.42 E 9.0 2.6 Mesh motion In the experimental setup, the outer casing of the shock absorber was stationary while the piston was forced through the fluid in a specified forcing function. In order to simulate the forcing functions, the mesh motion capabilities of STAR-CD were utilized. Computationally, it is easier to move the vertices where the cell density is less, rather than where the cell density is high. Therefore, instead of moving the piston, the outer casing was moved. The Change Grid (CG) commands in PROSTAR and the user subroutine to specify the vertex positions with time NEWXYZ, were used to specify the mesh motion. It was advised to maintain cell size constant by changing the number of cells in the mesh for moving mesh simulations. In this case, keeping a cell size constant was not so critical; l3 hence the number of cells in the mesh was kept constant, while the size of the cells near the top and bottom surfaces changed to accommodate the mesh motion. This was possible since the minimum length scale and maximum velocity for the critical Courant number C, which is defined as, _ MAI _ 1 C (2.1) where M local velocity 1 mesh dimensions were obtained from the gap where the cell density is relatively high. 2.7 Incompressible flow In the shock absorber, as the piston moves, the liquid volume changed, to accommodate the connecting shaft volume. The change in the liquid volume was compensated by changing the volume in a pressurized Nitrogen chamber located at the t0p portion of the shock absorber, see Figure 1. The change in the volume was accomplished by having a movable wall, separating the liquid and Nitrogen. The motion of the movable wall was defined based on the conservation of mass, as a ratio of the displacement of the top movable wall to the outer moving wall, D, was given by, D = [Elw— ] = 0.076 (2.2) cylinder / where dshafl is the diameter of the shaft and dcynnder is the diameter of the cylinder. 2.8 Forcing functions Several forcing functions were considered for this study. First a constant speed forcing function, with discrete values of piston speeds 2, 1.6, 1.2, 0.8 and 0.4 m/s was considered. Assuming that the flow is incompressible and the mass is conserved, the speed of the top wall was calculated using an equation similar to Eq. 2.2. The boundary region assigned for the outer wall was defined as a moving wall with the same speed as the forcing function. The speeds of the top and the outer walls for a given piston forcing speed are tabulated in Table 3. Table 3 Top and Outer boundary speeds for various constant speed piston forcing functions Forcing function (m/s) Top wall (m/s) Outer wall (m/s) 2 1.848 2 1.6 1.4784 1.6 1.2 1.1088 1.2 0.8 0.7392 0.8 0.4 0.3696 0.4 Next a sinusoidal forcing function with an amplitude of 101.6 mm and a frequency of 6.26 Hz was considered. This function was used for most of the cases considered. Using the conservation of mass, the motion profiles for the top and bottom boundaries were given by 15 am, = (273.64 + 46.939 Sin 8):) mm (2.3) (180,“,M = (50.8 + 50.8 Sin a») mm (2.4) respectively, where 63 = 39.33274002 rad/s. The outer wall speed was described by VOUT = 1.9981 Cosw t m/s (2.5) 2.9 Non-Newtonian model As mentioned earlier, the current study includes both Newtonian and Non- Newtonian fluids. The Non-Newtonian fluid behavior was restricted only to the region of the straight portion of the gap 1, which has a length of 23.82 mm, while the Newtonian model was used elsewhere. The stress-strain rate variations for the Non-Newtonian fluids are shown in Figure 5. The power law and modified Bingham plastic models were used to specify the Non-Newtonian behavior in the study. The power law was used to evaluate STAR-CD capabilities while the modified Bingham plastic model was of primary interest. l6 STRESS +NEWTONIAN ‘ -I— POWER LAW ; +POWER LAW . +B|NGHAM STRAIN RATE pseudoplastic L dilatant Figure 5 Behavior of stress-strain rate variation for Non—Newtonian fluid models 17 2.9.1 The power law model This two-parameter model is also known as Ostwald—de Waele model. The viscosity is calculated by n—l y = mll2 (2-6) where p. effective viscosity 1; second invarient of rate of deformation tensor m, n constants Pseudoplastic and dilatant fluids correspond to n < 1 and n > 1 respectively and the Newtonian fluids corresponds to n = 1. The constants m and n used in the study with the power law model are given in Table 4. The comparisons of effective viscosity as a function of the second invariant of the rate of deformation tensor with those of the modified Bingham plastic model for the real fluid are shown in Figure 6. Table 4 Parameters used in the power law model case m n 1 1350 0.999 2 87 0.85 3 24E-3 0.5 2.9.2 The Bingham plastic model The Bingham model is another two parameter model as shown in Figure 5, and if a substance follows this model it is called a Bingham plastic; it remains rigid when shear stress is smaller than the yield stress but flows like a Newtonian fluid when the shear stress is exceeded. In terms of the second invariant shear stress can be expressed as )> 2'02 (2.7) "1 i=[uO+————TO :lg f0r(_g_: :77: where g stress tensor 3 rate of deformation tensor (Q = O for (39> 2'02) I 2 second invarient of g To yield stress ,uo viscosity 2.9.3 The Modified Bingham plastic model When 12 approaches zero, computationally the above model creates a problem because of the discontinuity. By combining the Bingham fluid model with shear thickening power law the problem can be avoided [Ahamed, 1994]. ,U + foil —exp(— mIJI—Zi) O W II"! 2 (2.8) 19 where [10,10 and m are constants. This is a continuous model. The parameters used with the modified Bingham plastic model in this study are tabled in Table 5. Table 5 Parameters used in the modified Bingham plastic model Case To / Pa 110 / (Pas) m/(s) 1 34.5 E 3 31.97 E —3 4.0 E -2 2 32.0 E-3 31.97 E —3 7.0 E -5 — 1600.00 M 140000 1200.00 A I , . ,, - f / t-I—Mod. Bingham1 ‘ / 100000 *—O—Mod. Bingham 2 800.00 ‘ + Power 13w 1 / + Power law 2 j X 60000 +P_owe_rlaw3 / —~— — / 400.00 ’ 200.00 W» 0.00 1.00E+07 1.00E+02 1.00E-03 1.00E-08 Figure 6 The effective viscosity vs. second invariant of the rate of deformation tensor for the Non-Newtonian models considered in the study 20 The transition from a Newtonian to a Non-Newtonian fluid was defined by increasing the yield rate from zero to a given value in the transitional region. The variations of effective viscosity with the second invariant of the rate of deformation tensor for different models with various constants are shown in Figure 6. It indicates that the two models have different characteristics and that it would be difficult, but not impossible to model the modified Bingham plastic behavior with the power law and vice versa. 2.10 Solution control As the problem was inherently transient the PISO algorithm was used for the calculations. The load step was defined with a constant time step value. As discussed in detail later, it was found out that the time step size has a very important role in the solution process. Upwind Differencing (UD) was the primary choice for the discretization. SFCD, a second order scheme and QUICK, a third order scheme were also used only for verifying but there were no significant improvements in the results. The under relaxation factor for pressure was in the range of 1.0 to 0.8. The maximum number of PISO correctors used was 40. Increasing this any further did not improve the results. Number of inner sweeps for velocities and pressure, were given values of 100 and 1000 with solver tolerances at 0.01 and 0.001 respectively. 2.11 Post Processing The post processing was done using the PROSTAR capabilities. Velocity vectors and magnitudes and pressure contours were obtained at discrete time intervals. With the animation capabilities of STAR-CD it was possible to animate the fluid motion during 21 the forcing processes. The main interest was to obtain the applied force on the piston during the forcing. A subset of wall cells was defined to cover the surface of the piston and the rod. The total force, consisting of pressure and shear forces for each cell in this set, was calculated. The z direction component of the force for all the cells was summed up to give the total force on the piston. Together with the piston motion profile, the force vs. the piston displacement and the force vs. the piston speed are constructed and thereafter compared with those of the experimental results. Chapter 3 NUMERICAL TECHNIQUES AND FORMULATIONS 3.1.1 Basic conservation equations The mass and momentum conservation equations solved by STAR-CD for general incompressible fluid flow and a moving coordinate frame (the Navier Stokes equations) are, in Cartesian tensor notation: l a 0 ~ _ TEE(\/—g-p)+—a—;(puj )— Sm (3-1) fi%(\/gpui)+£7(pfijui —T,j)=——aa?pi+si (32) where t time x,- Cartesian spatial coordinates (i = 1,2,3) u,- absolute fluid velocity component in 1' direction uj ui—uq, relative velocity between fluid and local (moving) coordinate frame that moves with velocity up,- p piezometric pressure 22 23 p density 1;,- stress tensor components sm mass source s, momentum source component Vg determinant of the metric tensor The specialization of the above equations to a particular class of flow involves: 0 Application of ensemble or time averaging if the flow is turbulent. 0 Specification of a constitutive relationship, connecting the components of stress tensor Tij to the velocity gradient. 0 Specification of the source 5, which represents the sum of the body and other external forces, if present. 3.1.2 Constitutive relations In the case of laminar flow, for Newtonian and Non-Newtonian fluids, following constitutive relation exists: 2 Bu Tl} =2/JSU —§' 5;];5” (3.3) k where 11 viscosity 8g Kronecker delta sij rate of strain tensor 24 The rate of strain tensor, is given by: l Bu, Bu- ..2— —— ’ 3.4 Su2axj+axi) () Non-Newtonian fluids fall into shear thinning or shear thickening class, in which the fluid viscosity is dependent on the local value of the second invariant of the rate of strain tensor [2, which is defined as 1 12 E 5(51'131)’ — Sijsij) (35) There are several different viscosity models for Non-Newtonian fluids, such as power law, Bingham model, as discussed earlier. 3.1.3 Newtonian turbulent For turbulent flows 14,, p and other dependent variables, including Tjj, assume their ensemble averaged values (equivalent to time average for steady state situations) for Newtonian fluids 2 0L4 7—,? 2"]:211 sq. —§ fidU—puiuj (3.6) k where u’ is the fluctuations about the ensemble average velocity. The overbar denotes the ensemble averaging process. The rightmost term in the above equation represents the additional Reynolds stresses due to turbulent motion. These are linked to the mean velocity field via turbulence models. 25 3.2.1 Discretisation practice The differential equations governing the conservation of mass, momentum, within the fluid, are discretised by the finite volume (FV) method. Thus, they are first integrated over the individual computational cells and over a finite time increment in the case of a transient problem and then approximated in terms of the cell centered nodal values of the dependent variables. This method has the merit of ensuring that the discretised forms preserve the conservation properties of the parent differential equations. For the purposes of the FV discretisation, it is convenient to work with the following general coordinate free form of the conservation equations: —\/l=%(\/_g—p¢)+ div(pii,¢—Fograd¢)= 5‘, (3.7) 8 where ii, 2 ii - iic realative velocity between fluid and local coordinate ¢ any dependent variable I" diffusion coefficient 3 source coefficient. An exact form of the above equation, valid for arbitrary time varying volume V bounded by a moving closed surface S can be written as follows %Ip¢dV+I(pfi,¢—F¢grad¢).d3 =Is¢dV (3-8) V S V 26 where 5 surface vector ii, relative velocity between fluid and surface. If V and S are, respectively, taken to be the volume VP and the discrete faces 5] (j = 1,Nf) of a computational cell (see Figure 7), the above equation becomes %Jp¢dV +2 [(pE,¢—F0grad¢)dS= [sodV (3.9) VP 1' S, V l i l Ti 72 T3 From here onwards, approximations are introduced. Thus the first term, T1, of Eq.3.9 is discretized as T 2 (prov); —(p¢v ;: (3.10) ' 6t where the superscripts 0 and n refer to ‘old’ and ‘new’ time levels, respectively, separated by a time interval St. The second term of the Eq. 3.9, T2, is split into the separate contributions C} and Dj due to convection and diffusion, respectively, and each is expressed in terms of average values over cell faces, denoted by 0} thus: T2 z 209145)]. —2(F¢grad¢.3)j E 2C1_sz (3.11) ,- ,- j ,- 27 o N‘ \‘w": .P faccj Figure 7 Cell with centered node P and neighboring cell with centered Node N used in the finite volume method The diffusion terms D,- are approximated by face centered expressions of the form 07 '7 Bait-Wat -¢,, Fifi-WM] (3.12) where the first term in the brackets represents the gradient between the cell-centered node P and its neighboring node N (see Figure 7), and the second (summation) term is over all vertex pairs on face j, the vertex (i) values being interpolated from surrounding nodal values. Termfj are geometrical factors and F9]- is the (interpolated) face diffusivity. The approximation of the convective flux terms is discussed later under spatial flux discretisation. The third term of Eq. 3.9, T3 may in general contain components representing sources or sinks of the transported property, as well as additional flux terms; the exact form depends on the circumstances and the dependent variable. Fluxes and 28 other gradient containing terms are approximated in a similar fashion to the Cj and D, while non-gradient quantities are evaluated using the cell-centered nodal quantities. For convenience, the results of this process is written in the general quasi-linear form T3 = S. - 32¢ (3.13) P 3.2.2 Temporal discretization As already noted the FV equation applies over an arbitrary time increment spanning the “old” and “new” time levels. The fluxes and sources prevailing over this interval are taken to have their new time level values; i.e. a fully implicit formulation is used. This avoids a stability related time step restriction inherent in the otherwise attractive explicit approach, which becomes onerous in regions of small mesh spacing and high velocity or diffusion rate. In principle, the fully implicit formulation allows any magnitude of time step to be used, but in practice other consideration imposes limits, i.e. for transient problems, must be small enough to limit the temporal approximation errors to acceptable levels. 3.2.3 Spatial flux discretization The manner in which the convective and diffusive fluxes are expressed in terms of nodal values is one of the key factors determining accuracy and stability, for both steady state and transient calculations. At high Reynolds numbers, the choice of convective flux approximation is particularly important. There are three main classes of convective flux approximation in widespread use, namely: 29 0 “Low order” schemes, which characteristically generate discretised equation forms that are easy to solve, produce solutions, which obey the expected physical bounds, but sometimes give rise to smearing of gradients. i.e. numerical diffusion. This is a form of truncation error that diminishes, as the grid is refined, but at an increased cost of calculation. Upwind Differencing (UD) is a first order scheme. 0 “Higher order” schemes, which better preserve steep gradients, but may result in equations that, are more difficult to solve and / or have solutions exhibiting non-physical spatial oscillations, i.e. numerical dispersion. This too can be diminished by grid refinement. Linear Upwind Differencing (LUD) and Central Differencing (CD) are two second order schemes. 0 “Filtered” or “Blended” schemes, which combine a higher order formulation with a strategy for suppressing the spatial oscillations, usually by detecting their onset and then locally modifying the discretisation in a suitable way. The local modification is almost invariably a form of reversion to a lower order scheme, either completely or partially. Because these schemes must base their filtering/blending practice on the evolving solution, additional non-linearity and coupling are introduced which can adversely effect the performance in a calculation. Self-Filtered Central Differencing (SFCD), Gamma, Blended Differencing (BD) and Quadratic Upstream Interpolation of Convective Kinematics (QUICK) are some of these schemes. The Cj factors defined in Eq. 3.1 l are written as 30 C]. E [71¢]. (3.14) F]. a (pas). (3.15) where F,- is the mass flux through the face j and (1)], the average value at the face, is effectively interpolated from selected nodal values in accordance with the scheme used. The face values of the auxiliary properties such as p are T also obtained by interpolation, usually linear. 3.2.3.1 Upwind Differencing (UD) The main choice for differencing schemes for this study was the “Upwind Differencing (UD)”. This lower order scheme selects the nearest upwind neighbor value for (see Figure 8): (3.16) This form of interpolation preserves the correct physical bounds under all conditions, but can lead to numerical diffusion. 31 O O P 9ch - Figure 8 Node labeling convention for flux discretization for the finite volume method 3.2.4 Final FV equations The final form of the discrete finite volume equation is obtained by substituting the various approximated terms back into Eq. 3.9 and then invoking the discretised continuity equation, which can be written as 006—0012,]. 2 The result, in its most compact form, is Air); = ZAmtz); + s1+ 8P0}? where Am effects of convection and/or diffusion 0 B, (£2). 51 AP 22A", +s2 +B,, (3.17) (3.18) 32 There exists an equation such as Eq. 3.18 for every computational cell (suitably modified, where necessary, to incorporate boundary conditions). There are as many such equation sets as dependent variables, when the continuity set is taken into account. As described later, the solution strategy in STAR-CD involves iterative solution of these sets. It is pertinent in this regard to note that the choice of convective differencing scheme can have a strong bearing on the reliability and speed of iterative method. 3.3.1 Solution algorithm As mentioned earlier implicit methods are employed to solve the algebraic finite volume equations resulting from the discretising. Implicit method involves the solution of simultaneous algebraic equation sets, which sometime held to be a disadvantage. However very efficient, almost invariably, iterative methods have been developed for solving these equations with their characteristically sparse matrices. STAR-CD currently incorporates three different implicit algorithms, SIMPLE, SIMPISO, which are basically for steady state calculations, in an iterative mode. PISO is applicable to both transient and steady state calculations and particularly suitable for transient. It was the only solution algorithm used in this study. 0 Algorithms employ a form of predictor-corrector strategy, enabled by the use of operator splitting, to temporarily decouple the flow equations from each other so that they can be solved sequentially. 33 0 Continuity is enforced with the aid of an equation set for pressure, derived by combining the FV momentum and mass conservation equations. 0 The solution sequence involves a predictor stage, which produce a provisional velocity field derived from the momentum equations and a provisional pressure distribution. The provisional fields are then refined in the corrector stages by demanding simultaneous satisfaction, to some approximation, of both momentum and continuity balances. 0 Within the sequence, operator-split equation sets involved at any stage, only one of the field variables, i.e. the vector set of unknowns is split into a sequence of scalar sets. 0 The algebraic equation sets are solved by iterative means, though this is not an essential practice. 3.3.2 The FV momentum equation The FV momentum equations are extracted from the general transport Eq. 3.18 in the form of 4.4:. = H(u:i.. )+ BEqup + s. + 0.027.. — pi-) (3.19) where ”(aim )E 2 Amuim The last term in Eq. 3.19 is the FV approximation to the pressure gradient, DP being a geometric co-efficient (see Figure 9). —l l l [— Figure 9 Arrangement of variables and notations for PISO algorithm The FV continuity equation is written here as B; —B;; +2(pm;'sj)=0 (3.20) I where u j. velocities normal to the cell faces S j face areas In the collocated variable arrangement, the face velocity 14,- need to be expressed in terms of the nodal velocities (in order to calculate the mass fluxes) and neighboring pressures (to allow formulation of a suitable pressure equation). Assembling a cell face momentum equation of the following form does this. 35 XpuyzHufm +8214, +Si+DP(p;, —p;,) (3.21) where the overbars denote a form of averaging on the nodal momentum co-efficient appearing under them. This equation, when substituted into the continuity relation, Eq. 3.20, yields a pressure equation of the form App; =2Amp; +sl (3.22) where the “source” term is a function of the nodal velocity fields and other quantities. Thus, an equation set is produced from which the pressure may be calculated. 3.3.3 Solution Sequence Starting from initial values (1)" of the variable fields, PISO advances through a time increment 61 in the following sequence of steps. Predictor stage. Eq. 3.19 are assembled and solved in the following operator-split form for the provisional nodal velocity field 14,-”) A it“) =H(u,{',3,)+ Bgugf, +5, +1) ,(p<°3— pg?!) (3.23) p iP where p (0) correspond to the pressure field at the start of the time step. Following the solution of this equation, obtained by an iterative method, the provisional face m are calculated via Eq. 3.21, with the u,” and p" replaced by ut (I) and PO, velocities uj respectively. First corrector stage. The operative nodal momentum equation is now taken to be A um: H(u i‘))+ 8314”,, +5l +D P(p(')— pm) (3.24) p LP rm 36 and the face momentum equations are likewise approximated by replacing u,” and p" in Eq. 3.20 by u,- (2) and pm, respectively. Correspondingly, the pressure equation in its approximated form reads App; = Zamp; + s, (3.25) where s. is now a function of the (known) nodal velocities u; (I) and uio. Accordingly, these equations can be solved (here again by iteration) for the p (1) field, following which “'12) and 14,- (2) can be obtained from equation and its face counterparts. The resulting solution is an approximation to that of the original Eqs. 3.19 and 3.20. 0 Additional corrector stages. Further correctors are performed in the same manner as the first one, using the generalized equations 4.4.3:" = H(u§.:.?)+ Bzuzip + s. + D. (p532 — p133) (3.26) 71.16;”: EA," pig“ +5, (3.27) where q =1,2,.. is the corrector level. Note that the coefficients Ap are held constant. The solutions from the successive stages represent increasingly better approximation (q+l) (q+l) to the solution to the original equations; i.e. u,- and p tend towards u," and p" with increasing q. After completion of the required number of correctors, judged according to practices described below, the solution produce is taken as the starting field for the next time step and sequence is repeated from stage 1. If the calculation of scalar fields such as turbulence parameters and temperature is required, it is performed in further steps executed after the final flow corrector. 37 3.3.4 Completion tests For the unsteady calculation with PISO, the monitoring information is the global rate of change, C2, defined as C: = 2 (3;); where the k now denotes the number of time steps performed and the summation is over -|B;:¢;:l) (3.28) all cells. In the case of the mass conservation equation, the quantity summed is actually 3",» - B0,». The magnitudes of the C O quantities evidently provide a global measure of the rate of change of mass, momentum, energy etc. within the calculation domain and can therefore be used for variety of purposes. Chapter 4 RESULTS AND DISCUSSION 4.1 Flow Visualization In order to study the effects of the geometry of the shock absorber on the internal flow field, various geometries were considered for a given sinusoidal forcing function. Velocities and pressure fields were calculated for a transient flow simulation. The simulations were done for 3600 of piston motion. The piston position is defined as 00 when it is closest to the Nitrogen chamber. The results of the velocity and pressure fields are presented for different geometries with values of r/g = 0, 1, 5, and 10 at various piston positions, 900, 180°, 270°, and 3560. A detailed study of the flow field was carried out when the piston changes its direction (1800). The velocity vectors and magnitudes at the entrance to the gap, with r/g = 0,1,5 an at I: . - SCC. ares own 1n lures a-. 1meste , IOI'ICEIOVC d10 90°( 40152 ) h 'F'g 10dT' pAf hb cases is 8.87E-4 seconds. It was obvious that for the case of r/g = 0 a non-smooth passage for the fluid exists due to the sharp edges. It could also be observed that the boundary layer on the inner side was thicker than that of the outer side for all these cases (see Figure 2 for the definition of inner and outer sides). As the ratio r/g increases, flow through the entrance becomes smoother. 38 Figure 10b Velocity field at the entrance to the gap for r/ g = 1 at 90° Figure 10c Velocity field at the entrance to the gap for r/g = 5 at 900 Figure 10d Velocity field at the entrance to the gap for r/g = 10 at 90° 41 Figure 11a Pressure contours at the entrance to the gap for r/g = 0 at 90° Figure 11b pressure contours at the entrance to the gap for r/g = 1 at 90° 42 The absolute pressure contours for r/g = 0 and 1, at time t = 4.0 E-2 are given by Figureslla and b. Due to the effect of the sharp edges, a pressure wake can be seen at the entrance for the case with r/g = 0. The pressure contours in the gap for the case of r/g = 5 are given in Figures 12. The pressure drop decreased because of the smooth entrance. The pressure wake has disappeared too. Figure12 Pressure contours in the gap for r/ g = 5 at 900 The exit region, velocity vectors and magnitudes, at 90°, for the cases with r/g = 0,1,5, and 10 are given in the Figuresl3a-d. It can be seen in Figure 13a that the jet is pushed toward the outer wall by a strong swirling motion. Two vortices formed at the exit, one rotates in the counterclockwise direction and the other rotates in the clockwise direction. As the ratio of r/ g increases, the jet is straighten up, and for the case of r/g = 10 43 Figure 13b Velocity field at the exit from the gap for r/g = 1 at 900 Figure 13c Velocity field at the exit from the gap for r/g = 5 at 90° Figure 13d Velocity field at the exit from the gap for r/ g = 10 at 90° 45 the jet is almost straight. It can be assumed that for this ratio the jet is stronger than that of other cases. It can also, be observed that as the jet gets straighten up the fluid is getting entrained into the jet. The maximum value for the velocity magnitude was around 50 m/s for all these cases. These two vortices mentioned for the case of r/g = 0 could also be identified in the pressure contours (see Figure 14a). The two low—pressure regions next to the exit, one on the inner side and the other on the outer side can be seen in Figure 14a. The pressure contours for the case of r/g = 1, are plotted in Figure 14b. It also shows the two vortices at the exit, but with a lesser intensity. Figure 14a Pressure contours at the exit from the gap for r/g = 0 at 90° Figure 14b Pressure contours at the exit from the gap for r/g = 1 at 90° In the exit region, the velocity vectors and magnitudes, at the time t = 8.00E-2 seconds (180°) are given in Figures 15a-d for the cases of r/g = 0,1,5, and 10. At this moment the direction of piston motion is about to change and therefore the magnitude of the velocity has a smaller value. Still the large clockwise vortex motion can be observed at the center for all these cases. The smaller counter clockwise vortex next to the exit at the inner side for the case of r/g = 0, can be seen in Figure 15a. As the ratio r/g increases this small vortex moves into the exit region and for the case of r/g = 10 this counter clockwise vortex becomes significant. 47 Figure 15b Velocity field at the exit from the gap for r/g = 1 at 180° Figure 15d Velocity field at the exit from the gap for r/g = 10 at 1800 49 Figure 16 Pressure contours at the exit from the gap for r/g = 0 at 1800 Figure 17 Velocity field at the entrance to the gap at 180° 50 The pressure contours at the exit, for the case of r/g = 0 at 180° are given by Figure 16. A high-pressure region is formed next to the piston and this can be attributed to possible flow stagnation. In the entrance region, at the same piston position (180°), flow is almost stationary, see Figure 17. The velocity vectors and contours just before and after the direction change of the piston motion (at 178° and 182°) are shown in Figures 183 and b while the pressure contours are shown in Figures 19a and b. These figures show a detaibd flow development when a change of direction of motion of the piston takes place. Figurel8a Velocity field near the gap just before the direction change (178°) 51 Figure18b Velocity field near the gap just after the direction change (182°) Figure 19a Pressure contours near the gap just before the direction change (178°) 52 Figure 19b Pressure contours near the gap just after the direction change (182°) At the time t = 0.12 seconds (which corresponds to the piston position 270°) the velocity vectors and magnitudes are plotted in Figures 20a-d for the cases of r/g = 0,1,5 and 10, while the velocity vectors on pressure contours for the case of r/g = 0 is shown in Figure 21. As mentioned earlier, because of the existence of sharp comers for the case of r/g = 0, a vena contractor has formed (see Figure 20a). =,rirtt1 A .)}{arllllfivilof_i rpprljlrfirlvoxr‘.’ 1:11'rllrrlifcl' Figure20b Velocity field at the entrance to the gap for rl g = 1 at 270° Figure20c Velocity field at the entrance to the gap for rlg = 5 at 270° Figure20d Velocity field at the entrance to the gap for r/g = 10 at 2700 55 Figure21 Velocity and pressure fields at the entrance to the gap for rlg = 0 at 270° The velocity vectors and magnitudes in the exit region for the cases of r/g = 0 and 10, are given in Figures 22a and b respectively. The jet is deflected toward the outer side and the deflection is greater for the case of rlg = 0. This may be due to the fact that as the ratio decreases, the jet becomes weaker as mentioned earlier. The entrainment of the fluid into the jet in the inner side and a swirl development on the outer side can also been observed for the case of r/g = 10. Figure 22b Velocity field at the exit from the gap for r/g = 10 at 2700 57 In the exit region, the velocity field at the time t = 0.16 seconds (356°), are plotted in Figures 23a-d. The counter clockwise vortex motion is dominant at this time. The maximum velocity magnitude has reduced to about 5 m/s. The two vortices in the exit region are developing as the ratio r/g increases. They can be seen clearly in Figure 23d. The pressure contours for the case of r/g = 0 are shown in Figure 24. The low-pressure region at the center is due to the dominant vortex motion. Figure 23a Velocity field at the exit from the gap for r/ g = 0 at 3560 58 Figure 23b Velocity field at the exit from the gap for r/g = 1 at 356° Figure 23c Velocity field at the exit from the gap for r/ g = 5 at 356° 59 Figure 24 Pressure contours at the exit from the gap for r/g = 0 at 356° 60 4.2 Damping force diagrams and curves The shock absorber characteristics under various conditions are studied using damping force diagrams. The simulations for such studies were done with constant speed and sinusoidal forcing functions. Both Change Grid (CG) and user subroutine NEWXYZ methods (see section 2.6) were used to define the mesh motion based on the conservation of mass. For Newtonian fluids both laminar and turbulent flows were considered. The Power law and the modified Bingham models were used for the Non-Newtonian fluid calculations. The shortcomings of STAR-CD and means of overcoming them were also discussed. The following results were obtained using the standard geometry of r/g = 6.8, except for the cases in which the geometry effects were considered. 4.2.1. Constant speed forcing Initially a constant speed forcing function was used to simulate the forcing conditions. Each forcing speed would provide a single data point for the damping force vs. piston speed curves. Though it was a very inefficient way to generate the necessary data, this method provided means to identify the shortcomings of STAR-CD as discussed later. The mesh motion was defined, by using the Change Grid (CG) command, based on the conservation of mass. The damping forces for the Newtonian incompressible fluids with the constant speed forcing of 2, 1.6, 1.2, 0.8 and 0.4 m/s are shown in Figures 25. The expected behavior for the damping forces is to remain constant as time increases. But, as can be seen from the figure, after an initial period of time, the forces started to oscillate. The mean of the oscillation remains approximately a constant, which is close to the initial converged damping force. This was particularly true for the constant 61 + 0.4 mls + 0.8 m/s -- A 7 1.2 mls + 1.6 m/s + 2.0 mls Figure 25 The damping force vs. time for different constant speed forcing functions with the mesh motion defined using the Change Grid (CG) commands 62 speed forcing with speeds of 1.6, 1.2, 0.8 and 0.4 m/s. For the case of 2 m/s, the force never stayed constant but oscillated at a smaller amplitude and after a certain period of time, the amplitude of oscillation increased. These observations indicated that as time elapsed the damping force started to oscillate between two values, while keeping the mean approximately a constant. The time steps used for the above constant speed forcing functions are tabulated in Table 6. The time steps were of order 10 ’ 4 to 10 ’ 5. It seems that the amplitude of oscillations depended on the time step; smaller the time steps larger the amplitudes of oscillation. Table 6 Time steps used with various constant speed piston forcing functions Speed /(m/s) time step size /(s) 2.0 7.06E-5 1.6 8.82E-5 1.2 1.18E-4 0.8 1.76E-4 0.4 3.53E-4 The default values for the solution controls were used in the above simulations; the number of PISO correctors was 20; error tolerances were 0.1 and 0.05 for velocities and pressure respectively; number of sweeps was 100 for velocities and 1000 for pressure; discretisation was done by Upwind Differencing scheme (UD). The maximum Courant number was in the range of 56 to 58 and the mean Courant number was in the 63 range of 16 to 18. It was found that the number of PISO correctors used, exceeded the limit when the oscillation started. This indicates that the solution might not have been fully converged. The time step, therefore, was reduced by a half. The reduced time step, however, did not produce a converged solution. The time step was, then, further reduced to At = 1.47E-5 which is 1/8th of the original time step. The damping force vs. time for the time step of 1.47E-5, for the case of 1.2 m/s is given in Figure 26. (The result was obtained for a different geometry, hence the mean value is different.) 0.00 0.01 0.01 0.02 Figure 26 The damping force vs. time for 1.2 m/s constant speed forcing function with a time step of 1.47E-53. There was no improvement on the results by the reduction of the time step. Figure 27 indicates that the PISO correctors did not exceed the limit as before, and hence the solution was fully converged. 0.00 0.01 0.01 0.02 0.02 Figure 27 The number of PISO correctors vs. time for 1.2 m/s constant speed forcing function with the time step of 1.47E-5s. Next the error tolerances (0.001 for velocities and 0.0005 for pressure) and number of sweeps (1000 for velocities and 10000 for pressure) were adjusted. Various discretising schemes, such as Linear Upwind Differencing (LUD), Central Differencing (CD), Self-Filtered Central Differencing (SFCD), were also tested. None of these measures achieved the non-oscillatory results either. Later, it was found out, instead of a smaller time step, a larger time step (10 ' 3 to 10 ' 4), is needed to avoid oscillations. Once the cause for the oscillations in the damping force calculations had been identified, the constant speed forcing was not pursued further, as mentioned earlier, because of its inefficient way to generate the necessary data to study the damping characteristics. 65 4.2.2. Sinusoidal forcing 4.2.2.1 Mesh motion by Change Grid commands The oscillations of the damping forces in the simulation were not only restricted to the constant speed forcing functions but also for the sinusoidal forcing functions. When the Change Grid (CG) commands for mesh motion is used with a small time step (10 ' 4 to 10 ’ 5) the damping force as a function of the piston displacement exhibits the oscillatory behavior as shown in Figure 28. The time step used in this calculation was At = 2.22E-4. It can be seen that the mean values of the oscillations follow an ellipsoidal path. Note that for the last (approximate) 1200 of the motion of the piston, the amplitude of the oscillation has increased. By using a larger time step of At = 4.44E-4, as discussed previously, the amplitude of oscillations are minimized as indicated by Figure 29. As the time step increases, the amplitude of oscillations decreases. 225 245 265 285 305 325 Figure 28 The damping forces vs. piston displacement for a time step of 2.22E-4s. with the mesh motion defined by the CG commands 66 225 245 265 285 305 325 Figure 29 The damping forces vs. piston displacement for a time step of 4.44E-4s. with the mesh motion defined by the CG commands. 225 245 265 285 305 325 Figure 30 The damping forces vs. piston displacement for a time step of 8.87E-4s. with the mesh motion defined by the CG commands. 67 225 245 265 285 305 325 Figure 31 The damping forces vs. piston displacement for a time step of 1.13E-3s. with the mesh motion defined by the CG commands. This can be seen in Figure 30 for a time step of At = 8.87E-4 and in Figure 31 for a time step of At = 1.33E-3. It can be concluded from the above investigation, when the mesh motion is defined using the Change Grid (CG) commands, the time step has a major influence on the calculations. When the time step is small (At 5 10 ' 4 to 10 ' 5) the damping force shows oscillations and the oscillations, however, can be eliminated by increasing the time step (At 5 10 ' 3). The amplitude of oscillations, as mentioned for the constant speed forcing functions, increases when the time step decreases. Hence, CG commands are not suitable for mesh motion definition if a smaller time step has to be used to capture the transient flow details. The exact cause for this phenomenon is not known but may be due to the fact that CG commands are defined in PROSTAR, where single precision is used. 68 4.2.2.2 Mesh motion by user subroutine NEWXYZ For all the cases considered so far, the mesh motion was defined, using the Change Grid (CG) commands in PROSTAR based on the conservation of mass. It was concluded that the oscillations of forces, or the pressure of the flow field, exist when the CG commands are used with a small time step. It was suggested that the oscillations could be eliminated, by defining the mesh motion in the user subroutine NEWXYZ, which is double precision. The following results were obtained using the user subroutine NEWXYZ to define the mesh motion based on the conservation of mass. For a time step of At = 1.77 E4, the damping force vs. piston speed is shown in Figure 32. Though in this case piston speed is plotted rather than the piston displacement, it could be observed that the damping force shows no oscillatory behavior even with a smaller time step. Further reduction of time step did not produce oscillation either. This can be seen in Figure 33, damping force vs. piston speed for a time step of 1.77E-5. It can be concluded that using the user subroutine NEWXYZ avoids the oscillatory behavior of the damping forces even with a smaller time step (At = 104) for the transient analysis. 69 20 , , 15 i A ; I 10”— ifi l . _ l ' l , -1O -15 -2.00 -1.50 -1.00 ~0.50 0.00 0.50 1.00 1.50 2.00 Figure 32 The damping forces vs. piston speed for a time step of 1.77E-4s. with the mesh motion defined by the user-subroutine NEWXYZ. 15 10 0.00 0.50 1.00 1.50 2.00 Figure 33 The damping forces vs. piston speed for a time step of 1.77E-5s. with the mesh motion defined by the user-subroutine NEWXYZ. 70 4.2.3 Comparison of numerical and experimental results Numerical results need to be validated by the experimental results. It will be the basis for the subsequent numerical simulations under similar conditions. The numerical results of Case 1 of Table 1 are compared with that of the experimental data in Figure 34. Figure 34 Comparison of experimental and numerical results of damping force vs. piston speed for Case 1 of Table 1. It can be seen that the experimental and the numerical results agree very well with each other. Next, the results of Cases 2 and 3 of Table l are compared with those of experimental data. It was mentioned earlier that the flow inside the gap for these two cases, is in the transition stage from laminar to turbulent. First the numerical results from a laminar flow assumption are compared with the corresponding experimental data. This is shown in Figure 35. It can be seen from Figure 35 that the numerical results over predicts the damping force by about 10 to 15% at the maximum forcing speed. Second the 71 + CaseZ-CFD + Case3-CFD + CaseZ-TEST -l- Case3-TEST 0.00 -0.50 -1.00 -1.50 -2.00 -2.50 Figure 35 Comparison of experimental and numerical results of damping force vs. piston speed for Cases 2 and 3 of Table 1 under a laminar assumption. numerical results from a turbulent flow assumption are compared with the experimental data (see Figure 36). The turbulent calculations were done using the standard k-e model with a characteristic length of 0.735mm (gap width g) and a ’5’ value of 9.0 for the log law walL The time step for the turbulent calculation was At = 3.0E-5. The maximum Courant number was in the range of 35 to 55. The agreement between the experimental data and numerical data is better under the turbulent assumption than under the laminar assumption. In fact the Reynolds numbers for Cases 2 and 3 are 2600 and 2400, and the critical Reynolds number for an annular region is 2100 as mentioned in Chapter 2. Hence the turbulent assumption is more appropriate. 72 + CaseZ-CFD + CaseS-CFD + Case2-TEST -I— Cases-TEST 0.00 -0.50 -1.00 -1.50 -2.00 -2.50 Figure 36 Comparison of experimental and numerical results of damping force vs. piston speed for Cases 2 and 3 of Table 1 under a turbulent assumption. 4.2.4 Dimensional analysis From the dimensional analysis, it can be shown that the dimensionless damping force should be a function of the Reynolds number only for a given geometry. In order to study the Reynolds number effect in the numerical simulation, the results of Case 1 and Case 4 of Table 1 are compared. In Figure 37, damping forces vs. piston speed for the two cases are compared first in a dimensional form. Obviously different damping curves are shown. It the damping force is non-dimentionalized by 1/2 mezA, where A is a characteristic area and Vm is the maximum piston speed, the dimensionless damping curves overlap for the two cases (see Figure 38). This indicates that for a shock absorber of a given geometry, the damping curve is only a function of Reynolds number. The numerical results substantiate that of the dimensional analysis. .2, +case 4 +case 1 Figure 37 Damping forces vs. piston speed for Cases 1 and 4 of Table 1 -0.002 A‘:‘_ a (W251, -0. 6 ff I —l— ca: 1 -0. ’JJI‘S':‘II- -001 -1 .50 -1 . 1150 0- Figure 38 Non-dimensionalized damping force vs. non-dimensionalized piston speed for Cases 1 and 4 of Table l 74 4.2.5 Effect of the Entrance / Exit geometry The damping force vs. piston speed for different cases of r/g ratios is plotted in Figure 39. A time step of 8.87E-4 was used for these calculations. The force curve for the square edge (rlg = 0) has a higher force value at a given forcing speed than that of rounded edges, and at the maximum speed the force for the square edge is about 15% higher. The force curves for the rounded edges overlap each other indicating that the magnitude of the radius of curvature of the entrance/exit has a lesser influence on the forces than the shape of the edge. 2) ~15 / / I] -10 x" . --n’g=10 " +rrg=5 + rlg: 1 ‘5 ' +trg=0 O I 000 0.5) -1.(D -1 H) -2.(D Figure 39 Comparison of the damping forces vs. piston speed for different cases of r/ g. r4' 75 4.2.6 Non-Newtonian 4.2.6.1 The power law For the Non-Newtonian fluid calculations, the power law and the modified Bingham models were used. As mentioned earlier, the power law was used to evaluate the capabilities of the software, while the desired Non-Newtonian model is the modified Bingham plastic. The simulation was done by assuming a Non-Newtonian fluid behavior in the gap and a Newtonian one elsewhere. The constants used with the power law model are given in Table 4. The Case 3 of Table 4 represents the lowest effective viscosity (hypothetical). The parameters for this case are given by m = 24 E —3 and n = 0.5. The damping force vs. piston speed calculated using these parameters with a time step of 4.0 E —4 is plotted in Figure 40. The comparison of the effective viscosity calculated with these parameters, to that of the fluid used in the experiments, (Case 2, Table 5) is shown in Figure 6. The lower limit for the second invariant during the simulation of Case 3 of Table 4 lies in the region of 10 ‘ 4 to 10 ‘ 6. The effective viscosity does not reach the desired region with these parameters at all and therefore the parameters were changed to m = 87 and n = 0.85 in order to obtain an effective viscosity close to that of the real fluid. With a time step of 3.0 E-4 the damping force vs. piston speed was obtained and is shown in Figure 41. 76 -2.50 -2.00 -1.50 -1.00 -0.50 0.00 4. -2.50 -2.00 -1.50 -1.00 -0.50 0.00 Figure 41 The damping force vs. piston speed for the power law model of Case 2 of Table 4. Though the effective viscosity calculated with these parameters has reached the lower limit of the effective viscosity of the real fluid, it has not reached the upper limit yet 77 and hence the parameters for the power law model were changed again to m = 1350 and n = 0.999. This case is listed under Case 1 of Table 4 and has the highest effective viscosity. The effective viscosity calculated using these parameters lies in the upper limit for the real fluid. As 11 is close to unity, it simulates a fluid whose behavior could be approximated by a Newtonian model. The damping curve with a time step of 3.0 E-4 is given by Figure 42. It could be observed that as the effective viscosity increases the solution tends to be partially converged. 300000 ’ -100000 ~400000 ~500000 -2.50 -2.00 -1.50 -1 .00 -0.50 0.00 Figure 42 The damping force vs. piston speed for the power law model of Case 1 of Table 4 for a time step of 3.0 E-4s. It was found out later, that STAR-CD could not resolve for higher effective viscosities with a time step of order 104. and a converged solution could only be obtained with an extremely small time step. The calculation was repeated with a time step of 5.33 E-7 and the results are plotted in Figure 43. This very small time step is impracticable and 78 the calculation had to be terminated. For these calculations there is no test results to validate the numerical results. But the basic trend of the damping force vs. piston speed plots should be the same as that of the stress vs. strain rate plots. The trend should be parabolic for Cases 2 and 3, and a straight line for Case 1. Hence the above results are acceptable. STAR-CD could handle the power law model for Non—Newtonian fluid behavior only for a relatively low effective viscosity. The simulation becomes impractical for high effective viscosity. 340000 1 330000 320000 E==T 310000 300000.. 290000 ‘1 280000 270000 1.9977 1.9978 'L9979 1.9980 1.9981 1.9982 Figure 43 The damping force vs. piston speed for the power law model of Case 1 of Table 4 for a time step of 5.33 E-7s. 4.4.2 The Modified Bingham Next the modified Bingham model was used in the gap to simulate the Non- Newtonian behavior. The parameters used with this model are given in Table 5. The behavior of the second invariant for Case 2 indicated in Table 5 is similar to that of case 3 79 of Table 4 (the Power law). As discussed earlier this case has the lowest effective viscosity. The damping force vs. piston speed for a time step of 1.0 E —4 is given in Figure 49. Note that the magnitude of the force is different from that of the power law because of the different densities used for the base fluid. However, the trend is parabolic indicating that the modified Bingham model is capable of producing acceptable results when the effective viscosity is low. -2.00 -1.50 -1.00 -0.50 0.00 Figure 44 The damping force vs. piston speed for the modified Bingham model of Case 2 of Table 5 for a time step of 1.0 E-4s. After verifying that the modified Bingham model, the parameters in the model were changed to that of the fluid used in the experiments. These parameters are listed in Case 1 of Table 5. They are to = 34.5 E+3 Pa, 11 0 = 31.97 E-3 Pas, and y: 25 s'I. Again an extremely small time step had to be used in order for STAR~CD to resolve the flow field. 80 Figure 45 shows the damping force vs. piston speed for the time step of 5.33 E -7. 0.000 -0.020 -0.040 -0.060 -0.080 -0.100 -O.120 Figure 45 The damping force vs. piston speed for the modified Bingham model of Case 1 of Table 4 for a time step of 5.33 E78. It could be concluded that the STAR-CD capabilities of handling Non-Newtonian fluid are limited; for higher effective viscosities an impracticably small time step has to be used. Chapter 7 SUMMARY AND CONCLUSIONS The internal fluid flows of a controllable shock absorber have been studied by numerical simulation. The effects of the geometry at the entrance/exit regions on the fluid motion and the damping force were investigated. Both laminar and turbulent flow assumptions were used to calculate the damping forces for the standard geometry. Finally a combination of Newtonian and Non-Newtonian fluids were considered for the calculation of the damping forces. All the above calculations were done under an incompressible assumption. A summary of the results is presented below: 0 The case with the ratio of r/g = 0 does not provide a smooth entrance for the flow into the gap and creates a large pressure drop. As the ratio rlg increases, a smooth entrance to the flow is achieved and the pressure drop decreases, however the exit jet becomes stronger. 0 It has been shown that the geometric ratio r/g has not only affected the flow field inside the shock absorber but also the damping characteristics. The damping force for the straight edge (r/g = 0) has a higher value at a given piston speed and is about 15% higher than that of the rounded edges at the maximum forcing speed. 81 82 For the rounded edges there is no appreciable variation among the force curves, indicating that magnitude of the radius of curvatures is not a significant factor contributing to the forces. Using Change Grid (CG) commands to define the mesh motion is not recommended with a small time step. The user subroutine NEWXYZ overcomes the shortcomings of the earlier method. The dimensionless damping characteristic is a function of the Reynolds number only, for a shock absorber of a given geometry. Both laminar and turbulent models for Newtonian fluids produce acceptable results. The damping force calculated using the standard K-e model, with the given parameters and the gap width (g) as the characteristic length, is less than that of the corresponding laminar flow calculations and agrees better with the experimental results. The Non-Newtonian calculation is only successful when the effective viscosity is small and for higher values of effective viscosity (with large gradients) the solution fails to converge with a reasonable time step. An extremely small time step is required but this might not be a viable solution as the intense use of CPU resources. The model selected for the definition of the Non-Newtonian behavior does not influence the convergence. A modification of the current version of STAR—CD is necessary to provide a converged result for Non-Newtonian fluids with a reasonable time step. APPENDICES APPENDIX A APPENDIX A USER SUBROUTINE FOR BOUNDARY DEFINITION C************‘k'k******'k************************'k************************* SUBROUTINE BCDEFW(U,V,W,TORHF,SCALAR,RESWT,RSTSC) C Boundary conditions definition for walls C*********************************************************************** C _______________________________________________________________________ C STAR VERSION 2.300 C _______________________________________________________________________ INCLUDE ’comdb.inc' COMMON/USROOl/INTFLG(100) DIMENSION SCALAR(SO),RSTSC(SO) DIMENSION SCALC(50) INCLUDE ’usrdat.inc' EQUIVALENCE( UDAT12(OOl), ICTID ) EQUIVALENCE( UDAT01(OO6), ELOG ) EQUIVALENCE( UDAT02(O70), X ) EQUIVALENCE( UDAT02(071), Y ) EQUIVALENCE( UDAT02(O72), Z ) EQUIVALENCE( UDATO4(002), DENC ) EQUIVALENCE( UDATO4(003), EDC ) EQUIVALENCE( UDATO4(OOS), PRC ) EQUIVALENCE( UDATO4(OO9), SCALC(Ol) ) EQUIVALENCE( UDAT04(OO7), TC ) EQUIVALENCE( UDATO4(008), TEC ) EQUIVALENCE( UDAT04(059), UC ) EQUIVALENCE( UDATO4(060), VC ) EQUIVALENCE( UDATO4(061), WC ) EQUIVALENCE( UDATO4(O64), UCL ) EQUIVALENCE( UDATO4(O65), VCL ) EQUIVALENCE( UDATO4(O66), WCL ) This subroutine enables the user to specify WALL boundary conditions for U,V,W,TORHF,SCALAR,RESWT and RSTSC(IS). RESWT, RSTSC(IS) C C C C C ** Parameters to be returned to STAR: U,V,W,TORHF,SCALAR, C C C IF(IREG.EQ.2) THEN U=O. V=O. W: -l.998103193*SIN(39.33274002*TIME) ELSE END IF C _______________________________________________________________________ RETURN END C 83 APPENDIX B APPENDIX B USER SUBROUTINE FOR MODIFIED BINGHAM MODEL C*************************t*****‘k'k‘k'k‘k'k********************************** SUBROUTINE VISMOLiVISM) C Viscosity (molecular) C**********************'k*iri't******************************************** C _______________________________________________________________________ C STAR RELEASE 3.000 C _______________________________________________________________________ INCLUDE ’Comdb.inc’ COMMON/USROOl/INTFLG(100) INCLUDE ’usrdat.inc' DIMENSION SCALAR(50) EQUIVALENCE( UDAT12(OOl), ICTID ) EQUIVALENCE( UDATO3(OO9), DUDX EQUIVALENCE( UDATO3(010), DVDX EQUIVALENCE( UDATO3(Oll), DWDX EQUIVALENCE( UDATO3(012), DUDY EQUIVALENCE( UDATO3(013), DVDY EQUIVALENCE( UDATO3(Ol4), DWDY EQUIVALENCE( UDATO3(015), DUDZ EQUIVALENCE( UDATO3(016), DVDZ EQUIVALENCE( UDATO3(Ol7), DWDZ EQUIVALENCE( UDATO3(018), SECINV ) EQUIVALENCE( UDAT11(001), CP ) EQUIVALENCE( UDAT11(002), DEN ) EQUIVALENCE( UDATll(OO6), P ) EQUIVALENCE( UDAT11(OO7), T ) EQUIVALENCE( UDAT11(OO9), SCALAR(Ol) ) EQUIVALENCE( UDAT11(OS9), U EQUIVALENCE( UDATll(O60), V EQUIVALENCE( UDATll(O6l), W X Y Z vvvvvvvvv EQUIVALENCE( UDATll(O67), EQUIVALENCE( UDATll