Simulation of Hydrodynamics of the Jet Impingement using Arbitrary Lagrangian Eulerian Formulation by Hamid Maghzian B.Sc., Sharif University of Technology, Tehran, Iran, 1998 M.Sc., Sharif University of Technology,Tehran, Iran, 2000 A THESIS IN PARTIAL FULFILLMENT OF THE REQUIREMENT FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES (Mechanical Engineering) THE UNIVERSITY OF BRITISH COLUMBIA October 2007 © Hamid Maghzian, 2007 Abstract Controlled cooling is an important part of steel production industry that affects the properties of the outcome steel. Many of the researches done in controlled cooling are experimental. Due to progress in the numerical techniques and high cost of experimental works in this field the numerical work seems more feasible. Heat transfer analysis is the necessary element of successful controlled cooling and ultimately achievement of novel properties in steel. Heat transfer on the surface of the plate normally contains different regimes such as film boiling, nucleate boiling, transition boiling and radiation heat transfer. This makes the analysis more complicated. In order to perform the heat transfer analysis often empirical correlations are being used. In these correlations the velocity and pressure within the fluid domain is involved. Therefore in order to obtain a better understanding of heat transfer process, study of hydrodynamics of the fluid becomes necessary. Circular jet due to its high efficiency has been used vastly in the industry. Although some experimental studies of round jet arrays have been done, yet the characteristics of a single jet with industrial geometric and flow parameters on the surface of a flat plate is not fully understood. Study of hydrodynamics of the jet impingement is the first step to achieve better undrstanding of heat transfer process. Finite element method as a popular numerical method has been used vastly to simulate different domains. Traditional approaches of finite element method, Lagrangian and Eulerian, each has its own benefits and drawbacks. Lagrangian approach has been used widely in solid domains and Eulerian approach has been widely used in fluid fields. Jet impingement problem, due to its unknown free surface and the change in the boundary, falls in the category of special problems and none of the traditional approaches is suitable for this application. The Arbitrary Lagrangian Eulerian (ALE) formulation has emerged as a technique that can alleviate many of the shortcomings of the traditional Lagrangian and Eulerian formulations in handling these types of problems. Using the ALE formulation the computational grid need not adhere to the material (Lagrangian) nor be fixed in space (Eulerian) but can be moved arbitrarily. Two distinct techniques are ii being used to implement the ALE formulation, namely the operator split approach and the fully coupled approach. This thesis presents a fully coupled ALE formulation for the simulation of flow field. ALE form of Navier-Stokes equations are derived from the basic principles of continuum mechanics and conservation laws in the fluid. These formulations are then converted in to ALE finite element equations for the fluid flow. The axi-symmetric form of these equations are then derived in order to be used for jet impingement application. In the ALE Formulation as the mesh or the computational grid can move independent of the material and space, an additional set of unknowns representing mesh movement appears in the equations. Prescribing a mesh motion scheme in order to define these unknowns is problem-dependent and has not been yet generalized for all applications. After investigating different methods, the Winslow method is chosen for jet impingement application. This method is based on adding a specific set of partial differential Equations (Laplace equations) to the existing equations in order to obtain enough equations for the unknowns. Then these set of PDEs are converted to finite element equations and derived in axi-symmetric form to be used in jet impingement application. These equations together with the field equations are then applied to jet impingement problem. Due to the number of equations and nonlinearity of the field equations the solution of the problem faces some challenges in terms of convergence characteristics and modeling strategies. Some suggestions are made to deal with these challenges and convergence problems. Finally the numerical treatment and results of analyzing hydrodynamics of the Jet Impingement is presented. The work in this thesis is confined to the numerical simulation of the jet impingement and the specifications of an industrial test setup only have been used in order to obtain the parameters of the numerical model. iii Table of Contents Abstract^ ii Table of Contents^ iv List of Figures vi Nomenclature^ ix Acknowledgment xi 1 Introduction^ 1 1.1 Motivation 1 1.2 Steel Strip Production^ 1 1.3 Controlled Cooling 2 1.4 Cooling Systems^ 5 1.5 Circular Jets 7 1.6 Hydrodynamics of the jet impingement^ 8 1.7 Traditional Numerical Approaches 11 1.8 Scope of work^ 12 2 Literature Review 15 2.1 Overview of Jet Impingement Research^ 15 2.1.1^Heat Transfer and Hydrodynamics of the Jet^ 17 2.2 Numerical Techniques and Formulation^ 23 2.2.1^Numerical Techniques^ 23 2.2.2^Traditional Numerical Formulations^ 25 2.2.3^The Arbitrary Lagrangian Eulerian (ALE) Formulation^ 26 3 Finite Element ALE Equations for Jet Impingement^ 30 3.1 Flow Equations^ 30 3.1.1 Constitutive Model^ 30 3.1.2 Conservation Laws 33 iv 3.1.3 Field Equations^ 34 3.2 Boundary Conditions 35 3.3 ALE (Arbitrary Lagrangian Eulerian) Kinematics^ 36 3.3.1 Material and Spatial Configurations 36 3.3.2 Referential configuration^ 38 3.3.3 Time Derivatives in Different Configurations and their Relations^40 3.4 Derivation of ALE Finite Element Equations in the Flow Field^43 3.4.1 Transformation to Referential Domain^ 43 3.4.2 Conversion to FE Equations^ 44 3.4.3 ALE Finite Element Equations of the Flow in Axi-Symmetric Form^47 4 Mesh Motion Schemes^ 51 4.1 Treatment of Mesh Movement 51 4.2 Smoothing Technique^ 53 4.3 Implementation of Winslow method to Finite Element^ 55 4.4 Boundary Conditions for Mesh Motion Equations 56 5 Application and Results: Jet impingement problem^58 5.1 Overview of Experimental Program^ 58 5.2 Analysis setup and Results^ 59 5.2.1 Free Jet Region 60 5.2.1.1 Geometry and Boundary Conditions and Meshing^60 5.2.1.2 Approximation and Integration of Field Equation 63 5.2.1.3 Detailed Solution Procedures and Results^64 5.2.2 Impingement Region^ 79 5.2.2.1 Geometry and Boundary Conditions and Meshing^ 80 5.2.2.2 Detailed Solution Procedures and Results^ 83 5.2.2.3 Analytical verification of the solution 107 6 Conclusions and Future works^ 111 Bibliography^ 117 v List of Figures Figure 1.1 Figure 1.2 Figure 1.3 Figure 1.4 Figure 1.5 Figure 5.1 Figure 5.2 Figure 5.3 Figure 5.4 Figure 5.5 Figure 5.6 Figure 5.7 Figure 5.8 Figure 5.9 Figure 5.10 Figure 5.11 Figure 5.12 Figure 5.13 Steel Strip Production^ 2 Various types of cooling systems^ 5 Jet impingement on hot steel plate 8 Jet Configurations^ 9 Inviscid pressure and velocity profile and flow regions^ 10 Schematic layout of the Run Out Table at UBC 59 Free Jet Region (Initial Domain)^ 61 Velocity-Steady State Solution of Free Jet Region (phase one)^ 65 Velocity-Transient ALE Solution of Free Jet Region Problem (Time= 0, 0.5, 0.1, 0.2 s)^ 68 Velocity-Transient ALE Solution of Free Jet Region Problem (Time= 0.4, 0.6, 0.8, 2 s)^ 69 Pressure -ALE Solution of free Jet Region Problem (Time=2 s)^70 Displacement dr- ALE Solution of Free Jet Region Problem (Time=0.1, 0.2, 0.3, 0.4 s)^ 71 Displacement dr- ALE Solution of Free Jet Region Problem (Time=0.5, 1, 2 s)^ 72 Velocity-Jet End- ALE Solution of Free Jet Region Problem (Time=0, 0.1, 0.2, 0.3 s)^ 74 Velocity-Jet End- ALE Solution of Free Jet Region Problem (Time=0.4, 0.45, 0.47, 0.5 s)^ 75 Velocity-Jet End- ALE Solution of Free Jet Region Problem (Time=0.6, 0.7, 1 s)^ 76 Deformed Mesh, End Jet, ALE Solution of Free Jet Region Problem (Time=0, 0.4, 0.45, 0.5 s)^ 77 Deformed Mesh, End Jet, ALE Solution of Free Jet Region Problem (Time=0.55, 0.6, 0.7, 1 s) ^ 78 vi Figure 5.14 Figure 5.15 Figure 5.16 Figure 5.17 Figure 5.18 Figure 5.19 Figure 5.20 Figure 5.21 Figure 5.22 Figure 5.23 Figure 5.24 Figure 5.25 Figure 5.26 Figure 5.27 Figure 5.28 Figure 5.29 Figure 5.30 Figure 5.31 Figure 5.32 Figure 5.33 Figure 5.34 Figure 5.35 Figure 5.36 Figure 5.37 Figure 5.38 Figure 5.39 Figure 5.40 Figure 5.41 Figure 5.42 Jet Impingement Region^ 80 Boundary Mode, Jet Impingement Region^ 81 Meshing, Jet Impingement Region 82 Velocity, Phase one, steady state Jet Impingement Region^85 Pressure, Phase one, steady state Jet Impingement Region^ 86 Velocity, Second step, ALE Solution, Time=0.1 s^88 Velocity, Second step, ALE Solution, Time=0.2 s 89 Velocity, Second step, ALE Solution, Timc-0.5 s^89 Velocity, Second step, ALE Solution, Time=2 s 90 Velocity, Second step, ALE Solution, Time=4 s^90 Velocity, Second step, ALE Solution, Timc-6 s 91 Velocity, Second step, ALE Solution, Time=8 s^91 Velocity, Second step, ALE Solution, Time=10 s 92 Velocity, Second step, ALE Solution, Time-11 s^92 Velocity, Second step, ALE Solution, Time=13 s 93 Pressure, Second step, ALE Solution, Time=0.1 s^93 Pressure, Second step, ALE Solution, Time=0.2 s 94 Pressure, Second step, ALE Solution, Time=0.5 s^94 Pressure, Second step, ALE Solution, Time=2 s 95 Pressure, Second step, ALE Solution, Time=4 s^95 Pressure, Second step, ALE Solution, Time=6 s 96 Pressure, Second step, ALE Solution, Time=8 s^96 Pressure, Second step, ALE Solution, Time=10 s 97 Pressure, Second step, ALE Solution, Time=11 s^97 Pressure, Second step, ALE Solution, Time=13 s 98 Free Surface Movement, ALE Solution, Time=0.01 s^99 Free Surface Movement, ALE Solution, Time=0.02 s 99 Free Surface Movement, ALE Solution, Time=0.05 s^ 100 Free Surface Movement, ALE Solution, Time=0.07 s 100 vii Figure 5.43 Figure 5.44 Figure 5.45 Figure 5.46 Figure 5.47 Figure 5.48 Figure 5.49 Figure 5.50 Figure 5.51 Figure 5.52 Figure 5.53 Figure 5.54 Figure 5.55 Figure 5.56 Free Surface Movement, ALE Solution, Time=0.15 s^ Free Surface Movement, ALE Solution, Time=0.2 s Free Surface Movement, ALE Solution, Time=0.3 s^ Mesh Deformation, Curve, ALE Solution, Time=0 s Mesh Deformation, Curve, ALE Solution, Mesh Deformation, Curve, ALE Solution, Mesh Deformation, Curve, ALE Solution, Mesh Deformation, Curve, ALE Solution, Mesh Deformation, Curve, ALE Solution, Mesh Deformation, Curve, ALE Solution, Mesh Deformation, Curve, ALE Solution, Mesh Deformation, Curve, ALE Solution, Surface Pressure, Impingement Region^ Velocity Profile, Film Flow^ Time=0.01 s^ Time=0.02 s Time=0.04 s^ Time=0.06 s Time=0.08 s^ Time=0.1 s Time=0.2 s^ Time=0.3 s 101 101 102 102 103 103 104 104 105 105 106 106 110 110 viii Nomenclature 1^unit or identity tensor Ti,'^deviator stress tensor bulk viscosity D'y^deviator rate of deformation tensor mean pressure D^rate of deformation tensor D ij^components of rate of deformation tensor C^4th order isotropic tensor Cod components of 4 th order isotropic tensor v^velocity vector v,^velocity component P^pressure total stress tensor component T^total stress tensor kronecker delta function p^density b. ^acceleration component u g^mesh(grid) velocity in horizontal direction vg^mesh(grid) velocity in vertical direction vg^mesh(grid) velocity vector F^non-linear tensoric function velocity component time derivative t^small time step v^kinematic viscosity ix 1-1^dynamic viscosity X^material coordinates x^spatial coordiantes n unit normal vector to the surface(boundary) b^body acceleration vector J^jacobian matrix of mapping between material and spatial domains 3^jacobian matrix of mapping between referential and spatial domains c^convective velocity relative to mesh X^referential coordinates 51„,o,^moving domain O ft,^fixed or original domain D^determinant of jacobian matrix dx^horizontal incremental value in spatial (2D) domain dY^vertical incremental value in spatial (2D) domain dX^horizontal incremental value in material (2D) domain dY^vertical incremental value in material (2D) domain dr^mesh displacement in horizontal direction(axi-symmetric) dz^mesh displacement in vertical direction (axi-symmetric) R^radial coordinate in material domain Z^axial coordinate in material domain r^radial coordinate in spatial domain z^axial coordinate in spatial domain Acknowledgement Sincere thanks and appreciation are expressed to my supervisor, Professor Mohamed Gadala, for his invaluable advice and continuous guidance. His knowledge and experience greatly influenced this thesis. I am very glad for having the opportunity to study at the University of British Columbia. I would like to acknowledge the funding of my supervisor as well as of the Faculty of Graduate Studies. I am honored to be a recipient of the University Graduate Fellowship for two consecutive years. I would also like to thank my professors both in The University of British Columbia and Sharif University of Technology. I am grateful to my brother, Alireza, for his support and encouragement and finally a very special thanks and dedication goes to my parents. Words are incapable of their appreciation. xi Chapter 1 INTRODUCTION 1.1 Motivation The history of steel strip industry is based on a highly empirical form of technology. Perhaps one of the most expensive and time-consuming stages of this process is experimental and prototype development. This process is mainly based on trial and error. This is in a large part due to the fact that this process is too complex to model using the available tools. Due to advancements in the numerical methods it is now more feasible to model different parts of such a process. Today finite element method with its powerful capabilities and flexibilities could be used to simulate different types of complex problems. The vast increase in computational capabilities together with the advancements in numerical techniques made it possible to simulate different stages of the controlled cooling process in the production of steel plates and strips. Currently, the use of finite element method in the steel industry is practiced widely. As a by-product of the simulation, a better understanding of the process can be acquired. 1.2 Steel strip Production In the recent decades, great progress has been made in steel plate production technologies that include refining mechanical properties, and the advancement of controlled cooling. Steel strip, which is one of the most important steel products, is widely used in many industries, for example car and can industry. For optimum use of the material, it is required that steel strip has high strength and good toughness as well as better formability and weldability. 1 04, 0 414) al• /-1 4 I I 0 id V%9/V AtiAh Qehoot ingurnace Roughing ^Finishini)Mitts (Run—Ou) C.)Maio Mitts Most of steel products (i.e., plates, strips and tubes) are produced by hot rolling. The complete steel strip production (Figure 1.1) is composed of reheating furnace, roughing mills, finishing mills, run-out table (ROT) and coiler. Figure 1.1 Steel Strip Production [61] In the initial stage, steel strip is produced as continuous-cast slab with thickness of about 250 mm. In the reheating stage, the slab is normally heated to about 1200 °C. In roughing mills stage, the slab is rolled several passes and the thickness is reduced to about 50 mm Finishing mills usually consist of 5-7 stands and the final thickness of the strip is between 1 mm to 25 mm with a temperature of about 850-950 °C. After the finishing rolls, the strip is cooled down to the coiling temperature which is about 600 °C. All the micro-structural phase changes of the steel will occur in this cooling process. The properties of the outcome steel would be highly dependent on the phase change and the cooling characteristics and the only way to control this process is to control the cooling rates on the run out table. Therefore, controlled cooling has a significant contribution in the improvement of the outcome properties. After the controlled cooling on run out table, the strip is coiled and left for air-cooling. 1.3 Controlled Cooling Typically 50-100 meter in length, a "run out table" is made up of a set of cooling banks placed above and under the strips, and an array of motorized rollers. Cooling systems are required to 2 have a large cooling efficiency, and to be capable of producing the desired microstructures and mechanical properties. Controlled cooling of steel refers to the cooling process on run out table after finishing hot rolling where the cooling start and stop temperatures and the cooling rate are exercised accurately according to the expected property. While the geometry and surface quality are influenced by the rolling deformation, microstructure and mechanical properties highly depend on the cooling process applied immediately after the finishing stand. The main benefits of controlled cooling may include grain refinement, reduction of carbon content and controlling mechanical properties. Finer grain size can increase the strength of steel as well as notch toughness and resistance to brittle fracture. The reduction of carbon content improves formability and contributes to superior notch toughness and weldability. The controlled cooling process entails, however, some difficulties. Uniform cooling must be maintained in every location of the steel plate including top and bottom surfaces, edge and center, and leading and trailing ends. Also plate distortion and resulting residual stress due to the cooling process must be eliminated or minimized. In addition to that accurate control of cooling rate, cooling start and finish temperatures should be achieved. The above requirements are not simple ones because, for example, there always exist a temperature difference between the center and the surfaces of a plate even immediately after the finishing roll, and when water is applied to the top (or top and bottom) surface of the plate, the problem is exaggerated. Non-uniform cooling of the plate will cause deflection, edge wave and residual stress as well as non-uniform mechanical properties. In order to achieve uniform cooling, various measures have been taken in industry. These include: balancing water for the top and bottom surfaces at an optimum ratio, setting appropriate 3 transverse water amount gradient from center to edge along with edge masking devices to prevent transverse temperature gradient and buckling leading to edge waves, and controlling water amount for leading and trailing ends to prevent excessive temperature drop. All of these measures are highly experimental and the only way to apply them is by trial and error. Small changes in the geometry and conditions can cause significant change in the results. In Controlled Cooling, the starting temperature, the finishing temperature and the cooling rate for the steel must be controlled in order to achieve the balanced and appropriate properties required for the industry. Therefore Heat Transfer process on the "run out table" becomes significantly important. The amount of heat extracted from the plate is mostly achieved during run out table. Since the steel on the run out table is very hot (more than 900 °C), different heat transfer regimes happens on the run out table. These regimes include film boiling, transition and contact heat transfer. This makes the heat transfer analysis more complicated. That is why the majority of works done in the area of heat transfer are empirical. Empirical correlations are used to obtain heat transfer rates. These correlations normally contain parameters such as velocity and pressure of the cooling water, which entails the study of Hydrodynamics of the jet impingement. However there is no adequate and accurate research done on the hydrodynamics of the jet. In analyzing the heat transfer, often the hydrodynamical data being used are based on the very simplified model of the cooling jet, i.e. assumption of inviscid flow. Therefore an accurate study of hydrodynamics of the jet impingement with the realistic assumptions could achieve the velocity and pressure profile of the flow. The result of such a study could be used in future studies as a mean to obtain better correlation and as a result more accurate analysis of the heat transfer in its different regimes 4 1.4 Cooling Systems There are three main types of cooling systems used for the controlled cooling process of steel strip on ROT: water spray, laminar flow (circular jet) and water curtain (planar jet), as shown in Figure 1.2 [62]. In a water spray system, the water as coolant impinges from a row of specially designed nozzles on to the strip. In the water curtain system, the strip is cooled by a planar jet, which spans the entire width of the strip. One of the advantages of this system is its improvement of the uniformity of the cooling in the strip's transverse direction. I-- 10 a --I^Pyrometers Cross spray nozzles 71^T2 i^ Th i 6^T3II VIAAAV, llIrrTS-r, r•rvi-v-)^r-r-r&rir, rrirh-ri I{ ^T4 e 42.71....„"f6 4 i s■-11..44.1,1..aii.., 0 .6**171m .^; 1^ a^11 II ° /^ \ L........—.........,,..,...............v.-..................., \\*\Laminar-^Water curtain- \\\ V‘s\^\\\ \\*1/4 V *\\ .‘!// Vo l*\'‘*//qq\\■ , "itiolv‘divt\ \‘‘0 iiIii\\\ /"/11 1 1"\‘u:;..htitwyithiiii\‘‘‘‘ Figure 1.2 Various types of cooling systems [62] The laminar flow is normally prefered because it can penetrate the vapour film on the cooled surface and increase the coolant residence time. While each of the laminar cooling apparatus consists of two or four rows of U-pipes, this system creates low-pressure, laminar flow streams, and, therefore, has a higher cooling efficiency. 5 A combination between high cooling performance and uniform strip cooling is obtained by an array of round jets impinging on the steel sheet. The jets are arranged on headers at a height of 1-3 m above the plate. Controlled cooling is achieved by switching off and on various headers along the runout table and by controlling the water temperature. Similar arrangement is normally done for cooling the bottom of the plate. Although this cooling process has been used in the steel industry for decades, it is still mainly a trial and error procedure when new products are developed. This is because the thermal events during the cooling process are very complex, and a fundamental understanding of this process is not yet complete. In a typical hot strip rolling mill, a steel strip leaves the last finishing stand at temperature ranging from 750 to 1000 °C and is rapidly transported along a run out table (ROT); the strip is subsequently subjected to air-cooling and water-cooling. Due to the high temperatures of steel strip during the water-cooling periods, the nucleate boiling is typically confined to a small region beneath the jet, while the film boiling exists over most of the surface at locations downstreamfrom the stagnation point. There is also conductive heat transfer in three directions; width, thickness and length directions within the strip.. Typical cooling conditions on a ROT for the strip hot rolling are [63]: Cooling rate : 5-150^°C/s Cooling efficiency: 10-400 kJ/kg Heat flux: 105 -107 W/m2 Heat transfer coefficient (HTC): 103 -105 w/m2 . 0C The basic requirements for the cooling system are summarized as follows. First, cooling must be uniform in every location of the top and bottom surface, the edge and center, and the leading and trailing ends. Second, any distortion and residual stress due to water cooling must be minimized. 6 Third, the cooling rate and the cooling start and finish temperatures must be accurately controlled. These requirements which are, generally, not simple to achieve and in order to pave the way for the simulation of the process, detailed heat transfer and hydrodynamic analyses for the jet impingement are required. 1.5 Circular Jets Although some experimental studies of round jet arrays have been done [35], yet the heat transfer characteristics of single jet on a hot plate is not fully understood. As discussed above, cooling of high temperature steel strip (more than 800 °C), results in various phases of heat transfer. In stationary strip jet impingement, these regimes are present simultaneously at differing distances from the jet impingement center-line. (See Fig 1.3) Convective single-phase heat transfer is the mode of heat transfer when there is no boiling. The wall temperature has to be lower than or equal to the saturation temperature of the liquid at given pressure, for single-phase convection to take place. In jet impingement heat transfer, the single- phase convective heat transfer coefficient varies over the surface as the boundary layer develops and also due to hydrodynamic variations in the flow. If the wall temperature exceeds the saturation temperature of the liquid, nucleate boiling occurs. At the onset of nucleate boiling, bubbles begin to form and detach from the hot surface, hence enhancing fluid circulation and increasing the heat transfer. In transition boiling the vapor blankets formed by the coalescence of bubbles are unstable and collapse so that the surface is intermittently wetted resulting in large fluctuations in surface temperature and heat flux. Transition boiling is partially nucleate boiling and partially unstable film boiling and is recognized as the least understood of all boiling regimes. 7 free-jet region stagnation region wall jet region Water tronsitionfnucleate bat" film boiling 1.4 cc"' h""•maarao-i^film boiling (insulated)^Xis^(insulated) Steel Plate us In the film boiling region the strip surface is separated from the jet by a vapor layer which reduces cooling effectiveness considerably. In this region, the heat transfer from the surface to the liquid is across a vapor film. First the primary mode of heat transfer is thought to be forced convection within the vapor film but as the temperature is increased, radiation becomes a more dominant mode of heat transfer. Figure 1.3 Jet impingement on hot steel plate [21] 1.6 Hydrodynamics of Jet Impingement Jet impingement can be divided into five categories from a hydrodynamic point of view: free surface jets, plunging jets, submerged jets, confined jets and wall jets. These configurations are all shown in Figure 1.4. Furthermore, impinging jets can be classified by their shape, whether the jet is oriented normally or obliquely with respect to the impingement surface or weather the impingement surface is flat, convex, or concave. A free surface jet, Figure 1-4 (a), is where the liquid travels relatively undisturbed from the nozzle to the impingement surface. This is the case for a water jet injected into air impinging on an unconstrained flat surface. Unconstrained means that the water can flow off the edges and hence, does not pool on the surface. 8 I Nozzle place NNN.V.A.V.47V,ANNV,ANSNV,AANN a. Free-surface UqUid - NANNNWV00.,00‘,000....NANNN, b. Plum Nozzle Tmcvccromoccommccevc‘ c. Submersed SVCCRTMCVMCVMCCSNSC d. Confined Gas Liquid 100.X.0"NNWSNANNANANNNN "Vs e. Wall (free-surface) Figure 1.4 Jet configurations [31] On an industrial run out table there can be either free surface or plunging jets depending on the configuration. Usually the moving metal strip on a steel mill run out table is cooled by an array of water jets. The first row of jets are free surface jets but further down stream a layer of water is formed on top of the plate and the jets act as plunging jet. As the thickness of the water layer increases the effect of impingement decreases. For maximum heat transfer, all the jets should be free surface types. Figure 1.5 shows the inviscid pressure and velocity distributions for a free surface planar jet with uniform velocity along with the varying flow regimes. 9 11 ujx) 1.0 0.6 0.6 0.4 0.2 0.0 POO-P. Nozzle Free surface Stagnation Point C.,`, ....,,N•■■■■••■■■•■■""‘,^X A. Starastion !lesion O. 4earlaratio• **Owe C. Paralisl-Flor Iteslost Figure 1.5 Inviscid pressure and velocity profile and flow regions [31] Stagnation region, acceleration region and parallel-flow region are three main regions within the domain. The stagnation region has the same size as the jet. The combination of the stagnation and acceleration regions is referred to as impingement zone. In an inviscid flow at the stagnation point the stream wise velocity is zero and then increases rapidly within the impingement zone. In a viscous flow this assumption is not valid and the velocity of the fluid is zero along the surface. In the parallel laminar flow region, the velocity approaches the velocity of the jet as approaching the surface. Further downstream the flow becomes turbulent. Unlike the velocity, pressure is at its maximum at the stagnation point and decreases to the ambient pressure with the distance from the stagnation point. The pressure distribution on the surface of the plate also determines the local saturation conditions of the liquid along the surface which becomes important in heat transfer analysis. 10 The effect of gravity on the jet velocity at the impingement surface is significant if nozzle to surface distance is relatively large. The increase in velocity of the jet with distance from the nozzle exit causes the thickness of the jet to decrease. This is due to the continuity equation. This phenomenon is called "thinning". Often inviscid models and equations (such as Bernoulli equation) have been used to find the velocity and pressure in the impingement region. Water which often is used in cooling process is a viscous flow and the inviscid models are not suitable to simulate the jet impingement region. Therefore a need for numerical analysis of the water jet impingement is realized. 1.7 Traditional Numerical Approaches There are three main approaches to numerically model the above problem; finite element, finite difference and finite volume. Traditionally finite volume and finite difference methods have been used for a long time to model fluid flow problems whereas finite element is used to model solid domains. In the formulation of the problem, there are two main formulation methods; the Eulerian and the Lagrangian formulations. In the Lagrangian approach, the computational mesh is fixed to the material points of the deforming body and the grid point moves with the same velocity as the material point. While in Eulerian approach the mesh is fixed in space and the material is passing through the mesh. Since "finite difference" and "finite volume" methods are usually based on fixed grids, they would, generally, be categorized as Eulerian Approaches and since the mesh is fixed in this approach, the large deformation and fast movements of the fluid body would be easier to handle. 11 However, the use of Eulerian type of formulations in any of the above methods; finite volume, finite difference or finite element, entails some problems. First it is generally more difficult to account for material point quantities as we have to introduce the effect of convective terms. This is especially important if history dependent behavior is involved. Another important problem is the tracking of free surface and fluid-structure interaction. The need to keep track of the evolving boundary is particularly difficult to handle with a mesh fixed in space. This is obvious since tracking an evolving boundary means tracking material points which is a feature of the Lagrangian type of formulation. Generally, The above problems may be easier to handle in a Lagrangian type of formulation and this has been commonly used in the finite element formulation of nonlinear and large deformation problems. This benefit along with the flexibility and accuracy of the finite element method in dealing with curved boundaries makes the FEM a good candidate to deal with jet impingement problem. However due to the fast movements of the fluid, using Lagrangian FE formulation would entail the serious problems of mesh distortion and entanglement. This leads to the choice of the newly evolving "Arbitrary Lagrangian Eulerian" (ALE) formulation which may combine the advantages of both the Eulerian and the Lagrangian approaches. 1.8 Scope of Work The main objective of this work is to establish an effective and accurate simulation procedure to solve the hydrodynamics of jet impingement. In this work, we will only consider circular jet impingement on a stationary plate at room temperature. Although it is not the intention of this thesis to study the full thermal and hydrodynamic behavior of the jet, the results could be used in improving the heat transfer simulation on the plate surface in the future researches. From the 12 application side, the outcome of a detailed pressure and velocity distribution on the surface of the plate would be used to give better empirical correlations for the film coefficient on the plate surface, especially in the impingement region. In the impingement region due to high inertia of the flow the vapor layer does not exist and the flow is in direct contact with the steel plate (Figure 1.3), therefore it is hoped that the knowledge of the detailed hydrodynamics on the plate surface may also help in producing new and more accurate empirical correlations for the heat transfer analysis of the hot plate. This work would also provide the base for the further research on moving plate under impinging jets. In order to study the hydrodynamics of the jet, the following stages of research are considered: a) An Arbitrary Lagrangian Eulerian (ALE) formulation based on the fundamentals of continuum mechanics and fluid dynamics will be derived. This formulation is more general than the Eulerian and the Lagrangian approaches and could capture the benefits of both for fluid domains. b) The ALE continuum equations will be cast into finite element form. This would be obtained using Galerkin weighted-residual method. c) In order to model the stationary circular jet impingement problem, an axi-symmetric model is required. Therefore, the formulation would be described in a 2D axi-symmetric domain. d) In ALE formulation, one of the main problems is the extra set of unknowns representing grid or mesh velocities. Therefore based on the physics of the domain, a suitable approach in order to obtain extra equations for these unknowns is presented. This approach is based on solving an extra set of partial differential equations (PDEs). The solution of these PDEs is put into a finite element form through the use of Galerkin method. 13 e) Due to the continuous evolution of the free boundary in the jet impingement application, it is expected that convergence problems would be important and critical. Practical procedures would be thought to handle the convergence problems and special procedure will be established for the solution of the impingement problem. f) Due to the extra set of equations for the mesh displacement, the practical dimensions of the jet geometry and the fast velocity involved in the problem, special separation of the problem into two different regions of solution is proposed. This is thought to improve the handling of the convergence problem which is very much dependent on the mesh size and shape as well as mesh configuration with respect to the boundaries and initial conditions. g) In dealing with the free surface and in order to apply the free surface boundary condition, a Lagrange multiplier approach is proposed and used to maintain free-surface status of the fluid with the moving mesh. 14 Chapter 2 LITERATURE REVIEW 2.1 Overview of Jet Impingement Research In the last two decades, jet impingement heat transfer has received noticeable attention because of its many applications for high heat flux cooling or heating. There are numerous publications and review articles on the subject. The majority of this work is, however, experimental [see for example Jambunathan et al. (1992), Viskanta (1993), and Webb & Ma (1995) as well as references [1-26]]. Existing research in the literature may be generally categorized in three groups; Heat Transfer [1-26, 31-38], Micro-structural Evolution [64, 65], and Hydrodynamics [21-24, 27-28]. In the controlled cooling process as mentioned in chapter-1, the goal is to achieve novel properties for the steel. The mechanical properties are mainly controlled by the cooling process and are dependent on the dominant micro-structural component of the steel. The cooling rates of the steel strip on the runout table are exercised accurately according to the expected properties in terms of phase composition, size, and distribution. In Micro-structural Evolution research [64, 65] the focus is on the modeling of phase changes from austenite to ferrite. Several process parameters such as, coefficient of heat transfer, temperature at the exit of finishing stand, thickness and the speed of strip have been varied to determine their influence on the extent of phase transformation on the runout table. Another category is the study of heat transfer on the surface of the cooled strip [1-26]. The temperature of the steel at any point depends on the local heat transfer rate or heat transfer 15 coefficient. Due to the different heat transfer regimes existing simultaneously on the plate surface, heat transfer during water jet impingement on a hot steel plate is one of the most complicated industrial processes. The majority of the works done in this category is experimental and, therefore, many researchers aim at introducing empirical correlations to obtain the heat transfer coefficient on the plate surface for different heat transfer regimes. During the past three decades researchers have taken significant efforts attempting to quantify the thermal phenomena during run out table cooling operation and trying to control the cooling path to optimize the microstructure of the steel strip [29]. In order to study the problem under real life industrial. conditions, some researchers [4, 22, 35-37] focus on experiments with moving plate. The coupled heat transfer and the hydrodynamics problem is, actually, more complicated and many of the studies aiming at basic understanding still utilize stationary plate and mainly rely on experimental findings [4, 29-30]. The main focus of this experimental work has been to evaluate and modify existing heat transfer empirical correlations by analyzing the temperature measurements on the surface of the plate under impinging water flow [29-30]. The last category is the study of the hydrodynamics of the jet [21-24]. The heat transfer aspect of the jet is still the ruling element of these studies. The complications of heat transfer process has prevented most researchers in this area of research from studying the hydrodynamics of a single jet with heat transfer effects and only modest or very few publications focus on the hydrodynamic problem. While water which is a viscous fluid is being used in the industry for the cooling purpose, researchers [29, 30] in order to study the heat transfer correlations often use inviscid fluid relations in order to obtain velocity on the surface of the plate in the impingement region, or to find the pressure at the stagnation point as well as to find the jet diameter at the end 16 of the free jet. This is normally due to the lack of numerical analysis for the hydrodynamics of a single jet. The study of hydrodynamics of a viscous model is needed in order to obtain more accurate velocity and pressure profiles and to modify the empirical correlations existing for various modes of heat transfer on the existing on the plate surface and including convection, nucleate, transition boiling, and film boiling regimes. These correlations often contain hydro dynamical parameters such as velocity and pressure especially for correlations in the impingement region. As discussed in Chapter-1, the focus of this thesis is to study the hydrodynamics of a single jet impingement. Due to the fact that there are very few works done in the this area and the fact that these works are, in many cases, inseparable from the heat transfer researches, the review below will focus on both the heat transfer and the hydrodynamics aspects of the problem. 2.1.1 Heat Transfer and Hydrodynamics of the Jet Heat transfer analysis of impinging jets on hot plate has been the subject of many researches [1- 26]. In order to quantify how much heat is extracted from a hot plate, a local heat transfer coefficient has to be obtained. There are number of parameters that affect the heat transfer rate in a jet impingement configuration. These may include jet to target distance, jet exit velocity, jet confinement, jet angle of attack, water saturation temperature, plate geometry, velocity, surface condition and surface temperature [1-20, 29]. The flow and heat transfer conditions involve turbulent, free boundary flow, and boiling heat transfer conditions. 17 With respect to the heat transfer part, the problem is extensively handled both numerically and experimentally. Studying the hydrodynamic aspect, various numerical approaches exist including finite difference, finite volume, Eulerian finite element, and Eulerian-Lagrangian finite elements. From the experimental side, experimental testing included study of various parameters, verification of different turbulence models, study of various heat transfer simulations for boiling heat transfer as well as distribution of heat transfer coefficient. The Center for Turbulence Research (CTR) [7-16] has done extensive research on various turbulence models. Their focus has been on gas flow, which is not normally used in the steel industry cooling systems. Turbulent impinging jets have complex features due to entrainment, stagnation, and high streamline curvature. These features are proved to be incompatible with most of the popular turbulence models, which are essentially developed and tested for flows parallel to a wall. One of the most important problems of these models is a substantial over-prediction of heat transfer in the stagnation region. Some researchers tried to improve the widely used k-c model for the jet impingement problem [25]. Knowles [25] introduced corrections to the k-c turbulence model and used finite volume code to perform numerical modeling of an axi-symmetric jet impingement. He investigated the effect of pressure ratio, height and nozzle exit velocity profile and concluded that there are still significant errors in the predicted wall jet growth. He also performed numerical experiments and showed that there is great effect from impingement on the wall jet [25]. In 1991 [11] Durbin suggested a new turbulence model called `v2f for the jet impingement. The v2f model makes use of the standard k—c model, but extends it by incorporating near-wall 18 turbulence anisotropy and non-local pressure-strain effects, while retaining a linear eddy viscosity assumption. In the work done in the CTR [7-10, 12-13] data obtained in a fully developed impinging jet was used to show that the v2f model successfully and economically predicts the rate of heat transfer. The reports indicated the effects of parameters such as jet-to- target distance, geometry, and Reynolds number, as well as jet confinement. Behnia et al. [10, 13] used an elliptic relaxation turbulence model (v2-f ) that has been suggested by Durbin [11,14] to simulate the flow and heat transfer in circular confined and unconfined impinging jet. The model was validated by the available experimental data sets. They showed that flow characteristics in the nozzle strongly affect the heat transfer rates, especially in the stagnation region and that different nozzle velocity profiles can have different heat transfer rates. They studied the effect of several parameters such as nozzle to plate distance and Reynolds number on the heat transfer rates. In some previous studies these effects have been addressed, and results of the experiments performed by different investigators sometimes are contradictory because of the differences in the experimental conditions [I]. Jambunathan et al. (1992) pointed out this problem and noted that for a better understanding of the jet impingement heat transfer process, flow details, geometry and turbulence conditions are required so that a comparison between different experimental data can be made. Behnia et al. compared v2f turbulent model with the k-E model and showed that k-E model over predicts the rate of heat transfer whereas the v2f model provides more accurate predictions. In conclusion most predictions of jet impingement heat transfer in industry involve the use of standard or modified version of k-E turbulence model, which is widely available in existing CFD 19 packages. It seems that with all improvements to the well known turbulence models, they are still not suitable for jet impingement simulation, and that the v2f model which is specifically proposed for jet impingement provides a promising alternative. A considerable amount of work has been done on the heat transfer of the jet impingement during last few decades [1-20, 31-32]. Numerical models have been developed to model the heat transfer in boiling and during jet impingement [33-37]. The modes of heat transfer depend on the velocity, temperatures of the liquid and impinged surface as well as the physical properties of both liquid and solid. If the initial wall super heat is high enough, one should expect conditions of film boiling, transition boiling, nucleate boiling and finally convective heat transfer. One of the important issues in this area is investigation of heat flux and its maximum value on the plate surface. The critical heat flux at the stagnation point has been proven to be dependent on the velocity of the jet and the velocity distribution in the impingement zone. This indicates that accurate simulation of the flow is necessary to analyze the heat transfer [29, 30]. On the numerical side, the finite difference method is probably the most widely used method of analysis [7-10, 12, 13, 38]. For the spatial discretization of convective terms, upwind scheme is normally used and a central difference scheme is used for the diffusion terms. In the CTR work [10], finite difference with fine, non-uniform, orthogonal, cylindrical grids were used with a high resolution near all solid boundaries. Mesh sensitivity was carried out by varying the density of the mesh until the changes in the Nussult number in the impingement zone was less than one percent and, therefore, the solution was assumed to be grid-independent. 20 The research in CTR is basically focused on presenting better turbulent model for gas impinging flow. Free surface feature and hydro dynamical data have not been exclusively studied in their work as it is not in direct connection to the application of jet impingement in the steel industry cooling systems. Bierbrauer et al. [21, 22] used Eulerian — Lagrangian Method (ELM) to investigate the heat transfer on a hot moving steel strip impinged by a water jet. Their ELM approach effectively decouples the advective terms by solving a purely Lagrangian step followed by a purely Eulerian one. This is shown to improve the stability of traditional Lagrangian and Eulerian approaches. Their study show that under the impingement of the water jet the strip temperature at the stagnation point is initially depressed and subsequently rebounds towards a steady state. The heat flux increases with increasing convective acceleration of the jet. Bierbrauer et al. defined three thermal regions; free jet, stagnation, and wall-jet regions. In the free jet region the velocity and temperature distribution of the jet is not affected by the presence of the surface. In the stagnation region, single-phase forced convection takes place over several jet widths and cooling effectiveness is high. The wall jet region is separated into a small region of transition and nucleate boiling before entering the parallel film boiling region separating the surface from the jet by a vapor layer which reduces cooling effectiveness considerably. Phares et al. [28] used an analytical method to find solution for arbitrary velocity profile of the jet in the axi-symmetric and 2-D impingement flow. They obtained their analytical solution for inviscid flow outside the boundary layer to provide the boundary condition for the boundary layer near the stagnation point. 21 As mentioned above, most of the work done in the area of heat transfer is experimental. While this is essential in providing physical insight and understanding of various phenomena occurring during jet impingement, it entails many drawbacks. This includes the cost of experimentation, the long time needed for experiment preparation and setup, problems and efficiency of surface temperature measurements, errors in data measurements and acquisition, and the wide range of parameters that are required to be investigated. These problems are magnified if the plate is moving under the jet which is the actual case in practice. Because of these difficulties, many simulations exist for various regions of heat transfer on the plate surface and each is restricted to a limited range of parameters and conditions. In most of these simulations, the velocity profile and the hydrodynamics of the flow on the surface of the plate are required. Most existing work use approximate and/or empirical methods to provide the necessary information for the flow conditions. The above arguments indicate the need for an accurate independent numerical simulation to the flow conditions of jet impingement on stationary and moving plates. The numerical simulation of the problem will enable an efficient and inexpensive parametric study of the various parameters affecting the flow and pressure conditions. The aim of this thesis is to use a more accurate numerical method with a suitable mesh motion scheme to solve for the hydrodynamics of a single stationary jet impingement problem. The dimensions and specifications used in this work are taken from the actual experimental setup built and used in the run out table group in UBC for the multi—disciplinary project "Controlled Cooling of Steel Strip on the Run-out Table for Novel Properties ". 22 In the following section we will briefly review the numerical methods and techniques that may be used in the project. 2.2 Numerical Techniques and Formulations 2.2.1 Numerical Techniques The main numerical techniques used to solve a continuum problem of the above nature are the finite difference, the finite volume and the finite element. The Finite element method has been traditionally used in solids while finite difference and finite volume method in fluids. In both finite difference and finite volume techniques, the grid or computational points are fixed in space and the material is passing through these grid points. In other words, we are basically dealing with an Eulerian description of the domain. In finite element method, however, it is possible to use any of the Eulerian, Lagrangian or ALE (Arbitrary Lagrangian Eulerian) description. The ALE description has the advantage of the possibility to convert to pure Eulerian or pure Lagrangian one. Finite Element Method has been proven to be more capable of solving complicated boundary conditions. Free boundary, unknown interface between vapor and water on the surface and interaction of solid and fluid are some of the features of the jet impingement on the hot steel strip. Finite element method and specifically ALE approach seems to be more appropriate to deal with these features. 23 Due to its capabilities the Finite element method is continually gaining popularity in solving fluid problems. Although the Free surface and solid-fluid interaction are some of the main features of the jet impingement problem, and these are the strong reasons to use ALE method, there is basically very little work done in the hydrodynamics of jet impingement using this method. This is probably due to the fact that the main focus of the majority of numerical researches has been the heat transfer aspect of jet impingement rather than its Hydrodynamics In order to deal with the unknown position of free surface two different approaches have been used in the literature. For a simple case of a single-valued function, in which the position of the free surface could be described by one variable, a hyperbolic equation as the kinematic equation of the surface has been used. This approach has been used by Ramswamy and Kawahara (1987), Huerta and Liu (1988) and Souli and Zolesio (2001). A second approach can be obtained by imposing the condition that no particle can cross the free surface (because it is a material surface). This can be imposed by using a Lagrangian description in which the velocity of the computational mesh is equal to material velocity along the boundary of the free surface. This approach has been used by Huerta and Liu (1989) and Braess and Wriggers (2000). However using the Lagrangian boundary condition could create severe mesh distortion. Therefore this condition can be relaxed by imposing the constraint only in the normal direction. This will allow the mesh points to move tangential to the free surface. In the work of this thesis, the relaxed condition has been used in conjunction with the Lagrange Multiplier method. 24 2.2.2 Traditional Numerical Formulations Traditionally there are two main formulation approaches to deal with a continuum; the Eulerian and the Lagrangian approaches. The Eulerian approach is similar to studying the material particles from a constant position while the Lagrangian approach, is more like sitting on the particles and follow them while moving. In the Lagrangian approach, the computational mesh is fixed to the material points of the deforming body and the grid point moves with the same velocity as the material point. In the case of problems with large deformation and fluid problems, this leads to excessive distortion of the mesh and presents a major disadvantage of the formulation. Distorted mesh affects the solution accuracy and may ultimately result in a singular stiffness matrix. Another drawback of the Lagrangian approach is the difficulty to model non-material-associated boundary conditions. The only remedy for the mesh distortion problem in Lagrangian formulation is to interactively re-mesh the material domain whenever needed. Re-meshing, besides being very much time consuming, requires user intervention. The simulation process may need frequent manual re- meshing to reach the required deformation or flow level. As well as the fact that due to re- meshing of the domain the material properties and variables must be transferred from previous mesh to the new generated mesh. The analyst must either run the code interactively or have enough knowledge of the solution to write a detailed set of re-meshing instructions ahead of time. 25 In the Eulerian approach the computational mesh is fixed in the space, while the material is passing through it. When the material is fast moving, it is more suited to use this approach. As this approach is more suited for fluids, it is generally the base of the Finite Difference and Finite Volume methods. In specific fluid problems such as free surface and fluid-solid interactions, this approach has its shortcomings and could not be used in its original form. 2.2.3 The Arbitrary Lagrangian Eulerian (ALE) Formulation The above discussion indicates the need for a formulation that is more suited for the simulation of large deformation problems with boundary motions. The Arbitrary Lagrangian Eulerian (ALE) formulation has emerged as a technique that can alleviate many of the drawbacks of the traditional Lagrangian and Eulerian formulations. Hughes (1981) first proposed ALE in fluid mechanics Some other studies in this field are by Belytschko (1978), Liu (1981), Liu and Ma (1982), Belytschko and Liu(1985), Liu et al. (1991), Hu and Liu (1993) and Braess et al. ( 2000). Also, ALE has been recently introduced in solid media [40-47, 54]. The large distortion/deformation that characterizes non-linear solid problems clearly undermines the utility of the Lagrangian approach traditionally used in problems involving materials with path-dependent constitutive relations [54, 55]. Liu et.al. presented an ALE formulation for path-dependent materials using explicit time integration. They also have presented a new method for the treatment of convective terms [73]. Gadala et.al . have used a fully coupled implicit time stepping method, in order to derive an ALE formulation for the dynamic problems. They have used the derived formulation in large deformation problems such as metal cutting and crack propagation [44-47, 59-60]. 26 In the ALE formulation, the finite element mesh, or reference configuration, need not adhere to the material nor be fixed in space but can move arbitrarily. As the material deforms, the finite element mesh is continuously moved to meet any preset criterion (e.g. optimize elements shape) and the simulation should be completed without user intervention. Combining the merits of both the Lagrangian and Eulerian formulations, ALE can easily describe all types of boundary conditions and prevent mesh distortion. Certain types of boundary conditions such as sudden changes in material point may be easily dealt with in ALE method and they are more difficult to model with Lagrangian or Eulerian method. Fluid free surface may be tracked in Lagrangian formulation but serious mesh problems would be evident. Also, due to velocity of the fluid flow and the geometry of real life impingement problem, studying the whole body of fluid using Lagrangian formulation is not possible. Different approaches have been adopted to develop the ALE formulation. The differences among the available ALE formulations depend on the intended application, assumptions made in deriving the ALE equations and details of implementation. These factors determine the limitations and accuracy of the resulting formulation. ALE is usually termed as coupled formulation since material deformation and convective effects are coupled in the same equations. However, two distinct techniques are being used to implement the ALE equations. The first technique is referred to as an 'operator split' or a `fractional step' approach. Most of existing ALE analyses, especially in solid mechanics, are based on this strategy [41 ]. In this approach, material deformation and convective effects are 27 treated separately. Thus each time step may be divided into two steps: a regular Lagrangian step followed by an Eulerian step. In the Lagrangian step, the grid moves with the material, whereas in the Eulerian step, the Lagrangian solution is mapped to the reference grid and stresses are updated using convective effects. Time advances only during the Lagrangian step and there is no time associated with the Eulerian step. Thus it is not necessary to perform an Eulerian step in every time step. The alternative ALE approach, which has been used by fewer researchers, is known as the 'fully coupled' or the 'un-split' approach. In this approach, the governing ALE equations are implemented and solved without disruption. A solution algorithm that can handle convective effects simultaneously with material deformation must be used. The main advantage of the operator split over the fully coupled approach is the ease of implementation of ALE into current Lagrangian codes as the Lagrangian step is unchanged and only the Eulerian step algorithm needs to be added. Moreover, the decoupling of the Lagrangian and Eulerian steps results in simpler equations to be solved and simpler algorithms to be used. However, from the theoretical point of view, the fully coupled ALE approach represents a true kinematic description that employs a more rigorous scheme in considering equilibrium in each step relative to a moving reference configuration. Therefore, the operator split solution, which stems from computational convenience, is expected to be less accurate than the fully coupled solution. Using the ALE formulation, the finite element mesh can be moved arbitrarily to maintain a homogeneous mesh and properly represent boundary conditions throughout the deformation process. A mesh motion scheme is necessary to continuously adapt the positions of internal 28 nodes. On the boundaries, however, mesh motion must satisfy the boundary constraint that prevents the relative motion between the material points and mesh points in the direction normal to the boundary. This constraint ensures that the material and mesh configurations have the same boundary at all times. The main drawback of the ALE technique lies in the fact that the number of unknowns is increased, having both the material and mesh velocities (or displacements) and in the fact that a separate mesh motion scheme is required. Fully coupled ALE computations are expected to be more time consuming than the corresponding Lagrangian ones. The other critical point in ALE formulation is the proper choice and adaptation of a suitable mesh motion for the problem. Many of the new pieces of research focus on addressing these points but until now, the existing work pertains to specific application and lacks the generality of application to general large deformation and moving boundary problems. 29 Chapter 3 Finite Element ALE Equations for Jet Impingement In this chapter first the flow equations are derived based on the fundamental principles of continuum mechanics using conservation laws of fluid, then the ALE fundamentals are introduced and applied to the flow equations and finally these equations are derived in axi- symmetric form in order to be used for jet impingement application in chapter 5. 3.1 Flow Equations The following section discusses various constitutive models suitable for the jet impingement problem and presents stress and strain measures used in these models. The notation which is used is indicial notation as described by Malvern [39]. Vectors and tensors are shown in boldface and scalar variables are shown in regular letters. Tensors could also be shown in matrix form and by their scalar components. In that case the first index in the subscript would represent the row and the second index would represent the column in the matrix. A comma followed by a subscript (,i) represents partial derivative with respect to id' direction. A dot sign on top of a variable represents the absolute time derivative of that variable. A hat sign on top of a variable shows the weighting function based on the shape functions of that variable. Derivation of the flow equations is based on fundamentals of continuum mechanics and fluid conservation laws [39, 51-52]. 3.1.1 Constitutive Model Experiments show that a fluid at rest or in uniform flow can not maintain shear stress. 30 Therefore in a fluid at rest or in uniform flow the stress is purely hydrostatic and the stress may be given by: (3.1) Where p is the hydrostatic pressure, l is stress tensor and Sy is Kronecker Delta with the following definition: {0 i j = 1^i = j (3.2) In general, the thermodynamic pressure is a function of position and, therefore, it is one of the unknowns of the fluid problem. An ideal frictionless (non viscous) fluid is the one that can sustain no shear stress even in motion. This assumption may be justified if the pressure and body forces are much larger than the viscous effects. In such case, the fluid may be assumed ideal and we can assume fluid to be non-viscous except in the boundary layer (thin layer near a solid boundary). The boundary condition that should be imposed for the flow right at the interface with a solid boundary is vn= 0 (normal velocity at the interface equal to zero). The constitutive equation for an ideal non viscous fluid is given by Equation (3.1). In viscous fluid however shear stresses are developed. Stokes (1845) assumed that the difference between the stress in a deforming fluid and the static equilibrium stress is a function of rate of deformation tensor (D), so the total stress may be given by: T = -pl +F(D) (3.3) In which 1 is unit or identity tensor and D is rate of deformation tensor given by: 31 D = 2 -1 (Vv + (Vv) T ) = (vy + vfi ) ^ (3.4) Although Stokes only considered the case of linear viscosity, a nonlinear relation of F with respect to the deformation gradient D in equation (3.3) would be called stokesian or non- Newtonian fluid. When F is linear, the fluid is normally called Newtonian fluid or linear viscous fluid. For linear viscous fluid we may write the general form of the constitutive equation as: Tii = — p + C yid D ^ (3.5) In which Cijkl are constants representing the components of the 4 th order tensor C. For isotropic fluids, C should be an isotropic fourth order tensor, and since T and D are symmetric tensors, C would have the following general form which is completely defined by two independent constants: Cijk/ = ^'8 (aik gfi^5jk ^(3.6) Substituting equation (3.6) into equation (3.5), we get: ^Tij = -Pgi; ±ilDkkgu + 214 ^ (3.7) In which p and A are two independent parameters characterizing the viscosity of the fluid. If we adopt deviatoric stress tensor measures and rate of deformation tensor ( , DO, the Newtonian fluid law will take the following special form: (1-3 P)8./ ± (2 + —32 P)Dkkg y 2 ,1-11) ^ (3.8) In which ji is the mean pressure. Noting that T1,̂ 0 and D;1 0 , and letting i = j in equation (3.8), we get: 32 (i) 15')+(2± 2 /40k =03^k This relation may be used to reduce equation (3.8) to the following relations: (3.9) (3.10)T = 2 iuD 1-70 = p^k _ p+^dp p dt In which K is also known as bulk viscosity, K = A + —2 p 3 and due to continuity (as will be discussed below) in conservation laws: (3.11) (3.12) D kk = div (v) — 1 d p p dt (3.13) Assumption of stokes condition states: K^+ 0 3 (3.14) Which leads to the following equation: Tij =^ (3.15) 3.1.2 Conservation Laws The conservation laws adopted for fluid media are presented in the following: Conservation of linear momentum: pb, = pvl^(3.16) In which p is the density and pb, represents body force component. 33 or dt d p + pV • v = 0 (3.18) The next law is conservation of angular momentum which causes the stress tensor to be symmetric and the last law is conservation of mass: dp^av i+ p^= 0 dt^ax (3.17) 3.1.3 Field Equations Eliminating the stress T from equations (3.15-3.18), we obtain the Navier-Stokes equations for incompressible fluid in the form: - 1—Vp + vv 2 v+ b dt)9^dt (3.19) In which v is the kinematic viscosity () . For an incompressible fluid, the conservation of mass or the continuity equation reduces to the following simple form: V • v = 0^ (3.20) The vector equation (3.19) along with the scalar continuity equation (3.20) gives rise to four scalar equations in four unknowns; three velocity components and the pressure at a point. 34 3.2 Boundary Conditions In viscous flow any solid object applies shear stress on the flow and this causes the flow immediately in contact with the solid to become stationary. This condition is called "no slip" condition and basically is defined as v = 0 for the flow at any rigid wall. In the jet impingement case on a stationary plate, on the boundary of the water domain with the solid plate this condition applies. The pressure in Equation (3.21) is a function of position and has to be known at least at one point as the boundary condition. Generally, the solution obtained for pressure would be relative pressure with respect to the known boundary value. In case of jet impingement problem, this condition could be applied by setting the outflow pressure to zero, both in free jet region and in the impingement region, as will be shown in Chapter 5. Free surface is another boundary type that exists in the jet impingement problem (both in free jet and impingement region) and we may apply two types of boundary conditions on such surface; kinematical and dynamic. The dynamic condition is explained by stress-free situation, (n • T = 0) in which T is stress tensor as described by Eqn. (3.15) and n is the outward unit normal to the fluid domain. This condition after replacing the deviatoric tensor in Eqn. (3.15), could be shown to reduce to [—p1+ p * (Vv + (Vv) T )1 • n =0 . This condition is also known as the neutral boundary condition. The other type of boundary condition for the free surface is kinematical boundary condition. The base of this condition is the fact that the material can not pass the mesh boundary. This condition which is Lagrangian in nature is expressed by equating the material velocity to mesh velocity. The condition couples the mesh motion equation (to be presented in chapter-4) with the fluid FE equations presented in this chapter. Generally, this condition may cause undesirable distortions 35 in the mesh and it is better to relax the condition and only use it in the normal direction. In this case the mesh velocity is equal to the material velocity only in the direction normal to the free surface. Since this condition will create a kind of constraint on the flow equation on the free boundary which may still violate the free-stress condition, it is best to apply the condition using Lagrange Multiplier on the free surface. This approach will be explained in Chapter-4. 3.3 ALE (Arbitrary Lagrangian Eulerian) Kinematics ALE description as mentioned in the first chapter is a generalization of the Lagrangian and Eulerian descriptions. In the general description of the ALE formulation [see for example: 48- 50, 53], fundamentally three configurations are prescribed: Material Configuration, Spatial Configuration and Referential Configuration. Mathematically, the referential configuration is related to spatial and material configurations through two different mappings. In order to derive the fundamentals of ALE, first the relation between the spatial and material configurations is derived. These configurations are often used in Lagrangian and Eulerian Descriptions. Next the Referential configuration is explained and finally the ALE fundamental relations are derived. 3.3.1 Material and Spatial Configurations Two domains are commonly used in continuum mechanics: the material domain made up of material particles defined by X and the spatial domain consisting of spatial points defined by x. The Lagrangian viewpoint consists of following the material particles of the continuum in their motion. In order to perform the computation for the variables of a domain we need to define a computational domain introduced by grid and grid nodes. 36 In the Lagrangian description grid nodes are permanently connected to the same material points, so the reference configuration (identified by grid) will be defined by material coordinates. In that case the motion of the material points relates the material coordinates, to the spatial ones x. Therefore, at a given time instant, the spatial coordinate could be described as function of material coordinate and time: x = x(X,t)^ (3.21) In a two-dimensional case, the Jacobian mapping the two coordinates for the above transformation is defined by: J = ax,^ax, ax, ax2 ax2 ax2 ax, ax2 _ (3.23) The mapping should be objective and one-to-one, which means x -1 = X or in other words: for any point in the material domain there is one and only one point in the spatial domain at each time. In the Lagrangian description since the computational grid (reference grid) matches the material points during the whole motion, there are no convective effect in the calculation. However, when large material deformations do occur, for instance vortices in fluids, Lagrangian algorithms undergo a significant loss in accuracy and may even be unable to conclude a calculation due to excessive distortions of computational mesh linked to the material. These difficulties caused by excessive distortion of the finite element grid are alleviated by the use of the Eulerian formulation. The basic idea of the Eulerian formulation which is very popular in fluid mechanics, lies in examining, as time evolves, the physical quantities associated with the 37 fluid particles passing through a fixed region of space. In an Eulerian description, the finite element mesh is fixed and the continuum moves and deforms with respect to computational grid. The material velocity v at a given mesh node corresponds to the velocity of the material point coincident at the considered time t with the considered node. The velocity v is consequently expressed with respect to fixed-element mesh without any "reference" to the initial configuration of the continuum and the material coordinate X. In other words, the reference coordinates in the Eulerian formulation is assumed to be the same as spatial coordinates. Since the Eulerian formulation dissociates the mesh nodes from the material particles, convective effects appear because of relative motion between the deforming material and the computational grid. Eulerian algorithms present numerical difficulties due to non-symmetric character of the convection operators, but permit an easy treatment of complex material motion. By contrast to the Lagrangian description, serious difficulties are found in following the deforming material interfaces and mobile boundaries. 3.3.2 Referential configuration The above discussion of the classical Lagrangian and Eulerian descriptions shows the advantages and drawbacks of each approach. It also shows the interest in a generalized description capable of combining the interesting aspects of each approach while minimizing their drawbacks. This generalized approach is called ALE (Arbitrary Lagrangian Eulerian) method. In this approach the reference coordinates is neither attached to the material points nor to the spatial coordinates. Therefore we have a reference grid or "computational mesh" independent of material points. In this generalized formulation if we set the mesh movement to follow material point, the result would be Lagrangian and if we set it to have no movement the result would be 38 Eulerian. It is also possible that the "reference mesh" has its own independent velocity. The approach is an attempt to avoid drawbacks of each of the Eulerian and the Lagrangian approaches. In ALE, we define a Referential Coordinates x (independent of material and spatial Coordinates X, x , respectively) which is introduced to identify the grid points. The mesh or computational domain described by x , moves with grid velocity vg , relative to spatial configuration x. In the same way that we derived the Jacobian for the mapping between material and spatial coordinates, we can derive the jacobian for mapping to the referential coordinates. The jacobian for the referential domain in 2D could be obtained as: s = ax,^ax, aX1 axe axe ax2 ax, axe (3.23) The reference system locates the nodes and elements of the mesh. Therefore yi shows the velocity of the mesh. The arbitrariness of the mesh movement could be described in these three points: • If 2, = X then vg = v which means the velocity of the mesh is the same as velocity vector of the material, and the approach reduces to a Lagrangian formulation. • If x = x then vg = 0 which means that the velocity of the mesh becomes zero or the mesh is fixed in space and the approach reduces to an Eulerian formulation. • If vg {0, v} then the mesh is neither attached to the material nor fixed to the space. In that case mesh-motion has no physical interpretation and special scheme to identify the mesh motion is needed. 39 3.3.3 Time Derivatives in Different Configurations and Their Relations Since our calculations are all in reference coordinates or computational domain (mesh), we need to find the relation between the time derivatives of the variables in spatial and material coordinates (available in the field equations) to the time derivative in the referential system. We can define these derivatives in different systems as follow: Material time derivative: m off; = ,^X : constantat Spatial time derivative: (3.24) fiS aft^x : constant^ (3.25) Referential time derivative: f R at' ^: constant^ (3.26) Material time derivative of the spatial coordinate system is the material velocity (real velocity) while the Referential time derivative of the spatial coordinate system is mesh velocity. So we can write: and axv.—,at axvg = —at , X : constant^ (3.27) 2, : constant^ (3.28) If we define c as convective velocity relative to the mesh: 40 c = v — vg^ (3.29) By definition, if c becomes zero we would have a Lagrangian approach and if it becomes equal to v, then we would have an Eulerian approach. If the mesh velocity is not zero and it is independent of material velocity then c is an independent vector. Now if we assume a small time step, and assume that mesh-velocity doesn't vary within this time step. We may state the following: x(x, t) = x(x, t) — vg (x)t^(3.30) From (3.30) we can derive the Jacobian of the referential domain: axi^ay.g -5=^= (Sr + '̂ t ax;^ax; (3.31) Material time derivative of referential domain x based on (3.30) is: ax,^ ax; ^ t= x"at^ems^at^at Consequently we get: fax. at ax j • /X vonst = vi vi gaz; (3.32) (3.33) and with the definition of the velocity difference c and 3 (convective velocity and Jacobian of the referential system) we conclude: aX = Cat (3.34) From the proceeding equations we can now establish the material time derivative of any variable f and relate it to the referential system: 41 (Df Dt ) _(ai)^± (af axi)^ at ) zconst^X:const Therefore: (ptl _ rai at ) l^ax+ Vf • 3 — ^;e:const^at Consequently using (3.34), we obtain: (3.35) (3.36) Df -= (la^+ (Vf). c^(3.37) Dt at ) xxonst The above equation (3.37) is the ALE expression for a material time derivative for a scalar quantity. Since ALE is the generalized form of Lagrangian and Eulerian formulation we have to be able to obtain Eulerian and Lagrangian formulation from equation (3.37). If we set c = v we will find the Eulerian relation of time derivatives and if we set c = 0 , it is clear that the material time derivative is the same as referential one, which is exactly the case in Lagrangian approach in which the material and referential systems match. As a concluding remark for the above section, we state the following: If we pay attention to equation (3.37), we will find out that the only change from spatial domain in the referential domain (ALE description) is replacing v in convective terms with c = v — vg (convective velocity). This means that we should be able to change any equation in to referential domain by just replacing v , in convective terms , by c= v—vg or the convective velocity. We will use this important conclusion in the derivation of ALE form of flow equation. 42 3.4 Derivation of ALE Finite Element Equations in the Flow Field In this section first the suitable form of Navier-Stokes equations, based on the total stress is derived. Then this form is transformed in to Referential domain based on the discussion of Equation (3.37), then these equations are transformed into FE equations using Galerkin method by integration over the fixed domain and finally these equations are presented for an axi- symmetric domain to be applied to jet impingement application [53-58]. 3.4.1 Transformation to Referential Domain The total stress in a flow consists of pressure and viscous parts. This could be shown by substituting D' in (3.15) then derive the total stress tensor: T. [–pi tu*(Vv + (Vv) T )]^ (3.38) By rephrasing the flow equation (3.19) based on the total stress tensor, we obtain: p — dv — V.[–pI + ,u(Vv + (Vv) T )] = F dt (3.39) In which F is body force (body force). This form of the flow equation is more desirable, since in the finite element weak form the total force should be zero on the free surface and we would generally try to minimize it. dvThe absolute acceleration — in the first term of (3.39) can be decomposed into a local and dt convective acceleration. dv av + v • (Vv) dt at (3.40) 43 It will be shown below that one source of the non-linearity in the Navier-Stokes equations comes from the second term in the above equation which represents the convective acceleration. If we input the above equation into the flow equation, we may have: p av + pv.(Vv)— V.[---pI + ,u(Vv + (Vv) T )] = F at (3.41) Using the conclusion of Section (3.3.3) we transfer the above equation into the referential system with the coordinate system velocity of vg . Since the only convective term in (3.41) is the second term on the left hand side of the equation, we can replace the velocity v with the convective velocity c=v–vg. In this case we get: p— Ot + p(v – vg).(Vv) – V .[–pl + ,u(Vv + (Vv)T )1 F (3.42) Equation (3.42) is the ALE form of the flow equation in Referential domain. In the section to follow, this equation will be transformed into a Finite Element form. 3.4.2 Conversion to FE Equations In order to convert Equation (3.42) into a finite element form, we use the Galerkin Residual Method. In this method we use weighting functions which are the same as the shape functions. The weighting functions are multiplied by the residuals (in the domain of analysis and on the boundary) and then integrated to get the total residual for the system. Normally we use the "integration by parts" or the "Gauss divergence theorem" in order to decrease the order of the derivatives in the "total residual form", and the finite element equations are obtained after identifying the specific shape functions or weighting functions: 44 f[(2 p u ,x — p)i ,x + p(u ,y + v ,x )fi y + p((u — u g )u x + (v — v g )u ,y )fi + ,u(vx + u y )i) x + (2 itv y — p)i + p((u — u g )1, + (v — vg )v + (u x + v y ) + au „^avpu + at In f^ L mov = 0at (3.43) which (u, v) are the components of v and (ug, vg ) are the components of mesh velocity and (fx , fy ) are components of body force in the 2D domain. In (3.43), terms (2 ,uu — p)fi u(u y + v x )14 y u(v + t y )i) , (2 ,uv y — p)i) y are the result of integration by parts which results in decreasing the derivative order of the weighting function. This also makes the final finite equations to be symmetric. The term (u x + v y ):/3 arises from the continuity equation (3.20) and the shape function associated with this term is a pressure shape function. The neutrality of the free surface is naturally taken care of in this weak form. The boundary terms are automatically eliminated as this is the case in traditional Galerkin method. Since this integration is performed throughout the whole domain and in ALE the domain is changing with time, we need special treatment to deal with this integration over a changing domain A mapping from the original domain C2 0, with original coordinates (X, Y) to a deformed domain 0.„„ with moving coordinates (x, y) is defined: X,YEC2 fix^x(X ,Y ,t), y(X ,Y ,t) e^(3.44) To obtain the spatial relationship between the original and deformed coordinates incremental values, we can use the chain rule of differentiation in the form: 45 adYxdx --dXax + — ax^aY (3.45) and dy dX k-dY+ (3.46)= ax^ay In the matrix form we have: [ dX X X^YX drill CIX1 LdY] [Y,X ,Y1dY Ji_dyi (3.47) In which J is the Jacobian of the transformation. To formulate the inverse transformation, we use the inverse of the Jacobian matrix in the form: In which rdxi^rdxi Ldyi -Ji Ldyi 1[ Y,y^- x0,1 Ix, IXy l-1 D - Y ,x x X^I IYx^Yy (3.48) (3.49) Where D is the determinant of the Jacobian defined given by: D = x,x y,y - xyy,x^(3.50) In order to replace the incremental volume in Equation (3.43) in terms of fixed original incremental volume, we use: dC2,nc,„ = Ddf2 fix^(3.51) The shape function derivatives could now be given by: fi ,X1Xx +11 ,YIYx Zt =fi I +^Iy^,Y Yy^X Xy V =V I +X Xx ,Y I Yx (3.52) 46 ) I + ) Iy^,Y Yy^X Xy By substituting (3.52) and (3.51) in equation (3.43) we obtain: xIx + u y1 yx )(2 pu x – p) + (1; I)5, + ,x xy )(2^– p) -1) + (II Ylyy + 14 x' xy +1) ,x1xx +i), y yx )p(11,y -FV, x ) + p((u ug)v x + (v – vg )1? y )i) + p((u – ug)u + (v – vg )u y )fi + (u,x +vo,)13+ au^avpu + p -1) – fit – fy i))DdS2 fix = 0at^at Equation (3.53) is the transformed form of the equation (3.43) from the moving domain into the fixed domain. 3.4.3 ALE Finite Element Equations of the Flow in Axi-Symmetric Form In order to obtain the finite element formulation for axi-Symmetric domain, we integrate the axi- symmetric form of the ALE Navier-Stokes equations in moving domain using the Galerkin method. Applying the normal procedure of using the "integration by parts" or "Gauss theorem" we obtain: f[(2,tru„–rp)14„+ ,tr(u z +v)i + (2,uu I r– p)14+ cr((u–ug )14, + (v–vg )u + kr(v„+u)v„.+(2ritivz –rp)i lz + p((u–ug)v,+(v–vg)v,z )fi+ ^ [r(u„+v z )+urf3+^ (3.54) „rp—au u + rp av—v –rfru –rfzvjannio,= 0at^at In which (u, v) are components of v and (ug,vg) are components of mesh velocity and (fr , fz ) are components of body force in the axi-symmetric domain. (3.53) 47 Equation (3.54) is integrated in the moving domain. In order to transfer the equation to the original (fixed) domain we can use the same measures used for the 2D equations. A mapping is defined by: R,Z E C2 fix^r(R,Z,t),z(R,Z ,t) E C2,„0,^(3.55) To obtain the spatial relationship between the original and deformed incremental values, we use: dr = ar dR+— ar dZaR dz = — az dR+— az dZaR^az In the matrix format we can say: [dri = [rR rZ ][dRi = j[dR1 dz z R z Z dZ^dZ (3.56) (3.57) (3.58) In which J is called Jacobian of the transformation. To formulate the inverse transformation, we can use the Jacobian inverse matrix J -1 and we obtain: In which dR1=^dr] LdZ^Ldz (3.59) i[z z=- ' D r,z1 [IRr R^IZr IZz ^ (3.60) In which D is determinant of the Jacobian given by: D = r:R z ,z rz z^(3.61) And the shape function derivatives are give by: 48 =u ,RI +12,,IzrRr 14,z =^ (3.62) 1)„ = RIRr +1),zIzr 1;,z = 1;,R/Rz + 1),zizz After replacing the shape function derivatives using Equation (3.62) and the incremental volume in terms of fixed domain we will finally obtain: ,RI Rr^,z1 zr)(2Pru -rp)+ ,ur(14 ,zI zz +it R )(ti,z V,r (2/1U r - p)14 + pr((u - ug )u , r +(v - vg)u z )fi + 1,(V,RIRr^zr)(1)r + u z )+ (1) RI Rz +i) z1zz )(2r,uv,z - rp) + p((u -ug)v, + (v -vg)v z )i) + [r(ur +v z )+u]fr + av^„rp—fi + rp—i) - rfrfi - rfzviDc12 fix =0at^at (3.63) After factorization we will get: Rr + u,zI zr )(2,uru, - rp) +RIRz +i) , z1 zz )(2r,uv ,z - rp) + (2, I r - p)14 pr(fi zIzz + 1,1RIRz +V,RIRr +^zr )(u z + v r ) + pr((u -ug)u „ + (v -vg)u z )ii + p((u -ug)v + (v -vg)v z )f) + ^ (3.64) [r(u r +v z )+14]19 + rp at u + rp—at v - rfr u - rfzv]Ddl fix =0 Equation (3.64) is the full ALE finite element equation transferred to the fixed domain. The terms ug and vg are mesh velocity components and they represent an extra set of unknowns. 49 In chapter 4 a mesh motion scheme is presented to deal with the extra set of unknowns of the grid motion. Coupling equation (3.64) and the FE mesh motion scheme gives enough equations to solve the problem. The above set of equations has been programmed in the Comsol Multiphysics program. This program has versatile capabilities and high degree of flexibility that enables the programming and solving coupled problems containing two or more physics. One of the attractive features of this program is that the finite element formulation could be scripted into the program in the weak form. The user is able to define the variables as well as shape functions and the type of elements which are used for discretization of the domain. One of the limitations of this program is that this program does not have adaptive mesh capability. In ALE application there are two sets of equations, flow equations and mesh motion equations. Basically we can assume one physic for the flow equations and another one for the mesh motion equations. Therefore ALE formulation is amendable to this multiphysics capability and the special boundary conditions required in the jet impingement problem makes this programming option in Comsol an attractive one for the solution of the problem. In chapter 5 the written program has been used to solve the application problem of the jet impingement. 50 Chapter 4 Mesh Motion Schemes As discussed in chapter 3, one of the main differences between the ALE and other FE classical formulation approaches is the existence of mesh velocity as extra unknowns in the final equilibrium equations. In order to deal with these new unknowns we need to develop a scheme for mesh motion. In this chapter, first the main existing methods for mesh motion in the literature are described with discussion of their advantages, disadvantages and suitability for the jet impingement application. A new approach which seems to be more suitable for the jet impingement application is, then, proposed and presented in a FE context. Finally the boundary condition for the free surface problem is addressed and proposed methods of applications are discussed [53-58]. 4.1 Treatment of Mesh Movement In using the Lagrangian type of formulations for large deformation and moving boundary problems, the problem of mesh distortion and element entanglement poses a serious concern [46]. Although some automatic mesh rezoning methods are developed [70], these methods are not yet robust and efficient enough to remedy the mesh distortion for general large deformation and free surface problems [71]. On the other hand, Eulerian formulation methods are normally used to remedy the mesh distortion problems but it is mostly suitable for steady state cases [72]. In ALE formulation the reference system (normally called computational mesh) is not a priori fixed in space nor attached to the body and the finite element mesh may move arbitrarily relative to the material body. This feature provides the potential for efficiently handling mesh distortion 51 as well as contact boundary conditions. In general, the stiffness matrix of ALE may be rectangular because there are two unknowns for each degree of freedom; the material velocity (or displacement) and the mesh velocity (or displacement), while only one equation may be derived from the equilibrium condition. There are generally two ways to get a solvable set of equations in an ALE formulation. One is to specify the mesh velocities before solving the equilibrium equations. This method is straight forward and has been applied by most researchers [59, 70, 72]. In general, researchers have opted for the use of this method when there is no explicit formula for the mesh movement. Also, when the boundary has free movement or the suitable mesh motion as function of time is not predictable beforehand. The method, however, is not a general one and may not be capable of considering the effect of the current material motion. Another method is to set up supplementary constraint equations, i.e., the relations between material displacements or velocities and mesh displacements or velocities. This method is a general one, but it will double the number of unknowns and increase the calculation time significantly. In the case of free boundary, the second method may be easily applied if the boundary motion can be specified as function of time. This is, however, very restrictive and cannot be generally realized in most of the practical applications. Other techniques utilize algebraic interpolation. These are generally simpler to apply, but they introduce curve fitting errors in the description of regions whose boundary may not be exactly described by polynomials of the same order as those appearing in the interpolation functions [69]. The mesh motion algorithm is a critical aspect in the ALE formulation and it affects the computation efficiency significantly. Schreurs [67] introduced a mesh motion technique for an ALE formulation that requires the solution of a simultaneous set of equations. Generating mesh by solving differential equations will generally create high quality mesh [68]. The main 52 drawback is, however, the time and cost of solving the set of differential equations and the errors that may be introduced in the process. This technique is particularly attractive in the jet impingement application because of the existence of free surface and the fact that the free surface tends to correct itself very rapidly. This entails a constant and rapid regulation of the new mesh which is practically possible with this method. This technique is the one that will be used here and will be discussed in more detail. 4.2 Smoothing Technique This technique is normally applied by defining some Partial Differential Equations (PDEs) as governing equations for the mesh displacements or mesh velocity. The choice of the PDEs is critical in the solution process. They are normally chosen based on the physics of the problem. In fluids the PDE that is normally used is the Laplace equation. The main reason for such choice is that Laplace equation is the governing equation for the inviscid fluid velocity, which is the special case for the general Navier-Stokes equations and, therefore, it seems reasonable to choose Laplace equations as governing equation for mesh motion. The Laplacian PDE for a scalar quantityf is given by: V 2 f =^= 0^ (4.1) For a vector field F, this equation could be written as: 0 2 F= f, =0^ (4.2) Equation (4.2) may be used in two different ways for the mesh motion scheme; Laplace Smoothing Method and Winslow Smoothing Method. 53 Laplace Smoothing Method: In this method the term F in Equation (4.2) is assumed to be the mesh (grid) velocity vector vg . In such a case and with the assumption of (ug ,vg) as components of mesh velocity and (x , y) as reference coordinates and (x, y ) as spatial coordinates, mesh velocity has to satisfy the following two equations: a2ug a2ug ^ +^ = 0 a2^ay2 a2 vg a2vg +ax 2 ay2̂ = 0 with ax^g aygu = -1, = at^and^at Winslow Smoothing Method: In this method first the mesh displacements dx=x _x and dy , y _y are defined as degrees of freedom (DOF), and equation (4.6) and (4.7) are solved for these DOFs. a2^a2 X 2 (dx)+ 2 (dx) = 0 a2 axe ( a2 )+ —T(dY) = C/ If the mesh displacement becomes large, the mesh elements eventually have a bad quality or even become inverted. (4.3) (4.4) (4.5) (4.6) (4.7) 54 Depending on the specific model, either smoothing technique may do a better job avoiding inverted mesh elements. The choice of smoothing method and the way boundary conditions are introduced are important in preventing numerical and convergence problems. In the case of jet impingement, and due to the rapid change of boundary, the winslow method which seems to work more smoothly has been applied. In an axi-symmetric domain, with the assumption of R and Z as reference coordinates and r and z as spatial coordinates, the following equations (4.8 and 4.9) must be solved: a 2dr a2dr^_0 art^a ▪ z2 a2dZ a2dZ ^= art^a ▪ z2 4.3 Implementation of Winslow Method in Finite Element We use the regular Galerkin weighted residual method to obtain the finite element form of the winslow smoothing method and we choose a weighting function to be the same as the shape functions. Following the regular Galerkin procedure for Equations (4.8 and 4.9) we get: ( 432 dr + a2 dr )(dP)dCI =0 art^az2 f( 82 dz + dz )(c1I)c1C2 =0ar t az 2 (4.10) (4.11) After applying the Gauss divergence theorem to the above equations, we obtain the FE form of winslow method in axi-symmertic form as: (4.8) (4.9) 55 adr^adr ar (`^',r)+) --- (cfz))an= 0az „(—(dadz z)+ —adz (di )c/r2 = 0ar^r^az^z (4.12) (4.13) Equations (4.12) and (4.13) are two extra FE equations, with degrees of freedom (dr) and (dz). The above equations along with the axi-symmetric ALE FE equations obtained in chapter 3, constitute the proper equations to solve for the unknown material and mesh motions. In the transient process of the solution and in each step, the program solves simultaneously for mesh displacements, material velocities and pressures in the domain. 4.4 Boundary Conditions for Mesh Motion Equations Boundary conditions for the mesh motion are essential to finish the ALE solution. The choice of grid motion boundary condition can cause the mesh to be distorted during the process, and the solution not to converge. In order to understand the exterior boundary condition for the mesh, it should be noted that the mesh cannot exceed the material boundary in the spatial domain. This is an obvious condition since the computation domain must be within the material domain. Therefore, the boundary conditions must be defined in a way to keep the mesh points on the material boundary in space. As discussed above if the constraint is applied in both the tangential and normal directions, it may cause the mesh to be distorted and convergence problems may occur. For this reason, it would be more appropriate to apply the constraint only in the normal direction, i.e., normal to the boundary and let the mesh be regulated in the tangential direction by the Winslow Laplace equations. 56 In a transient problem the above condition is satisfied by the following equation: vg•n=v•n^ (4.14) In which vg is the mesh velocity, v is the material velocity and n is the unit normal vector to the deformed mesh in the spatial domain. This allows the mesh to slip relative to the material along the tangent of the boundary and, generally, improves mesh quality. The above form of boundary condition couples the mesh motion variables (mesh velocities) and the material velocities that would be obtained from Navier-Stokes equation. This coupling may create an unavoidable constraint on the flow with reaction forces on the flow equation on the free surface. The constraint given by Equation (4.14) is generally applied through a Lagrange Multiplier technique. This type of implementation would add another set of unknowns (normally called Lagrange multipliers) and would generally generate zero diagonal terms in the global assembled matrix which requires special pivot searching treatment in the solution process. The process is achieved by adding the following expression to the weak form of the FE equation: A(v — v g).n + 2(ar).n (4.15) where A represents the Lagrange multiplier, dr is the mesh displacement vector and the hat (A) denotes respective weighting function. Based on weighted residual theory the first term in (4.15) will enforce the ALE boundary condition (no material passes the mesh boundary) while the second term prevent the coupling from imposing any undesired constraint on the free surface. The accuracy and effectiveness of this method in case of jet impingement is shown in chapter 5. 57 Chapter 5 Application and Results: Jet Impingement Problem 5.1 Overview of experimental program An overview of the UBC run out table facility is presented in Figure 5.1. The facility consists of an upper tank, a pump, a header, a furnace and necessary pipe and pipe fittings. The overall height of the facility is approximately 6.5 m. The upper tank is located on a 6.5 m tower. The size of the tank is 1.5 x 1.5 x 1 m. The tank is equipped with inlet and outlet pipes and valves, as well as the overflow and drain pipes. The maximum capacity of the upper tank is 1.35 m 3 . Water is usually filled up to a depth of 0.55 m and heated. A 30-kW heater lies inside the upper tank which can heat water up to 65 °C. In the current setup the water temperature is the room temperature. An orifice fitting is installed upstream of the header and connected with a pressure gage. Flow rate is measured by reading the gage. The header could hold up to three nozzles, with the spacing which could be adjusted from 0.05 to 0.09 m. The height of the nozzles from the plate surface could be changed from 0.6 to 2.0 m which is set to 1.5 m in the current experiments. The diameter of the nozzle is 0.019 m. The flow rate can be adjusted up to 1.5 x10 -3 m3/s and is set to 2.5 x10 .-4 m3/s in the current experimental setup. The gravitational acceleration for water is 9.8 m/s 2 and its viscosity is assumed to be 10 -3 (N.s/m2) at the room temperature. In the performed experiments, a sample steel plate with eight thermocouples embedded in the center and different distances from the center of the plate is heated in the furnace, when the plate is around 860 °C it is placed underneath the nozzle (in a way that the center of the plate matches the center of the jet) and the cooling water is turned on. The information from thermocouples is collected and converted to temperature values by a data acquisition system. Using these 58 temperature measurements and existing correlations available for different heat transfer regimes, various heat transfer film coefficients may be calculated on the plate surface. These correlations often contain hydrodynamical parameters such as velocity and pressure within the water flow and on the surface of the plate. Therefore, a hydrodynamic analysis of the jet impingement is necessary in order to study different heat transfer regimes and to analyze and examine the existing correlations. Figure 5.1 Schematic layout of the run out table at UBC The described setup is only used for supplying parameters to the numerical model. Furthur study on the implementation of the outcome of the numerical analysis in to emperical correlations are subject for future researches. 5.2 Analysis setup and Results The elevation of the water head is about 1.5 m above the plate and the impingement action occurs on the surface of the plate. The flow from the nozzle exit down to just above impingement is a free falling jet which would be quite different from the hydrodynamics of 59 impingement. In addition, the dimensions and scale of the experimental and industrial setup lends itself to dividing the problem into two parts: First is a free jet region starting from the nozzle outlet down to just above the impingement and the second is the impingement region starting just few centimeters above the impingement centre and covering the flow on the plate surface. 5.2.1 Free Jet Region In the free jet region the dominant force is garvity. The jet velocity will increase as the water jet moves downward and the jet diameter will decrease to maintain continuity causing what is termed "thinning". This behavior is difficult to capture by Lagrangian type of formulations because of the immense progression of the mesh and the very long distance between the nozzle and the plate. Using Eulerian formulation may also be a problem because of the possible creation of empty elements due to the thinning phenomenon. It is thought that the ALE formulation with the capability of continuous mesh adaptation would be a suitable and logical one to use in this case. 5.2.1.1 Geometry and Boundary Conditions and meshing For the analysis of the axisymmetric problem, we choose a rectangular domain with the left side being the axis of symmetry, the right being the free moving surface, the upper side defines the inflow from the nozzle outlet, and the lower side is the outflow. The analysis domain is shown schematically in Figure 5.2. It should be noted that in Figure 5.2, the scale of the vertical y-axis is quite different from that of the horizontal x-axis since the height is about 1.46 m and the jet radius (domain width) is only about 9.5 x 10 -3 m. The vertical line shown by (R=0) is the axial 60 1 ^4 Free SurfaceAxis ofSym. Inflow Outflow Z-z0- - -0.005^0.005^0.015 symmetry line which is the centerline of the jet. The boundaries are numbered for future references. 111 1.5 1.3 1.1 0.9 0.7 0.5 0.3 0.1 -°7-1).015 Figure 5.2 Free Jet Region (Initial Domain) The boundary conditions used for the analysis are as follows: a) Navier-Stokes Boundary Conditions (based on the numbering scheme in Fig. 5.2): 1:at R=0 (Axis of symmetry), u=0^ (5.1) 2: at Z=0 (Outflow), p=0^ (5.2) 3: at Z=1.46 (Inflow), v=v(r) and u=0^ (5.3) (form of v(r) will be discussed later) 4: at R= 9.5 x 10 -3 (free surface), Total force = 0 (Neutral) or [—pI+,0(Vu+(Vu) -r )iii = 0^ (5.4) 61 The ambient pressure is assumed to be zero and the FE pressure results should be adjusted according to the actual ambient pressure value. b) Mesh Motion Boundary Condition ( for winslow method) 1: at R=0 (axis of symmetry), (dr=r-R=0)^(5.5) Mesh displacement in the horizontal direction is zero 2: at Z=0 (outflow), (dz=z-Z=0)^ (5.6) Mesh displacement in the vertical direction is zero 3: at Z=1.46 (inflow) ,(dz=z-Z=0, dr=r-R=0)^(5.7) Mesh displacement in horizontal and vertical direction is zero 4: at R= 9.5 x10 -3 , vg •n = v-n^ (5.8) The last boundary condition (Equation 5.8) insures that the mesh velocity normal to the free surface is equal to fluid velocity in the direction normal to the boundary. This boundary condition is applied using Lagrange Multiplier as described in chapter 4. It is worth noting that on boundaries number 2 and number 1, the other boundary condition which is not specified is a dirichlet boundary condition. The rectangular domain of the analysis is meshed with uniform mesh of rectangular elements two (along the width) and 400 (along the Length). Three different types of input velocity profile at z =1.46m for the inflow may be used in the analysis. These are: (i) Uniform velocity with a value equal to the mean velocity of the nozzle. 62 (ii) Velocity profile for a fully developed laminar flow in a pipe: It is assumed that the nozzle is long enough and that the flow inside the nozzle is fully developed. This profile may be given by the following expression: v(r)=2v,n (1—( r ) 2 ) ^ (5.10) Where vm is the mean velocity at the nozzle exit. (iii) Using (1/7)th power law profile in a pipe [66]: In this method the velocity is decreased from the center according to the following equation: v(r)=v,n /.817(1— (r I D)) 117 ^ (5.11) The nozzle may not be long enough to assume the inflow velocity as fully developed, however due to the long height of the jet and the fact that the flow velocity approaches uniform profile, second method above is the one used in all analyses to follow. The water flow rate through the nozzle is 2.5 x 104 m3/sec and the diameter of the nozzle is 0.019m and, therefore, the velocity profile may be given by: v(r) = —1.78*(1— (2r/D) 2 )^ (5.12) 5.2.1.2 Approximation and integration of field equation In the mixed finite element procedure used in this analysis, pressure and velocity are the primary unknowns in the analysis. Shape functions (approximation functions) used for this analysis are Lagrange shape functions. In the analysis of free jet region the approximation functions for velocity components (u, v) is assumed second order (quadratic) while those for pressure (p) are assumed to be linear (first order). The shape functions for the displacements (dr, dz ) in the mesh motion equations is also assumed to be of the second order. 63 5.2.1.3 Detailed Solution Procedures and Results: The axi-symmetric form of Navier-Stokes ALE equations and the winslow method are used to solve the free jet problem. Numerous numerical experiments showed that the solution of the problem is highly sensitive to different parameters. The important ones are the choice of initial conditions, mesh density and distribution, mesh aspect ratio and mesh configuration with respect to the free boundary. In many cases the solution diverges and no results could be obtained. In the following a brief account of the procedures used to handle convergence problems is presented. i) Choice of Initial Condition: Due to the nature of the ALE equations and the highly nonlinear and transient form of the problem, the choice of initial conditions highly affect the convergence characteristics of the solution. Numerical experiments showed that starting with zero initial velocities would always lead to convergence problems. In a nonlinear problem, using itterative methods such as Newton- Raphson starting from an initial value which is far from the solution may not lead to the solution. This problem is more significant in ALE problem due to the number of equations and the magnitude of non-linearity. In order to obtain better initial values for velocity components and pressure, the problem is solved in two steps: First; only the Navier-Stokes equations for the domain are solved with fixed boundary (i.e., boundary edge number 4 in Figure 5.2 is fixed). Second; the solution obtained from the fixed boundary case is used as the initial values for velocity components and pressure in the second phase. In this phase, the fully coupled ALE-Navier-Stokes equations with free boundary (boundary edge number 4 in Figure 5.2) are solved. The following boundary conditions are used in the first solution phase: 64 w. 1.781 Ib .4 1 1 0.4 0^001 Figure 5.3 Steady State Velocity Solution (m/s) numerical trials. 1,5 1.4 1.3 12 11 1 0.9 Os as a5 a4 03 0.2 -0.1 ^ .0.01 Navier-Stokes Boundary conditions: 1: at R=0 (Axis of symmetry), u=0^ (5.13) 2: at Z=0 (Outflow), p=0^ (5.14) 3: at Z=1.46 (Inflow), v=v(r) and u=0^ (5.15) (form of v(r) is given by equation 5.12) But for the last boundary: 4: at R= 9.5 x 10 -3 Slip/symmetry boundary condition. (u = 0 , t.[—pI + ,u(Vu + (Vu) T in = 0 )^ (5.16) The reason for choosing slip boundary condition for boundary number 4, is that solving this initial problem, will create zero horizontal velocity on this boundary, which will provide consistent initial velocity for the second phase. This conclusion was only realized after many 65 The solution of free jet region or phase one is shown in Figure (5.3) and it indicates that the velocity doesn't have a uniform distribution at the end of the jet, if the free boundary is assumed fixed in space. This solution is used as initial value for velocity components in the fully coupled ALE problem, or the second phase of the problem. (ii) Choice of mesh density and configurartion: ALE method is very sesitive to the mesh size. In the numerical experiments performed, it was found that very fine mesh often leads to numerical problems and the solution will, generally, not converge. This is normally due to the fact that very fine mesh is prone to be distorted and/or inverted faster. Larger mesh sizes may also create problems in the process of adjusting the mesh regularity, in the mesh motion scheme and may lead to convergence problems. It is unfortunate that there is no systematic procedure for the choice of mesh density in both regions of the above analysis and that the successful mesh size and arrangement are to be obtained by trial and error process. (iii) Choice of mesh aspect ratio: Numerical experiments showed that element aspect ratios close to unity may enhance the convergence characteristics of the problem. Accordingly, quadrilateral mesh with elements close to be square or having aspect ratios close to unity is use in the analysis. (iv) Choice of mesh configurarion with respect to the free boundary: In the jet impingement problem both in free jet region and in the impingement region, the mesh configuration is chosen in a way which is perpendicular to the free boundary. This helps the solution process due to the rapid change of free surface in the first few second of the process and the boundary condition of the free surface which allows the mesh to move perpendicular to the free surface. 66 These choices are taken in to account for both free jet region and impingement region. In the second phase of solution for the free jet region, we utilize the fully coupled ALE Navier- Stokes equations with the original boundary conditions (equations 5.1 to 5.8) for mesh motion and fluid dynamics to solve the problem. The solution to the first step (steady state Navier-Stokes solution) is used as initial conditions for the velocity components. Figures (5.4)-(5.5) show snap shots of the solution results and gives a clear picture of the development of the flow. As shown in the Figures (5.4, 5.5) the free surface starts changing from the beginning of the free jet and the thinning developes throughout the process and continues until it reaches the end of the jet. This process happens so fast with total time of about 2 seconds and the solution reaches the steady state condition and the velocity and the jet thickness doesn't change after this initial time period. The final end diameter of the jet due to thinning is found to be 0.0076 m. The velocity increases as the jet diameter decreases and the more we get close to the end of the jet, the flow takes almost uniform velocity. 67 b) Maw 2.259 115 : 05 d) Atm 3.693 35 3 -125 OS 0 Mn: 0 Mn.: 0.01^0.02 1.5 1,4 1.3 1.2 1,1 1 0.9 0.8 ? 0.7 0.6 as 0,4 - 0.3 0.2 0.1 0 -0.1 ^ -0.01 12 . 04 0.2 .4 Fint0 1.5 1.4 1.3 1.2 1.1 1 0.9 as ? al 0.6 0,5 0.4 0.3 - 0.2 01 -0.1 ^ .0 01 a) 0.01^0 . 02 (m) t5 1.4 1.3 1.2 1 0.9 aB • ? 0.7 • as a4 0.3 - 0.2 01 Mu: 2.73 0 1.5 2 .5 2 1.5 1.4 1.3 1.2 1.1 1 oo 0.8 ? 07 06 • 0.5 0.4 - 0.3 0.2 0 1 0 -01 ^ -0.01 0^001^0 . 02 (m) Figure 5.4 ALE Velocity Solution of Free Jet Region Problem 0^0.01 (m) 0.02 Mn: 0 -0.1 ^ -0.01 Times (a) 0 (b) 0.5 s (c) 0.1 s (d) 0.2 s 68 i2 °Mn: 0 ^1.5^a) 1.4 1.3 - 1.2 - 1.1 1 0.9 - 0.8 ? 0.7 0.6 0.5 0 4 0.3 02 01 - 0 ^ -01 ^ -O. 1 Mau 5.57 3 0 Mn: 0 .5^b) 1.4 1.3 - 1,2 11 1 0.9 0.8 0.7 0 .6 0.5 04 0 3 - 0.2 - OA a ^. .1 ^ -0.01 Mu: 5427 -13 az kin: 0 0.020.01^0.02 0 01 (ml (m1 Miss:5596 Mu: 5.5% Is1.51.41.3 1.2 1.1 1 0.9 0,6 ? 0.7 0.6 0.5 0.4 0.3 0.2 01 0 c) ' 3 Mn: 0 -0.1 ^ -0.01 0 3^0,8 "-E' 0.7 0.6 as 0.4 0.3 0.2 0.1 0 -0.1 0.01^0.02^-0.01^0^0.01 0.02 (m)^ (m) Figure 5.5 ALE Velocity Solution of Free Jet Region Problem Times (a) 0.4 s (b) 0.6 s (c) 0.8 s (d) 2 s 69 It should be noted that in Figures (5.4, 5.5), the scale of the vertical y-axis is quite different from that of the horizontal x-axis since the height is about 1.5 m and the jet radius is only about 9.5 x10-3 m. Figure (5.5) clearly shows the "thinning phenomenon" and shows that the velocity has become almost uniform at the end of the jet. The transient nature of the flow and the development of the steady state pattern are captured in the snap shots given by Figures. Fig (5.5-d) shows the velocity in the domain at time = 2s. At this stage the flow has reached the steady solution. The thinning phenomenon is clear in this Figure. As expected due to gravity and the length of the jet, the velocity increases and the jet diameter decreases. Also the velocity becomes almost uniform at the end of the jet. The other important point is how fast the free surface changes in time. The change in the pressure throughout the domain is negligible due to the fact that in the free jet region the flow is free and the pressure within the flow is equal to ambient pressure. Figure (5.6) shows pressure distribution at t = 2 second. 1.5 14 1.3 1.2 1.1 1 as as "E.' 07 0.6 0.5 0.4 0.3 0.2 01 0 Mn:-39.679 -0.1 ^ -001 0 001^002 (m) Figure 5.6 ALE-Pressure Solution of free Jet Region Problem (time=2 s) 70 0-0 . 01 0 . 01 0.02 r 1.595.5 z10- ' 1.5 b) 0 14 -02 1.3 12 1.1 -04 1 -0.11 0.9 0.8 1 12 06 0.5 0.4 0.3 0.2 0.1 0 Mtn, -2.015e-3 -0.1 0.01^0.02 Matt 2.042*-5 s10-. 0 '25 -35 Mint -3.516..3 Mae: 0 In order to see how mesh motion degrees of freedom change during the transient process, mesh displacement in r direction (dr = r - R) is shown as function of time Figures (5.7, 5.8). The mesh displacement in z direction although goes through some small local variations, remains close to zero within the domain during the process and therefore is not shown. 1.5-^a) 1.4 - 1.3 12 - 1.1 1 0.9 0.8 - 07 - 06 0.5 g4 0.3 0.2 01 0 -0.1 ^ -001 (m) Most 2.33205 ^1.5 •10' 0 ^ 1,4 1.3 1.2 1 . 1 1 0.9 C I-2^ 0.8 -E-* 0. 7 0.6 OS 0.4 0.3 0.2 01 0 0 0.01 0.132 -0.01 0 -01 ^ Mn: -52114e-3 0.01^0.02 15 1.4 1.3 1.2 1.1 1 09 0.8 E 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 c) -0.1 ^ -0.01 d) I .44.5 1.4nt t4.551.-3 (m)^ (m) Figure 5.7 Displacement dr - ALE Solution of Free Jet Region Problem Times: a) 0.1 s b) 0.2 s c) 0.3 s d) 0.4 s 71 Km 0 t5^a) 1.4 1.3 1.2 1.1 1 0.9 - 0.8 ?0.7 0.6 0.5 0.4 0.3 02 0.1 0 -0.1 ^ -001 Max: 0 c10'-o ^ .5^b) 1.4 - 1.3 - 1.2 1.1- 1 0.9 0.8 `E‘ 0 7 0.6 a5 - 0.4 0.3 as 0.1 - 0 - ^ 0.1 ^ 1301 sae'a z CI Mn: -5.676.3 0 . 020^0.01 (rn) a02 0.01 (m) 1,6 1.4 1.3 12 1.1 1 0.9 0.8 ? 0.7 as 0.5 0.4 as a2 o 0 c) .3 him -5.-675.-.1 -01 ^ -0 01 0 0.01^002 (m) Figure 5.8 Displacement dr- ALE Solution of Free Jet Region Problem Times : a) 0.5s b) 1 s c) 2s 72 Ma, 2.732a) b) 0.035 0.035 2.0 .4 1 1.2 0.025 5 0.015 0.025 - 1,0 E 0 015 104 l it 0.4 0.0050.005 0.0 0 144, 0 10^15 x 10 3 -0.005 ^ -5 -0.005 ^ .55^10 (m) 15 x 104 PlioF1 4.650 .5 0,5 Mi, 0 15 x10 3 d)c) 0.025 0.015 0.005 KI.2.601 o 0.025 E 0.015 0.005 046r: 0 ^-0005^ 15^.5 x 10-3 0^5^10 (m) Figure 5.9 Velocity at the Jet End - ALE Solution of Free Jet Region Times a) 0 b) 0.1 s c) 0.2 s d) 0.3 s -0005 ^ -5 5 ^ 10 (m) 74 b)a) MI= 5,09 0.0350035 H 4 40 0250 025 EE 0015 2 0.015 0.105 Mkt 0 t^ I 10^15^-5 x104 Mint 9 10^15 X 10 5 (m) d)c) 0 035 0235 0.0250.625 5 (m) 0.105 E 00150,015 z 0 105an MT: 6 ^ 10^15 fra) x 10-3 Figure 5.10 Velocity-Jet End- ALE Solution of Free Jet Region Problem Times: a) 0.4 s b) 0.45 s c) 0.47 s^d) 0.5 s ^-0 0135 ^ -0 0e6 ^ ^ -5 0 ^ 5 ^ 10 ^ 15^-5 (m) x103 75 0 035 0E125 E 0.015 E05 a) Mat: $45,27 b) 0035 0.025 E 0.015 0 0)5 -0915 5^ 5^10 15 x^ (m) x 104 -0105 5 ^10 (m) c) him S.,5915 0.035 E 0.025 0015 0.005 0 ^ 5 ^ 10^15 (m) x 104 Figure 5.11 Velocity-Jet End- ALE Solution of Free Jet Region Problem Times a) 0.6 s b) 0.7 s c) 1 s -0.005 ^ -5 76 a) Thu-4 Scle& Vitizilo Otlf IW 0,025 E 0.015 -0 005 5 5 (m) 10 c) TO4440.45^Sodom: VoiccOo^Pt Max; OA d) 0.035 0.035 4 0.025 0.025 E 0015 0,015 0105 0E06 0 -0.005 ' -0105 0^5^10 15 0-5 -5 (1n) X 10.3 Max: 5540 5^10^15 01) X 10 .3 Mar CM Mar 1.7E1 lb 1.2 0.035 : 0.025 Tinoto0.4^Sodom Velocity OW fro b) E 0,11 0.015 0.4 0.405 0.2 0 000-: -0.035 5^10-515 (m) X10 3 xe Fig.5.12 Deformed Mesh., End Jet, ALE Solution of Free Jet Region Problem Times a) 0 s b) 0.4 s c) 0.45 s d) 0.5 s 77 E0^5 ft) 10^15 x10.3 Mu: 5141 b) 0.035 0.025 S 0.015 0005 e -01155^10 15 -5 (m) x 10-3 -0;035 -6 x10 3(m) Mar E.595 d) 0.035 0.025 E 0.015 ran: 0.005 -0 CO55^10 15 -5 V7llaree. *maw new tfty 0 5^10^15 (m) e c) 0.035 0.025 0.015 o 0 6 Mira Fig.5.13 Deformed Mesh., End Jet, ALE Solution of Free Jet Region Problem Times a) 0.55 s b) 0.6 s c) 0.7 s d) 1s 78 Verification of the Solution: In order to verify this solution, the continuity equation is checked for the free jet region. After the free jet reaches the steady state solution, due to the fact that no flow is going in or out of the free boundary (boundary number 4), the inflow rate from the boundary number 3 must be equal with the outflow rate from boundary number 2. The following surface integrals could verify the accuracy of ALE method. As expected, accuracy is one of the strong benefits of ALE formulation. The automatic nature of mesh movement and fully coupling of mesh motion equations and flow equations are the contributing factor to this advantage. Inflow: Value of surface integral, Boundary number 3: -2.523406e-4[m 3/s] Outflow:Value of surface integral, Boundary number 2: -2.5230e-4[m 3/s] error: 1.609e-2 % 5.2.2 Impingement Region This section presents the analysis setup and results for impingement region .In this domain the flow of water starts from a few centimeters from the surface and impinges on the surface and radially spread over the surface. Same as free jet region we can model this part with 2D axi- symmetric model. Fig.(5.14) shows a schematic form of the impingement domain. 79 Free Surface Axi- Sym Inflow Out flow No-slip Figure 5.14 Jet Impingement Region For this domain, we can obtain the inflow velocity profile from the solution of the outflow of "free jet region" discussed in setion (5.2.1). The free surface has a Neutral boundary condition which means the force applied at the free surface is zero. Since we don't know the position of the free boundary beforehand, therefore using ALE method as presented in chapters 3 and 4 is the reasonable method to solve this domain. 5.2.2.1 Geometry and Boundary Conditions and meshing After the water impinges on the surface of the plate and spread over it, the depth of the flow is called "film thickness". In the impingement domain the film thickness is unknown before the analysis. This fact and the availability of a curved boundary, makes the choices mentioned in section 5.2.1.1 to avoid convergance problems even more important and critical. In this case in order to prevent the numerical problems, for the initial geometry a reasonable film thickness (0.001 m) is chosen. As the inflow of this region is the outflow of the "free jet region", 80 0.006 E 0.012 0.009 0.003 0.000 the velocity profile and jet diameter of the outflow of the previous region is assumed as the velocity profile and jet diameter of the inflow for this region. The velocity profile at the outflow of the free jet region is uniform velocity ( v = —5 .5m / s) and the jet diameter is 0.00765 m. The axi-symmetric domain of the jet is assumed to have 0.012 height from the surface and 0.025 span from the centerline of the jet. In order to obtain a suitable mesh density and configuration with respect to the free boundary according to choices (ii), (iii) and (iv) of section (5.2.1.1), this domain is first devided in to few subdomains with four boundaries for each subdomain (Fig.(5.15)). Then these subdomains are meshed using mapped mesh (quadrilateral structured mesh) Fig(5.16). This will create a mesh configuration perpendicular to the free surface along the length of this boundary. Also this allows a mesh size proportional to the thickness of the domain at all locations as well as maintaining the aspect ratio of the mesh close to unity. -0.003 ^ -0.001 0.002 0.005 0.008 0.011 0.014 0.017 0.020 0.023 0.026 Fig.5.15 Boundary Mociem, )Jet Impingement Region 81 0.012 0.009 0.006 E 0.003 0.000 -0.003 ^ -0.001 0.002 0.005 0.008 0.011 0.01 4 0.017 0.020 0.023 0.026 (m) Figure 5.16 Meshing , Jet Impingement Region Based on the boundary numbers introduced in Fig (5.15), the Boundary Conditions for the ALE analysis of the jet impingement domain is described for flow equations and mesh motion equations: a) Navier-Stokes Equations: Boundaries: 1, 4, 6^:Axial Symmetry (u = 0) (5.17) 3, 11, 13 :no slip boundary condition (u = v = 0) (5.18) 8 :Inflow (u = 0,v =-5.5m1 s) (5.19) 9,16,17,18,14 :Neutral (Total Force=0) [—pi + p* (Vu + (Vu) T )ln = 0 (5.20) 15 :Outflow (p=0) (5.21) 82 b) Mesh Motion Boundary conditions (for winslow method): Boundaries: 1,4,6^:Mesh displacement in horizontal and vertical direction is zero. (dr = r — R = dz = z —Z =0)^(5.22) 3, 11, 13^:Mesh displacement in horizontal and vertical direction is zero. (dr=r—R=dz=z—Z=0)^(5.23) 8^:Mesh displacement in horizontal and vertical direction is zero.( dr = dz = 0).^ (5.24) 15^:Mesh displacement in horizontal directon is zero (dr =0)^ (5.25) 9,16,17,18,14 :vg .n.v•n (5.26) Mesh velocity normal to the free surface is equal to fluid velocity in the normal direction. This boundary condition is applied using Lagrange Multiplier as described in chapter 4. It is worth noting that for the outflow (boundary 15), only the horizontal mesh displacement is set to zero. This is due to the fact that the film thickness is supposed to change freely from the initial value in order to capture the final position. Therefore the mesh must be allowed to regulate itself in vertical direction. 5.2.2.2 Detailed Solution Procedures and Results: The axi-symmetric form of Navier-Stokes ALE equations and the winslow method are used to solve the impingement region. Numerous numerical experiments showed that the solution of the problem is highly sensitive to different parameters. These parameters as described in section (5.2.1.3) are also applied to the impingement region and most of them are taken in to account in 83 the meshing of the domain. The reasonable initial film thickness has been obtained based on numerous numerical experiments. However due to unknown nature of the final value of "film thickness", higher inflow velocity and more complicated geometry a new treatment for convergance of the problem is suggested and applied. In order to obtain a better initial values for velocity components and pressure ,the problem is solved in two phases. In the first one only the Navier-Stokes equations are solved for the fixed boundaries and much lower inflow velocity as an steady state problem. In the next phase, the obtained solution for velocity and pressure throughout the domain is used as the initial value for velocity and pressure to solve the transient ALE problem with the free bundary. The inflow velocity during this transient process is increased gradually to reach the final inflow velocity. Next this treatment is described in details. a) Phase one :Only solve the Navier-Stokes equations for the above domain with fixed boundary (Fig( 5.15): boundaries 9,16,17,18,19,14) . Boundaries: 1, 4, 6^:Axial Symmetry (u=0)^(5.27) 3, 11, 13^:no slip boundary condition (u = v = 0)^(5.28) 8^:Inflow ( u = 0, v = —0.5)^(5.29) 9,16,17,18,14^: Slip/symmetry boundary condition. v.n = 0 , t.[—pI + ,u(Vu + (Vu)nn = 0)^(5.30) 15^:Outflow (p=0) ^ (5.31) 84 14n: 0 0.012 0.009 -E- 0.006 0.003 0.000 Maxi 0.%456 ! 4 D 10 oa -az 1 As described in section (5.2.1.3) the reason for chosing slip/symmetry boundary condition on the free surface, is the fact that the solution to the above problem will create a better and more consistent initial condition for the fully coupled equations. In this step the inflow velocity is chosen to be (v=-0.5). Fig(5.17) shows the velocity profile for steady state solution in the phase one. Fig(5.18) shows the pressure profile for steady state solution in phase one. Theese solutions for velocity and pressure are then used as initial values for the transient problem. -0.001 0.002 0.005 0.008 0.011 0.014 0.017 0.020 0.023 0.026 (m) Figure 5.17 Velocity, phase one, steady state Jet Impingement Region 85 0.012 0.009 -E` 0.006 0x703 0:000 -0.001 0.002 0.006 0.006 0.011 0.014 0.017 0.020 0.023 0.026 (m) Figure 5.18 Pressure, Phase one, steady state Jet Impingement Region b) Phase Two: Solve the Navier-Stokes equations with mesh motion equations . i) Navier-Stokes Boundary condition: Boundaries: 1, 4, 6^:Axial Symmetry (u = 0)^(5.32) 3, 11, 13^:no slip boundary condition (u = v = 0)^(5.33) 8^:inflow u=0, (v as described in (5.42))^(5.34) 9,16,17,18,14^:Neutral (Total Force=0) {—pI+,u *(Vu + (Vu) T )in = 0^(5.35) 15^:Outflow (p=0)^ (5.36) 86 $4661 0.702 10.6 1 +b,: 0.7 Although the inflow velocity based on this function reaches the desired value at 10 second, the solution process is continued up to 13 second. This is to make sure that we have reached the steady state solution. In the following pictures the development of the transient solution at different times during 13 seconds is presented. Figures (5.19)-(5.28) show the change in velocity and Figures (5.29)-(5.38) show the change in pressure. In order to see how the velocity and pressure are changing with the transient velocity inflow ,we have to pay attention to the maximum and minimum values of the legend in each picture. -0.001 0.002 0.005 0.008 0.011 0.014 0.017 0.020 0.023 0.026 (m) 0.012 0.009 Es 0 006 0.003 0.000 Figure 5.19 Velocity, Phase Two ,ALE Solution, Time=0.1 s 88 0.012 0.009 0.006 0 003 0.000 -0.001 0.002 0.005 0.008 0.011 0.014 0.017 0.020 0.023 (m) 14t4: 0.741 0. 0 Mint 0 0.026 103 .2 0 Figure 5.20 Velocity, Second step ,ALE Solution, Time=0.2 s 0.012 0.009 'ff` 0.006 0.003 0.000 -0.001 Mn: 0 '^ t. 0.002 0.005 0.008 0.011 0.014 0.017 0.0 '20 0.023 0.026 (m) Figure 5.21 Velocity, Second step, ALE Solution, Time=0.5 s 89 0.012 0.009 -E" 0.006 0.003 0.000 14ac 1561 1.4 iat 0. 0 -0.001 0.002 0.005 0.008 0.011 0.014 0.017 0.020 0.023 0.026 (m) Figure 5.22 Velocity, Second step, ALE Solution, Time=2 s Mac 2$3S 2.5 a 0.012 0.009 0.006 0.003 0.000 -0.001 0.002 0.005 0.008 0.011 0.014 0.017 0.020 0.023 0.026 (m) Fig.5.23 Velocity, Second step ,ALE Solution, Time=4 s 90 0.012 0.009 0.006 0.003 0.000 Max: 3.53 .5 ZS 0.4 0 141,; 0 -0.001 0.002 0.005 0.008 0.011 0.014 0.017 0.020 0.023 0.026 (m) Figure 5.24 Velocity, Second step ,ALE Solution, Time=6 s Mop: 4526 4.5 0.012 0.009 0.006 0.000 0.003 2.5 1 0.50 Rm. 0 -0.001 0.002 0.005 0.008 0.011 0.014 0.017 0.020 0.023 0.026 (m) Figure 5.25 Velocity, Second step, ALE Solution, Time=8 s 91 0.5 12 Min: 0 1 O 012 0 009 -E 0.006 0 003 0.000 -0.001 0.002 0.005 0.008 0.011 0.014 0.017 0.020 0.023 0.026 (m) Figure 5.26 Velocity, Second step, ALE Solution, Time=10 s -0.001 0.002 0.005 0.008 0.011 0.014 0.017 0.020 0.023 0.026 (m) Figure 5.27 Velocity, Second step, ALE Solution, Time=1 I s 0.012 0.009 -E's 0.006 0.003 0.000 92 M 5 -5 27 3 Yin: 0 -0.001 0.002 0.005 0.008 0.011 0.014 0.017 0.020 0.023 0.026 (m) Figure 5.28 Velocity, Second step, ALE Solution, Time=13 s 0.012 0.009 0.006 0.003 0.000 Kw 247.3U 0.012 0.009 0.006 0.003 0.000 Kw. 47.496 -0.001 0.002 0.005 0.008 0.011 0.014 0.017 0.020 0.023 0.026 (m) Figure 5.29 Pressure, Second step, ALE Solution, Time=0.1 s 93 I. ^3701-442 !150 I-^1150 0.012 0.009 -E 0.006 0.003 0.000 Min- -16'!65 0 026-0.001 0.002 0.005 0.008 0.011 0.014 0.017 0 020 0.023 (m) Figure 5.30 Pressure, Second step, ALE Solution, Tirne=0.2 s Mint: 276.262 0.012 0.009 0.006 0.003 0.000 1^f^t^t in: 47.314 -0.001 0.002 0.005 0.008 0.011 0.014 0.017 0.020 0.023 0.026 (m) Figure 5.31 Pressure, Second step, ALE Solution, Time=0.5 s 94 0.012 0.009 0.006 0.003 0.000 -0.001 0.002 0.005 0.008 0.011 0.014 0.017 0 020 0.023 0.026 (m) Figure 5.32 Pressure, Second step, ALE Solution, Time-2 s Max: 3245.541 0.012 0.009 0.006 0.003 0.000 tin; -3.0S2 -0.001 0.002 0.005 0.008 0.011 0.014 0.017 0.020 0.023 0.026 (m) Figure 5.33 Pressure, Second step, ALE Solution, Time-4 s 95 Figure 5.34 Pressure, Second step, ALE Solution, Time=6 s '0.002 0.005 0.008 0.011 0.014 0.017 0.020 0.023 (m) Figure 5.35 Pressure, Second step, ALE Solution, Time=8 s O.012 0 009 0 006 0 003 0 000 -0 001 96 Max, 6269.993 Mip: -3,962 '0.002 0.005 0.0080.011 0.014 0.017 0.020 0.023 0.026 (m) 0.012 0.009 "E` 0 006 0 003 0.000 -0 001 P4a.: 15340 404 1.4 11.2 02 Ilex, 1.47544 404 0 0.4 02 0.012 0.009 0.006 0.003 0.000 -0.001 0.002 0.005 0.008 0.011 0.014 0.017 0.020 0.023 0.026 (m) Figure 5.36 Pressure, Second step, ALE Solution, Time=10 s 0.012 0.009 0.006 0.003 0.000 tin: -7007 -0.001 0.002 0.005 0.008 0.011 0.014 0.017 0.020 0.023 0.026 (m) Figure 5.37 Pressure, Second step, ALE Solution, Time=11 s 97 14ax: 1534.4 x1.04 .4 0.4 Od 0 Km,-7007 -0.001 0.002 0.005 0.008 0.011 0.014 0.017 0.020 0.023 0.026 (m) Figure 5.38 Pressure, Second step, ALE Solution, Time=13 s What we can see in the pictures is that the free surface takes almost its final form very fast in less than one tenth of a second and based on the inflow velocity function, velocity gradually increases troughout the domain ,when we reach 10 seconds ,we have reached the maximum velocity for the inflow ,and due to the gradualness of increase in velocity we reach the steady state at almost the same time ,and after that the solution does not change much. In order to show how the free surface takes its approximate form in the very beginning of the first second ,Figs.(5.39)-(5.45) are obtained in very small time steps in the first second.You can see the change in the position of the free surface in the mesh deformation form in these figures. In order to see a more clear image of the change in the curved surface of the impingement region during the preocess some snapshots (Figs(5.46)-(5.54)) are presented.These figures show how the impingement curve originally chosen to be quarter of a circle changes in time. It also shows how in this ALE transient process the mesh smoothly regulates itself and contracts in the normal direction to the free surface. 0.012 0.009 0.006 0.003 0.000 98 Mi.= 0.012 0.009 -Es 0.006 11111111111110•111 1111111,1111111•11 11111111111111111111 111111111111111111111111 1111111111111111111111 1•111111111111111111 1111111•111MUM11 11111111011•1111 111111111111•111111 11111.1111111•111 1111111111111MMII 111111111111111111111111 INIUME1111111 11111111liw 11100.1OW" 0.5 1w 0.2 0.003 0 000 Obn: 0 -0.001 0.002 0.005 0.008 0.011 0.014 0.017 0.020 0.023 0.026 (m) Figure 5.39 Free Surface Movement, ALE Solution, Time=0.01 s Max. 0 593 0.012 0.009 0.006 0 003 0,4 .03 0 000 41, APTdr Min: 0 -0.001 0.002 0.005 0.008 0.011 0.014 0.017 0.020 0.023 0.026 (m) Figure 5.40 Free Surface Movement ,ALE Solution, Time=0.02 s 99 101 0 114i,. 0 Ku; 0 ,605 0.012 0.5 monso■ monssassaumma■ mussam 11111111111111111111 11111111111111111, 0.009 RIMis SR 11111•111111111111 'E"' 0.006 0.3■ •a. 0.003 0 0.000 0 6441:. -0.001 0.002 0.005 0.008 0.011 0.014 0.017 0.020 0.023 0.026 (m) Figure 5.41 Free Surface Movement, ALE Solution, Time=0.05 s 0.012 IN111111111111111111■um= 11111111111111111111111aaaaaaa 111111111111111111 11.1111111110 0.009 as MOM 0.006 0.003 0.000 11 0"Alr -0.001 0.002 0.005 0.008 0.011 0.014 0.017 0.020 0.023 0.026 (m) Figure 5.42 Free Surface Movement, ALE Solution, Time-0.07 s 100 0.1 11111111111111111111111111 1111111•111111111111 111111111111111111101111 111111•11111111111111. 11•111111681111111 111111111•111111.1111111 11111111111111•111111 11111111111115111111111 11111•1111111111111111esanumeaNaimes 11111010 :PPP"• 0.012 0.009 `E- 0.006 0 003 0.000 1464 0.741 0.7 0.6 1^1 0.4 03 2 0 1 0 Air 0.012 0.009 0 006 0.003 0.000 1^ I^I^I^I^I^1 -0.001 0.002 0.005 0.008 0.011 0.014 0.017 0.020 0.023 0.026 (m) Figure 5.43 Free Surface Movement, ALE Solution, Time=0.15 s -0.001 0.002 0.006 0 008 0.011 0.014 0.017 0.020 0.023 0.026 (m) Figure 5.44 Free Surface Movement, ALE Solution, Time=0.2 s 101 0.012 0.009 --E" 0.006 0.003 0.000 0,2 OJ 46n: 0 -0.001 0.002 0.005 0.006 0.011 0.014 0.017 0.020 0.023 0.026 (m) Figure 5.45 Free Surface Movement, ALE Solution, Time=0.3 s 0.007 714,4n=0 Surface: Vekcity Add WO ^ 1144,Q 0.564 0.006 0.005 - 0.004 - 0.003 - 0.002 - 0.001 - 0.000 - -0.001 -0.001 1^I^1^1^I^1^1^I^I^I Mn. 0 0.000 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 (m) Figure 5.46 Mesh Deformation, ALE Solution, Time=0 102 0.006 0.005 0.004 ." 0.003 0.002 0.001 0.000 0.007 Mxz 0572 Os .4 0.3 0. 1 Kw 0 -0.001 -0.001 0.000 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 (m) Figure 5.47 Mesh Deformation, ALE Solution, Time=0.01 s 0.007 0.006 0.005 0.004 0.003 0.002 0.001 0.000 Mix: 0.380 0.1 01, 0 -0.001 ^ -0.001 0.000 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 (m) Figure 5.48 Mesh Deformation, ALE Solution, Time-0.02 s 103 Tr44-.06 Surfrxe: velecty Odd 1 ,4/41 41ax, 0.547 0.S .4 0.3 0 0 0.007 0.006 0.005 0.004 0.003 0.002 0.001 0.000 0 Wn: 0 -0.001 '^' -0.001 0.000 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 (m) Figure 5.49 Mesh Deformation, ALE Solution, Time=0.04 s 0.007 0.006 0.005 0.004 0.003 0.002 0.001 0.000 0 SW 0 -0.001 ^ -0.001 0.000 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 (m) Figure 5.50 Mesh Deformation, ALE Solution, Time=0.06 s Ku: 01113 0.6 0.5 0.4 03 02 0.1 104 0.4 10.1 0,2 0 No: 0 1444.. 0.636 0.6 0.4 0 0.2 0.1 a Kw: 0.007 0.006 0.005 0.004 0.003 0.002 0.001 0.000 -0.001^ I^ 1 -0.001 0.000 0.001 0.002 0 003 0.004 0.005 0.006 0.007 0.008 0.009 (m) Figure 5.51 Mesh Deformation, ALE Solution, Time=0.08 s Mac 0.702 0.7 0.000 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 (m) 0.6 0.007 0.006 0.005 - 0.004 - -E- 0.003 - 0.002 - 0.001 - 0.000 - -0.001 ^ -0.001 Figure 5.52 Mesh Deformation, ALE Solution, Time=0.1 s 105 0.007 0.006 0.005 0.004 "Ls 0.003 0.002 0.001 0.000 M 0.743 0.7 0.6 0 4 i 0.3 0.2 0 .1 0 0 -0.001 ^ -0.001 0.000 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 (m) Figure 5.53 Mesh Deformation, ALE Solution, Time=0.2 s Mu. 0.782 0 03 0.1 0 0.007 0.006 0.005 0.004 0.003 0.002 0.001 0.000 -0.001 ^ -0.001 0.000 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 (m) Figure 5.54 Mesh Deformation, ALE Solution, Time=0.3 s 106 Verification of Solution: In order to verify this solution, continuity equation is checked. So after the jet impingement reached steady state condition the inflow rate (Boundary number 8 in Fig.(5.15)) must be equal with the outflow rate(Boundary number 15). This is due the fact that after we reach steady state condition there is no flow going in and out of the free surface.The following values for surface integrals show the accuracy of the ALE method. Value of surface integral , Boundary number 8: -2.5279e-4[m3/s] Value of surface integral , Boundary number 15: -2.5280e-4[m3/s] error: 2.452e-3 % In order to check the accuracy of the Lagrange Multiplier application on free surface boundary condition, the total force in horizontal and vertical direction on the free surface was integrated. The values for total forces show the effectiveness of the suggested method for free surface boundary condition. Value of total force in horizontal direction: 1.715e-4 [N] Value of total force in vertical direction: -6.820e-6 [N] 5.2.2.3 Analytical verification of the solution: As mentioned in chapter 2, often researchers have used simplified inviscid models in order to study the jet impingement. In order to verify the solutions obtained in chapters 5.2.2.2 and 5.2.1.3 the inviscid model can be used. In the free jet region, Bernoulli's equation could be used between the nuzzle exit and free jet end in order to check the answers. The assumption is the pressure in the free jet region is uniform and equal to the ambient pressure. This assumption is 107 valid for free jet region due to the fact that based on definition in the free jet region the flow is unaware of the plate existence. In that case the velocity of jet at the end of free jet could be obtained by equation (5.43). V = \IV 2 + 2 .g.H (5.43) in which v1 is the is the velocity at the end of free jet, v, is the velocity input from the nuzzle, g is the gravitational acceleration and H the free jet height. Using this equation with the industrial setup parameters used in presented numerical analysis, with the assumption of v, = 0.89 /2-1 the velocity at the end of free jet will be v i = 5.42 11 . In comparison to the numerical results with the average velocity of 5.5 —m , the error of 1.47 % will be obtained. This small margin shows a very close result from viscous numerical model and non-viscous analytical model. Using the analytical result, the jet diameter could be estimated by continuity equation. The flow rate should be equal in the nuzzle exit and the end of the jet. Therefore the diameter of the jet at the end of the jet could be obtained with this formulation: =^ (5.44) v . In which d, is the diameter at the nuzzle and d is the diameter at the end of the jet. After inputting di =19e – 3m in equation (5.44) we would obtain di = 7.7e – 3m . The numerical results show the jet diameter at the end of the jet is 7.66e-3m . The comparison 108 gives the error of 0.52 % which shows a very close estimation. In order to verify the solution in the impingement region the best measure is the pressure at the stagnation point. In order to do that we can use the Bernoulli's eqaution between the nuzzle exit and the stagnation point. Due to the fact that the velocity at the stagnation point is zero, we can derive the following from Bernoulli's equation: 1PS .- 2 pv, 2 + p.g.H (5.45) In which v, is the velocity at the nuzzle exit, H is the distance from the nuzzle exit and the plate, g is the gravitational acceleration p is water density and p, is the pressure at the stagnation point. The pressure value obtained from (5.45) is the relative pressure with respect to ambient pressure. Using (5.45) after substituting the setup appropriate properties the pressure will be obtained as p5 =14900.05Pa . The pressure on the surface of the plate in the impingement region obtained in the numerical analysis is shown in Figure (5.55) . The pressure at the stagnation point is 15343.71 [Pa]. In comparison to analytical solution, the error would be obtained 2.89 % percent. As for the velocity profile the previous research [29] indicates the velocity (outside the boundary layer) of the flow reaches the velocity at the end of the free jet (in which the flow is not aware of the existence of the plate). Figure (5.56) shows the velocity profile at the end of the film flow. The figure shows a boundary layer as expected for a viscous flow and the velocity reaches 5.5 —m . Which is the same velocity as the velocity at the end of free jet region. This verifies the numerical solutions for velocity in the impingement region. 109 0 1 4 (m) 8 9 x10-3 ity fi nVsl Pressure (Pal 4 Figure (5.55) Surface pressure, Impingement region Figure (5.56) Velocity profile, film flow 110 1.6 1.4 1.2 0.8 0.6 0,4 0.2 Chapter 6 Conclusions and Future works In this work an ALE formulation for the flow equations is obtained. This formulation is then transferred into finite element domain. In order to apply this formulation to jet impingement problem, it is then transformed to an axi-symmetric domain Mesh motion is an essential and challenging aspect of ALE method. The way mesh moves can be critical in the solution process and it has a significant impact on the convergence of the problem. Each problem based on specifications, geometry and dimensions is unique in this aspect. Therefore the suitable mesh motion scheme is very much problem dependent. Some of the most important features of jet impingement problem are free surface and fluid-solid interaction. In order to deal with these features, it is suitable to choose a set of partial deferential equations (PDEs) to govern the mesh motion (displacement or velocity). The choice of PDEs is dependent on the physics of the problem. In case of flow the Laplace equations are appropriate since they are the governing equations for the special case of inviscid flow. Therefore the Winslow method based on Laplace equations has been developed in 2D finite element domain and then in axi-symmetric mode. The Winslow method for the case of jet impingement due to the fast moving of free surface has been proven to be more smooth and gradual. This is very critical in the jet impingement problem. If the mesh does not move smoothly, it can cause mesh distortion very fast and the solution will not converge. The applied ALE method in this work is numerically fully coupled. There is a set of equations for the flow equations in which the mesh velocity or displacement appears as extra unknowns. Mesh motion scheme provides another set of equations with mesh displacements or velocities as unknowns. These two sets of equations are then solved simultaneously for the hydrodynamical 111 variables, velocity and pressure, as well as mesh motion variables, mesh velocity or displacements. This coupling can cause complications and create numerical and convergence problems. To avoid these problems some treatments and matures are suggested for the application of jet impingement. Dividing the problem in to two separate problems of free jet region in which the flow is not aware of the solid surface and jet impingement region in which flow changes direction is one of these measures. The outflow velocity profile of the free jet region is used as the inflow for the jet impingement region. In the free jet region due to long distance from the nuzzle exit to the surface the gravitational force is the dominant force. This force causes the velocity to increase considerably while a phenomenon known as "thinning phenomenon" happens and the jet diameter decreases to less than half of the initial diameter. This phenomenon is completely and clearly predicted using ALE finite element method. In ordinary finite element this phenomenon could not be easily predicted using classical formulation methods. The change in the free surface happens very fast, this could cause the mesh to be distorted and the solution process to stop. Some treatments are suggested to prevent this from happening. The solution process of ALE method is extremely sensitive to the choice of initial conditions, mesh density and distribution, mesh aspect ratio and mesh configuration. In order to achieve a suitable initial condition, the problem is divided in to two phases. First the steady state problem with fixed boundary and slip/symmetry boundary condition. The velocity solution of this phase is then used as an initial condition to solve the ALE domain with moving boundary and neutral boundary condition. The mesh size is also critical. If the mesh density is very high, mesh becomes distorted soon and, on the other hand, if the mesh density is low, the accuracy will be affected and the solution may divert and the run may be aborted. 112 The mesh aspect ratio close to unity can help the convergence of the ALE solution. The suitable mesh density is best obtained by trial and error. The mesh configuration is suitable to be chosen normal to the free surface. This helps the surface motion to happen smoothly and avoid convergence problems. In case of impingement region due to its curved boundary, further treatment is necessary. A reasonable position for the initial surface is chosen. The problem is divided in to two phases. First the domain with the fixed chosen boundary and slip/symmetry boundary condition is solved in the steady state mode. In order to avoid convergence problems the inflow velocity in this phase is initially chosen as a portion of the true velocity (less than one tenth). The solution of this phase is then used as initial condition for the ALE problem with moving boundary and neutral free boundary condition. In the transient process of the solution the inflow velocity is gradually increased based on a suggested function until it reaches the true value of the velocity profile. This treatment helped significantly in maintaining stability to the solution process and achieving ALE solution for the velocity and pressure. The other treatments such as mesh density and configuration are also important in this region. The same concept of perpendicular mesh with respect to free boundary is used for impingement region. Since the presented model is viscous model a boundary layer will be formed on the surface of the plate. This is not possible to predict in inviscid models previously used in some researches. The suggested Lagrange Multiplier method in application of ALE free boundary condition is proven to be effective and prevent the coupling of two sets of equations (flow equations and mesh motion equations) to impose external force on the free boundary. 113 In summary the following main points have been achieved and solved in the research project: - ALE finite element formulation for 2D flow has been obtained. - ALE finite element formulation for axi-symmetric domain has been obtained. - Suitable mesh motion scheme has been chosen and converted to axi-symmetric finite element domain. A Lagrange multiplier approach has been proposed for the treatment of ALE kinematic free surface boundary condition and successfully applied to jet impingement free surface. In order to apply the formulation to the industrial setup of the jet impingement, the problem is divided into two separate numerical domains, free jet region and impingement region. - Suitable choices regarding mesh size, mesh density and mesh configuration are suggested to avoid numerical problems of solution process in each region. Special numerical treatment has been suggested to achieve suitable initial condition for the ALE solution process in free jet region and impingement region. Special numerical treatment has been suggested to achieve stability and convergence of the solution process in the impingement region. - Two sets of ALE coupled equations (flow equations and mesh motion equations) with appropriate boundary conditions have been solved simultaneously for hydrodynamical variables (velocity and pressure) and mesh motion variables (mesh displacements) in each region. The results of this analysis are presented. This research can show in applying the ALE formulation to the jet impingement application, generally we face various numerical challenges. These challenges include the high sensitivity of 114 the convergence of the solution process to parameters such as mesh density, mesh configuration with respect to free surface, mesh aspect ratio and initial values. These challenges could be obstacles for the robustness of the simulation code. Future work in this field could be in two basic areas. First the hydrodynamical aspect of jet impingement. Based on the presented formulation extended consideration could be given to ALE solution of jet impingement on a moving plate. In that case the interaction between the velocity of the plate and the impinging flow would require special treatments in defining boundary conditions for the flow. This could create new challenges in dealing with mesh motion on the surface which would require new treatments. On the next level the hydrodynamics of interaction between two or more jets could be studied. In order to achieve more robust ALE formulation for this application, different mesh motion schemes such as Laplace scheme could be applied to the jet impingement problem. In this scheme mesh velocity instead of mesh displacement is taken as mesh motion variable. The new scheme could be adopted in to finite element method. The fully coupled systems of equations for the flow and mesh motion scheme could be solved simultaneously for the mesh velocity as well as hydrodynamical variables. In another effort, VOF (Volume of Fluid) method could be used in order to capture the movement of the free surface of the jet impingement. This method which is based in Finite volume numerical technique has been used for discrete two phase flow problems such as "droplet breakup" and sloshing. 115 Second area for the future work is the heat transfer on the surface of the plate. Based on the presented simulation, the hydrodynamical variables (velocity and pressure) could be used in existing heat transfer correlations to obtain improved and more accurate heat transfer coefficient. Heat transfer coefficient is the necessary factor in any heat transfer analysis on the surface of the steel plate. In addition to that, the velocity and pressure profiles could be used in obtaining new and more accurate correlations for heat transfer analysis. 116 Bibliography [1] Jambunathan, K., Lai, E., Moss, M. & Button, B. 1992 A review of heat transfer data for single circular jet impingement. Int. J. Heat and Fluid Flow. 13(2), 106-115. [2] Viskanta, R. 1993 Heat transfer to impinging isothermal gas and flame jets. Exp. Thermal Fluid Sci. 6, 111-134. [3] Webb, B. & Ma, C.-F. 1995 Single-phase liquid jet impingement heat transfer. Advances in heat transfer. 26, 105-217. [4] D. A. Zumbrunnen, F. P. Incropera and R. Viskanta,1990, Method and Apparatus for measuring heat transfer distributions on moving and stationary plates cooled by a planar liquid jet, Experimental Thermal and Fluid Science, Vol. 3 , pp. 202-213. [5] Baughn, J., Hechanova, A. & Yan, X. 1991 An experimental study of en- trainment effects on the heat transfer from a flat surface to a heated circular impinging jet. J. Heat Transfer. 113, 1023-1025. [6] Baughn, J. & Shimizu, S. 1989 Heat transfer measurements from a surface with uniform heat flux and an impinging jet. J. Heat Transfer. 111, 1096-1098. [7] Behnia, M., Parneix, S. & Durbin, P., 1996 Simulation of jet impingement heat transfer with the k-0- v 2 model.Annual Research Briefs, Center for Turbulence Research, Nasa Ames/Stanford Univ., pp .3 -16. [8] Behnia, M., Parneix, S. & Durbin, P., 1997 Accurate predictions of jet impingement heat transfer. HTD-Vol. 343 National Heat Transfer Conference. 5, Book No H0190, 111-118. [9] Behnia, M., Parneix, S. & Durbin, P.,1998, prediction of heat transfer in an axisymmetric turbulent jet impinging on a flat plate". International Journal of heat and mass transfer ,v.41,pp.1845-1855 117 [10] Behnia, M., Parneix, S. & Durbin, P.,1997 , Accurate modeling of impinging jet heat transfer, Annual Research Briefs , Center for Turbulence Research [11] Durbin, P. 1991 Near-wall turbulence closure without damping functions. Theoretical and Computational Fluid Dynamics 3(1), 1-13. [12]Behnia, M., Parneix ,S. ,Durbin, P. ,1999, "Predictions of turbulent heat transfer in an axisymmetric jet impinging on a heated pedestal " Journal of Heat Transfer ,v.121, pp 43-49, [13] Behnia, M. ,Parneix, S. ,Shabany,Y. ,Durbin ,P. ,1999 "Numerical study of turbulent heat transfer in confined and unconfined impinging jets", International Journal of Heat and Fluid Flow ,v.20 ,pp 1-9 [14] Durbin, P. 1993a A Reynolds-stress model for near-wall turbulence. J. Fluid Mech. 249, 465-498. [15] Durbin, P. 1993b Application of a near-wall turbulence model to boundary layers and heat transfer. Int. J. Heat and Fluid Flow. 14(4), 316-323. [16] Durbin, P. 1995 Separated flow computations with the k-0-v 2 model. AIAA J. 33(4), 659-664 [17]Colucci, D. & Viskanta, R. 1996 Effect of nozzle geometry on local convective heat transfer to a confined impinging air jet. Experimental Thermal and Fluid Science. 13, 71-80. [18] Craft, T., Graham, L. & Launder, B. 1993 Impinging jet studies for turbulence model assessment--II. An examination of the performance of four turbulence models. Int. J. Heat Mass Transfer. 36(10), 2685-2697. [19]Mesbah, M. 1996 An experimental study of local heat transfer to an impinging jet on non-flat surfaces: a cylindrical pedestal and a hemispherically concave 118 surface. PhD Thesis, University of California, Davis. [20]Yan, X. 1993 A preheated-wall transient method using liquid crystals for the measurement of heat transfer on external surfaces and in ducts. PhD Thesis, University of California, Davis. [21] F. Bierbrauer, W.-K. Soh and W.Y.D. Yuen,2001, An Eulerian-Lagrangian Immersed Interface Method for the Cooling of a Hot Steel Strip by an Axisymmetric Water-Jet, Proceedings of the ASME 2001 International Mechanical Engineering Congress and Exposition. [22] F. Bierbrauer, W.-K. Soh and W.Y.D. Yuen, 2001 ,Application of the Eulerian-Lagrangian Method to Water-Jet Cooling of a Hot Moving Strip, ICHMT International Symposium on Advances in Computational Heat Transfer . [23] S.H.Chuang "Numerical Simulation of an impinging jet on a flat plate",Internatinal Journal for numerical Methods in Fluids, v.9 ,pp1413-1426 ,1989 [24] G.Mejak "Finite Element Analysis of Axisymmetric free jet impingement",International Journal for Numerical Methods in Fluids, v.13 ,pp.491-505 ,1991 [25] K.Knowles ,"Comutational studies of impinging jets using k-U turbulence models", International Journal for Numerical Methods in Fluids ,v.22 , pp.799-810 ,1996 [26] Lytle, D. & Webb, B. 1994 Air jet impingement heat transfer at low nozzle-plate spacings. Int. J. Heat Mass Transfer. 37, 1687-1697. [27] H.Braess , P.Wriggers ,2000,"Arbitrary lagrangian eulerian finite element analysis of free surface flow", Computational methods in applied mechanical engineering" , vol.190,pp 95-109 [28] D.J.Phares , G.T.Smedley ,R.C.Flagan ,2000, "The inviscid impingement of a jet with arbitrary velocity profile" Physics of Fluids ,v.12 , 2046-2055 119 [29]A.T.Hauksson, 2001, Experimental study of boiling heat transfer during water jet impingement on a hot steel plate, Master Thesis ,UBC [30] Q.Meng , 2002, Experimental study of transient cooling of a hot steel plate by an impinging circular jet , Master Thesis ,UBC [31] D.H.Wolf, F.P.Incropr, R.Vikanta,1993,"jet impingement boiling" ,Advances in heat transfer ,V.23,1-131 [32] L.Haung and M.S.EI-Genk,1994, "Heat transfer of an impinging jet on a flat surface" ,Int.J.Heat and Mass Transfer, 1915-1923 [33]C.Pan, J.Y.Hwang and T.L.Lin: "The mechanism of heat transfer in transition boiling", Int.J. Heat and Mass Transfer, 32-7, 1989, pp1357-1349 [34]H.Fujimoto, H.Takuda, N.Hatta and R.Viskanta: "Numerical simulation of transient cooling of a hot solid by an impinging free surface jet", Numerical Heat Transfer, Part A, Vol. 36,1999, pp767-780 [35]J.Filipovic, R.Viskanta, F.P.Incropera and T.A.Veslocki: "Thermal behavior of a moving steel strip cooled by an array of planar water jets", Steel Research (Germany), vol. 63, no.10, 1992, pp438-446 [36]N.Hatta, Y.Tanaka, H.Takuda and J.Kokado: "A numerical study on cooling process of hot steel plates by a water curtain", ISJI International, Vol. 29, No. 8, 1989, pp673-679 [37]N.Hatta and H.Osakabe: "Numerical modeling for cooling process of a moving hot plate by a laminar water curtain", ISJI International, Vol. 29, No. 11, 1989, pp. 919-925 [38]H.Fujimoto , N.Hatta ,R. Viskanta ,1999 ,"Numerical simulation of convective heat transfer to a radial free surface jet impinging on a hot solid" ,Heat and Mass Transfer ,V.35,266-272 120 [39]L.E.Malvern ,1969 , Introduction to the mechanics of continuous medium ,Prentice-Hall ,Inc. [40] J.Wang, Arbitrary Lagrangian-Eulerian Method and Its Applications in Solid Mechanics, PhD Thesis, The University of British Columbia, Jan. 1998. [41] H.N.Bayoumi, Arbitrary Lagrangian-Eulerian Formulation for Quasi-Static and Dynamic Metal Forming Simulation, PhD Thsis, The University of British Columbia, May 2001. [42] M. R. Movahhedy, ALE Simulation of Chip Formation in Orthogonal Metal Cutting Process, PhD Thesis, The University of British Columbia, Jan. 2001. [43] A.I.Abdelgalil ,Modeling of Dynamic Fracture Problems Using ALE Finite Element Formulation ,PhD Thesis ,The University of British Columbia ,Feb. 2002 [44] M.S. Gadala and J. Wang, "A practical procedure for mesh motion in arbitrary Lagrangian Eulerian formulation", Engineering with Computers, 14 (1998) 223-234 [45] M.S. Gadala and J. Wang, "ALE formulation and its application in solid mechanics", Computer Methods in Applied Mechanics and Engineering, 167 (1998) 33-55. [46] J. Wang and M.S. Gadala, "Formulation and survey of ALE method in nonlinear solid mechanics", Finite Elements in Analysis and Design, 24 (1997) 253-269. [47] M.S. Gadala "Recent trends in ALE formulation and its applications in solid mechanics" , Computer Methods in Applied Mechanics and Engineering, 193 (2004) 4247-4275. [48] Donea J, Giuliani S, and Halleux JP. An Arbitrary Lagrangian-Eulerian Finite element method for transient dynamic fluid-structure interactions. Comput. Meth. Appl. Mech. Eng. 1982;33(1-3):689-723. 121 [49]Hughes TJR, Liu WK and Zimmermann TK. Lagrangian-Eulerian Finite element formulation for incompressible viscous flows. Comput. Meth. Appl. Mech. Eng. 1981; 29(3):329-349. [50]B. Ramaswamy and M. Kawahara. _Arbitrary Lagrangian Eulerian Finite element method for unsteady,convective, incompressible viscous free surface fluid flows. Int. J. Num. Meth. Vol. 7, p.1053-1075, 1987 [51] John H. Heinbockel , Introduction to Tensor Calculus and Continuum Mechanics [52] Frank M. White ,Fluid Mechanics [53] J. Donea, A. Huerta, J.-Ph. Ponthot, A. Rodrigues-Ferran, "Arbitrary Langrangian-Eulerian Methods", Encyclopedia of Computational Mechanics, Edited by Erwin Stein, Rene de Borst and Thomas J.R. Hughes, John Wiley & Sons, 2004. [54] Belytschko T, Liu WK and Moran B. Nonlinear Finite Elements for Continua and Structures.Wiley: Chichester, 2000 [55] Donea J. Arbitrary Lagrangian-Eulerian Finite element methods. In Computational Methods for Transient Analysis, Belytschko T, Hughes TJR (eds). North-Holland: Amsterdam,1983; 474- 516 . [56] Donea J and Huerta A. Finite Element Methods for Flow Problems. Wiley: Chichester, 2003. [57] Eriksson LE. Practical three-dimensional mesh generation using transfinite interpolation. SIAM J. Sci. Stat. Comput. 1985; 6(3):712-741. [58] Comsol Multiphysics Manual 2005 [59] Gadala MS and Wang J. Simulation of metal forming processes with Finite element methods. Int. J. Numer. Methods Eng. 1999; 44(10):1397-1428. 122 [60] Gadala MS, Movahhedy MR and Wang J. On the mesh motion for ALE modeling of metal forming processes. Finite Elem. Anal. Des. 2002; 38(5):435-459. [61] M.korchynsky, Development of "controlled cooling" practice, Accelerated cooling of steel, P. D. Southwick (ed.),1986,pp.3-14. [62] G. Tacke, H. Litzke and E. Raquet, Investigations into the efficiency of cooling systems for wide-strip hot rolling mills and computer-aided control of strip cooling, Accelerated cooling of steel,P .D. Southwick (ed.), 1986, pp.35-54 [63] S. J. Chen, A. A. Tseng and F. Han, Spray and jet cooling in steel rolling ,Heat Transfer in metals and containersless processing and manufacturing American society of mechanical engineers, Heat-Transfer ,Vol. 162, 1991,pp. 1-11 [64] N. Nakata and M. Militzer , Modelling of microstructure evolution during hot rolling of a 780 MPa high strength steel" ISIJ International vol 45 , pp 82-90 (2005). [65] P. Aran , J. Sudhakar , N. S. Mishra , Modelling of microstructural evolution during accelerated cooling of hot strip on the runout table , steel research 1995, vol. 66, n 10, pp. 416- 423 [66] H. Fujimoto, N. Hatta, R. Viskanta, Numerical simulation of convective heat transfer to a radial free surface jet impinging on a hot solid, Heat and Mass Transfer,1999 vol.35, pp.266-272 [67] Schreurs, P. J. G.; Veldpaus, F. E.; Brekelmans, W. A. M. (1988) Simulation of forming processes, using the arbitrary Eulerian-Lagrangian formulation, Computer Methods Applications in Mechanical Engineering 58, 19-36 [68] Dabke, P. D.; Hague, I.; Srikrishna, M.; Jackson, J. E. (1992) A knowledge-based system for generation and control of finite-element meshes in forging simulation, Journal of Materials Engineering Performance 1, 415-428 123 [69] Harber, R.; Shephard, M. S.; Abel, J. F.; Gallagjer, R. H.; Greenger, D. P. (1981) A general two dimensional graphic finite element preprocessor utilizing discrete transfinite mapping, International Journal of Numerical Methods Engineering, 17,1015-1044 [70] J. Cheng, N. Kikuchi, A mesh re-zoning technique for finite element simulations of metal forming process, Int. J. Numer. Methods Eng. 23 (1986) 219-288. [71] R.B. Haber, A mixed Eulerian—Lagrangian displacement model for large deformation analysis in solid mechanics, Comput. Methods Appl. Mech. Eng. 43 (1984) 277-292. [72] Z.C. Zienkiewicz, Flow formulation for numerical solution of forming processes, in: Pitmann et al. (Eds.),Numerical Analysis of Forming Processes, Wiley, New york, 1984. [73] W. K. Liu, T. Blytschko and H. Chang ,"An Arbitrary Lagrangian-Eulerian Finite Element Method for Path-Dependent materials", Computer Methods in Applied Mechanics and Engineering, Vol.58, pp. 227-245,1986. 124
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Simulation of hydrodynamics of the jet impingement...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Simulation of hydrodynamics of the jet impingement using Arbitrary Lagrangian Eulerian formulation Maghzian, Hamid 2007
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Simulation of hydrodynamics of the jet impingement using Arbitrary Lagrangian Eulerian formulation |
Creator |
Maghzian, Hamid |
Publisher | University of British Columbia |
Date Issued | 2007 |
Description | Controlled cooling is an important part of steel production industry that affects the properties of the outcome steel. Many of the researches done in controlled cooling are experimental. Due to progress in the numerical techniques and high cost of experimental works in this field the numerical work seems more feasible. Heat transfer analysis is the necessary element of successful controlled cooling and ultimately achievement of novel properties in steel. Heat transfer on the surface of the plate normally contains different regimes such as film boiling, nucleate boiling, transition boiling and radiation heat transfer. This makes the analysis more complicated. In order to perform the heat transfer analysis often empirical correlations are being used. In these correlations the velocity and pressure within the fluid domain is involved. Therefore in order to obtain a better understanding of heat transfer process, study of hydrodynamics of the fluid becomes necessary. Circular jet due to its high efficiency has been used vastly in the industry. Although some experimental studies of round jet arrays have been done, yet the characteristics of a single jet with industrial geometric and flow parameters on the surface of a flat plate is not fully understood. Study of hydrodynamics of the jet impingement is the first step to achieve better understanding of heat transfer process. Finite element method as a popular numerical method has been used vastly to simulate different domains. Traditional approaches of finite element method, Lagrangian and Eulerian, each has its own benefits and drawbacks. Lagrangian approach has been used widely in solid domains and Eulerian approach has been widely used in fluid fields. Jet impingement problem, due to its unknown free surface and the change in the boundary, falls in the category of special problems and none of the traditional approaches is suitable for this application. The Arbitrary Lagrangian Eulerian (ALE) formulation has emerged as a technique that can alleviate many of the shortcomings of the traditional Lagrangian and Eulerian formulations in handling these types of problems. Using the ALE formulation the computational grid need not adhere to the material (Lagrangian) nor be fixed in space (Eulerian) but can be moved arbitrarily. Two distinct techniques are being used to implement the ALE formulation, namely the operator split approach and the fully coupled approach. This thesis presents a fully coupled ALE formulation for the simulation of flow field. ALE form of Navier-Stokes equations are derived from the basic principles of continuum mechanics and conservation laws in the fluid. These formulations are then converted in to ALE finite element equations for the fluid flow. The axi-symmetric form of these equations are then derived in order to be used for jet impingement application. In the ALE Formulation as the mesh or the computational grid can move independent of the material and space, an additional set of unknowns representing mesh movement appears in the equations. Prescribing a mesh motion scheme in order to define these unknowns is problem-dependent and has not been yet generalized for all applications. After investigating different methods, the Winslow method is chosen for jet impingement application. This method is based on adding a specific set of partial differential Equations(Laplace equations) to the existing equations in order to obtain enough equations for the unknowns. Then these set of PDEs are converted to finite element equations and derived in axi-symmetric form to be used in jet impingement application. These equations together with the field equations are then applied to jet impingement problem. Due to the number of equations and nonlinearity of the field equations the solution of the problem faces some challenges in terms of convergence characteristics and modeling strategies. Some suggestions are made to deal with these challenges and convergence problems. Finally the numerical treatment and results of analyzing hydrodynamics of the Jet Impingement is presented. The work in this thesis is confined to the numerical simulation of the jet impingement and the specifications of an industrial test setup only have been used in order to obtain the parameters of the numerical model. |
Extent | 7518240 bytes |
Subject |
Lagrangian Eulerian jet impingement hydrodynamics numerical simulation |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2008-02-21 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0066270 |
URI | http://hdl.handle.net/2429/421 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mechanical Engineering |
Affiliation |
Applied Science, Faculty of Mechanical Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2008-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2008_spring_maghzian_hamid.pdf [ 7.17MB ]
- Metadata
- JSON: 24-1.0066270.json
- JSON-LD: 24-1.0066270-ld.json
- RDF/XML (Pretty): 24-1.0066270-rdf.xml
- RDF/JSON: 24-1.0066270-rdf.json
- Turtle: 24-1.0066270-turtle.txt
- N-Triples: 24-1.0066270-rdf-ntriples.txt
- Original Record: 24-1.0066270-source.json
- Full Text
- 24-1.0066270-fulltext.txt
- Citation
- 24-1.0066270.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0066270/manifest