**1.3 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;

Yabe 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 (<sup>≈</sup> 3 g cm−3) differs from that of the gas disk

Hydrodynamics of a Droplet in Space 387

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

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

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

The governing equations are the equation of continuity and the Navier-Stokes equation as

where *ρ* is the density of fluid, *u* is the velocity, *p* is the pressure, and *μ* is the viscosity. The

where *n*<sup>i</sup> is the unit normal vector of the surface of the droplet, *n*<sup>g</sup> is the unit vector pointing the direction in which the gas flows, and *r*<sup>i</sup> 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

*<sup>∂</sup><sup>t</sup>* <sup>+</sup> ∇ · (*<sup>ρ</sup><sup>u</sup>*) = 0, (1)

*F*g, is exerted on the surface of the droplet and

+*g*, (2)

*F*<sup>s</sup> + *F*g

*F*<sup>g</sup> = −*p*fm(*n*<sup>i</sup> ·*n*g)*n*g*δ*(*r* −*r*i) for *n*<sup>i</sup> ·*n*<sup>g</sup> ≤ 0, (3)

*ρ*

three-dimensional shapes of chondrules based on the hydrodynamics simulations.

chondrules. We conclude our hydrodynamics simulations in section 6.

*∂u*

ram pressure of the high-velocity gas flow,

*∂ρ*

*<sup>∂</sup><sup>t</sup>* + (*<sup>u</sup>* · <sup>∇</sup> )*<sup>u</sup>* <sup>=</sup> −∇ *<sup>p</sup>* <sup>+</sup> *<sup>μ</sup>*∇<sup>2</sup>*<sup>u</sup>* <sup>+</sup>

(<sup>≈</sup> <sup>10</sup>−<sup>8</sup> g cm−<sup>3</sup> or smaller) by many orders of magnitude.

section 3.2.2.

**1.4 Aim of this chapter**

**2. Governing equations**

given by (Sekiya et al., 2003)

follows;

Fig. 4. Three-dimensional shapes of chondrules (Tsuchiyama et al., 2003, and their unpublished data). *a*, *b*, and *c* 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 *<sup>r</sup>*<sup>∗</sup> <sup>≡</sup> (*abc*)1/3; the largest circle corresponds to *<sup>r</sup>*<sup>∗</sup> <sup>=</sup> 1129 ¯m. For the data of crypto-crystalline, *r*<sup>∗</sup> = 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.

6 Will-be-set-by-IN-TECH

Fig. 4. Three-dimensional shapes of chondrules (Tsuchiyama et al., 2003, and their unpublished data). *a*, *b*, and *c* 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

prolate with axial ratio of *b*/*a* ≈ 0.7 − 0.8.

crypto-crystalline (filled square). The radius of each symbol is proportional to the effective radius of each chondrule *<sup>r</sup>*<sup>∗</sup> <sup>≡</sup> (*abc*)1/3; the largest circle corresponds to *<sup>r</sup>*<sup>∗</sup> <sup>=</sup> 1129 ¯m. For the data of crypto-crystalline, *r*<sup>∗</sup> = 231 ¯m. Chondrule shapes are classified into two groups: group-A shows the relatively small deformation from the perfect sphere, and group-B is

Yabe 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 (<sup>≈</sup> 3 g cm−3) differs from that of the gas disk (<sup>≈</sup> <sup>10</sup>−<sup>8</sup> g cm−<sup>3</sup> 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.

#### **1.4 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.

#### **2. Governing equations**

The governing equations are the equation of continuity and the Navier-Stokes equation as follows;

$$\frac{\partial \rho}{\partial t} + \vec{\nabla} \cdot (\rho \vec{u}) = 0,\tag{1}$$

$$\frac{\partial \vec{u}}{\partial t} + (\vec{u} \cdot \vec{\nabla})\vec{u} = \frac{-\vec{\nabla}p + \mu \nabla^2 \vec{u} + \vec{F}\_\mathbf{\hat{s}} + \vec{F}\_\mathbf{\hat{g}}}{\rho} + \vec{\mathbf{g}}\_\mathbf{\hat{\epsilon}} \tag{2}$$

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)

$$\vec{F}\_{\rm g} = -p\_{\rm fm}(\vec{n}\_{\rm i} \cdot \vec{n}\_{\rm g})\vec{n}\_{\rm g}\delta(\vec{r} - \vec{r}\_{\rm i}) \quad \text{for } \vec{n}\_{\rm i} \cdot \vec{n}\_{\rm g} \le 0,\tag{3}$$

where *n*<sup>i</sup> is the unit normal vector of the surface of the droplet, *n*<sup>g</sup> is the unit vector pointing the direction in which the gas flows, and *r*<sup>i</sup> 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

Parameter Sign Value

Droplet radius *r*<sup>0</sup> 500 *μ*m

for advection equation (section 3.1.1). The non-advection phases can be written as

*∂u <sup>∂</sup><sup>t</sup>* <sup>=</sup> <sup>−</sup>

*∂p <sup>∂</sup><sup>t</sup>* <sup>=</sup> <sup>−</sup>*ρ<sup>c</sup>*

The input parameters adopted in this chapter are listed in Table 2.

description.

and 3.2.3, respectively.

**3.1 Advection phase 3.1.1 CIP method**

written as

velocity *u*.

given at the time step *n* and denoted by *f <sup>n</sup>*

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

We solve above equations using the R-CIP method, which is the oscillation preventing method

∇*p ρ*

2

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

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

*∂ f*

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

Let us consider that the values of *<sup>f</sup>* on the computational grid points *xi*−1, *xi*, and *xi*<sup>+</sup><sup>1</sup> are

*<sup>i</sup>* , and *<sup>f</sup> <sup>n</sup>*

*<sup>i</sup>*−1, *<sup>f</sup> <sup>n</sup>*

*∂ f <sup>∂</sup><sup>t</sup>* <sup>+</sup> *<sup>u</sup>* + *Q ρ* ,

<sup>s</sup>∇ · *<sup>u</sup>*, (9)

*<sup>∂</sup><sup>x</sup>* <sup>=</sup> 0, (10)

*<sup>i</sup>*+1, respectively. In Fig. 5, *<sup>f</sup> <sup>n</sup>* are shown

*f*(*x*;*t*) = *f*(*x* − *ut*; 0), when *u* is constant, (11)

Momentum of gas flow *p*fm 4000 dyn cm−<sup>2</sup> Surface tension *γ*<sup>s</sup> 400 dyn cm−<sup>1</sup> Viscosity of droplet *μ*<sup>d</sup> 1.3 g cm−<sup>1</sup> s−<sup>1</sup> Density of droplet *ρ*<sup>d</sup> 3 g cm−<sup>3</sup> Sound speed of droplet *<sup>c</sup>*s,d <sup>2</sup>×105 cm s−<sup>1</sup> Density of ambient *ρ*<sup>a</sup> 10−<sup>6</sup> g cm−<sup>3</sup> Sound speed of ambient *c*s,a 10−<sup>5</sup> cm s−<sup>1</sup> Viscosity of ambient *μ*<sup>a</sup> 10−<sup>2</sup> g cm−<sup>1</sup> s−<sup>1</sup>

Hydrodynamics of a Droplet in Space 389

pressure does not work for �*n*<sup>i</sup> ·�*n*<sup>g</sup> > 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)

$$\vec{F}\_{\mathbf{s}} = -\gamma\_{\mathbf{s}} \kappa \vec{n}\_{\mathbf{i}} \delta(\vec{r} - \vec{r}\_{\mathbf{i}}) , \tag{4}$$

where *γ*s is the fluid surface tension and *κ* is the local surface curvature. Finally, we consider the equation of state given by

$$\frac{dp}{d\rho} = c\_{\text{s}}^2 \tag{5}$$

where *c*s is the sound speed.

#### **3. 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; *φ* = 1 for fluid 1, *φ* = 0 for fluid 2, and a region where 0 < *φ* < 1 for the interface. The density of a fluid element is given by

$$
\rho = \phi \rho\_1 + (1 - \phi)\rho\_2. \tag{6}
$$

where *ρ*<sup>1</sup> and *ρ*<sup>2</sup> are the inherent densities for fluid 1 and fluid 2, respectively. The governing equation for *φ* is given by

$$
\frac{
\partial \phi
}{
\partial t
} + \vec{\nabla} \cdot (\phi \vec{u}) = 0.
\tag{7}
$$

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.2 and 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 *ρ*<sup>1</sup> = *ρ*<sup>d</sup> and *ρ*<sup>2</sup> = *ρ*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 *ρ*, namely, *μ* = *φμ*<sup>d</sup> + (1 − *φ*)*μ*<sup>a</sup> and *c*<sup>s</sup> = *φc*s,d + (1 − *φ*)*c*s,a, respectively.

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

$$\begin{aligned} \frac{\partial \vec{u}}{\partial t} + (\vec{u} \cdot \vec{\nabla}) \vec{u} &= 0, \\ \frac{\partial p}{\partial t} + (\vec{u} \cdot \vec{\nabla}) p &= 0. \end{aligned} \tag{8}$$

8 Will-be-set-by-IN-TECH

pressure does not work for �*n*<sup>i</sup> ·�*n*<sup>g</sup> > 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, �

where *γ*s is the fluid surface tension and *κ* is the local surface curvature. Finally, we consider

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; *φ* = 1 for fluid 1, *φ* = 0 for fluid 2, and a region where 0 < *φ* < 1 for the

where *ρ*<sup>1</sup> and *ρ*<sup>2</sup> are the inherent densities for fluid 1 and fluid 2, respectively. The governing

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

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 *ρ*<sup>1</sup> = *ρ*<sup>d</sup> and *ρ*<sup>2</sup> = *ρ*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

*<sup>∂</sup><sup>t</sup>* + (�*<sup>u</sup>* · <sup>∇</sup>� )�*<sup>u</sup>* <sup>=</sup> 0,

*∂φ <sup>∂</sup><sup>t</sup>* <sup>+</sup> ∇ ·

R-CIP-CSL2 method with anti-diffusion technique (sections 3.1.2 and 3.1.3).

namely, *μ* = *φμ*<sup>d</sup> + (1 − *φ*)*μ*<sup>a</sup> and *c*<sup>s</sup> = *φc*s,d + (1 − *φ*)*c*s,a, respectively.

*∂*�*u*

*∂p*

*dp <sup>d</sup><sup>ρ</sup>* <sup>=</sup> *<sup>c</sup>*<sup>2</sup>

*F*<sup>s</sup> = −*γ*s*κ*�*n*i*δ*(�*r* −�*r*i), (4)

*ρ* = *φρ*<sup>1</sup> + (1 − *φ*)*ρ*2, (6)

� (*φ*�*u*) = 0. (7)

*<sup>∂</sup><sup>t</sup>* + (�*<sup>u</sup>* · <sup>∇</sup>� )*<sup>p</sup>* <sup>=</sup> 0. (8)

<sup>s</sup> , (5)

�

is given by (Brackbill et al., 1992)

the equation of state given by

where *c*s is the sound speed.

equation for *φ* is given by

as

**3. Numerical methods in hydrodynamics**

interface. The density of a fluid element is given by

*F*s,


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

$$\begin{aligned} \frac{\partial \vec{u}}{\partial t} &= -\frac{\vec{\nabla}p}{\rho} + \frac{\vec{Q}}{\rho},\\ \frac{\partial p}{\partial t} &= -\rho c\_s^2 \vec{\nabla} \cdot \vec{u}, \end{aligned} \tag{9}$$

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.

#### **3.1 Advection phase**

#### **3.1.1 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

$$\frac{\partial f}{\partial t} + \mu \frac{\partial f}{\partial x} = 0,\tag{10}$$

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

$$f(\mathbf{x};t) = f(\mathbf{x} - \boldsymbol{\mu}t; \mathbf{0}), \quad \text{when } \boldsymbol{\mu} \text{ is constant}, \tag{11}$$

which indicates a simple translational motion of the spatial profile of *f* with the constant velocity *u*.

Let us consider that the values of *<sup>f</sup>* on the computational grid points *xi*−1, *xi*, and *xi*<sup>+</sup><sup>1</sup> are given at the time step *n* and denoted by *f <sup>n</sup> <sup>i</sup>*−1, *<sup>f</sup> <sup>n</sup> <sup>i</sup>* , and *<sup>f</sup> <sup>n</sup> <sup>i</sup>*+1, respectively. In Fig. 5, *<sup>f</sup> <sup>n</sup>* are shown

In the CIP method, *fx* is treated as an independent variable and updated independently from

Hydrodynamics of a Droplet in Space 391

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 *fx* is given by *∂Fi*/*∂x*. The non-advection term can be solved analytically by considering that *∂u*/*∂x* is

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

*Fi*(*x*) = *ai*(*<sup>x</sup>* <sup>−</sup> *xi*)<sup>3</sup> <sup>+</sup> *bi*(*<sup>x</sup>* <sup>−</sup> *xi*)<sup>2</sup> <sup>+</sup> *ci*(*<sup>x</sup>* <sup>−</sup> *xi*) + *di*

where *α<sup>i</sup>* and *β<sup>i</sup>* are coefficients. The expressions of these coefficients are shown in (Xiao et al., 1996). Usually, we adopt *α<sup>i</sup>* = 1 to prevent oscillation. This method is called as the R-CIP

The CIP-CSL2 method is one of the numerical methods for solving the conservative equation.

*∂*(*u f*)

*xi*<sup>+</sup><sup>1</sup>

*f dx*. For *f* being density, *σi*+1/2 corresponds to the mass contained in a

*u f dt* is the transported value of *σ* from a region of *x* < *xi* to that of *x* > *xi*

*∂u*

*<sup>∂</sup><sup>x</sup>* , (13)

<sup>1</sup> <sup>+</sup> *<sup>α</sup>iβi*(*<sup>x</sup>* <sup>−</sup> *xi*) , (14)

*<sup>∂</sup><sup>x</sup>* <sup>=</sup> 0. (15)

*xi* = 0, (16)

*<sup>i</sup>*+1/2 − *Ji*<sup>+</sup><sup>1</sup> + *Ji*, (17)

 *x xi*−<sup>1</sup>

*Fi*(*x*)*dx* for the

*∂ fx <sup>∂</sup><sup>x</sup>* <sup>=</sup> <sup>−</sup> *fx*

*f* as follows. Differentiating Eq. (10) with respect to *x*, we obtain

constant.

as (Xiao et al., 1996)

**3.1.2 CIP-CSL2 method**

where *σi*+1/2 ≡

where *Ji* ≡

*∂ fx <sup>∂</sup><sup>t</sup>* <sup>+</sup> *<sup>u</sup>*

method. The model with *α<sup>i</sup>* = 0 corresponds to the normal CIP method.

*∂ f <sup>∂</sup><sup>t</sup>* <sup>+</sup>

*σi*+1/2

unit area and per unit time, the time evolution of *σ* is determined by *σn*+<sup>1</sup> *<sup>i</sup>*+1/2 <sup>=</sup> *<sup>σ</sup><sup>n</sup>*

within Δ*t*. The CIP-CSL2 method uses the integrated function *Di*(*x*) ≡

Eq. (13), we can obtain the updated value *f <sup>n</sup>*+<sup>1</sup> as well as *f <sup>n</sup>*+<sup>1</sup> *<sup>x</sup>* in the CIP method.

*<sup>∂</sup><sup>t</sup>* <sup>+</sup> [*u f* ]

computational cell between *i* and *i* + 1, so it should be conserved during the time integration. Since the physical meaning of *u f* in the second term of the left-hand side is the flux of *σ* per

interpolation when *ui* > 0. The function *Di*(*x*) is a cubic polynomial satisfying following four constraints; *Di*(*xi*−1) = 0, *Di*(*xi*) = *<sup>σ</sup>i*−1/2, *<sup>∂</sup>Di*/*∂x*(*xi*−1) = *Fi*(*xi*−1) = *fi*−1, and *∂Di*/*∂x*(*xi*) = *Fi*(*xi*) = *fi*. Moreover, since Eq. (15) can be rewritten as the same form of

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 *Di*(*x*) (Nakamura et al., 2001). This

In one-dimension, the conservative equation is written as

Integrating Eq. (15) over *x* from *xi* to *xi*+1, we obtain

 *xi*<sup>+</sup><sup>1</sup> *xi*

method is called as the R-CIP-CSL2 method.

 *<sup>t</sup>n*+<sup>1</sup> *tn*

Fig. 5. Interpolate functions with various methods: CIP (solid), Lax-Wendroff (dashed), and first-order upwind (dotted). The filled circles indicate the values of *f* defined on the digitized grid points *xi*−1, *xi*, and *xi*<sup>+</sup><sup>1</sup> before updated.

by filled circles. From Eq. (11), we can obtain the values of *fi* at the next time step *n* + 1 by just obtaining *f <sup>n</sup> <sup>i</sup>* at the upstream point *x* = *xi* − *u*Δ*t*, where Δ*t* is the time interval between *t<sup>n</sup>* and *tn*<sup>+</sup>1. If the upstream point is not exactly on the grid points, which is a very usual case, we have to interpolate *f <sup>n</sup> <sup>i</sup>* with an appropriate mathematical function composed of *<sup>f</sup> <sup>n</sup> i*−1, *f n <sup>i</sup>* , and so forth. There are some variations of the numerical solvers by the difference of the interpolate function *Fi*(*x*). One of them is the first-order upwind method, which interpolates *f n <sup>i</sup>* by a linear function and satisfies following two constraints; *Fi*(*xi*−1) = *<sup>f</sup> <sup>n</sup> <sup>i</sup>*−<sup>1</sup> and *Fi*(*xi*) = *<sup>f</sup> <sup>n</sup> i* (here we assume that *u* > 0 and the upstream point for *f <sup>n</sup> <sup>i</sup>* locates left-side of *xi*). The other variation is the Lax-Wendroff method, which uses a quadratic polynomial satisfying three constraints; *Fi*(*xi*−1) = *<sup>f</sup> <sup>n</sup> <sup>i</sup>*−1, *Fi*(*xi*) = *<sup>f</sup> <sup>n</sup> <sup>i</sup>* , and *Fi*(*xi*<sup>+</sup>1) = *<sup>f</sup> <sup>n</sup> <sup>i</sup>*+1. We show these interpolation functions in Fig. 5.

On the contrary, the CIP method interpolates using a cubic polynomial, which satisfies following four constraints; *Fi*(*xi*−1) = *<sup>f</sup> <sup>n</sup> <sup>i</sup>*−1, *Fi*(*xi*) = *<sup>f</sup> <sup>n</sup> <sup>i</sup>* , *<sup>∂</sup>Fi*/*∂x*(*xi*−1) = *<sup>f</sup> <sup>n</sup> <sup>x</sup>*,*i*−1, and *∂Fi*/*∂x*(*xi*) = *f <sup>n</sup> x*,*i* , where *fx* ≡ *∂ f* /*∂x* is the spatial gradient of *f* . The interpolation function is given by

$$F\_i(\mathbf{x}) = a\_i(\mathbf{x} - \mathbf{x}\_i)^3 + b\_i(\mathbf{x} - \mathbf{x}\_i)^2 + c\_i(\mathbf{x} - \mathbf{x}\_i) + d\_i. \tag{12}$$

where *ai*, *bi*, *ci*, and *di* are the coefficients determined from *f <sup>n</sup> <sup>i</sup>*−1, *<sup>f</sup> <sup>n</sup> <sup>i</sup>* , *<sup>f</sup> <sup>n</sup> <sup>x</sup>*,*i*−1, and *<sup>f</sup> <sup>n</sup> x*,*i* . The expressions of these coefficients are shown in (Yabe & Aoki, 1991). We show the profile of *Fi*(*x*) in Fig. 5 with *f <sup>n</sup> <sup>x</sup>*,*i*−<sup>1</sup> <sup>=</sup> *<sup>f</sup> <sup>n</sup> <sup>x</sup>*,*<sup>i</sup>* <sup>=</sup> 0. In the CIP method, therefore, we need the values of *<sup>f</sup> <sup>n</sup> <sup>x</sup>* in addition of *f <sup>n</sup>* for solving the advection phase.

10 Will-be-set-by-IN-TECH

Fig. 5. Interpolate functions with various methods: CIP (solid), Lax-Wendroff (dashed), and first-order upwind (dotted). The filled circles indicate the values of *f* defined on the digitized

by filled circles. From Eq. (11), we can obtain the values of *fi* at the next time step *n* + 1 by

*t<sup>n</sup>* and *tn*<sup>+</sup>1. If the upstream point is not exactly on the grid points, which is a very usual

*<sup>i</sup>* , and so forth. There are some variations of the numerical solvers by the difference of the interpolate function *Fi*(*x*). One of them is the first-order upwind method, which interpolates

variation is the Lax-Wendroff method, which uses a quadratic polynomial satisfying three

On the contrary, the CIP method interpolates using a cubic polynomial, which satisfies

expressions of these coefficients are shown in (Yabe & Aoki, 1991). We show the profile of

*<sup>i</sup>* , and *Fi*(*xi*<sup>+</sup>1) = *<sup>f</sup> <sup>n</sup>*

*<sup>i</sup>*−1, *Fi*(*xi*) = *<sup>f</sup> <sup>n</sup>*

, where *fx* ≡ *∂ f* /*∂x* is the spatial gradient of *f* . The interpolation function is

*Fi*(*x*) = *ai*(*<sup>x</sup>* <sup>−</sup> *xi*)<sup>3</sup> <sup>+</sup> *bi*(*<sup>x</sup>* <sup>−</sup> *xi*)<sup>2</sup> <sup>+</sup> *ci*(*<sup>x</sup>* <sup>−</sup> *xi*) + *di*, (12)

*<sup>x</sup>*,*<sup>i</sup>* <sup>=</sup> 0. In the CIP method, therefore, we need the values of *<sup>f</sup> <sup>n</sup>*

*<sup>i</sup>* by a linear function and satisfies following two constraints; *Fi*(*xi*−1) = *<sup>f</sup> <sup>n</sup>*

(here we assume that *u* > 0 and the upstream point for *f <sup>n</sup>*

*<sup>i</sup>*−1, *Fi*(*xi*) = *<sup>f</sup> <sup>n</sup>*

where *ai*, *bi*, *ci*, and *di* are the coefficients determined from *f <sup>n</sup>*

*<sup>i</sup>* at the upstream point *x* = *xi* − *u*Δ*t*, where Δ*t* is the time interval between

*<sup>i</sup>* with an appropriate mathematical function composed of *<sup>f</sup> <sup>n</sup>*

*i*−1,

*i*

*<sup>x</sup>*,*i*−1, and

*x*,*i* . The

*<sup>x</sup>* in

*<sup>i</sup>*−<sup>1</sup> and *Fi*(*xi*) = *<sup>f</sup> <sup>n</sup>*

*<sup>i</sup>* locates left-side of *xi*). The other

*<sup>i</sup>*+1. We show these interpolation

*<sup>i</sup>* , *<sup>∂</sup>Fi*/*∂x*(*xi*−1) = *<sup>f</sup> <sup>n</sup>*

*<sup>i</sup>* , *<sup>f</sup> <sup>n</sup>*

*<sup>x</sup>*,*i*−1, and *<sup>f</sup> <sup>n</sup>*

*<sup>i</sup>*−1, *<sup>f</sup> <sup>n</sup>*

grid points *xi*−1, *xi*, and *xi*<sup>+</sup><sup>1</sup> before updated.

following four constraints; *Fi*(*xi*−1) = *<sup>f</sup> <sup>n</sup>*

*<sup>x</sup>*,*i*−<sup>1</sup> <sup>=</sup> *<sup>f</sup> <sup>n</sup>*

addition of *f <sup>n</sup>* for solving the advection phase.

just obtaining *f <sup>n</sup>*

*f n*

*f n*

case, we have to interpolate *f <sup>n</sup>*

constraints; *Fi*(*xi*−1) = *<sup>f</sup> <sup>n</sup>*

*x*,*i*

functions in Fig. 5.

*∂Fi*/*∂x*(*xi*) = *f <sup>n</sup>*

*Fi*(*x*) in Fig. 5 with *f <sup>n</sup>*

given by

In the CIP method, *fx* is treated as an independent variable and updated independently from *f* as follows. Differentiating Eq. (10) with respect to *x*, we obtain

$$\frac{\partial f\_{\mathbf{x}}}{\partial t} + \mu \frac{\partial f\_{\mathbf{x}}}{\partial \mathbf{x}} = -f\_{\mathbf{x}} \frac{\partial \mu}{\partial \mathbf{x}'} \tag{13}$$

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 *fx* is given by *∂Fi*/*∂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)

$$F\_{l}(\mathbf{x}) = \frac{a\_{l}(\mathbf{x} - \mathbf{x}\_{l})^{3} + b\_{l}(\mathbf{x} - \mathbf{x}\_{l})^{2} + c\_{l}(\mathbf{x} - \mathbf{x}\_{l}) + d\_{l}}{1 + a\_{l}\beta\_{l}(\mathbf{x} - \mathbf{x}\_{l})},\tag{14}$$

where *α<sup>i</sup>* and *β<sup>i</sup>* are coefficients. The expressions of these coefficients are shown in (Xiao et al., 1996). Usually, we adopt *α<sup>i</sup>* = 1 to prevent oscillation. This method is called as the R-CIP method. The model with *α<sup>i</sup>* = 0 corresponds to the normal CIP method.

#### **3.1.2 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

$$\frac{\partial f}{\partial t} + \frac{\partial (uf)}{\partial x} = 0.\tag{15}$$

Integrating Eq. (15) over *x* from *xi* to *xi*+1, we obtain

$$\frac{\partial \sigma\_{i+1/2}}{\partial t} + [\![\mu f] \!]\_{\mathbf{x}\_i}^{\mathbf{x}\_{i+1}} = \mathbf{0},\tag{16}$$

where *σi*+1/2 ≡ *xi*<sup>+</sup><sup>1</sup> *xi f dx*. For *f* 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 *u f* 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

$$
\sigma\_{i+1/2}^{n+1} = \sigma\_{i+1/2}^{n} - f\_{i+1} + f\_{i\nu} \tag{17}
$$

where *Ji* ≡ *<sup>t</sup>n*+<sup>1</sup> *tn u f dt* is the transported value of *σ* from a region of *x* < *xi* to that of *x* > *xi* within Δ*t*. The CIP-CSL2 method uses the integrated function *Di*(*x*) ≡ *x xi*−<sup>1</sup> *Fi*(*x*)*dx* for the interpolation when *ui* > 0. The function *Di*(*x*) is a cubic polynomial satisfying following four constraints; *Di*(*xi*−1) = 0, *Di*(*xi*) = *<sup>σ</sup>i*−1/2, *<sup>∂</sup>Di*/*∂x*(*xi*−1) = *Fi*(*xi*−1) = *fi*−1, and *∂Di*/*∂x*(*xi*) = *Fi*(*xi*) = *fi*. Moreover, since Eq. (15) can be rewritten as the same form of Eq. (13), we can obtain the updated value *f <sup>n</sup>*+<sup>1</sup> as well as *f <sup>n</sup>*+<sup>1</sup> *<sup>x</sup>* 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 *Di*(*x*) (Nakamura et al., 2001). This method is called as the R-CIP-CSL2 method.

Fig. 6. Numerical solutions of the one-dimensional advection or conservation equation solved by various methods: (a) first-order upwind, (b) Lax-Wendroff, (c) CIP, (d) R-CIP, (e)

or 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

Hydrodynamics of a Droplet in Space 393

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

R-CIP-CSL2 without anti-diffusion, and (f) R-CIP-CSL2 with anti-diffusion.

solution in which the numerical diffusion is prevented effectively.

comparing with the above methods.

#### **3.1.3 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

$$\frac{\partial \sigma}{\partial t} = \frac{\partial}{\partial x} \left( u \frac{\partial \sigma}{\partial x} \right) . \tag{18}$$

Eq. (18) can be separated into two equations as

$$\frac{\partial \sigma}{\partial t} = -\frac{\partial f'}{\partial x}'\tag{19}$$

$$\mathbf{J}' = -\alpha \frac{\partial \sigma}{\partial \mathbf{x}'} \tag{20}$$

where *J*� indicates the anti-diffusion flux per unit area and per unit time. Using the finite difference method, we obtain

$$
\sigma\_{i+1/2}^{\*\*} = \sigma\_{i+1/2}^{\*} - (\hat{f}\_{i+1}^{\prime} - \hat{f}\_{i}^{\prime})\_{\prime} \tag{21}
$$

$$\mathbf{j}'\_{i} = -\mathbb{A}\_{i} \times \text{minmod}(\mathbf{S}\_{i-1\prime}\mathbf{S}\_{i\prime}\mathbf{S}\_{i+1})\_{\prime} \tag{22}$$

where <sup>ˆ</sup>*<sup>J</sup>* <sup>≡</sup> *<sup>J</sup>*� /(Δ*x*/Δ*t*) is the mass flux which has the same dimension of *<sup>σ</sup>*, *<sup>α</sup>*<sup>ˆ</sup> <sup>≡</sup> *<sup>α</sup>*/(Δ*x*2/Δ*t*) is the dimensionless diffusion coefficient, and *Si* ≡ *<sup>σ</sup>i*+1/2 − *<sup>σ</sup>i*−1/2. 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*, and *c* have the same sign. The value of the diffusion coefficient *α*ˆ is also important. Basically, we take *α*ˆ = −0.1 for the anti-diffusion. Here, it should be noted that *σ* takes the limited value as 0 ≤ *σ* ≤ *σ*m, where *σ*<sup>m</sup> is the initial value for inside of the droplet. The undershoot (*σ* < 0) or overshoot (*σ* > *σ*m) are physically incorrect solutions. To avoid that, we replace *<sup>α</sup>*ˆ*<sup>i</sup>* = 0.1 only when *<sup>σ</sup>i*−1/2 or *<sup>σ</sup>i*+1/2 are out of the appropriate range. We insert the anti-diffusion calculation after the CIP-CSL2 method is completed.

#### **3.1.4 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) 12 Will-be-set-by-IN-TECH

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,

. (18)

*<sup>∂</sup><sup>x</sup>* , (19)

*<sup>∂</sup><sup>x</sup>* , (20)

*<sup>i</sup>*), (21)

2007). In our model, we have an additional diffusion equation about *σ* as

Eq. (18) can be separated into two equations as

difference method, we obtain

where <sup>ˆ</sup>*<sup>J</sup>* <sup>≡</sup> *<sup>J</sup>*�

is completed.

**3.1.4 Test calculation**

*∂σ <sup>∂</sup><sup>t</sup>* <sup>=</sup> *<sup>∂</sup> ∂x α ∂σ ∂x* 

> *∂σ <sup>∂</sup><sup>t</sup>* <sup>=</sup> <sup>−</sup> *<sup>∂</sup>J*�

*J* � = −*α*

*σ*∗∗ *<sup>i</sup>*+1/2 = *σ*<sup>∗</sup>

ˆ*J* �

filled circles indicate the numerical results after 300 time steps.

*∂σ*

� *<sup>i</sup>*+<sup>1</sup> <sup>−</sup> <sup>ˆ</sup>*<sup>J</sup>* �

/(Δ*x*/Δ*t*) is the mass flux which has the same dimension of *<sup>σ</sup>*, *<sup>α</sup>*<sup>ˆ</sup> <sup>≡</sup> *<sup>α</sup>*/(Δ*x*2/Δ*t*)

*<sup>i</sup>* = −*α*ˆ*<sup>i</sup>* × minmod(*Si*−1, *Si*, *Si*+1), (22)

where *J*� indicates the anti-diffusion flux per unit area and per unit time. Using the finite

*<sup>i</sup>*+1/2 <sup>−</sup> (ˆ*<sup>J</sup>*

is the dimensionless diffusion coefficient, and *Si* ≡ *<sup>σ</sup>i*+1/2 − *<sup>σ</sup>i*−1/2. 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*, and *c* have the same sign. The value of the diffusion coefficient *α*ˆ is also important. Basically, we take *α*ˆ = −0.1 for the anti-diffusion. Here, it should be noted that *σ* takes the limited value as 0 ≤ *σ* ≤ *σ*m, where *σ*<sup>m</sup> is the initial value for inside of the droplet. The undershoot (*σ* < 0) or overshoot (*σ* > *σ*m) are physically incorrect solutions. To avoid that, we replace *<sup>α</sup>*ˆ*<sup>i</sup>* = 0.1 only when *<sup>σ</sup>i*−1/2 or *<sup>σ</sup>i*+1/2 are out of the appropriate range. We insert the anti-diffusion calculation after the CIP-CSL2 method

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

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)

**3.1.3 Anti-diffusion**

Fig. 6. Numerical solutions of the one-dimensional advection or conservation equation solved by various methods: (a) first-order upwind, (b) Lax-Wendroff, (c) CIP, (d) R-CIP, (e) R-CIP-CSL2 without anti-diffusion, and (f) R-CIP-CSL2 with anti-diffusion.

or 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.

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−<sup>3</sup> for the panel (b), respectively. We adopt *p*fm = 5000 dyn cm−<sup>2</sup> in this figure.

Hydrodynamics of a Droplet in Space 395

where *φ*¯ is the smoothed profile of *φ* (see section 3.2.4), and *Mi*<sup>+</sup><sup>1</sup> = *Mi* for *φ*¯*i*+<sup>1</sup> < *φ*¯*<sup>i</sup>* because the momentum flux does not increase when the molecular flow goes outward from inside of

*<sup>F</sup>*g,*xi* <sup>=</sup> <sup>−</sup> *Mi* <sup>−</sup> *Mi*<sup>+</sup><sup>1</sup>

from the left equation in Eq. (26). The momentum flux at upstream is *M*<sup>0</sup> = *p*fm. First, we solve Eq. (27) and obtain the spatial distribution of the molecular gas flow in all computational

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

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.

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

*<sup>i</sup>*) for *<sup>φ</sup>*¯*i*+<sup>1</sup> <sup>≥</sup> *<sup>φ</sup>*¯*i*, (27)

<sup>Δ</sup>*<sup>x</sup>* , (28)

Using the finite difference method to the right equation in Eq. (26), we obtain *Mi*<sup>+</sup><sup>1</sup> <sup>=</sup> *Mi* <sup>−</sup> *<sup>p</sup>*fm(*φ*¯*i*+<sup>1</sup> <sup>−</sup> *<sup>φ</sup>*¯

the droplet. Similarly, we obtain

flow.

**3.2.3 Surface tension**

domain. Then, we calculate the ram pressure by Eq. (28).

#### **3.2 Non-advection phase**

#### **3.2.1 C-CUP method**

Using the finite difference method to Eq. (9), we obtain (Yabe & Wang, 1991)

$$\frac{\vec{u}^{\ast \ast \ast} - \vec{u}^{\ast}}{\Delta t} = -\frac{\vec{\nabla} p^{\ast \ast}}{\rho^{\ast}} + \frac{\vec{Q}}{\rho^{\ast}} \quad \frac{p^{\ast \ast} - p^{\ast}}{\Delta t} = -\rho^{\ast} c\_{\text{s}}^{2} \vec{\nabla} \cdot \vec{u}^{\ast \ast} \, \tag{23}$$

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

$$\vec{\nabla} \cdot \left( \frac{\vec{\nabla} p^{\*\*}}{\rho^\*} \right) = \frac{p^{\*\*} - p^\*}{\rho^\* c\_s^2 \Delta t^2} + \frac{\vec{\nabla} \cdot \vec{u^\*}}{\Delta t} + \vec{\nabla} \cdot \left( \frac{\vec{Q}}{\rho^\*} \right) . \tag{24}$$

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).

#### **3.2.2 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*), where *s* is the collisional cross section of gas molecules and *n* is the number density of the nebula gas. Typically, we adopt *n* ≈ 1014 cm−<sup>3</sup> 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 *<sup>s</sup>* <sup>≈</sup> <sup>10</sup>−<sup>16</sup> cm−<sup>2</sup> 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

$$F\_{\mathbf{g},\mathbf{x}} = p\_{\mathbf{fm}} \delta(\mathbf{x} - \mathbf{x}\_{\mathbf{i}})\_{\prime} \tag{25}$$

where *x*<sup>i</sup> is the position of the droplet surface. This equation can be separated into two equations as

$$F\_{\mathbf{g},\mathbf{x}} = -\frac{\partial M}{\partial \mathbf{x}} \; \prime \quad \frac{\partial M}{\partial \mathbf{x}} = -p\_{\mathbf{f}\mathbf{m}} \delta(\mathbf{x} - \mathbf{x}\_{\mathbf{i}}) \, \tag{26}$$

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.

14 Will-be-set-by-IN-TECH

*<sup>ρ</sup>*<sup>∗</sup> , *<sup>p</sup>*∗∗ <sup>−</sup> *<sup>p</sup>*<sup>∗</sup>

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

<sup>s</sup>Δ*t*<sup>2</sup> <sup>+</sup>

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

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*), where *s* is the collisional cross section of gas molecules and *n* is the number density of the nebula gas. Typically, we adopt *n* ≈ 1014 cm−<sup>3</sup> 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 *<sup>s</sup>* <sup>≈</sup> <sup>10</sup>−<sup>16</sup> cm−<sup>2</sup> 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

Consider that the molecular gas flows for the positive direction of the *x*-axis. The *x*-component

where *x*<sup>i</sup> is the position of the droplet surface. This equation can be separated into two

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

*<sup>∂</sup><sup>x</sup>* , *<sup>∂</sup><sup>M</sup>*

∇ · *u*<sup>∗</sup>

<sup>Δ</sup>*<sup>t</sup>* <sup>=</sup> <sup>−</sup>*ρ*∗*<sup>c</sup>*

<sup>Δ</sup>*<sup>t</sup>* <sup>+</sup> ∇ · *Q ρ*∗

2 s∇ ·

*F*g,*<sup>x</sup>* = *p*fm*δ*(*x* − *x*i), (25)

*<sup>∂</sup><sup>x</sup>* <sup>=</sup> <sup>−</sup>*p*fm*δ*(*<sup>x</sup>* <sup>−</sup> *<sup>x</sup>*i), (26)

*u*∗∗, (23)

. (24)

Using the finite difference method to Eq. (9), we obtain (Yabe & Wang, 1991)

∇*p*∗∗ *<sup>ρ</sup>*<sup>∗</sup> <sup>+</sup> *<sup>Q</sup>*

<sup>=</sup> *<sup>p</sup>*∗∗ <sup>−</sup> *<sup>p</sup>*<sup>∗</sup> *ρ*∗*c*<sup>2</sup>

disturbed except of the direct collision with the droplet (free molecular flow).

*<sup>F</sup>*g,*<sup>x</sup>* <sup>=</sup> <sup>−</sup> *<sup>∂</sup><sup>M</sup>*

**3.2 Non-advection phase**

*u*∗∗ −*u*<sup>∗</sup>

the right equation. Then we obtain an equation

∇ · <sup>∇</sup> *<sup>p</sup>*∗∗ *ρ*∗

calculate *u*∗∗ by solving the left equation in Eq. (23).

**3.2.2 Ram pressure of free molecular flow**

of the ram pressure *F*g,*x* is given by

equations as

per unit area.

<sup>Δ</sup>*<sup>t</sup>* <sup>=</sup> <sup>−</sup>

**3.2.1 C-CUP method**

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−<sup>3</sup> for the panel (b), respectively. We adopt *p*fm = 5000 dyn cm−<sup>2</sup> in this figure.

Using the finite difference method to the right equation in Eq. (26), we obtain

$$M\_{i+1} = M\_i - p\_{\rm fm}(\vec{\phi}\_{i+1} - \vec{\phi}\_i) \quad \text{for } \vec{\phi}\_{i+1} \ge \vec{\phi}\_{i\prime} \tag{27}$$

where *φ*¯ is the smoothed profile of *φ* (see section 3.2.4), and *Mi*<sup>+</sup><sup>1</sup> = *Mi* for *φ*¯*i*+<sup>1</sup> < *φ*¯*<sup>i</sup>* because the momentum flux does not increase when the molecular flow goes outward from inside of the droplet. Similarly, we obtain

$$F\_{\mathbf{g}, \mathbf{x}\_i} = -\frac{M\_i - M\_{i+1}}{\Delta \mathbf{x}} ,\tag{28}$$

from the left equation in Eq. (26). The momentum flux at upstream is *M*<sup>0</sup> = *p*fm. First, we 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).

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.

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.

#### **3.2.3 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

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 = 104 dyn cm−2, *r*<sup>0</sup> = 500 ¯m, and *μ*<sup>d</sup> = 1.3 poise for

Hydrodynamics of a Droplet in Space 397

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

Fig. 9 shows 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

calculations.

droplet.

discontinuous interface to the transition region which has some width. According to them, the surface tension is expressed as

$$
\vec{F}\_{\text{s}} = \gamma\_{\text{s}} \kappa \vec{\nabla} \phi / [\phi]\_{\text{\textquotedblleft}} \tag{29}
$$

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

$$\kappa = -(\vec{\nabla} \cdot \vec{n}),\tag{30}$$

where

$$
\vec{n} = \vec{\nabla}\phi / |\vec{\nabla}\phi|.\tag{31}
$$

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).

#### **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

$$\bar{\phi} = \frac{1}{2}\phi\_{i,j,k} + \frac{1}{2}\frac{\phi\_{i,j,k} + \mathcal{C}\_1 \sum\_{L\_1}^6 \phi\_{L\_1} + \mathcal{C}\_2 \sum\_{L\_2}^{12} \phi\_{L\_2} + \mathcal{C}\_3 \sum\_{L\_3}^8 \phi\_{L\_3}}{1 + 6\mathcal{C}\_1 + 12\mathcal{C}\_2 + 8\mathcal{C}\_3},\tag{32}$$

where *L*1, *L*2, and *L*<sup>3</sup> indicate grid indexes of the nearest, second nearest, and third nearest from the grid point (*i*, *j*, *k*), for example, *L*<sup>1</sup> = (*i* + 1, *j*, *k*), *L*<sup>2</sup> = (*i* + 1, *j* + 1, *k*), *L*<sup>3</sup> = (*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, twelve for *L*2, and eight for *L*3, respectively. The coefficients are set as

$$\mathbb{C}\_1 = 1/(6 + 12/\sqrt{2} + 8/\sqrt{3}), \quad \mathbb{C}\_2 = \mathbb{C}\_1/\sqrt{2}, \quad \mathbb{C}\_3 = \mathbb{C}\_1/\sqrt{3}. \tag{33}$$

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.

#### **4. 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*- and *y*-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 16 Will-be-set-by-IN-TECH

discontinuous interface to the transition region which has some width. According to them,

where [*φ*] is the jump in color function at the interface between the droplet and the ambient

*κ* = −(∇ ·

The finite difference method of Eq. (31) is shown in (Brackbill et al., 1992). When we calculate

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

where *L*1, *L*2, and *L*<sup>3</sup> indicate grid indexes of the nearest, second nearest, and third nearest from the grid point (*i*, *j*, *k*), for example, *L*<sup>1</sup> = (*i* + 1, *j*, *k*), *L*<sup>2</sup> = (*i* + 1, *j* + 1, *k*), *L*<sup>3</sup> = (*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, twelve for *L*2, and eight for *L*3, respectively. The coefficients are

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

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*- and *y*-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

*<sup>L</sup>*<sup>1</sup> *<sup>φ</sup>L*<sup>1</sup> <sup>+</sup> *<sup>C</sup>*<sup>2</sup> <sup>∑</sup><sup>12</sup>

3), *C*<sup>2</sup> = *C*1/

1 + 6*C*<sup>1</sup> + 12*C*<sup>2</sup> + 8*C*<sup>3</sup>

√

*<sup>F</sup>*<sup>s</sup> <sup>=</sup> *<sup>γ</sup>*s*κ*<sup>∇</sup> *<sup>φ</sup>*/[*φ*], (29)

*<sup>n</sup>* <sup>=</sup> <sup>∇</sup> *<sup>φ</sup>*/|∇ *<sup>φ</sup>*|. (31)

*<sup>L</sup>*<sup>2</sup> *<sup>φ</sup>L*<sup>2</sup> <sup>+</sup> *<sup>C</sup>*<sup>3</sup> <sup>∑</sup><sup>8</sup>

2, *C*<sup>3</sup> = *C*1/

*<sup>L</sup>*<sup>3</sup> *φL*<sup>3</sup>

√

, (32)

3. (33)

*n*), (30)

gas. In our definition, we obtain [*φ*] = 1. The curvature is given by

the surface tension, we use the smoothed profile of *φ* (see section 3.2.4).

1 2 *φi*,*j*,*<sup>k</sup>* + *C*<sup>1</sup> ∑<sup>6</sup>

<sup>2</sup> <sup>+</sup> 8/<sup>√</sup>

the surface tension is expressed as

where

**3.2.4 Smoothing**

calculated by

set as

*<sup>φ</sup>*¯ <sup>=</sup> <sup>1</sup> 2 *φi*,*j*,*<sup>k</sup>* +

interface is kept unchanged.

**4.1 Vibrational motion**

*<sup>C</sup>*<sup>1</sup> <sup>=</sup> 1/(<sup>6</sup> <sup>+</sup> 12/<sup>√</sup>

**4. Deformation of droplet by gas flow**

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 = 104 dyn cm−2, *r*<sup>0</sup> = 500 ¯m, and *μ*<sup>d</sup> = 1.3 poise for calculations.

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.

Fig. 9 shows 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. 10. Same as Fig. 9 except of *μ*<sup>d</sup> = 100 poise.

angular velocity.

*<sup>N</sup>*, we obtain *<sup>N</sup>* <sup>≈</sup> *<sup>f</sup>πr*<sup>3</sup>

(Miura, Nakamoto & Doi, 2008)

*ω<sup>f</sup>* ≈

<sup>=</sup> <sup>140</sup> *<sup>f</sup>*

the angular momentum is maintained during melting.

15 *fπp*fm/8*r*<sup>2</sup>

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

Hydrodynamics of a Droplet in Space 399

The angular velocity *ω<sup>f</sup>* can be roughly estimated by *Iω<sup>f</sup>* ≈ *N*Δ*t*, where *I* 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

moment of inertia for a sphere with an uniform density *ρ*d, we obtain the angular velocity

104 dyn cm−<sup>2</sup>

Therefore, in the shock-wave heating model, the droplet should be rotating rapidly if most of

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 *μ*<sup>d</sup> =

sign of *N* would change after half-rotation. Substituting *I* = (8/15)*πr*<sup>5</sup>

<sup>0</sup>*ρ*<sup>d</sup>

0.01 1/2 *<sup>p</sup>*fm

<sup>0</sup> *p*fm. We can set Δ*t* ≈ *π*/*ω<sup>f</sup>* (a half-rotation period) because the

1/2 *<sup>r</sup>*<sup>0</sup>

1 mm

−<sup>1</sup>

<sup>0</sup>*ρ*d, which is the

rad s<sup>−</sup>1. (34)

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*. Each curve shows the calculation result for the different value of the ram pressure *p*fm. We use *r*<sup>0</sup> = 500 ¯m and *μ*<sup>d</sup> = 1.3 poise for calculations.

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 *P*vib = 2*π ρ*d*r*<sup>3</sup> <sup>0</sup>/8*γ*<sup>s</sup> ≈ 2.15 msec (Landau & Lifshitz, 1987).

#### **4.2 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.

#### **4.3 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, 18 Will-be-set-by-IN-TECH

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*. Each curve shows the calculation result for the different value of the ram pressure *p*fm. We use

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

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.

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,

<sup>0</sup>/8*γ*<sup>s</sup> ≈ 2.15 msec (Landau & Lifshitz, 1987).

*r*<sup>0</sup> = 500 ¯m and *μ*<sup>d</sup> = 1.3 poise for calculations.

*P*vib = 2*π*

 *ρ*d*r*<sup>3</sup>

**4.3 Effect of droplet rotation**

**4.2 Overdamping**

Fig. 10. Same as Fig. 9 except of *μ*<sup>d</sup> = 100 poise.

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 *ω<sup>f</sup>* can be roughly estimated by *Iω<sup>f</sup>* ≈ *N*Δ*t*, where *I* 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 *<sup>N</sup>*, we obtain *<sup>N</sup>* <sup>≈</sup> *<sup>f</sup>πr*<sup>3</sup> <sup>0</sup> *p*fm. We can set Δ*t* ≈ *π*/*ω<sup>f</sup>* (a half-rotation period) because the sign of *N* would change after half-rotation. Substituting *I* = (8/15)*πr*<sup>5</sup> <sup>0</sup>*ρ*d, which is the moment of inertia for a sphere with an uniform density *ρ*d, we obtain the angular velocity (Miura, Nakamoto & Doi, 2008)

$$\begin{split} \omega\_f &\approx \sqrt{15f \pi p\_{\rm fm}/8r\_0^2 \rho\_{\rm d}} \\ &= 140 \left( \frac{f}{0.01} \right)^{1/2} \left( \frac{p\_{\rm fm}}{10^4 \text{ dyn } \text{cm}^{-2}} \right)^{1/2} \left( \frac{r\_0}{1 \text{ mm}} \right)^{-1} \text{ rad s}^{-1}. \end{split} \tag{34}$$

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 *μ*<sup>d</sup> =

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

rotation axis.

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

Hydrodynamics of a Droplet in Space 401

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 *μ*<sup>d</sup> = 103 poise, *p*fm = 104 dyn cm−2, *ω* = 100 rad s−1, and *r*<sup>0</sup> = 1 mm.

103 poise, *p*fm = 104 dyn cm−2, *ω* = 100 rad s−1, and *r*<sup>0</sup> = 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.

#### **4.4 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*<sup>0</sup> + *r*1(*θ*), where *r*<sup>0</sup> is the unperturbed droplet radius and *r*<sup>1</sup> 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

$$\frac{r\_1(\theta)}{r\_0} = \frac{W\_\ell}{12} \left(\frac{19}{20} - R\right) P\_2(\cos \theta),\tag{35}$$

20 Will-be-set-by-IN-TECH

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 *μ*<sup>d</sup> = 103 poise, *p*fm = 104 dyn cm−2, *ω* = 100 rad s−1, and *r*<sup>0</sup> = 1 mm.

103 poise, *p*fm = 104 dyn cm−2, *ω* = 100 rad s−1, and *r*<sup>0</sup> = 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

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*<sup>0</sup> + *r*1(*θ*), where *r*<sup>0</sup> is the unperturbed droplet radius and *r*<sup>1</sup> 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

*r*1(*θ*) *r*0

<sup>=</sup> *We* 12

19 <sup>20</sup> <sup>−</sup> *<sup>R</sup>*

*P*2(cos *θ*), (35)

words, a rugby-ball-like shape.

**4.4 Origin of prolate chondrule**

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

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. 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 *We*. *R* = 19/20 (vertical dashed line) is a critical value for the droplet shape to be prolate

Hydrodynamics of a Droplet in Space 403

(*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

the break-up droplet. The droplet radius is *r*<sup>0</sup> = 2 cm, which corresponds to *We* = 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

Susa & Nakamoto (Susa & Nakamoto, 2002) suggested that the fragmentation of the droplets in high-velocity gas flow limits the sizes of chondrules (upper limit). They considered the balance between the surface tension and the inhomogeneity of the ram pressure acting on the droplet surface, and derived the maximum size of molten silicate dust particles above which the droplet would be destroyed by the ram pressure of the gas flow using an order of magnitude estimation. In their estimation, they adopted the experimental data in which

shows a condition in which the droplet will be fragmented by rapid rotation.

parent droplet breaks up to many smaller pieces.

where *We* (Weber number) is the ratio of the ram pressure of the gas flow to the surface tension of the droplet defined as

$$\mathcal{W}\_{\ell} = \frac{p\_{\rm fm} r\_0}{\gamma s},$$

*R* is the ratio of the centrifugal force to the ram pressure defined as

$$R = \frac{\rho\_{\rm d} r\_0^2 \omega^2}{p\_{\rm fm}},$$

*ω* is the angular velocity of the rotation, and *Pl*(cos *θ*) is Legendre polynomials. This solution is applicable under the assumption of *r*<sup>1</sup> � *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 *We* 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 *We* > 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 *We* > 0, but *b*/*a* = 1. As *We* 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 a wide range of *We* 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 = 104 dyn cm−<sup>2</sup> and the radius of chondrule of *r*<sup>0</sup> = 1 mm, we obtain *We* = 2.5 for *γ*<sup>s</sup> = 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 *We* and *R* depend on many factors, e.g., *p*fm, *r*0, and *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.

#### **5. Fragmentation**

#### **5.1 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 22 Will-be-set-by-IN-TECH

where *We* (Weber number) is the ratio of the ram pressure of the gas flow to the surface tension

*We* <sup>=</sup> *<sup>p</sup>*fm*r*<sup>0</sup> *γ*s

*<sup>R</sup>* <sup>=</sup> *<sup>ρ</sup>*d*r*<sup>2</sup>

*ω* is the angular velocity of the rotation, and *Pl*(cos *θ*) is Legendre polynomials. This solution is applicable under the assumption of *r*<sup>1</sup> � *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

Fig. 14 shows the droplet shape as functions of the Weber number *We* 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 *We* > 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 *We* > 0, but *b*/*a* = 1. As *We* 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

Let us consider the shape of chondrule expected from the shock-wave heating model. Adopting ram pressure of the gas flow of *p*fm = 104 dyn cm−<sup>2</sup> and the radius of chondrule of *r*<sup>0</sup> = 1 mm, we obtain *We* = 2.5 for *γ*<sup>s</sup> = 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 *We* and *R* depend on many factors, e.g., *p*fm, *r*0, and *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.

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

0*ω*<sup>2</sup> *p*fm

*R* is the ratio of the centrifugal force to the ram pressure defined as

*R* = 19/20 because the ram pressure balances with the centrifugal force.

, (36)

, (37)

of the droplet defined as

a wide range of *We* and *R*.

**5. Fragmentation**

**5.1 Direct fragmentation**

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 *We*. *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.

the break-up droplet. The droplet radius is *r*<sup>0</sup> = 2 cm, which corresponds to *We* = 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.

Susa & Nakamoto (Susa & Nakamoto, 2002) suggested that the fragmentation of the droplets in high-velocity gas flow limits the sizes of chondrules (upper limit). They considered the balance between the surface tension and the inhomogeneity of the ram pressure acting on the droplet surface, and derived the maximum size of molten silicate dust particles above which the droplet would be destroyed by the ram pressure of the gas flow using an order of magnitude estimation. In their estimation, they adopted the experimental data in which

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

Hydrodynamics of a Droplet in Space 405

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,

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 *We* = *p*fm*r*0/*γ*<sup>s</sup> ≈ 6. Since the ram pressure of the gas flow is typically *<sup>p</sup>*fm <sup>≈</sup> 103 <sup>−</sup>105 dyncm−2, we obtain the upper limit of chondrule sizes as *<sup>r</sup>*max <sup>≈</sup> 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 *We* < 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

by the eddy. We use *μ*<sup>d</sup> = 1.3 poise and *p*fm = 4000 dyn cm−2.

Eq. (38) provide us an insight of the cavitation criterion qualitatively.

**5.3 Comparison with chondrule properties**

be investigated in the future.

Fig. 15. Three-dimensional view of the fragmentation of molten droplet. We use *μ*<sup>d</sup> = 1 poise, *p*fm = 4000 dyn cm−2, and *r*<sup>0</sup> = 2 cm. The calculation was performed on a 100 × 100 × 100 grid.

the droplets suddenly exposed to the gas flow fragment for *We* >∼ 6 (Bronshten, 1983, p.96). This results into the fragmentation of droplet for *r*<sup>0</sup> >∼ 6 mm if we adopt our calculation conditions: *p*fm = 4000 dyn cm−<sup>2</sup> and *γ*<sup>s</sup> = 400 dyn cm−1. Our hydrodynamics simulations agree with the criterion for fragmentation.

#### **5.2 Fragmentation via cavitation**

Fig. 16 shows the internal pressure inside the droplet for various droplet sizes: *r*<sup>0</sup> = 3, 4, and 5 mm from panels (a) to (c). We use *μ*<sup>d</sup> = 1.3 poise and *p*fm = 4000 dyn cm−2. These droplets 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*<sup>2</sup> circ/*r*eddy ≈ *p*/*r*eddy, where *v*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*<sup>0</sup> from the Young-Laplace equation and *v*circ ≈ *v*max = 0.112*p*fm*r*0/*μ*<sup>d</sup> (Sekiya et al., 2003), we obtain

$$r\_{0, \text{cav}} \approx \left(\frac{2\gamma\_s \mu\_{\text{d}}^2}{0.112^2 \rho\_{\text{d}} p\_{\text{fm}}^2}\right)^{1/3}.\tag{38}$$

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*<sup>0</sup> = 4 mm or larger. 24 Will-be-set-by-IN-TECH

Fig. 15. Three-dimensional view of the fragmentation of molten droplet. We use

100 × 100 × 100 grid.

droplet.

agree with the criterion for fragmentation.

**5.2 Fragmentation via cavitation**

gradient force around eddies as *ρ*d*v*<sup>2</sup>

(Sekiya et al., 2003), we obtain

*μ*<sup>d</sup> = 1 poise, *p*fm = 4000 dyn cm−2, and *r*<sup>0</sup> = 2 cm. The calculation was performed on a

the droplets suddenly exposed to the gas flow fragment for *We* >∼ 6 (Bronshten, 1983, p.96). This results into the fragmentation of droplet for *r*<sup>0</sup> >∼ 6 mm if we adopt our calculation conditions: *p*fm = 4000 dyn cm−<sup>2</sup> and *γ*<sup>s</sup> = 400 dyn cm−1. Our hydrodynamics simulations

Fig. 16 shows the internal pressure inside the droplet for various droplet sizes: *r*<sup>0</sup> = 3, 4, and 5 mm from panels (a) to (c). We use *μ*<sup>d</sup> = 1.3 poise and *p*fm = 4000 dyn cm−2. These droplets 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

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

around the eddy, *r*eddy is the radius of the eddy, and *p* is the pressure inside the droplet. Substituting *p* = 2*γ*s/*r*<sup>0</sup> from the Young-Laplace equation and *v*circ ≈ *v*max = 0.112*p*fm*r*0/*μ*<sup>d</sup>

<sup>2</sup>*γ*s*μ*<sup>2</sup>

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*<sup>0</sup> = 4 mm or larger.

d 0.1122*ρ*<sup>d</sup> *p*<sup>2</sup>

fm

1/3

*r*0,cav ≈

circ/*r*eddy ≈ *p*/*r*eddy, where *v*circ is the fluid velocity

. (38)

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 by the eddy. We use *μ*<sup>d</sup> = 1.3 poise and *p*fm = 4000 dyn cm−2.

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.

#### **5.3 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 *We* = *p*fm*r*0/*γ*<sup>s</sup> ≈ 6. Since the ram pressure of the gas flow is typically *<sup>p</sup>*fm <sup>≈</sup> 103 <sup>−</sup>105 dyncm−2, we obtain the upper limit of chondrule sizes as *<sup>r</sup>*max <sup>≈</sup> 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 *We* < 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.

URL: *http://www.sciencedirect.com/science/article/B6V66-48C8H7W-BK/2/bdd79a0a3820*

Brackbill, J. U., Kothe, D. B. & Zemach, C. (1992). A continuum method for modeling surface

Hydrodynamics of a Droplet in Space 407

Chandrasekhar, S. (1965). The stability of a rotating liquid drop, *Proceedings of the Royal Society*

Ciesla, F. J. & Hood, L. L. (2002). The nebula shock wave model for chondrule formation:

Ciesla, F. J., Hood, L. L. & Weidenschilling, S. J. (2004). Evaluating planetesimal bow shocks

Desch, S. J. & Jr., H. C. C. (2002). A model of the thermal processing of particles in solar

Doi, M. (2011). *Formation of cosmic spherule: chemical analysis and theory for shapes, compositions,*

Fredriksson, K. & Ringwood, A. (1963). Origin of meteoritic chondrules, *Geochimica et*

Gooding, J. L. & Keil, K. (1981). Relative abundances of chondrule primary textural types

Harold C. Connolly, J. & Hewins, R. H. (1995). Chondrules as products of dust collisions with

Hayashi, C. K., Nakazawa, K. & Nakagawa, Y. (1985). *Formation of the solar system*, Protostars

Hewins, R. H. & Radomsky, P. M. (1990). Temperature conditions for chondrule formation,

Hollenbach, D. & McKee, C. F. (1979). Molecular formation and infrared emission in fast

Hood, L. L. (1998). Thermal processing of chondrule precursors in planetesimal bow shocks,

Hood, L. L. & Horanyi, M. (1991). Gas dynamic heating of chondrule precursor grains in the

Hood, L. L. & Horanyi, M. (1993). The nebular shock wave model for chondrule formation -

Hughes, D. W. (1978). A disaggregation and thin section analysis of the size and mass

URL: *http://www.sciencedirect.com/science/article/pii/0012821X78901139* Iida, A., Nakamoto, T., Susa, H. & Nakagawa, Y. (2001). A shock heating model for chondrule

distribution of the chondrules in the bjurbîŠTe and chainpur meteorites, ˇ *Earth and*

nebula shocks: Application to the cooling rates of chondrules, *Meteorit. Planet. Sci.*

URL: *http://www.sciencedirect.com/science/article/B6V66-48C8FST-5/2/baf33fee9f7f0a8d0*

in ordinary chondrites and their bearing on conditions of chondrule formation,

totally molten droplets within a dust-rich nebular environment: An experimental

URL: *http://www.sciencedirect.com/science/article/pii/002199919290240Y*

Shock processing in a particle-gas suspension, *Icarus* 158: 281–293.

as sites for chondrule formation, *Meteorit. Planet. Sci.* 39: 1809–1821.

tension, *Journal of Computational Physics* 100(2): 335 – 354.

Bronshten, V. A. (1983). *Physics of Meteoric Phenomena*, Dordrecht: Reidel.

*and textures*, PhD thesis, Tokyo Institute of Technology.

investigation, *Geochimica et Cosmochimica Acta* 59: 3231–3246.

and Planets II, Univ. of Arizona Press, Tucson, pp. 1100–1153.

interstellar shocks. i. physical proceses, *Astrophys. J.* 41: 555–592.

*Cosmochimica Acta* 27(6): 639 – 641.

*of London. Ser. A, Mathematical and Physical Sciences* 286: 1–26.

*afc4d06ac02bdc7cfaa7*

37: 183–207.

*a3ef0a92ac38cf7*

*Meteoritics* 16: 17–43.

*Meteoritics* 25: 309–318.

*Meteorit. Planet. Sci.* 33: 97–108.

solar nebula, *Icarus* 93: 259–269.

*Planetary Science Letters* 38(2): 391 – 400.

one-dimensional calculations, *Icarus* 106: 179–189.

formation in a protoplanetary disk, *Icarus* 153: 430–450.

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.
