Hydrodynamics of a Droplet in Space

The phenomena related to the flow of fluids are generally complex, and difficult to quantify. New approaches - considering points of view still not explored - may introduce useful tools in the study of Hydrodynamics and the related transport phenomena. The details of the flows and the properties of the fluids must be considered on a very small scale perspective. Consequently, new concepts and tools are generated to better describe the fluids and their properties. This volume presents conclusions about advanced topics of calculated and observed flows. It contains eighteen chapters, organized in five sections: 1) Mathematical Models in Fluid Mechanics, 2) Biological Applications and Biohydrodynamics, 3) Detailed Experimental Analyses of Fluids and Flows, 4) Radiation-, Electro-, Magnetohydrodynamics, and Magnetorheology, 5) Special Topics on Simulations and Experimental Data. These chapters present new points of view about methods and tools used in Hydrodynamics.


Introduction 1.1 Droplet in space
It is considered that our solar system 4.6 billion years ago was composed of a proto-sun and the circum-sun gas disk.In the gas disk, originally micron-sized fine dust particles accumulated by mutual collisions to be 1000 km-sized objects like as planets.Therefore, to understand the planet formation, we have to know the evolution of the dust particles in the early solar gas disk.One of the key materials is a millimeter-sized and spherical-shaped grain termed as "chondrule" observed in chondritic meteorites.
Chondrules are considered to have been formed from molten droplets about 4.6 billion years ago in the solar gas disk (Amelin et al., 2002;Amelin & Krot, 2007).Fig. 1 is a schematic of the formation process of chondrules.In the early solar gas disk, aggregation of the micron-sized dust particles took place before planet formation (Nakagawa et al., 1986).When the dust aggregates grew up to about 1 mm in size (precursor), some astrophysical process heated them to the melting point of about 1600 − 2100 K (Hewins & Radomsky, 1990).The molten dust aggregate became a sphere by the surface tension (droplet), and then cooled again to solidify in a short period of time (chondrule).The formation conditions of chondrules, such as heating duration, maximum temperature, cooling rate, and so forth, have been investigated experimentally by many authors (Blander et al., 1976;Fredriksson & Ringwood, 1963;Harold C. Connolly & Hewins, 1995;Jones & Lofgren, 1993;Lofgren & Russell, 1986;Nagashima et al., 2006;Nelson et al., 1972;Radomsky & Hewins, 1990;Srivastava et al., 2010;Tsuchiyama & Nagahara, 1981-12;Tsuchiyama et al., 1980;2004;Tsukamoto et al., 1999).However, it has been controversial what kind of astronomical event could have produced chondrules in early solar system.The chondrule formation is one of the most serious unsolved problems in planetary science.
The most plausible model for chondrule formation is a shock-wave heating model, which has been tested by many theoreticians (Ciesla & Hood, 2002;Ciesla et al., 2004;Desch & Jr., 2002;Hood, 1998;Hood & Horanyi, 1991;1993;Iida et al., 2001;Miura & Nakamoto, 2006;Miura et al., 2002;Morris & Desch, 2010;Morris et al., 2009;Ruzmaikina & Ip, 1994;Wood, 1984).Fig. 2 is a schematic of dust heating mechanism by the shock-wave heating model.Initially, the chondrule precursors were floating in the gas disk without any large relative velocity against the ambient gas (panel (a)).When a shock wave was generated in the gas disk, the gas behind the shock front was accelerated suddenly.On the other hand, the chondrule Fig. 1.Schematic of formation process of a chondrule.The precursor of chondrule is an aggregate of µm-sized cosmic dusts.The precursor is heated and melted by some mechanism, becomes a sphere by the surface tension, then cools to solidify in a short period of time.
precursors remain un-accelerated because of their inertia.Therefore, after passage of the shock front, the large relative velocity arises between the gas and dust particles (panel (b)).The relative velocity can be considered as fast as about 10 km s −1 (Iida et al., 2001).When the gas molecule collides to the surface of chondrule precursors with such large velocity, its kinetic energy thermalizes at the surface and heats the chondrule precursors, as termed as a gas drag heating.The peak temperature of the precursor is determined by the balance between the gas drag heating and the radiative cooling at the precursor surface (Iida et al., 2001).The gas drag heating is capable to heat the chondrule precursors up to the melting point if we consider a standard model of the early solar gas disk (Iida et al., 2001).

Physical properties of chondrules
The chondrule formation models, including the shock-wave heating model, are required not only to heat the chondrule precursors up to the melting point but also to reproduce other physical and chemical properties of chondrules recognized by observations and experiments.These properties that should be reproduced are summarized as observational constraints (Jones et al., 2000).The reference listed 14 constraints for chondrule formation.To date, there is no chondrule formation model that can account for all of these constraints.
Here, we review two physical properties of chondrules; size distribution and three-dimensional shape.The latter was not listed as the observational constraints in the literature (Jones et al., 2000), however, we would like to include it as an important constraint for chondrule formation.As discussed in this chapter, these two properties strongly relate to the hydrodynamics of molten chondrule precursors in the gas flow behind the shock front.

Size distribution
Fig. 3 shows the size distribution of chondrules compiled from measurement data in some literatures (Nelson & Rubin, 2002;Rubin, 1989;Rubin & Grossman, 1987;Rubin & Keil, 1984).The horizontal axis is the diameter D and the vertical axis is the cumulative fraction of  The precursors of chondrules are in a gas disk around the proto-sun 4.6 billion years ago.The gas and precursors rotate around the proto-sun with almost the same angular velocity, so there is almost no relative velocity between the gas and precursors.(b) If a shock wave is generated in the gas disk by some mechanism, the gas behind the shock front is suddenly accelerated.In contrast, the precursor is not accelerated because of its large inertia.The difference of their behaviors against the shock front causes a large relative velocity between them.The precursors are heated by the gas friction in the high velocity gas flow.chondrules smaller than D in diameter.Table 1 shows the mean diameter and the standard deviation of each measurement.It is found that the chondrule sizes vary according to chondrite type.The mean diameters of chondrules in ordinary chondrites (LL3 and L3) are from 600 µm to 1000 µm.In contrast, ones in enstatite chondrite (EH3) and carbonaceous chondrite (CO3) are from 100 µm to 200 µm.
It should be noted that the true chondrule diameters are slightly larger than the data shown in Fig. 3 and Table 1 because of the following reason.This data was obtained by observations on thin-sections of chondritic meteorites.The chondrule diameter on the thin-section is not necessarily the same as the true one because the thin-section does not always intersect the center of the chondrule.Statistically, the mean and median diameters measured on the thin section are, respectively, √ 2/3 and √ 3/4 of the true diameters (Hughes, 1978).However, we do not take care the difference between true and measured diameters because it is not a substantial issue in this chapter.
It is considered that in the early solar gas disk the dust aggregates have the size distribution from ≈ µm (initial fine dust particles) to a few 1000 km (planets).In spite of the wide
size range of solid materials, sizes of chondrules distribute in a very narrow range of about 100 − 1000 µm.Two possibilities for the origin of chondrule size distribution can be considered; (i) size-sorting prior to chondrule formation, and (ii) size selection during chondrule formation.In the case of (i), we need some mechanism of size-sorting in the early solar gas disk (Teitler et al., 2010, and references therein).In the case of (ii), the chondrule formation model must account for the chondrule size distribution.The latter possibility is what we investigate in this chapter.

Deformation from a perfect sphere
It is considered that spherical chondrule shapes were due to the surface tension when they melted.However, their shapes deviate from a perfect sphere and the deviation is an important clue to identify the formation mechanism.Tsuchiyama et al. (Tsuchiyama et al., 2003) measured the three-dimensional shapes of chondrules using X-ray microtomography.They selected 20 chondrules with perfect shapes and smooth surfaces from 47 ones for further analysis.Their external shapes were approximated as three-axial ellipsoids with axial radii of a, b,andc (a ≥ b ≥ c), respectively.Fig. 4 shows results of the measurement.It is considered that the deviation from a perfect sphere results from the deformation of a molten chondrule before solidification.For example, if the molten chondrule rotates rapidly, the shape becomes oblate due to the centrifugal force (Chandrasekhar, 1965).However, the shapes of chondrules in group-B are prolate rather than oblate.Tsuchiyama et al. (Tsuchiyama et al., 2003) proposed that the prolate chondrules in group-B can be explained by spitted droplets due to the shape instability with high-speed rotation.However, it is not clear whether the transient process such as the shape instability accounts for the range of axial ratio of group-B chondrules or not.

Hydrodynamics of molten chondrule precursors
If chondrules were melted behind the shock front, the molten droplet ought to be exposed to the high-velocity gas flow.The gas flow causes many hydrodynamics phenomena on the molten chondrule droplet as follows.(i) Deformation: the ram pressure deforms the droplet shape from a sphere.(ii) Internal flow: the shearing stress at the droplet surface causes fluid flow inside the droplet.(iii) Fragmentation: a strong gas flow will break the droplet into many tiny fragments.Hydrodynamics of the droplet in high-velocity gas flow strongly relates to the physical properties of chondrules.However, these hydrodynamics behaviors have not been investigated in the framework of the chondrule formation except of a few examples that neglected non-linear effects of hydrodynamics (Kato et al., 2006;Sekiya et al., 2003;Uesugi et al., 2005;2003).
To investigate the hydrodynamics of a molten chondrule droplet in the high-velocity gas flow, we performed computational fluid dynamics (CFD) simulations based on cubic-interpolated propagation/constrained interpolation profile (CIP) method.The CIP method is one of the high-accurate numerical methods for solving the advection equation (Yabe & Aoki, 1991; 385 Hydrodynamics of a Droplet in Space www.intechopen.comYabe et al., 2001).It can treat both compressible and incompressible fluids with large density ratios simultaneously in one program (Yabe & Wang, 1991).The latter advantage is important for our purpose because the droplet density (≈ 3gc m −3 ) differs from that of the gas disk (≈ 10 −8 gcm −3 or smaller) by many orders of magnitude.
In addition, we should pay a special attention how to model the ram pressure of the gas flow.
The gas around the droplet is so rarefied that the mean free path of the gas molecules is an order of about 100 cm if we consider a standard gas disk model.The mean free path is much larger than the typical size of chondrules.This means that the gas flow around the droplet is a free molecular flow, so it does not follow the hydrodynamical equations.Therefore, in our model, the ram pressure acting on the droplet surface per unit area is explicitly given in the equation of motion for the droplet by adopting the momentum flux method as described in section 3.2.2.

Aim of this chapter
The hydrodynamical behaviors of molten chondrules in a high-velocity gas flow are important to elucidate the origin of physical properties of chondrules.However, it is difficult for experimental studies to simulate the high-velocity gas flow in the early solar gas disk, where the gas density is so rarefied that the gas flow around droplets does not follow the hydrodynamics equations.We developed the numerical code to simulate the droplet in a high-velocity rarefied gas flow.In this chapter, we describe the details of our hydrodynamics code and the results.We propose new possibilities for the origins of size distribution and three-dimensional shapes of chondrules based on the hydrodynamics simulations.
We describe the governing equations in section 2 and the numerical procedures in section 3.In section 4, we describe the results of the hydrodynamics simulations regarding the deformation of molten chondrules in the high-velocity rarefied gas flow and discuss the origin of rugby-ball-like shaped chondrules.In section 5, we describe the results regarding the fragmentation of molten chondrules and consider the relation to the size distribution of chondrules.We conclude our hydrodynamics simulations in section 6.

Governing equations
The governing equations are the equation of continuity and the Navier-Stokes equation as follows; where ρ is the density of fluid, u is the velocity, p is the pressure, and µ is the viscosity.The ram pressure of the high-velocity gas flow, F g , is exerted on the surface of the droplet and given by (Sekiya et al., 2003) where n i is the unit normal vector of the surface of the droplet, n g is the unit vector pointing the direction in which the gas flows, and r i is the position of the liquid-gas interface.The delta function δ( r − r i ) means that the ram pressure works only at the interface.The ram

387
Hydrodynamics of a Droplet in Space pressure does not work for n i • n g > 0 because it indicates the opposite surface which does not face the molecular flow.The ram pressure causes the deceleration of the center of mass of the droplet.In our coordinate system co-moving with the center of mass, the apparent gravitational acceleration g should appear in the equation of motion.The surface tension, F s , is given by (Brackbill et al., 1992) where γ s is the fluid surface tension and κ is the local surface curvature.Finally, we consider the equation of state given by where c s is the sound speed.

Numerical methods in hydrodynamics
To solve the equation of continuity (Eq.( 1)) numerically, we introduce a color function φ that changes from 0 to 1 continuously.For incompressible two fluids, a density in each fluid is uniform and has a sharp discontinuity at the interface between these two fluids if the density of a fluid is different from another one.By using the color function, we can distinguish these two fluids as follows; φ = 1forfluid1,φ = 0 for fluid 2, and a region where 0 < φ < 1forthe interface.The density of a fluid element is given by where ρ 1 and ρ 2 are the inherent densities for fluid 1 and fluid 2, respectively.The governing equation for φ is given by The conservation equation for φ (Eq.( 7)) is approximately equivalent to the original one (Eq.( 1)) through the relationship between ρ and φ given by Eq. ( 6) (Miura & Nakamoto, 2007).Therefore, the problem to solve Eq. ( 1) results in to solve Eq. ( 7).We solve Eq. ( 7) using R-CIP-CSL2 method with anti-diffusion technique (sections 3.1.2and 3.1.3).
In this study, the fluid 1 is the molten chondrule and the fluid 2 is the disk gas around the chondrule.The inherent densities are given by ρ 1 = ρ d and ρ 2 = ρ a , where subscripts "d" and "a" mean the droplet and ambient gas, respectively.The other physical values of the fluid element (viscosity µ and sound speed c s ) are given in the same manner as the density ρ, The Navier-Stokes equation (Eq.( 2)) and the equation of state (Eq.( 5)) are separated into two phases; the advection phase and the non-advection phase.The advection phases are written as 388 Hydrodynamics -Advanced Topics www.intechopen.com

Parameter
Sign Value Momentum of gas flow p fm 4000 dyn cm −2 Surface tension γ s 400 dyn cm −1 Viscosity of droplet Density of ambient ρ a 10 −6 gcm −3 Sound speed of ambient c s,a 10 −5 cm s −1 Viscosity of ambient µ a 10 −2 gcm −1 s −1 Droplet radius r 0 500 µm Table 2. Canonical input physical parameters for simulations of molten chondrules exposed to the high-velocity rarefied gas flow.We ought to use these parameters if there is no special description.
We solve above equations using the R-CIP method, which is the oscillation preventing method for advection equation (section 3.1.1).The non-advection phases can be written as where Q is the summation of forces except for the pressure gradient.The problem intrinsic in incompressible fluid is in the high sound speed in the pressure equation.Yabe and Wang (Yabe & Wang, 1991) introduced an excellent approach to avoid the problem (section 3.2.1).It is called as the C-CUP method (Yabe & Wang, 1991).The numerical methods to calculate ram pressure of the gas flow and the surface tension of droplet in Q are described in sections 3.2.2 and 3.2.3,respectively.
The input parameters adopted in this chapter are listed in Table 2.

CIP method
The CIP method is one of the high-accurate numerical methods for solving the advection equation (Yabe & Aoki, 1991;Yabe et al., 2001).In one-dimension, the advection equation is written as where f is a scaler variable of the fluid (e.g., density), u is the fluid velocity in the x-direction, and t is the time.When the velocity u is constant, the exact solution of Eq. ( 10) is given by which indicates a simple translational motion of the spatial profile of f with the constant velocity u.
Let us consider that the values of f on the computational grid points x i−1 , x i ,a n dx i+1 are given at the time step n and denoted by f n i−1 , f n i ,and f n i+1 , respectively.In  by filled circles.From Eq. ( 11), we can obtain the values of f i at the next time step n + 1b y just obtaining f n i at the upstream point x = x i − u∆t,w h er e∆t is the time interval between t n and t n+1 .If the upstream point is not exactly on the grid points, which is a very usual case, we have to interpolate f n i with an appropriate mathematical function composed of f n i−1 , f n i , and so forth.There are some variations of the numerical solvers by the difference of the interpolate function F i (x).One of them is the first-order upwind method, which interpolates f n i by a linear function and satisfies following two constraints; (here we assume that u > 0 and the upstream point for f n i locates left-side of x i ).The other variation is the Lax-Wendroff method, which uses a quadratic polynomial satisfying three constraints; We show these interpolation functions in Fig. 5.
On the contrary, the CIP method interpolates using a cubic polynomial, which satisfies following four constraints; x,i ,where f x ≡ ∂ f /∂x is the spatial gradient of f .The interpolation function is given by where a i , b i , c i ,a n dd i are the coefficients determined from T h e expressions of these coefficients are shown in (Yabe & Aoki, 1991).We show the profile of F i (x) in Fig. 5 with f n x,i−1 = f n x,i = 0.In the CIP method, therefore, we need the values of f n x in addition of f n for solving the advection phase.In the CIP method, f x is treated as an independent variable and updated independently from f as follows.Differentiating Eq. ( 10) with respect to x, we obtain where the second term of the left-hand side indicates the advection term and the right-hand side indicates the non-advection term.The interpolate function for the advection of f x is given by ∂F i /∂x.The non-advection term can be solved analytically by considering that ∂u/∂x is constant.
Additionally, there is an oscillation preventing method in the concept of the CIP method, in which the rational function is used as the interpolate function.The rational function is written as (Xiao et al., 1996) where α i and β i are coefficients.The expressions of these coefficients are shown in (Xiao et al., 1996).Usually, we adopt α i = 1 to prevent oscillation.This method is called as the R-CIP method.The model with α i = 0 corresponds to the normal CIP method.

CIP-CSL2 method
The CIP-CSL2 method is one of the numerical methods for solving the conservative equation.
In one-dimension, the conservative equation is written as Integrating Eq. (15) over x from x i to x i+1 , we obtain where σ i+1/2 ≡ x i+1 x i fdx.F o rf being density, σ i+1/2 corresponds to the mass contained in a computational cell between i and i + 1, so it should be conserved during the time integration.Since the physical meaning of uf in the second term of the left-hand side is the flux of σ per unit area and per unit time, the time evolution of σ is determined by where t n ufdt is the transported value of σ from a region of x < x i to that of x > x i within ∆t.The CIP-CSL2 method uses the integrated function Moreover, since Eq. ( 15) can be rewritten as the same form of Eq. ( 13), we can obtain the updated value f n+1 as well as f n+1 x in the CIP method.
Additionally, there is an oscillation preventing method in the concept of the CIP-CSL2 method, in which the rational function is used for the function D i (x) (Nakamura et al., 2001).This method is called as the R-CIP-CSL2 method.

391
Hydrodynamics of a Droplet in Space www.intechopen.com

Anti-diffusion
To keep the sharp discontinuity in the profile of φ, we explicitly add an diffusion term with a negative diffusion coefficient α (anti-diffusion) to the CIP-CSL2 method (Miura & Nakamoto, 2007).In our model, we have an additional diffusion equation about σ as Eq. ( 18) can be separated into two equations as where J ′ indicates the anti-diffusion flux per unit area and per unit time.Using the finite difference method, we obtain where Ĵ ≡ J ′ /(∆x/∆t) is the mass flux which has the same dimension of σ, α ≡ α/(∆x 2 /∆t) is the dimensionless diffusion coefficient, and Here, it should be noted that σ takes the limited value as 0 ≤ σ ≤ σ m ,w h e r eσ m is the initial value for inside of the droplet.The undershoot (σ < 0) or overshoot (σ > σ m ) are physically incorrect solutions.To avoid that, we replace αi = 0.1 only when σ i−1/2 or σ i+1/2 are out of the appropriate range.We insert the anti-diffusion calculation after the CIP-CSL2 method is completed.

Test calculation
In order to demonstrate the advantage of the CIP method, we carried out one-dimensional advection calculations with various numerical methods.Fig. 6 shows the spatial profiles of f of the test calculations.The horizontal axis is the spatial coordinate x.The initial profile is given by the solid line, which indicates a rectangle wave.We set the fluid velocity u = 1, the intervals of the grid points ∆x = 1, and the time step for the calculation ∆t = 0.2.These conditions give the CFL number ν ≡ u∆t/∆x = 0.2, which indicates that the profile of f moves 0.2 times the grid interval per time step.Therefore, the right side of the rectangle wave will reach x = 80 after 300 time steps and the dashed line indicates the exact solution.The filled circles indicate the numerical results after 300 time steps.
The upwind method does not keep the rectangle shape after 300 time steps and the profile becomes smooth by the numerical diffusion (panel a).In the Lax-Wendroff method, the numerical oscillation appears behind the real wave (panel b).Comparing with above two methods, the CIP method seems to show better solution, however, some undershoots ( f < 0)

392
Hydrodynamics -Advanced Topics www.intechopen.comor overshoots ( f > 1) are observed in the numerical result (panel c).In the R-CIP method, although the faint numerical diffusion has still remained, we obtain the excellent solution comparing with the above methods.
We also show the numerical results of the one-dimensional conservative equation.We use the same conditions with the one-dimensional advection equation.Note that Eq. ( 15) corresponds to Eq. ( 10) when the velocity u is constant.The panel (e) shows the result of the R-CIP-CSL2 method, which is similar to that of the R-CIP method.In the panel (f), we found that the combination of the R-CIP-CSL2 method and the anti-diffusion technique shows the excellent solution in which the numerical diffusion is prevented effectively.

C-CUP method
Using the finite difference method to Eq. ( 9), we obtain (Yabe & Wang, 1991) where the superscripts * and ** indicate the times before and after calculating the non-advection phase, respectively.Since the sound speed is very large in the incompressible fluid, the term related to the pressure should be solved implicitly.In order to obtain the implicit equation for p * * , we take the divergence of the left equation and substitute u * * into the right equation.Then we obtain an equation The problem to solve Eq. ( 24) resolves itself into to solve a set of linear algebraic equations in which the coefficients becomes an asymmetric sparse matrix.After p * * is solved, we can calculate u * * by solving the left equation in Eq. ( 23).

Ram pressure of free molecular flow
The ram pressure of the gas flow is acting on the droplet surface exposed to the high-velocity gas flow.It should be noted that the gas flow around a mm-sized droplet does not follow the hydrodynamical equations because the nebula gas is too rarefied.The mean free path of the nebula gas can be estimated by l = 1/(ns),w h e r es is the collisional cross section of gas molecules and n is the number density of the nebula gas.Typically, we adopt n ≈ 10 14 cm −3 based on the standard model of the early solar system at a distance from the sun of an astronomical unit (Hayashi et al., 1985).Substituting s ≈ 10 −16 cm −2 for the hydrogen molecule (Hollenbach & McKee, 1979), we obtain l ≈ 100 cm.On the other hand, the typical size of chondrules is about a few 100 µm (see Fig. 3).Since the object that disturbs the gas flow is much smaller than the mean free path of the gas, the free stream velocity field is not disturbed except of the direct collision with the droplet (free molecular flow).
Consider that the molecular gas flows for the positive direction of the x-axis.The x-component of the ram pressure F g,x is given by where x i is the position of the droplet surface.This equation can be separated into two equations as where M is the momentum flux of the molecular gas flow.The right equation in Eq. ( 26) means that the momentum flux terminates at the droplet surface.The left equation in Eq. ( 26) means that the decrease of the momentum flux per unit length corresponding to the ram pressure per unit area.

394
Hydrodynamics -Advanced Topics www.intechopen.comUsing the finite difference method to the right equation in Eq. ( 26), we obtain where φ is the smoothed profile of φ (see section 3.2.4),and M i+1 = M i for φi+1 < φi because the momentum flux does not increase when the molecular flow goes outward from inside of the droplet.Similarly, we obtain from the left equation in Eq. ( 26).The momentum flux at upstream is M 0 = p fm .F i r s t ,w e solve Eq. ( 27) and obtain the spatial distribution of the molecular gas flow in all computational domain.Then, we calculate the ram pressure by Eq. ( 28).We calculate the momentum flux M and the ram pressure F g at every time step in numerical simulations.Therefore, these spatial distributions are affected by droplet deformation.

Surface tension
The surface tension is the normal force per unit interfacial area.Brackbill et al. (Brackbill et al., 1992) introduced a method to treat the surface tension as a volume force by replacing the 395 Hydrodynamics of a Droplet in Space www.intechopen.comdiscontinuous interface to the transition region which has some width.According to them, the surface tension is expressed as where [φ] is the jump in color function at the interface between the droplet and the ambient gas.In our definition, we obtain [φ]=1.The curvature is given by where The finite difference method of Eq. ( 31) is shown in (Brackbill et al., 1992).When we calculate the surface tension, we use the smoothed profile of φ (see section 3.2.4).

Smoothing
We can obtain the numerical results keeping the sharp interface between the droplet and the ambient region.However, the smooth interface is suitable for calculating the smooth surface tension.We use the smoothed profile of φ only at the time to calculate the surface tension and the ram pressure acting on the droplet surface.In this study, the smoothed color function φ is calculated by where L 1 , L 2 ,a n dL 3 indicate grid indexes of the nearest, second nearest, and third nearest from the grid point (i, j, k), for example, L 1 =(i + 1, j, k), L 2 =(i + 1, j + 1, k), L 3 =(i + 1, j + 1, k + 1), and so forth.It is easily found that in the three-dimensional Cartesian coordinate system, there are six for L 1 ,t w e l v ef o rL 2 ,a n deigh tforL 3 , respectively.The coefficients are set as We iterate the smoothing five times.Then, we obtain the smooth transition region of about twice grid interval width.We use the smooth profile of φ only when calculating the surface tension and the ram pressure.It should be noted that the original profile φ with the sharp interface is kept unchanged.

Deformation of droplet by gas flow 4.1 Vibrational motion
We assume that the gas flow suddenly affects the initially spherical droplet.Fig. 8 shows the time sequence of the droplet shape and the internal velocity.The horizontal and vertical axes are the x-a n dy-axes, respectively.The solid line is the section of the droplet surface in xy-plane.Arrows show the velocity field inside the droplet.The gas flow comes from the left side of the panel.The panel (a) shows the initial condition for the calculation.The panel (b) shows a snapshot at t = 0.55 msec.The droplet begins to be deformed due to the gas ram pressure.The fluid elements at the surface layer, which is directly facing the gas flow, are blown to the downstream.In contrast, the velocity at the center of the droplet turns to upstream of the gas flow because the apparent gravitational acceleration takes place in our coordinate system.The droplet continues to be deformed further, and after t = 1.0 msec, the degree of deformation becomes maximum (see panel (c)).After that, the droplet begins to recover its shape to the sphere due to the surface tension.The recovery motion is not all but almost over at the panel (d).The droplet repeats the deformation by the ram pressure and the recovery motion by the surface tension until the viscosity dissipates the internal motion of the droplet.beginning because the initial droplet shape is a perfect sphere.The axial ratio decreases as time goes by because of the compression.After about 1 msec, c/b reaches minimum and then increases due to the surface tension.After this, the axial ratio vibrates with a constant frequency and finally the vibrational motion damps due to viscous dissipation.The calculated frequency of the vibrational motion is about 2 msec not depending on p fm .The calculated frequency is consistent with that of a capillary oscillations of a spherical droplet given by (Landau & Lifshitz, 1987).

Overdamping
Fig. 10 shows the time variation of the axial ratio c/b when the viscosity is 100 times larger than that in Fig. 9.It is found that the axial ratio converges on the value at steady state without any vibrational motion.This is an overdamping due to the strong viscous dissipation.

Effect of droplet rotation
We carried out the hydrodynamics simulations of non-rotating molten droplet in previous sections.However, the rotation of the droplet should be taken into consideration as the following reason.A chondrule before melting is an aggregate of numerous fine particles,  so the shape is irregular in general.The irregular shape causes a net torque in an uniform gas flow.Therefore, it is naturally expected that the molten chondrule also rotates at a certain angular velocity.
The angular velocity ω f can be roughly estimated by Iω f ≈ N∆t,whereI is the moment of inertia of chondrule and ∆t is the duration to receive the net torque N. Assuming that the small fraction f of the cross-section of the precursor contributes to produce the net torque N, we obtain N ≈ f πr 3 0 p fm .W e c a n s e t ∆t ≈ π/ω f (a half-rotation period) because the sign of N would change after half-rotation.Substituting I =( 8/15)πr 5 0 ρ d ,w h i c hi st h e moment of inertia for a sphere with an uniform density ρ d , we obtain the angular velocity (Miura, Nakamoto & Doi, 2008) Therefore, in the shock-wave heating model, the droplet should be rotating rapidly if most of the angular momentum is maintained during melting.
In addition, it should be noted that the rotation axis is likely to be perpendicular to the direction of the gas flow unless the chondrule before melting has a peculiar shape as windmill.
Fig. 11 shows the deformation of a rotating droplet in gas flow in a three-dimensional view.The rotation axis is set to be perpendicular to the direction of the gas flow.We use µ d = 399 Hydrodynamics of a Droplet in Space Fig. 11.Three-dimensional view of a rotating molten droplet exposed to a high-velocity gas flow.The object shows the external shape of the droplet (iso-surface of the color function of φ = 0.5).The gas flow comes from the left side (arrow).The rotation axis of the droplet is perpendicular to the direction of the gas flow.After t = 1.0 sec, the droplet shape becomes a prolate.We use µ d = 10 3 poise, p fm = 10 4 dyn cm −2 , ω = 100 rad s −1 ,andr 0 = 1 mm. 10 3 poise, p fm = 10 4 dyn cm −2 , ω = 100 rad s −1 ,andr 0 = 1 mm.It is found that the droplet elongates in a direction of the rotation axis as the time goes by.Fig. 12 shows the time variation of the axial ratios b/a (solid) and c/b (dashed).The major axis a corresponds to the droplet radius in a direction of the gas flow, so the decrease of b/a means the droplet elongation.The axial ratio b/a reaches a steady value of 0.76 after 1 sec.The axial ratio c/b is kept at a constant value of ≈ 0.95 during the calculation, which means that two droplet radius perpendicular to the rotation axis is almost uniform.The droplet shape at the steady state is prolate, in other words, a rugby-ball-like shape.

Origin of prolate chondrule
Why did the droplet shape become prolate?The reason, of course, is due to the droplet rotation.If there is no rotation on the droplet, its shape is only affected by the gas which comes from the fixed direction (see Fig. 13a).In this case, the droplet shape becomes disk-like (oblate) shape because only one axis, which corresponds to the direction of the gas flow, becomes shorter than the other two axes (Sekiya et al., 2003).In contrast, let us consider the case that the droplet is rotating.If the rotation period is much shorter than the viscous deformation timescale, the gas flow averaged during one rotation period can be considered to be axis-symmetrical about the rotation axis (see Fig. 13b).Therefore, the droplet shrinks due to the axis-symmetrical gas flow along directions perpendicular to the rotation axis and becomes prolate if the averaged gas ram pressure is strong enough to overcome the centrifugal force.
Doi (Doi, 2011) derived the analytic solution of deformation of a rotating droplet in gas flow in a case that the gas flow can be approximated as axis-symmetrical around the rotation axis as shown in Fig. 13(b).He considered that the droplet radius is given by r(θ)=r 0 + r 1 (θ), where r 0 is the unperturbed droplet radius and r 1 is the deviation from a perfect sphere.θ is the angle between the position (the origin is the center of the droplet) and the rotation axis.
According to his solution, the droplet deformation is given by 400 Hydrodynamics -Advanced Topics www.intechopen.comwhere W e (Weber number) is the ratio of the ram pressure of the gas flow to the surface tension of the droplet defined as R is the ratio of the centrifugal force to the ram pressure defined as ω is the angular velocity of the rotation, and P l (cos θ) is Legendre polynomials.This solution is applicable under the assumption of r 1 ≪ r 0 .Eq. ( 35) shows that the particle radius becomes the maximum at θ = 0, and minimum at θ = π/2.R = 19/20 is a critical value for the droplet shape to be prolate (R < 19/20) or oblate (R > 19/20).The droplet shape is sphere when R = 19/20 because the ram pressure balances with the centrifugal force.
Fig. 14 shows the droplet shape as functions of the Weber number W e and the normalized centrifugal force R using Eq. ( 35).R = 19/20 (vertical dashed line) is a critical value for the droplet shape to be prolate (R < 19/20) or oblate (R > 19/20).In the prolate region, the axial ratio b/a is less than unity for W e > 0 as shown by contours, but c/b = 1.On the other hand, in the oblate region, the axial ratio c/b is less than unity for W e > 0, but b/a = 1.As W e increases, the degree of deformation increases as shown in decrease of axial ratio b/a or c/b.The blue and red regions show ranges of axial ratios of group-A spherical chondrules and group-B prolate chondrules, respectively.We carried out the hydrodynamics simulations for a wide range of parameters and displayed on this diagram by symbols.It is found that the hydrodynamics simulation results show a good agreement with the analytic solution for awiderangeofW e and R.
Let us consider the shape of chondrule expected from the shock-wave heating model.Adopting ram pressure of the gas flow of p fm = 10 4 dyn cm −2 and the radius of chondrule of r 0 = 1 mm, we obtain W e = 2.5 for γ s = 400 erg cm −2 .According to Eq. ( 34), we evaluate R = 0.06 for f = 0.01.The evaluated value of R is smaller than the critical value of 19/20, so the expected droplet shape is prolate.In addition, the axial ratio b/a comes into a range of group-B prolate chondrules (see Fig. 14).This suggests that the origin of group-B prolate chondrules can be explained by the shock-wave heating model.Of course, it should be noted that the shock-wave heating model does not reproduce the group-B prolate chondrules for arbitrary conditions because W e and R depend on many factors, e.g., p fm , r 0 ,a n d f .Namely, it is possible that different shock conditions produce different chondrule shapes, even out of the range of group-A or -B.This fact, on the contrary, indicates that the chondrule shapes constrain shock conditions suitable for formation of these chondrules.The data of three-dimensional chondrule shapes measured by Tsuchiyama et al. (Tsuchiyama et al., 2003) is definitely valuable, however, the number of samples is twenty at most.We need more data to constrain the chondrule formation mechanism from their three-dimensional shapes.

Direct fragmentation
When the droplet size is too large for the surface tension to keep the droplet shape against the gas ram pressure, the fragmentation will occur.Fig. 15 shows the three-dimensional views of  the droplets suddenly exposed to the gas flow fragment for W e >∼ 6 (Bronshten, 1983, p.96).This results into the fragmentation of droplet for r 0 >∼ 6 mm if we adopt our calculation conditions: p fm = 4000 dyn cm −2 and γ s = 400 dyn cm −1 .Our hydrodynamics simulations agree with the criterion for fragmentation.

Fragmentation via cavitation
Fig. 16 shows the internal pressure inside the droplet for various droplet sizes: r 0 = 3, 4, and 5 mm from panels (a) to (c).We use µ d = 1.3 poise and p fm = 4000 dyn cm −2 .T h e s ed r o p l e t s reach steady states, so their hydrodynamics do not change significantly after these panels.We found a high pressure region at the front of the droplet, and low pressure regions at centers of eddies in all cases.The high pressure is due to the ram pressure of the gas flow.The low pressure in eddy is clearly due to the non-linear effect caused by the advection term in Eq. (2).Surprisingly, the pressure in eddy decreases to almost zero in panels (b) and (c).In the "zero"-pressure region, the vaporization (or boiling) of the liquid would take place because the vapor pressure of the liquid exceeds the internal pressure.This phenomenon is well known as cavitation.We did not take into account the cavitation in our simulations, so no vaporization occurred in the calculation.If the cavitation was taken into consideration, the eddies are no longer maintained because of the cavitation, which would cause the fragmentation of the droplet.
Miura & Nakamoto (Miura & Nakamoto, 2007) proposed the condition for the "zero"-pressure region to appear by considering the balance between the centrifugal force and the pressure gradient force around eddies as ρ d v 2 circ /r eddy ≈ p/r eddy ,w h e r ev circ is the fluid velocity around the eddy, r eddy is the radius of the eddy, and p is the pressure inside the droplet.Substituting p = 2γ s /r 0 from the Young-Laplace equation and v circ ≈ v max = 0.112p fm r 0 /µ d (Sekiya et al., 2003), we obtain This equation gives the critical radius of the droplet above which the cavitation takes place in the center of the eddy.We obtain r 0,cav = 1.3 mm for the calculation condition.In our hydrodynamic simulations, we observed the "zero"-pressure region for r 0 = 4mmorlarger .The inconsistency of cavitation criterion between hydrodynamics simulation and Eq. ( 38) might come from the fact that we substitute the linear solution into v circ .The Sekiya's solution did not take into account the non-linear term in the Navier-Stokes equation.On the other hand, the cavitation would be caused by the non-linear effect.The substitution of the linear solution into the non-linear phenomenon might be a reason for the inconsistency.However, Eq. ( 38) provide us an insight of the cavitation criterion qualitatively.

Comparison with chondrule properties
It was found from the chondrule size distribution (see Fig. 3) that chondrules larger than a few mm in radius are very rare.The origin of the chondrule size distribution has been considered as some size-sorting process prior to chondrule formation in the early solar gas disk (Teitler et al., 2010, and references therein).On the other hand, in the framework of the shock-wave heating model, the upper limit of chondrule sizes can be explained by the fragmentation of a molten chondrule in high-velocity gas flow.The criterion of fragmentation is given by W e = p fm r 0 /γ s ≈ 6.Since the ram pressure of the gas flow is typically p fm ≈ 10 3 − 10 5 dyncm −2 , we obtain the upper limit of chondrule sizes as r max ≈ 0.2 − 20 mm.This is consistent with the fact that chondrules larger than a few mm in radius are very rare.
In addition, our hydrodynamics simulations show a new pathway to the fragmentation by cavitation.The cavitation takes place for W e < 6 if viscosity of the molten chondrule is small.The viscosity decreases rapidly as temperature of the droplet increases.This suggests the following tendency: chondrules that experienced higher maximum temperature during melting have smaller sizes that that experienced lower maximum temperature.On the other hand, the data obtained by Nelson & Rubin (Nelson & Rubin, 2002) showed the tendency opposite from our prediction.They considered the reason of the difference in mean sizes among chondrule textural types being due mainly to parent-body chondrule-fragmentation events and not to chondrule-formation processes in the solar nebula.Therefore, to date, there is no evidence regarding the dependence of chondrule sizes on the maximum temperature.The relation between the chondrule sizes and the maximum temperature should be investigated in the future.

405
Hydrodynamics of a Droplet in Space www.intechopen.com How about the distribution of sizes smaller than the maximum one?Kadono and his colleagues carried out aerodynamic liquid dispersion experiments using shock tube (Kadono & Arakawa, 2005;Kadono et al., 2008).They showed that the size distributions of dispersed droplets are represented by an exponential form and similar form to that of chondrules.In their experimental setup, the gas pressure is too high to approximate the gas flow around the droplet as free molecular flow.We carried out the hydrodynamics simulations of droplet dispersion and showed that the size distribution of dispersed droplets is similar to the Kadono's experiments (Yasuda et al., 2009).These results suggest that the shock-wave heating model accounts for not only the maximum size of chondrules but also their size distribution below the maximum size.
In addition, we recognized a new interesting phenomenon relating to the chondrule formation: the droplets dispersed from the parent droplet collide each other.A set of droplets after collision will fuse together into one droplet if the viscosities are low.In contrary, if the set of droplets solidifies before complete fusion, it will have a strange morphology that is composed of two or more chondrules adhered together.This is known as compound chondrules and has been observed in chondritic meteorites in actuality.The abundance of compound chondrules relative to single chondrules is about a few percents at most (Akaki & Nakamura, 2005;Gooding & Keil, 1981;Wasson et al., 1995).The abundance sounds rare, however, this is much higher comparing with the collision probability of chondrules in the early solar gas disk, where number density of chondrules is quite low (Gooding & Keil, 1981;Sekiya & Nakamura, 1996).In the case of collisions among dispersed droplets, a high collision probability is expected because the local number density is high enough behind the parent droplet (Miura, Yasuda & Nakamoto, 2008;Yasuda et al., 2009).The fragmentation of a droplet in the shock-wave heating model might account for the origin of compound chondrules.

Conclusion
To conclude, hydrodynamics behaviors of a droplet in space environment are key processes to understand the formation of primitive materials in meteorites.
We modeled its three-dimensional hydrodynamics in a hypervelocity gas flow.Our numerical code based on the CIP method properly simulated the deformation, internal flow, and fragmentation of the droplet.We found that these hydrodynamics results accounted for many physical properties of chondrules.

Fig. 2 .
Fig. 2. Schematic of the shock-wave heating model for chondrule formation.(a)The precursors of chondrules are in a gas disk around the proto-sun 4.6 billion years ago.The gas and precursors rotate around the proto-sun with almost the same angular velocity, so there is almost no relative velocity between the gas and precursors.(b) If a shock wave is generated in the gas disk by some mechanism, the gas behind the shock front is suddenly accelerated.In contrast, the precursor is not accelerated because of its large inertia.The difference of their behaviors against the shock front causes a large relative velocity between them.The precursors are heated by the gas friction in the high velocity gas flow.

Fig. 4 .
Fig.4.Three-dimensional shapes of chondrules(Tsuchiyama et al., 2003, and their  unpublished data).a, b,andc are axial radii of chondrules when their shapes are approximated as three-axial ellipsoids (a ≥ b ≥ c), respectively.The textures of these chondrules are 16 porphyritic (open circle), 3 barred-olivine (filed circle), and 1 crypto-crystalline (filled square).The radius of each symbol is proportional to the effective radius of each chondrule r * ≡ (abc) 1/3 ; the largest circle corresponds to r * = 1129 ¯m.For the data of crypto-crystalline, r * = 231 ¯m.Chondrule shapes are classified into two groups: group-A shows the relatively small deformation from the perfect sphere, and group-B is prolate with axial ratio of b/a ≈ 0.7 − 0.8.
Fig. 5, f n are shown 389 Hydrodynamics of a Droplet in Space www.intechopen.com
The superscripts * and ** indicate the time step just before and after the anti-diffusion.The minimum modulus function (minmod) is often used in the concept of the flux limiter and has a non-zero value of sign(a) min(|a|, |b|, |c|) only when a, b,a n dc have the same sign.The value of the diffusion coefficient α is also important.Basically, we take α = −0.1 for the anti-diffusion.

Fig. 7 .
Fig. 7. Spatial distributions of the momentum flux M (a) and the ram pressure F g (b) of the free molecular gas flow around a spherical droplet in xy-plane.The dashed circles are sections of the droplet surfaces in xy-plane.Units of the gray scales are p fm for the panel (a) and dyn cm −3 for the panel (b), respectively.We adopt p fm = 5000 dyn cm −2 in this figure.

Fig. 7 (
Fig. 7(a)  shows the distribution of momentum flux M around two droplets in xy-plane.The dashed circles are the external shapes of large and small droplets.The gray scale is normalized by p fm , so unity (white region) means undisturbed molecular flow and zero (dark region) means no flux because the free molecular flow is obstructed by the droplet.It is found that the gas flow is obstructed only behind the droplets.Fig.7(b)shows the distribution of the ram pressure F g,x calculated from the momentum flux distribution.The ram pressure is acting at the droplet surface where M changes steeply.Note that no ram pressure acts at bottom half of the smaller droplet because the molecular flow is obstructed by the larger one.As shown in Fig.7, the model of ram pressure shown here well reproduces the property of free molecular flow.

Fig. 8 .
Fig. 8. Time evolution of molten droplet exposed to the gas flow.The gas flow comes from the left side of panels.We use p fm = 10 4 dyn cm −2 , r 0 = 500 ¯m, and µ d = 1.3 poise for calculations.

Fig. 9
Fig.9shows the time variation of axial ratio c/b of the droplet.Each curve shows the calculation result for the different value of the ram pressure p fm .The droplet is compressed unidirectionally by the gas flow, so the length of minor axis c corresponds to the half length of droplet axis in the direction of the gas flow.The axial ratio c/b is unity at the

Fig. 9 .
Fig. 9. Vibrational motions of molten droplet; the deformation by the ram pressure and the recovery motion by the surface tension.The horizontal axis is the time since the ram pressure begins to affect the droplet and the vertical axis is the axial ratio of the droplet c/b.E a c h curve shows the calculation result for the different value of the ram pressure p fm .W eu s e r 0 = 500 ¯m and µ d = 1.3 poise for calculations.

Fig. 12 .
Fig. 12.Time evolutions of axial ratios b/a and c/b in the case of Fig. 11.

Fig. 13 .
Fig. 13.The reason why the rotating droplet exposed to the gas flow is deformed to a prolate shape is illustrated.(a) If the droplet does not rotate, it is deformed only by the effect of the gas ram pressure.(b) If the droplet rotates much faster than the deformation due to the gas flow, the time-averaged gas flow can be approximated as axis-symmetrical around the rotation axis.

Fig. 14 .
Fig. 14.Shapes of rotating droplets in gas flow.The horizontal axis is the centrifugal force normalized by ram pressure of the gas flow R. The vertical axis is the Weber number W e .R = 19/20 (vertical dashed line) is a critical value for the droplet shape to be prolate (R < 19/20) or oblate (R > 19/20).Solid lines are contours of axial ratios of b/a (R < 19/20) or c/b (R > 19/20).A ranges of axial ratios of chondrules are shown by colored regions for group-A spherical chondrules (blue) and for group-B prolate chondrules (red), respectively.Symbols are results of hydrodynamics simulations (see legends in figure).Grayed region shows a condition in which the droplet will be fragmented by rapid rotation.thebreak-up droplet.The droplet radius is r 0 = 2 cm, which corresponds to W e = 20.The gas flow comes from the left side of the view along the x-axis.It is found that the droplet shape is deformed as the time goes by (panels (a) and (b)), and leads to fragmentation (panel (c)).The parent droplet breaks up to many smaller pieces.

Fig. 16 .
Fig. 16.Internal pressure inside droplet for different droplet radius r 0 : (a) 3 mm, (b) 4 mm, and (c) 5 mm.The pressure at a region surrounded by a white line decreases to almost zero bytheeddy .W euseµ d = 1.3 poise and p fm = 4000 dyn cm −2 .