spreading with solidification microstructure formation in plasma spraying

International Journal of Heat and Mass Transfer 150 (2020) 119267

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/hmt

Numerical modeling of YSZ droplet impact/spreading with solidification microstructure formation in plasma spraying Mingguang Shen a, Ben Q. Li b,∗, Yu Bai c a

State Key Laboratory for Manufacturing Systems Engineering, Xi’an Jiaotong University, Xi’an, Shaanxi 710049, PR China Department of Mechanical Engineering, University of Michigan, Dearborn, MI 48128, United States c State Key Laboratory for Mechanical Behavior of Materials, Xi’an Jiaotong University, Xi’an, Shaanxi 710049, PR China b

a r t i c l e

i n f o

Article history: Received 8 August 2019 Revised 1 December 2019 Accepted 21 December 2019

Keywords: Phase-field Multiphase flow Supersonic impact Polycrystalline growth

a b s t r a c t A comprehensive numerical model has been developed, which is capable of representing the impact and spreading of a high-velocity yttria-stabilized zirconia (YSZ) molten droplet with microstructure formation in plasma thermal spraying processes. The numerical model entails the explicit finite difference solution of the Navier-Stokes equations and the energy balance equations, coupled with the Cahn-Hilliard equation to track the liquid-gas two phase interface and with a phase field model for solidification microstructure formation involving the polycrystalline growth. Numerical model procedures are given. The model, after being checked for mesh-independence, was applied to study the spreading and solidification of a high-temperature high-velocity YSZ droplet impinging upon a preheated substrate surface under the conditions typical of supersonic plasma thermal spraying processes. Extensive numerical simulations were carried out and results reveal that the solidification microstructure formed in the spreading droplet consists primarily of columnar grains. Computed results are also compared with thermal spray experiments, and gratifying agreement is obtained. © 2019 Published by Elsevier Ltd.

1. Introduction Plasma thermal spraying offers an important means to fabricate metallic or ceramic thermal barrier coatings covering the blade surfaces of aircraft engines and gas turbines to prevent hot corrosion and reduce substrate temperature, thereby raising serving life of these components and bettering operating performance. During thermal spraying, a coating is formed by the pileup of a large number of solidified particles, which, having been accelerated to a high velocity and melted by the supersonic plasma gas jet, solidifies quickly after impinging and spreading onto a preheated surface. The process of coating formation involves rather complex coupled physical phenomena, including multiphase fluid flow, heat transfer and solidification microstructure evolution. From practical perspective, the mechanical and thermal performance of the coating under harsh combustion environment depends to a large extent on the final microstructures formed in the spread-out solidified droplets and yet the solidification microstructure formation is strongly affected by the fluid flow and heat transfer during spreading of a droplet along a solid surface.

∗

Corresponding author. E-mail address: [email protected] (B.Q. Li).

https://doi.org/10.1016/j.ijheatmasstransfer.2019.119267 0017-9310/© 2019 Published by Elsevier Ltd.

Because of the significant engineering and industrial importance, considerable numerical and experimental research work has been conducted, which have helped to gain good physical insight into the fundamental physics governing splat formation under plasma spraying conditions. With a commercial CFD Software Flow-3 D, Bobzin et al. [1] investigated YSZ impact on flat and rough surfaces under plasma spraying conditions, with the volume of fluid method to track the free surface between a liquid and a gas. The impacting velocities were quite high, but still subsonic. Tabbara and Gu [2] modelled the impact of partially molten zirconia droplets, and described in detail the freezing induced splashing, finding that it was the deceleration of the solidified that triggered splashing. A two-dimensional axisymmetric model based on a finite element method was developed by Oukach et al. [3] to simulate droplet spreading and solidification on a stainless steel substrate under plasma spraying conditions. They also derived a correlation for the spread factor with respect to the Reynolds number. Using a coupled Eulerian and Lagrangian method, Zhu et al. [4] examined both fully and partially molten YSZ impact under practical conditions, asserting that the smaller the solid core, the larger the flattening ratio. For hollow zirconia droplets, Safaei et al. [5], using a volume of fluid method for compressible flows, studied the effect of the trapped air on the impact dynamics, concluding that at velocities exceeding 100 m/s the trapped air played

2

M. Shen, B.Q. Li and Y. Bai / International Journal of Heat and Mass Transfer 150 (2020) 119267

a dominant role with the droplet splashing explosively. All these studies take the classical equilibrium solidification model, that is, the interface is defined at the location where the equilibrium freezing point resides. But such a model doesn’t take into account undercooling, which is of great importance to rapid solidification in splat quenching. Combining a volume of fluid method and a non-equilibrium rapid solidification model, Shukla et al. [6] probed into the flattening of a molten zirconia droplet with nucleation undercooling, observing that a larger nucleation undercooling produces a relatively lower interface temperature and a higher interface velocity. As discussed above, though significant amount of work has been carried out on the transport phenomena associated with droplet impingement, to the best knowledge of the authors, little, if any at all, appears to have been conducted to understand the fundamentals governing the coupled droplet spreading and solidification microstructure formation associated with a molten YSZ droplet in plasma spray processes. This paper aims to address this issue by presenting a comprehensive numerical model of transport phenomena and solidification microstructure formation associated with a molten YSZ droplet during plasma spraying. The model entails the numerical solution of the phase field equation of microstructure formation along with the multiphase flow and heat transfer equations enhanced with the Cahn-Hilliard equation for droplet-gas interface tracking. The paper is organized as follows. First, the mathematical model, including boundary conditions and numerical schemes, is explained. Then the model is applied to study the impingement and spreading of a high-velocity molten YSZ droplet on a preheated solid surface while undergoing solidification with columnar microstructure formation. The model predictions were compared with the plasma spray experiments, showing gratifying agreement. 2. Mathematical statement Let us consider a molten sphere typical in a plasma coating process. The droplet impinges on a planar substrate with a specified velocity and a given initial temperature, as shown in Fig. 1. Upon impacting, the droplet spreads outward along the surface of the substrate. Meantime, the temperature of the droplet decreases rapidly and the droplet solidifies while spreading, as the surface of the substrate is at a temperature far below the melting point of the droplet. This creates a thermal gradient inside the droplet, which favors the epitaxial growth of columnar microstructures upon solidification. The equations governing the complex fluid

flow, heat transfer and solidification microstructure formation associated with the spreading process are described below. 2.1. Mass, momentum and energy equations With incompressible Newtonian fluids, the equations governing mass, momentum and energy conservations are obtained as follows for the gas-liquid flow.

∇ ·u=0 ∂u ∇p ∇ · σ G∇ c + u · ∇u = − + +g+ +S ∂t ρ (c ) ρ (c ) ρ (c )   ρ L ∂φ ∂T ρ (c )cP (c ) + u · ∇ T = ∇ · [k(c )∇ T ] + l l ∂t 2 ∂t

(1) (2)

(3)

where u is velocity, t time, p mechanical pressure, σ = μ(c )[∇ u + (∇ u )T ] Newtonian stress tensor, and g local gravity acceleration. The fourth term on the right hand side of Eq. (2) represents the liquid-gas surface tension effect, with G being the chemical potential. While many numerical approaches may be used, the phase field model is employed in this paper to track the liquid-gas interface, which permits the use of fixed grids to treat a moving free surface as encountered in droplet spreading. This use of the phase field modeling introduces a phase field parameter, c, which renders the thermophysical properties as a function of c (e.g., ρ (c ) = (ρg − ρl )c + ρl where g denotes gas and l liquid). The term S in Eq. (2) is a source term to account for solid formation in the droplet, and is related to solidification microstructure formation. The evolution of phase field parameter c and microstructure formation is given in the following two sections. Eq. (3) shows energy conservation, with some parameters recast as functions of the order parameter c, such as density ρ (c), specific heat cP (c), and thermal conductivity k(c). The second term on the right hand side of Eq. (3) represents latent heat released, with ρ l and Ll denoting liquid density and latent heat per unit mass. Note that ∂ φ /∂ t is updated only where c ≤ 0.5 and will be discussed in detail in the following section. For the portion where c > 0.5, ∂ φ /∂ t = 0. 2.2. Phase field equation for an evolving gas-liquid interface In this paper, the motion of the gas-liquid interface is described by the Cahn-Hilliard (CH) equation,

ct + u · ∇ c − M∇ 2 G = 0

(4)

where M, the strength of diffusion, modifies the value of c through diffusion, and G = δ f /δ c is the chemical potential, with f = 12 ξ γ α|∇ c|2 + ξ −1 γ α · 14 c2 (1 − c )2 being free energy density per unit volume, where ξ stands for a measure of interfacial thick√ ness, γ surface tension, and α a constant equal to 6 2. The equation governs the evolution of the order parameter c, which is a measure of phase, or a phase indicator. Generally, two different integers will be assigned to c, each denoting a distinct phase. In this paper c = 1 is to denote a gas and c = 0 a liquid, with the value in between indicating the diffuse interface. 2.3. Solidification microstructure formation

Fig. 1. Schematic of the problem with boundary conditions defined.

To emulate polycrystalline growth with a phase field model, the key lies in describing a large number of grains with varied crystalline orientations. Recently, with the basic principles of phase field theory, the community has developed three major types of models: the so called multiphase or multiorder parameter models, and the orientational field model (OPF) that was championed by

M. Shen, B.Q. Li and Y. Bai / International Journal of Heat and Mass Transfer 150 (2020) 119267

3

Kobayashi and coworkers [7–9]. The starting point for an orientational phase field for a pure polycrystalline material is a free energy functional in terms of θ , φ , u, with θ representing the local orientation of crystals with respect to a fixed coordinate, φ the solid-liquid order parameter, and u = (T − Tm )/(L/c p ) the reduced temperature. Not that φ is another order parameter with φ = −1 denoting a liquid, or more broadly, a fluid, and φ = 1 a solid. The free energy functional is given by

 F=



dV f (φ , u ) +

εφ2 (ϕ ) 2

|∇φ| + sh(φ )|∇θ | + 2

εθ2 2



h(φ )|∇θ |

2

where f(φ , u) is the bulk free energy density, and the second term accounts for the energy penalty within the interface, with ɛφ (ϕ ) being the anisotropic gradient energy coefficient. It should be noted that the third and fourth terms are newly added, the former to ensure stable existence of grain boundaries and the latter possible boundary motion, and that h(φ ) is so chosen that it activates only in solid and vanishes in liquid. s and ɛθ specify the strength of the coupling between φ and ∇θ . ϕ is the inclination angle of solid-liquid interface (n = ∇ φ /|∇ φ|) with respect to the x axis, defined as ϕ = tan−1 (∂ φ /∂ y )/(∂ φ /∂ x ). With the bulk free energy density properly chosen [10, 11], the kinetic equations for φ and θ are given by

τφ ( ϕ )

∂φ = [φ − λu(1 − φ 2 )](1 − φ 2 ) + ∇ · ( · ∇φ ) ∂t 2 − 2s(1 + φ )|∇θ | − εθ2 (1 + φ )|∇θ |

    ∂θ s 2 2 P τθ ( 1 + φ ) = ∇ · (1 + φ ) + εθ ∇θ ∂t |∇θ | 2

(5)

(6)

where τ φ (ϕ ) and τ θ are interface attachment kinetic times, λ specifies the coupling between phase field φ and diffusion fields, and P (εθ |∇θ | ) = 1 + (μ/εθ − 1 )e−βεθ |∇θ| is the inverse mobility function. μ governs grain rotation and β grain boundary mobility.  is a second order tensor, with its components in two dimensions given by

⎡

εφ2

⎢ =⎣ ∂ε εφ φ ∂ϕ

−εφ

⎤ ∂εφ ∂ϕ ⎥ ⎦

εφ (ϕ ) = ε˜φ {1 + ε4 cos [4(ϕ − θ )]} ɛ4 is a constant, representing the strength of anisotropy. Likewise, τ φ can be made anisotropic also, taking the form

τφ (ϕ ) = τ˜φ {1 + ε4 cos [4(ϕ − θ )]}2 The solidification microstructure having been obtained, the source term S in the momentum balance equation can then be calculated. The S term in essence governs the flow in a mushy region, where solid and liquid coexist. S takes the form of

(1 − Fl )2 Fl 3 + b

u

motion as expected. While b is a quite small number to avoid division by zero. In this paper, d and b are chosen to be on the order of 108 ~ 109 and 10−3 , respectively. 2.4. Boundary conditions The solution of the above equations is subject to physical constraints, which are described by the boundary conditions. For the problem under consideration, the boundary conditions are given in Fig. 1. Note that symmetry is used to save computational resources. Some points should be made here. First, for computation of heat transfer, the two-phase flow region and the substrate is combined with all boundaries adiabatic; second, for computation of fluid motion and phase transformation, the substrate is excluded; third, ψ contains all field variables, except for the normal or all velocity components at the axis of symmetry and at the wall, respectively, and for the order parameter c at the substrate surface which satisfies ξ γ α n · ∇ c + fw (c ) = 0, where n is the unit normal pointing to the wall. fw (c ) = [−γ cos θS (4c3 − 6c2 + 1) + γw1 + γw2 ]/2 is the wall free energy density. fw (c = 0 ) gives the liquid–solid interfacial tension γ w1 and fw (c = 1 ) produces the gas-solid interfacial tension γ w2 , the two determining the static contact angle θ S through Young’s equation γw2 − γw1 = γ cos θS . 2.5. Numerical procedures

εφ2

If anisotropy is not considered, then the non-diagonal terms will vanish and ɛφ will be independent of ϕ , simply acting as a constant, and so will τ φ . If anisotropy does enter consideration, it takes the usual form with square symmetry

S = −d

Fig. 2. A modified staggered grid (schematic) for discretization.

(7)

where d is a rather large number to sweep out all other momentum sources when the liquid fraction Fl = (1 − φ )/2, expressing the volume percentage of the fluid component, is coming close to zero. Specifically, when a liquid is near the completion of solidification, its velocity is made to vanish. If Fl is unity, namely, when no solidification occurs, this term drops out, exerting no effect on fluid

The set of equations are explicitly discretized using the finite difference method, with upwind scheme to deal with convective terms. The flow equation is solved via the popular projection method. Besides, no special numerical techniques are involved. A modified staggered grid is shown in Fig. 2, where the square ABCD comprises of the computational domain. Note that pressure is stored at the cell center, marked by the filled circle in black, while velocity and other field variables, like temperature, are stored at the vertices of the cell, marked by the unfilled circle. Moreover, a series of ghost nodes, marked in brown and red, render boundary nodes inner ones. No modification is thus to be made to the discretized equations when considering boundary nodes. The no slip condition on the bottom is applied at j = 1 while row j = 0 represents the nodes in the substrate. The thermal contact resistance is thus applied between row j = 1 and j = 0. Some remarks should be made about the assignment of the two order parameters, namely c and φ . For better explanation, a complete, not halved, droplet is employed. Initially, no points with φ = 1 exist, except some nuclei randomly seeded at the substrate surface. As solidification progresses, the pair of c = 0 and φ = 1 emerges and gradually fills the interior of the droplet, as shown in Fig. 3. The time step is chosen empirically given the trade-off between the stability requirement and computational efficiency. Throughout

4

M. Shen, B.Q. Li and Y. Bai / International Journal of Heat and Mass Transfer 150 (2020) 119267

Fig. 3. Order parameter distribution (a) initial (b) intermediate.

the paper, the time step t is on the order of 10−12 s or smaller, which is small enough to satisfy the stability constraint. During a computational time interval, computation starts from the evolution of φ and θ , which is confined to the region with c ≤ 0.5, then proceeds to heat transfer and fluid motion, thus completing one loop. It will, however, not stop until the time duration set is reached. In addition, parallel programing based on OpenMP is applied to accelerate computation. 3. Results and discussion This section begins with simulation of a YSZ impact without solidification, followed by the analytical analysis of heat transfer in splat quenching and by the YSZ impact with microstructure formation. Finally, convection effect on solidification morphology is presented. The thermophysical properties are listed in Table 1. A fixed grid of uniform spatial step is adopted. ξ = ε˜φ = s = 1.25x, λ = a1 ε˜φ /d0 , and τφ = a2 λε˜φ2 /D, where a1 = 0.8839 and a2 = 0.6267. d0 is the thermal capillary length of YSZ, and D the thermal diffusivity. ε4 = 0.05, τθ = 20τφ , εθ = 10−3 ε˜φ , μ = 103 , and β = 105 . The choice of parameters starts from ξ and ε˜φ , both concerning interfacial thickness, the first of which, measuring liquid-gas interfacial thickness, often arises in the dimensionless constant Cn = ξ /Lc , where Lc is a characteristic length, for instance, droplet diameter. For cases without contact line movement, CH model convergence can be attained with Cn = O(10−2 ) in practice [15]; but for those with contact line motion, the picture becomes murkier, for the contact line mobility M matters now, and the debate is still open as to how to relate the two to realize the sharp interface limit [16]. Nonetheless, some empirical values of Cn, typically on the order of 0.01, along with a mobility roughly of interface thickness squared, produced reasonable results [17–20]. In the present paper, Cn = 0.0104 unless otherwise stated, and the mobility is so chosen that it satisfies a criteria proposed by Yue √ et al. [15], namely, Cn ≤ 4 Mμe /Lc , where μe = μg μl is the effective viscosity, with g denoting gas and l liquid. But, this is not the whole story. ε˜φ , defining solid-liquid interfacial width, should be as small as possible, or nearly vanished, to approach the socalled sharp interface limit for the OPF model. This is numerically intractable, however. Fortunately, the community was relieved by

Table 1 Thermophysical parameters used in simulations.

Density (kg/m3 ) Specific heat (J/kg · K) Thermal conductivity (W/m · K) Viscosity (mPa · s) Surface tension (N/m) Latent heat of fusion (kJ/kg) Melting point (K)

YSZ

Air

Substrate

5890 [12] 713 [12] 2.32 [13] 27.8 [12] 0.43 [14] 707 [13] 2923

1.18 1006 0.0263 0.0185

8400 575 18.8

the renowned thin interface limit by Karma and Rappel [10], which allowed for much thicker ε˜φ while maintaining the sharp interface limit. Besides, extensive benchmark tests have shown that simulation results become virtually independent of ε˜φ if it is about an order of magnitude smaller than the smallest radius of curvature of interest [21, 22]. In plasma spraying, solidification front is mostly planar, indicating a quite large radius of curvature. In turn, ε˜φ can be chosen much larger, and for convenience’s sake, is made equal to ξ . With the two parameters reasonably defined and a grid spacing no larger than ξ or ε˜φ [18, 21, 23], the following results are thought to be model-convergent and grid-independent. 3.1. Model convergence and grid independence test In the previous section, it is presumed that with Cn = O(10−2 ) and M∝Cn2 and a proper grid spacing to resolve the thin interface, model convergence and mesh independence are ensured. To validate this hypothesis, a supersonic molten YSZ impact is to be carried out on three different meshes, with Cn = 0.0104, Cn = 0.0078 and Cn = 0.0052, respectively, the first corresponding to λ = 145.36, the second λ = 109.02 and the last λ = 72.68. Our approach is to reduce mesh size and monitor some key factors, such as morphology of solidification microstructure, the maximum spread factor, and the substrate surface temperature at the impacting point. No further mesh refinement is made, since the total nodes on the mesh with Cn = 0.0052 reach up to 0.6 million, taking over 14 h to finish the computation on a common compute node, equipped with 24 processors, at the HPC Platform, Xi’an Jiaotong University. Besides, finer mesh means smaller time step and larger round-off errors, thus not necessarily producing better results. The numerical configurations are as follows. A YSZ droplet, 4.8 μm in diameter, is made to impinge at 518 m/s onto a solid surface that is preheated to 423 K. The initial temperature is 30 0 0 K and contact angle 120◦ . Interfacial heat transfer coefficient is set to 108 W/m2 · K. Vardelle et al. [24] asserted that when zirconia particles impinged on smooth steel substrates at 423 K, the intimate contact, between the disk-shaped splat and the substrate, would lead to a cooling rate greater than 400 K/μs, and that the interfacial thermal contact resistance varied between 10−7 and 10−8 m2 · K/W. The thermophysical properties of zirconia differs only slightly from that of YSZ so that we choose to use the smaller one given a higher impacting velocity 518 m/s compared with that studied in the paper, which is around 200 m/s. A higher impacting velocity may lead to more intimate contact and hence lower thermal contact resistance. The computational domain has a uniform spatial step of x = 4 × 10−8 m, x = 3 × 10−8 m and x = 2 × 10−8 m, corresponding to a coarse grid of 751 × 201, a medium grid of 1001 × 271 and to a fine grid of 1501 × 401, respectively. Besides, a number of supercritical nuclei of different orientations are randomly seeded at the surface. The supercritical nuclei are seeded randomly with an

M. Shen, B.Q. Li and Y. Bai / International Journal of Heat and Mass Transfer 150 (2020) 119267

5

Fig. 5. Substrate temperature at the impacting point with respect to time.

Fig. 4. Morphology of solidification microstructure at 0.3 μs, (a) Cn = 0.0052 (b) Cn = 0.0078 (c) Cn = 0.0104.

average size of R = 0.001 × x, where x is the grid spacing. These nuclei are also kept a few spatial steps apart to avoid the effect of mesh size as well as the overlapping of any two seeds. The number density is 65/750 in the coarse grid, 58/10 0 0 in the medium, and 65/1500 in the fine, where the denominator means the number of seeds and the nominator equals the horizontal length of the computational domain in dimensionless form. The morphology of solidification microstructure at 0.3 μs, the instant that is very close to the end of solidification, is presented in Fig. 4. Moreover, the phase field mobility M = 2.2 × 10−13 m3 · s/kg is calculated in response to Cn = 0.0104 and is used throughout the three cases. As can be seen in Fig. 4, reasonable agreement is reached among the three, all showing columnar grains growing perpendicular to the substrate surface. A slight difference, however, exists in the solidification thickness which is due to the discrepancy in the heat extraction rate caused partially by the interface thickness and by the flattening degree of the droplet. For row (c) in Fig. 4, the spread factor is about 6.86, relatively larger than that in row (a), 6.46, so that the steady state solidification thickness in row (c) is probably less than in row (a). Solidification begins a little earlier when Cn = 0.0052 because of a thinner interface, as demonstrated in Fig. 5, which shows the temperature history of the impacting point on the substrate surface. It is obvious that substrate temperature at the impacting point rises quite rapidly from the initial 423 K to about 1100 K within 0.1 μs. The steady state temperatures corresponding to the different interface thicknesses deviate from one another within 5%. The substrate surface temperature will keep constant while solidification is undergoing, as seen from 0.1 μs to 0.3 μs. Once the droplet freezes completely, it will drop accordingly. In addition, the last, but not the least, important criterion is the maximum spread

factor, which measures the flattening degree, turns out to be between 6.46 and 6.88 according to in Fig. 4, displaying a high level of agreement. To sum up, model convergence and mesh independence could be obtained in the limit of vanished interfacial thickness if the phase field mobility M stays unchanged. Here we mean both solidliquid interface and liquid-gas interface. Their importance is embodied in the dimensionless constants Cn and λ. It depends on practical purposes and the trade-off between desired numerical accuracy and available computational resources to select economical values for Cn and λ. As shown above, even with Cn = 0.0104 and λ = 145.36 it is still possible to get a fairly satisfactory result. Therefore, in the following subsections, we will use Cn = 0.0104 and λ = 145.36 wherever appropriate except when the very details of fluid field or solidification dynamics need to be unveiled. 3.2. A supersonic impact without solidification For one thing, a supersonic impact without phase change is examined to anatomy the fluid dynamics involved. With a diameter of 4.8 μm, a YSZ droplet impinges onto a solid surface perpendicularly at 518 m/s. The mesh size is 751 × 201, corresponding to a physical domain of 30 μm × 8 μm. The mesh is uniform and the Cahn number Cn = 0.0104, denoting the ratio of the capillary length ξ to droplet diameter. The phase field mobility M = 2.2 × 10−13 m3 · s/kg. The characteristic Reynolds number is Re = ρl V Dl /μl = 526.79, where Dl is the diameter of the droplet and Weber number is W e = ρl V 2 Dl /γ = 1.76 × 104 . Initially, the droplet perches quite close to the substrate. Besides, the contact angle is set to 120◦ and the temperature of the droplet is assume to be 30 0 0 K. As demonstrated in Fig. 6, the droplet, pushed by large pressure gradient along the solid surface, spreads rather swiftly, reaching the maximum spread around 0.08 μs. In the spreading process, the kinetic energy is mostly transformed into surface energy and the rest is consumed through, on the one hand, the violent interaction between the droplet and surrounding gas, and the other, its inner viscous dissipation. Upon reaching the maximum spread, the droplet is motionless temporarily so that the change of the horizontal velocity of the droplet is not monotonic, it first being accelerated by the pressure gradient and then being slowed down because primarily of the capillary force. Just like the external flow around a cylinder, the fluid velocity first increases due to a

6

M. Shen, B.Q. Li and Y. Bai / International Journal of Heat and Mass Transfer 150 (2020) 119267

3.3. Analytical analysis of solidification rate in splat quenching The case above anatomies fluid dynamics in high speed impact. This subsection is, however, devoted to the analytical analysis of heat transfer in splat quenching. In plasma spraying, spreading time is much smaller than solidification time, and splat thickness comes much thinner than its diameter. It is therefore appropriate to assume a one-dimensional heat conduction model to simplify the analysis. The schematic of the problem is shown below in Fig. 7, with a couple of temperature nodes identified. s represents solidification thickness, k thermal conductivity of the solidified, hc interfacial heat transfer coefficient, H substrate height, δ splat thickness, and ksub substrate thermal conductivity. TP , Tm , Ti , TB , Tss , and Tsub , stand for the temperature of splat top, equilibrium freezing temperature, solid-liquid interface, splat bottom, substrate surface, and substrate bottom, respectively. Note that the invisible gap between splat bottom and substrate surface  is artificially enlarged for better illustration and that Rk is the thermal resistance of kinetics. To find the solidification rate, or interface velocity, we proceed from the well-known Gibbs-Thomson condition and the energy balance equation at the interface, without considering the curvature effect.

Vi = μi (Tm − Ti ) 



(8)

ρ LVi = qs − ql = k

∂T ∂T | − k |i+ ∂ y i− ∂y

(9)

i− and i+ are signaling that the temperature gradient is sought from the solid and liquid side, respectively. Considering the flow of heat flux, we have

Ti − Tsub



qs = Fig. 6. A YSZ droplet impact at 518 m/s without phase change.

favorable pressure gradient, hitting a maximum, and then decreases due to an unfavorable one. But the pressure in the present case will play out its role after accelerating the droplet while surface tension force slows the velocity in the late stage of the spreading process. On the other hand, once the droplet extends to the limit, retracting begins, which is prominent at 2 μs while the surface tension force in this stage is rounding the tip and pulling the droplet backwards. The contact angle varies during spreading and retracting, but it will eventually settle at 120◦ . Note that the retracting is still ongoing. Hence the state at 2 μs is not final.



ql =

s k

1 hc

+

+

H ksub

TP − Tm δ −s k

TP is the temperature of splat top and can be seen unchanged during early stages of solidification, because the Biot number Bi = hc δ /k is far bigger than unity. It takes the value of 30 0 0 K, quite close to the equilibrium freezing temperature of YSZ. Besides, the   denominators in qs and ql are on the same order of magnitude 

while the nominators differ dozens of times so that ql can be safely neglected. Upon substitution, Eq. (9) becomes

Vi =

Ti − Tsub

ρL

s

k

+

1 hc

+

H ksub



Fig. 7. A simplified model for heat transfer in plasma spraying and the left is the thermal circuit.

(10)

M. Shen, B.Q. Li and Y. Bai / International Journal of Heat and Mass Transfer 150 (2020) 119267

Fig. 8. Interface velocity versus solidification thickness under different substrate heights.

Equating Eqs. (8) and (10) and rearranging gives

Ti − Tsub =

ρL

s

k

+



1 hc

+

H ksub



1 s 1 H μi + ρ L k + hc + ksub

 (Tm − Tsub )

Therefore the interface velocity can be expressed as

Vi =

ds = dt

Tm − Tsub



s 1 H μi + ρ L k + hc + ksub 1



(11)

It is obvious that the presence of a substrate always serves to reduce solidification rate. The faster the latent heat is dissipated through the substrate, the larger the solidification rate will be. Perfect contact, meaning bigger interfacial heat transfer coefficient, will also facilitate solidification. According to our simulation in subsection 3.5, the final splat thickness is about 0.5 μm. Fig. 8 displays the evolution of Vi with respect to s, which ranges from 0 to 0.5 μm. μi is taken to be 0.0012 m/s · K, according to Wang et al. [25]. Note that in this case Tsub is set to 423 K, hc = 108 W/m2 · K, and ksub = 18.8 W/m · K. According to Fig. 8, the interface velocity for various substrate heights H, though decreased with respect to solidification thickness, are all beyond the absolute stability velocity for planar solidification, about 0.3 m/s for YSZ by Wang et al. [25]. Hence planar solidification is more likely to occur. What’s more, the influence of substrate on interface velocity is evident and should be taken into account when the assumption of no isothermal delay fails or the Biot number Bi = hc H/ksub is large. 3.4. Validation of the full model Having examined pure fluid motion and analyzed heat transfer in splat formation, we are in a position to perform a case of real plasma spraying, with which experimental results can be compared. A YSZ droplet, 4.8 μm in diameter, is made to impinge at 200 m/s onto a solid surface of YSZ that is preheated to 1173 K. Note that the composition of the droplet and substrate differs. The substrate is made of 3YSZ (3% mol Y2 O3 –ZrO2 ). The thermal conductivity of the substrate is 3 W/m · K, the rest physical properties being nearly the same. The droplet and surrounding is at 30 0 0 K initially. Since the contact angle between the droplet and substrate

7

exerts no effect on fluid motion when the impacting velocity exceeds 1 m/s [26], we define it to be 120◦ . Interfacial heat transfer coefficient is set to 108 W/m2 · K. This value was in the range used by Wang et al. [25] in their work for ZrO2 droplet impact on ZrO2 substrate so that it is also utilized in YSZ droplet impact on YSZ substrate thanks to the slight differences in thermophysical properties of the two materials. The mesh size is 451 × 151, corresponding to a physical domain of 18 μm × 6 μm. Besides Cn = 0.0104 and λ = 145.36. The phase field mobility M = 2.2 × 10−13 m3 · s/kg. The characteristic Reynolds number is 203.40 and Weber number is 2629.95. Initially, a number of supercritical nuclei of different orientations are randomly seeded at the surface. The number density of the seed is 35/450. The numerical outcome is shown below (Fig. 9). At 0.1 μs, a few nuclei near the impacting point have started growing while the droplet is still spreading, but is close to an end. Heat transfer occurs immediately upon impacting, as collaborated by the temperature distribution at corresponding instants in Fig. 10. As time goes by, these crystals expand laterally and then impinge on each other, as seen at 0.2 μs. What is interesting is the secondary droplet at 0.2 μs, which may be generated even without solidification coming into play. Amazingly, the secondary droplet was later coalesced into the large one. Subsequently, the melt continues to be pulled backwards by surface tension, curling up as at 1 μs, after which time splat shape changes little. Besides, it takes about 1 μs to complete solidification, with much of it taken by solidification of the central part, which is relatively thicker. Fig. 11 shows the comparison of final splat microstructure between numerical and experimental results [27, 28], displaying reasonable agreement. Both show columnar grains that are perpendicular to the substrate surface and grow through the thickness of the splat. In other words, grain height is equal to splat thickness. When the droplet approaches the substrate inside a surrounding gas, the lubrication pressure inside the thin gas layer beneath the droplet will build up due to gas viscosity. The pressure force will further dimple the droplet, with the droplet above flowing outwards laterally because of mass conservation, transforming the point contact into one along a ring. A thin disk of air is thus trapped. Fig. 12 shows the dynamics of how the gas is trapped around the impacting point. The dimple is visible at 8 ns. The trapped gas layer, with the continued battering of the droplet impact, becomes thinner. Eventually at 20 ns, the amount of the entrapped gas is diminished due to the nature of the phase field model. The gas layer is too thin to be resolved by the numerical model so that it disappears a few time steps later. Such a phenomenon has been observed by Khatavkar et al. [29] and Yue et al. [30]. For finer meshes, the trapped gas will stay longer, but the physics will still be the same. Its effect on solidification is not as visible as in Fig. 16, where solidification on the gas-trapped region commences later than outside the region and the crystals around it are outgrown by those that first nucleate near where the thermal gradient, steered by convection, aligns with the flow. In the present case, even if solidification takes place later around the impacting point, the growth direction of these crystals will still be oriented vertically and will hardly be outgrown by the crystals nearby due partly to the tiny amount of gas entrapped and partly to the steep thermal gradient perpendicular to the substrate. 3.5. A supersonic impact with microstructure formation In this subsection, we investigate a supersonic impact with microstructure formation under practical conditions. A YSZ droplet, 4.8 μm in diameter, is made to impinge at 518 m/s onto a solid surface that is preheated to 423 K. The initial temperature is 30 0 0 K and contact angle 120◦ . Interfacial heat transfer coefficient is set to 108 W/m2 · K. The mesh size is 751 × 201, corresponding to a

8

M. Shen, B.Q. Li and Y. Bai / International Journal of Heat and Mass Transfer 150 (2020) 119267

Fig. 9. Snapshot sequences of a YSZ droplet impact at 200 m/s, with the red in the left column representing solid, the blue liquid, and the green gas, and with the different colors in the right column indicating grains of different orientations.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 10. Temperature distribution for the YSZ impact at 200 m/s.

M. Shen, B.Q. Li and Y. Bai / International Journal of Heat and Mass Transfer 150 (2020) 119267

9

Fig. 11. Comparison of microstructure with experimental results. (a) numerical (b) experimental.

Fig. 12. Gas entrapment for the YSZ impact at 200 m/s.

physical domain of 30 μm × 8 μm. The Cahn number Cn, the ratio of capillary length ξ to droplet diameter, is 0.0104, and λ = 145.36. The phase field mobility M = 2.2 × 10−13 m3 · s/kg. The characteristic Reynolds number is 526.79 and Weber number is 1.76 × 104 . Besides, a number of supercritical nuclei of different orientations are randomly seeded at the surface. The number density of the seed is 65/750. The numerical outcome is presented in Fig. 13. Not surprisingly, the spreading process is rather fast, and so is solidification. At 0.08 μs, the spreading process is close to an end, with a tiny portion on the bottom solidified. The droplet morphology at this instant, compared with that in Fig. 6, shows that the role of solidification may be neglected during this stage, a conclusion supporting the hypothesis of spreading first and solidification second. It is worth pointing out that at the end of spreading process, the velocity within the droplet becomes greatly reduced but still parallel to the substrate surface as at 0.08 μs, and that the droplet extends to a such a degree that heat transfer within the droplet can be assumed only in the vertical direction and, given the vanished vertical velocity component, mostly via conduction, as confirmed by the relatively flat isothermal lines at

corresponding instants in Fig. 14. In other words, fluid motion along the substrate surface does not influence heat transfer, which takes place mostly in the vertical direction where the vertical velocity component vanishes, only conduction surviving. Note that during this interval, from 0.08 μs to 0.2 μs, cooling rate is dropped because of increased solidification thickness. When the droplet reaches its maximum spread, its kinetic energy is transformed mostly into surface energy. Subsequently, under surface tension, the remaining molten part begins to be pulled backwards. It is after this time that solidification’s role becomes prominent, obstructing droplet retracting and further leaving the edge to curl up as at 0.2 μs. Heat transfer now concentrates mainly on the edge, but it progresses quite slowly because of increased edge thickness and of a low thermal conductivity of YSZ. With solidification occurring meantime, retracting gradually ceases, leading to a stable flattening ratio of about 7 as shown at 0.5 μs. On the other hand, solidification takes place first near the impacting point and then expands as the droplet spreads. On the rightmost column in Fig. 13, different colors represent different orientations. These crystals first expand laterally to cover the

10

M. Shen, B.Q. Li and Y. Bai / International Journal of Heat and Mass Transfer 150 (2020) 119267

Fig. 13. Snapshot sequences of a YSZ droplet impact at 518 m/s, with the red in the left column representing solid, the blue liquid, and the green gas, and with the different colors in the right column indicating grains of different orientations.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 14. Temperature distribution for the YSZ impact at 518 m/s.

M. Shen, B.Q. Li and Y. Bai / International Journal of Heat and Mass Transfer 150 (2020) 119267

11

Fig. 15. Steady state splat morphology and solidification microstructure for each impact, with the velocity increasing from 100 m/s to 400 m/s for the first four impact with an increment of 100 m/s.

substrate, then impinge on one another and grow upward vertically through droplet thickness and in alignment with the direction of thermal gradient. Besides, there is a competition between the driving force provided by the thermal gradient and the lower free energy along the axis of surface tension anisotropy. Dendritic structure is more likely to arise when the surface tension anisotropy gets an upper hand, seaweed when the two are neck and neck, and cellular or even planar when the thermodynamic driving forces hold the whip hand. The observed planar structure is due to the stronger thermodynamic driving forces, or equally, higher solidification rate, as elaborated on in subsection 3.3. The temperature distribution is given in Fig. 14. To summarize, in supersonic plasma spraying, solidification occurs mainly after spreading, fluid motion exerts no effect on solidification microstructure morphology, which consists of expected columnar grains with limited inter-grain competition. In addition, to investigate the effect of impacting velocity, a few simulations are run with the impacting velocity varying from 100 m/s to 400 m/s. The numerical configuration is the same as the impact at 518 m/s, except for the mesh size, which is now 601 × 151, corresponding to a physical domain of 24 μm × 6 μm. The domain has been demonstrated large enough to minimize its effect on the final result. The freezing time for each impact, along with the steady state splat morphology and solidification microstructure, is shown in Fig. 15. The last is the YSZ impact at 518 m/s. Intuitively speaking, the maximum spread is proportional to the impact velocity. The higher the impacting velocity, the larger the maximum spread. This is seen in Fig. 15. On the other hand,

as the impacting velocity diminishes, the splat thickness increases due to volume conservation and to a smaller flattening ratio, resulting in longer freezing times. The freezing time for each impact predicted by the numerical model is quite close to those anticipated by Vardelle et al. [24]. On the other hand, the undercooling of the droplet is not changed due to a fixed temperature of the droplet and that of the substrate at the beginning. But with the increase of solidification thickness, the growth direction of crystals is more likely to change their course, from perpendicular to the substrate to some certain direction determined by the competition with crystals nearby. In other words, anisotropy comes into play with the decrease of undercooling that is caused by the incremental solidified thickness. Such a phenomenon is most obvious for the impact at 100 m/s but it still produce a columnar structure on the whole. 3.6. Convection effect on solidification morphology In the last two sections, we have validated the full model and applied it to a practical case in plasma spraying, with the growth direction being perpendicular to the substrate, a phenomenon due largely to a strong vertical thermal gradient there. It seems that although there is strong fluid motion it exerts no significant effect on crystal growth. But this is not always the case. In this section, an ingeniously designed simulation, although experimentally difficult to realize, is performed to show the effect of convection. A YSZ droplet, 0.64 μm in diameter, is made to impinge at 500 m/s onto a solid surface. The initial temperature is 2100 K and

12

M. Shen, B.Q. Li and Y. Bai / International Journal of Heat and Mass Transfer 150 (2020) 119267

Fig. 16. Spreading and solidification evolution for the YSZ impact at 500 m/s, with the right column representing grains of different orientations, and with the red in the left column indicating solid and the blue liquid. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

contact angle 120◦ . The mesh size is 501 × 201, corresponding to a physical domain of 2 μm × 0.8 μm. Cn = 0.0078 and λ = 14.53. The phase field mobility M = 8 × 10−15 m3 · s/kg. The characteristic Reynolds number is 67.80 and Weber number is 2191.63. Besides, a number of supercritical nuclei of different orientations are seeded at the surface. The number density of the seed is 50/500 and the seeds are controlled equidistantly for better explanation of the convection effect. The numerical outcome is presented in Fig. 16. Note that the computation domain is filled with undercooled liquid at the beginning and that the boundary condition is symmetric for all variables except the order parameter c. Solidification process is fast due to the large undercooling. By 10 ns, the droplet has been almost solidified except for a little portion at the tip. Strikingly, solidification begins not from the impacting point but from a small distance apart. This may be caused by gas entrapment around the impact point. It is also obvious from Fig. 16 that the liquid around the impacting point starts solidifying later than the intermediate part, as seen at 5 ns, and that the crystals, which begins solidifying first, are also of the largest sizes, tending to grow upstream. Convection matters here. However, the grains behind the ones of large sizes are influenced much less by convection, showing a columnar structure. 4. Concluding remarks This paper has presented a numerical model for impingement and spreading of a high velocity molten YSZ droplet with solidification microstructure formation on a solid substrate. The model

is based on the explicit finite difference solution of mass, momentum and energy balance equations coupled with the Cahn-Hilliard equation to monitor the moving interface associated with droplet spreading and with the phase field equation for microscopic polycrystalline growth as the droplet is flowing along the preheated solid substrate surface. The explicit finite difference computational scheme was further enhanced via the shared-memory parallelism OpenMP. Extensive numerical experiments were conducted with YSZ droplets being of typical sizes in supersonic plasma spraying. The computed results are compared with experimental measurements. The major findings are summarized below: First, given the rapid spreading of YSZ droplets along a substrate surface, solidification may be safely assumed to take place mostly after spreading, even though there exists a large thermal gradient between the droplet and preheated substrate. Second, despite a supersonic impacting velocity, disk-like splats are predicted in the simulations. Disk-like splats are preferable in plasma spraying while fragmented ones are easy to deteriorate interlayer cohesion by increasing thermal contact resistance through absorbing the surrounding gas. Third, convection effect is negligible in columnar grain formation under practical plasma spraying conditions. The steep thermal gradient dominates the growth directions of multiple grains. Solidification takes places earliest not in the impacting center, but some distance apart, when gas entrapment occurs. Fourth, given the negligible role of fluid motion in the formation of columnar grains within a splat, it is desirable to increase the impacting velocity to achieve higher flattening ratios,

M. Shen, B.Q. Li and Y. Bai / International Journal of Heat and Mass Transfer 150 (2020) 119267

which, on one hand, facilitates the bonding between a splat and a substrate, alleviating gas entrapment, and on the other, reduces solidification time since more liquid is to be in contact with the substrate. Moreover, increasing the impact velocity may benefit splat/splat bonding with an even splat surface and prompt epitaxial growth of grains across splat surface. Declaration of Competing Interest There are no conflicts to declare. Acknowledgement This work is supported by HPC Platform, Xi’an Jiaotong University. Thanks also go to Professors Changjiu Li, Qingzhen Yang and Dr. Yu Wang for constructive discussions and helpful suggestions. This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors. References [1] K. Bobzin, N. Bagcivan, D. Parkot, I. Petkovic´ , Simulation of PYSZ particle impact and solidification in atmospheric plasma spraying coating process, Surf. Coat. Technol. 204 (8) (2010) 1211–1215. [2] H. Tabbara, S. Gu, Numerical study of semi-molten droplet impingement, Appl. Phys. A-Mater. Sci. Process. 104 (4) (2011) 1011–1019. [3] S. Oukach, H. Hamdi, M.E. Ganaoui, B. Pateyron, Numerical study of the spreading and solidification of a molten particle impacting onto a rigid substrate under plasma spraying conditions, Therm. Sci. 19 (1) (2015) 277–284. [4] Z. Zhu, S. Kamnis, S. Gu, Numerical study of molten and semi-molten ceramic impingement by using coupled Eulerian and Lagrangian method, Acta Mater. 90 (2015) 77–87. [5] H. Safaei, M.D. Emami, H.S. Jazi, J. Mostaghimi, Application of compressible volume of fluid model in simulating the impact and solidification of hollow spherical ZrO2 droplet on a surface, J. Therm. Spray Technol. 26 (8) (2017) 1959–1981. [6] R.K. Shukla, V. Patel, A. Kumar, Modeling of rapid solidification with undercooling effect during droplet flattening on a substrate in coating formation, J. Therm. Spray Technol. 27 (2018) 269–287. [7] R. Kobayashi, J.A. Warren, W.C. Carter, Vector-valued phase field model for crystallization and grain boundary formation, Physica D 119 (3–4) (1998) 415–423. [8] R. Kobayashi, J.A. Warren, W.C. Carter, A continuum model of grain boundaries, 140 (1–2) (20 0 0) 141–150. [9] J.A. Warren, R. Kobayashi, A.E. Lobkovsky, W.C. Carter, Extending phase field models of solidification to polycrystalline materials, Acta Mater. 51 (20) (2003) 6035–6058. [10] A. Karma, W.J. Rappel, Quantitative phase-field modeling of dendritic growth in two and three dimensions, Phys. Rev. E 57 (4) (1998) 4323–4349.

13

[11] Y. Shu, X. Ai, B.Q. Li, Phase-field modelling of solidification microstructure formation using the discontinuous finite element method, In. J. Num. Meth. Eng. 69 (6) (2006) 1194–1211. [12] C.W. Kang, J.K. Tan, L.S. Pan, C.Y. Low, A. Jaffar, Numerical and experimental investigations of splat geometric characteristics during oblique impact of plasma spraying, Appl. Surf. Sci. 257 (2011) 10363–10372. [13] C.L. Bot, S. Vincent, E. Meillot, F. Sarret, J. Caltagirone, L. Bianchi, Numerical simulation of several impacting ceramic droplets with liquid/solid phase change, Surf. Coat. Technol. 268 (2015) 272–277. [14] C.M. Wu, M. Bussmann, J. Mostaghimi, The impact of a partially molten YSZ particle, J. Therm. Spray Technol. 18 (5–6) (2009) 957–964. [15] P.T. Yue, C.F. Zhou, J.J. Feng, Sharp-interface limit of the Cahn-Hilliard model for moving contact lines, J. Fluid Mech. 645 (3) (2010) 16. [16] F. Magaletti, F. Picano, M. Chinappi, L. Marino, C.M. Casciola, The sharp-interface limit of the Cahn-Hilliard/Navier-Stokes model for binary fluids, J. Fluid Mech. 714 (2013) 95–126. [17] V.V. Khatavkar, P.D. Anderson, P.C. Duineveld, H.E.H. Meijier, Diffuse-interface modelling of droplet impact, J. Fluid Mech. 581 (2007) 97. [18] C.F. Zhou, P.T. Yue, J.J. Feng, C.F. Ollivier-Gooch, H.H. Hu, 3 D phase-field simulations of interfacial dynamics in Newtonian and viscoelastic fluids, J. Comput. Phys. 229 (2) (2010) 498–511. [19] Q. Zhang, T.Z. Qian, X.P. Wang, Phase field simulation of a droplet impacting a solid surface, Phys. Fluids 28 (2) (2016) 022103. [20] L. Luo, Q. Zhang, X. Wang, X.C. Cai, A parallel finite element method for 3 D two-phase moving contact line problems in complex domains, J. Sci. Comput. 3 (2017) 1–27. [21] M. Plapp, Three-dimensional phase-field simulations of directional solidification, J. Cryst. Growth 303 (1) (2007) 49–57. [22] S. Gurevich, A. Karma, M. Plapp, R. Trivedi, Phase-field study of three-dimensional steady-state growth shapes in directional solidification, Phys. Rev. E: Stat., Nonlinear Soft Matter. Phys. 81 (1) (2010) 011603. [23] B. Nestler, D. Danilov, P. Galenko, Crystal growth of pure substances: phase– field simulations in comparison with analytical and experimental results, J. Comput. Phys. 207 (1) (2005) 221–239. [24] M. Vardelle, A. Vardelle, A.C. Leger, P. Fauchais, D. Gobin, Influence of particle parameters at impact on splat formation and solidification in plasma spraying processes, J. Therm. Spray Technol. 4 (1) (1995) 50–58. [25] G.X. Wang, R. Goswami, S. Sampath, V. Prasad, Understanding the heat transfer and solidification of plasma-sprayed yttria-partially stabilized zirconia coatings, Mater. Manuf. Process. 19 (2) (2004) 259–272. [26] M. Pasandideh-Fard, J. Mostaghimi, On the spreading and solidification of molten particles in a plasma spray process effect of thermal contact resistance, Plasma Chem. Plasma Process. 16 (1995) S83–S98. [27] E.J. Yang, X.T. Luo, G.J. Yang, C.X. Li, C.J. Li, M. Takahashi, Epitaxial grain growth during 8YSZ splat formation on polycrystalline YSZ substrates by plasma spraying, Surf. Coat. Technol. 274 (2015) 37–43. [28] E.J. Yang, Study on the crystalline structure and microstructure formation of plasma sprayed ceramic coatings, PhD thesis, Xi’an Jiaotong University, Xian (2014). [29] V.V. Khatavkar, P.D. Anderson, P.C. Duineveld, H. Meijer, Diffuse interface modeling of droplet impact on a pre-patterned solid surface, Macromol. Rapid Commun. 26 (2005) 298–303. [30] P.T. Yue, C.F. Zhou, J.J. Feng, Spontaneous shrinkage of drops and mass conservation in phase-field simulations, J. Comput. Phys. 223 (2007) 1–9.