Numerical Simulation of Dendritic Pattern Formation in an Isotropic Crystal Growth Model on Curved Surfaces

: In this paper, we present several numerical simulation results of dendritic pattern formation using an isotropic crystal growth model, which is based on phase-ﬁeld modeling, on curved surfaces. An explicit time-stepping method is used and the direct computing method to the Laplace–Beltrami operator, which employs the point centered triangulation approximating Laplacian over the discretized surface with a triangular mesh, is adopted. Numerical simulations are performed not only on simple but also on complex surfaces with various curvatures, and the proposed method can simulate dendritic growth on complex surfaces. In particular, ice crystal growth simulation results on aircraft fuselage or metal bell-shaped curved surfaces are provided in order to demonstrate the practical relevance to our dendrite growth model. Furthermore, we perform several numerical parameter tests to obtain a best ﬁtted set of parameters on simple surfaces. Finally, we apply this set of parameters to numerical simulation on complex surfaces.


Introduction
The beauty of nature's crystals has led many researchers to investigate them over the past decades. From a physical point of view, crystal growth can be explained as a phase transition by heat transfer from a liquid phase to a solid phase. Accordingly, many researchers have experimented with forming dendritic patterns, a large amount of research results however are needed to find out the influence of the parameters of the explainable model so that, the analysis through computer simulation has been also actively conducted therefore. There are several methods that have been devised and developed to investigate dendritic pattern formation of crystals such as the front tracking method [1], Monte-Carlo simulation approach [2], the level set method [3], and the phase-field method [4][5][6]. In particular, the phase-field modeling, which captures the interface in explaining the dynamics when dendritic patterns are formed on surface, is a field that has been continuously studied. The phase-field method has been widely applied in numerical simulation of dendritic crystal growth model [7][8][9][10][11], however most literatures only consider the solidification of pure substances in undercooled melts in two-and three-dimensional spaces; not on curved surfaces. Using the phase-field method, for instance, the authors in [12] successfully simulate the crystal growth with an anisotropic diffusion scheme in two dimensional space. A highly nonlinear system combined with heat equation was proposed for numerical simulations of dendritic crystal growth model in [13] using a similar anisotropic scheme. Li et al. [14] propose the thermodynamic consistent model for the numerical approximation of energy and entropy conservation of crystal growth model, some numerical experiments of dendritic crystal growth are presented to demonstrate the accuracy of their proposed scheme. A phase-field crystallization model is proposed for the formation of dendrite morphology in the process of single component melt growth by Kobayashi [15]. In [5], a rescaling rule is presented for simulating dendritic growth in a two-dimensional domain at various undercoolings with an adaptive phase-field model and is also extended to the three-dimensional adaptive phase-field model. Lin et al. [16] compare the simulated equilibrium shapes of different anisotropic strengths with the exact solutions, and then further explore the dendrite growth of pure materials with different anisotropic strengths under different undercoolings. In the study of [17], robust and fast phase-field simulations of dendritic crystal growth are performed in both two-and three-dimensional spaces, their computational results are consistent with the numerical experiments.
Meanwhile, there are several researches that have focused on pattern formation and computation on surfaces. García et al. [18] focus on the mechanism of a crystal phase in two-dimensional curved space. Ren et al. [19] simulate the dendrite growth by introducing an orientation field which is the degree of freedom, the algorithm can generate arbitrary symmetric graphs with different degrees of asymmetric damage effects on arbitrarily curved surfaces.
In this paper, we present several numerical simulation results of an isotropic crystal growth model, which is based on phase-field approach, on curved surfaces with triangular discretizations. In order to compute numerically the surface diffusion operator, called the Laplace-Beltrami operator, we employ Xu's approach [20], while there are several other approaches to calculate this operator [21,22]. Furthermore, a comparative study is conducted on the effect of parameters affecting dendritic pattern formation on simple surfaces. Subsequently, a fitted parameter set is employed to simulate the dendrite growth on complex surfaces. Especially, we provide the dendritic crystal growth simulations on aircraft body or metal bell-shaped surfaces. As these surfaces are used in real life, the aircraft de/anti-icing system is a very important problem in aerodynamics [23], and the process of growing ice crystals on metal surfaces is a common phenomenon.
The manuscript is organized as follows. We introduce the mathematical model in Section 2. In Section 3, we present the numerical solution. In Section 4, we perform the numerical analysis and numerical experiments on surfaces. Finally, conclusions are drawn in Section 5.

Mathematical Model
The model of the crystal growth on a smooth closed surfaces S in R 3 is derived from a Lyapunov function [24]. As the crystal growth phase-field model which is consisted with the anisotropic function has been studied in [17], and the model performed fast, robust, and accurate results of crystal growth for phase-field simulation whether in two or three-dimensional space. In this study, we use phase-field model to simulate the crystal growth without anisotropic function on surfaces. The solidification process is modeled by using the standard dimensionless phase-field equations [25] and is represented as follows: In Equations (1) and (2), U is the temperature field, D is a constant which is defined by D = ατ 0 /E 0 2 , where τ 0 is the characteristic time, α is the thermal diffusivity, and E 0 is the characteristic length. λ is the coupling parameter and is defined by λ = D/0.6267 [25]. The order parameter φ is 1 or −1 in the solid phase or liquid phase, respectively. The interface of two phases is defined by φ = 0. Note that all the parameters are dimensionless. For more details, refer to [25] and the references therein.

Numerical Solution
In this section, we apply the explicit Euler method to solve the crystal growth systems (1) and (2) on a smooth closed surface. At first, we present a discretization of Laplace-Beltrami operator [20,26] on the curved surfaces. Figure 1a shows that M is a triangular surface mesh of given surfaces S. The set {x} N i=1 with N points means the surface vertex set of M. Considering each surface vertex x i with index n, we define that N 1 (i) = {i 1 , i 2 , · · · , i n } is a set of surface vertex indices representing one-ring neighbors with i 1 = i n (see Figure 1b). Here, T j is a triangle composed of three vertices x i , x j − and x j as shown in Figure 1c. We write simply for convenience as φ(x i ) = φ i and U(x i ) = U i . Now, we consider the sum of triangles T j sharing the vertex x i and denote it by A(x i ) (see Figure 1d). A(x i ) is expressed as follows [20]: Then, we get the following discretization systems of Laplace-Beltrami operator over given surface S [26].
where θ ij and θ ij + are angles in triangles T j and T j+ , respectively (see Figure 1c).
where ∆t is the time step size. Then the governing system is discretized as

Numerical Results
We perform computational experiments to show the efficiency and robustness of the proposed scheme on curved surfaces.

Crystal Growth on Sphere Surface
In this section, we test the crystal growth on a sphere. We take the condition for a sphere surface and the initial condition of sphere surface is defined in domain Ω = (−200, 200) 3 as follows.
where r is a given zero level set radius of the initial condition, U = 0 represents the melting temperature of the pure material, and U = ∆ corresponds to the initial undercooling.
A satisfied evolution effect depends on the appropriate choice of parameters. In following simulations, it is important to set the time step size ∆t and mesh size appropriately. Then we first test to know the effect of ∆t and mesh size and decide proper two parameters in our simulations. One of the factors influencing crystal growth can be the curvature of the sphere. We study the effect of the curvature through the below test. In addition, to verify the robustness of the proposed scheme, we investigate the impact of the zero level set radius of the initial condition r, and the initial undercooling ∆. Unless otherwise specified, we choose D = 7 and λ = 11.29. Then we perform various simulations on a set of ∆ and r.

Effect of Time Step Size
First, we confirm the effect of time step size ∆t and decide the appropriate ∆t for our simulations. Here, we set the initial radius r = 30, undercooling ∆ = −0.75 and use the same mesh size. Under these conditions, we select four time step size such as ∆t = 0.3, ∆t = 0.2, ∆t = 0.1, and ∆t = 0.05. As shown in Figure 2, the states of evolution look similar at same time. However, when ∆t = 0.2 and ∆t = 0.3, the results of crystal show a bit less-growing than other two conditions, and Figure 2c,d show a satisfied performance. Therefore, we use ∆t = 0.1 as proper time step size for our simulations, while a smaller ∆t will take more computational cost.

Effect of Mesh Size
In the following test, we examine the effect of mesh size. Because our proposed scheme is to solve the Equations (1) and (2) on surface directly, the mesh size is decided by the average of edge lengthh. We use DistMesh to represent the surface and decide the initial edge length h. Based on h, this algorithm finds optimized surface and each edge length varies widely after optimization. Then, we calculate the average of the varied edge length and denote byh.
We take R = 165, r = 30, ∆ = −0.75, ∆t = 0.1, D = 7 and λ = 11.29 and focus on the impact which the mesh size affects. We select four conditions for initial edge length such as h = 5.4, 4.7, 4, and 3.3. Then we can get the average edge lengthh = 4.8, 4.1, 3.5, and 2.9. Figure 3 shows snapshots of evolution states of crystal growth sources with different mesh size in surface at different time (t = 0, 500∆t, 1000∆t, and 1500∆t).
As depicted in Figure 3, various mesh sizes affect the growth of dendritic crystals. In other words, h must be kept at a certain level in order to be relatively more suitable for dendrite pattern formation and to have similar dendrite tip velocities to actual physical phenomena. In particular, results of Figure 3c,d have a little difference in growth of crystal. This phenomenon indicates that the average edge lengthh = 3.5 is sufficient for crystal growth simulation, which saves the computing resources. Based on the above simulations, we use ∆t = 0.1 andh = 3.5 for following tests in this paper.

Effect of Curvature
Next we consider the effect of the curvature of sphere. It is known that the curvature of sphere depends on the radius. As you know, there is an inverse relation between the curvature κ and the radius R. That is, κ becomes larger when R is smaller and smaller when R is larger. The radius of sphere used on the above simulations is R = 165. In this test, we use varied radii of the sphere by R = 130, 165, and 200, and compare each result. Here, we use r = 30 and ∆ = −0.75 without changing other parameters for the test.
As shown in Figure 4, from the perspective of visual effect, the states of dendritic crystals growth show almost no difference on the sphere surfaces. In other words, each crystal source on the sphere surface having different radius grows well under different curvatures. Therefore, we can conclude that the curvature has little effect on the growth of crystal on the object surface.

Effect of Undercooling Parameter ∆
The initial undercooling ∆ affects on growth of crystal on surface. Therefore, we study the effect of ∆ and set three conditions of ∆ = −0.55, −0.75 and −0.95. Setting the initial radius r = 30 on the sphere, we compare results from different values of ∆. Figure 5 shows the evolution of crystal growth on a sphere with different undercooling ∆ at different time (t = 0, 500∆t, 1000∆t, and 1500∆t).
As shown in Figure 5, we confirm the effect of ∆ easily. The bigger the value of ∆, the slower the growth of crystal appears. Conversely, results using ∆ = −0.95 show more fast growth. These results with ∆ = −0.55 and ∆ = −0.95 are hard to call crystal as shown in Figure 5a,c. However, crystal which grew up using ∆ = −0.75 show the shape we want and we decide proper the value of ∆.

Effect of Radius r
Now, we investigate the effect of radius r of the initial crystal growth source. Unchanging other parameters and conditions, we only set r differently. That is, we take three radius r = 15, 30, and 60, and simulate under ∆ = −0.75. Figure 6 shows states of crystal growth on sphere with different initial radius r over time (t = 0, 500∆t, 1000∆t, and 1500∆t). The difference among three results is clear. Figure 6a shows that crystal sources with r = 15 are less grown and branches are less stretched. However, results using r = 60 show that crystals grow and branches stretch more as shown in Figure 6c. From this test, we can see that when the initial radius r of the source is large, growth is fast and branches of crystal are extended further. Next, we observe the crystal growth states when sources are asymmetrically arranged on the sphere. As shown in Figure 7a, six crystal growth sources set asymmetrically on sphere surface. In addition, these sources are located the different positions and radius (r = 15, 20, 25, 30, 35, and 40). Parameters for this test are set as ∆t = 0.1, D = 7, λ = 11.29, and ∆ = −0.75. Figure 7b-d are snapshots of evolution of crystal growth sources which are arranged asymmetrically on sphere at 500∆t, 1000∆t, and 1500∆t, respectively.

Crystal Growth on Ellipsoid Surface
Asymmetrical arrangement of crystal sources is also applied on the ellipsoid surface. The initial condition of ellipsoid surface is defined in domain Ω = (−200, 200) 3 as follows.
where f = 0 denotes the ellipsoid surface. Three sources having different radius (r = 15, 30, and 45) are located at different positions on the ellipsoid surface. Parameters used for this test are the same with the above test and Figure 8 shows the initial condition of ellipsoid with three sources and snapshots of evolution of these crystal growth at t = 300∆t, 600∆t, and 900∆t.

Crystal Growth on Torus Surface
At last, we perform the crystal growth process on the torus surface. Three sources with radius r = 10, 20, and 30 randomly distributed on the torus surface. Parameters used for this test are the same with the above test and Figure 9 shows the initial condition of torus with three sources and snapshots of evolution of these crystal growth at t = 800∆t, 1600∆t, and 2400∆t. Through these simulations, our method can be applied not only to one crystal source on the surface, but many crystal sources that have different radius.

Crystal Growth on Complex Surface
In this section, we present the evolutions of Equations (4) and (5) on complex surface. For a concrete example, we adopt at first a point cloud data for an aircraft body to construct the discretized surface using a DistMesh function in MATLAB. Numerical simulation of the structure in which ice crystals form in the fuselage of airplane is a practically relevant experiment that can help in de/anti-icing it from aerodynamics [23]. The surface is bounded by [−114, 114] 39,39]. Similar to the above investigation cases, we manipulate the mean of spatial step sizeh as close as possible to 3.5 by multiplying the scaling factor. Figure 10 shows the evolutions of crystal growth on the surface of aircraft body and wings. Initial states of crystals on aircraft surface are set to be the most frequent parts in reality. Here, we use all the specific parameters are the same as above except for ∆t = 0.125.
For the next example, we perform the similar test to the metal bell-shaped surface. We fit the isosurface of the metal bell-shaped surface to [−69, 69] × [−69, 69] × [−69, 69]. The cubic interpolation is used to make the signed distance function of the metal bell-shaped surface in MATLAB. Figure 11 shows snapshots of the evolutions of crystal growth on the metal bell-shaped surface. Here, we usē h ≈ 3.5, r = 3.1883, and ∆t = 0.125. (c) (d) Figure 11. Snapshots of the evolutions of crystal growth on the metal bell-shaped surface. We usē h ≈ 3.5, r = 3.1883, and ∆t = 0.125. The evolution time at each row is (a) t = 0, (b) t = 180∆t, (c) t = 360∆t, and (d) t = 540∆t, respectively. Note that each column in each row shows the isosurface from different angle at the same time.

Conclusions
The dendritic pattern formation in crystal growth is an important topic and has a bunch of interesting applications in phase-field modelings. There are many studies related crystal growth simulation on two-and three-dimensional domains using phase-field modeling so far, we present the isotropic crystal growth simulation on curved surfaces in this paper. The Lyapunov function on smooth closed surfaces S is the base of our mathematical model and the Laplace-Beltrami operator is evaluated directly on triangularly discretized surfaces M by using Xu's approach [20]. The first-order time-stepping method is employed due to ease of implements. There are many parameters that affect the dendritic pattern formation in crystal growth, such as time step size ∆t, average mesh sizeh, and undercooling parameter ∆. We perform several numerical investigations to obtain a suitable set of parameters and choose ∆t = 0.125,h ≈ 3.5, ∆ = −0.75 as a default setting when a computational domain is contained in boxed region which is bigger than [−120, 120] 3 but smaller than [−180, 180] 3 approximately, because this set has shown that dendrite crystals are growing physically correctly over time and, further, it does not use excessive resources for computation. By using these parameters, we describe the snapshots of isotropic crystal growth on various curved surfaces. In particular, we provide dendrite growth simulations on physically and practically relevant shapes such as aircraft body and metal bell-shaped surfaces. In addition, simulation results that can be applied to the actual de/anti-icing problem are derived not only in the initial condition of the circular shape but also in the initial state of the curvilinear or rim shape. Because this model does not contain anisotropic cross derivatives, extending a mathematical model with anisotropic terms and computing the model directly on surfaces is left to one of our future works.