Stable and Anisotropic Freezing Framework With Interaction Between IISPH Fluids and Ice Particles

In this paper, we present a novel method that can stably express the directional ice form caused by freezing of flowing water. The key to the proposed framework is to reflect the flow of fluids with viscosity in the direction of ice growth. Water is simulated by applying a new viscous technique to the implicit incompressible fluid simulation, and the proposed anisotropic freezing solution is used to express directional ice and glaze effects. The conditions under which water particles turn into ice particles are calculated according to a new energy function based on humidity and water flow. The humidity is approximated based on the virtual water film on the surface of the object, and the flow of fluid is incorporated into our anisotropic freezing solution to guide the growth direction of the ice. As a result, the proposed technique reliably produces glaze and directional freezing effects according to the flow direction of viscous water.


I. INTRODUCTION
In various VR (Virtual Reality) contents such as movies and games, the freezing effect is used in various scenes expressing freezing scenes in cold regions or extreme environments: ''The Day after Tomorrow'' (2004), ''The Last Airbender'' (2010), ''The Huntsman: Winter's War'' (2016). VR content is often produced by bending ice in a certain direction due to factors such as the surface of an object, the flow of water, or wind direction, to express immersively the extreme winter, not just ice that grows downward by gravity. In the field of computer graphics, there have been many studies to show the creation and growth of ice. Im et al. [1] proposed a technique to express the interior of opaque ice, and Kharitonsky and Gonczarowski [12] proposed a technique to generate and grow icicles using physics-based simulation. Kim et al. [2] and Gagnon and Paquette [13] proposed a technique for expressing icicles and glaze effects appearing on the surface of an object using procedural methods. Since these techniques focus on expressing the icicle and glaze effects generated by slowly falling water droplets, it is difficult to express when The associate editor coordinating the review of this manuscript and approving it for publication was Songwen Pei .
considering only the growth of ice due to gravity without interacting with water.
In this paper, we propose a physics-based simulation technique that reflects the flow of water in the freezing process. Our method does not use a temperature-based model, but a crystallization-based state change model like Im et al.'s method [32]. Like their method, we also drew inspiration from the supercooled water phenomena, which stays liquid at temperatures below zero degrees and then freezes quickly when colliding with a cold surface(see Figure 1). In addition, in order to more accurately express the characteristics of freezing occurring as water flows, the growth direction of ice is modeled anisotropically according to the fluid flow. In the previous approaches, it was difficult to express the phenomenon of freezing as water flows along the surface of a solid, and because water particles are scattered by collisions between water and solid, the freezing modeling becomes unstable, resulting in lower freezing quality. In addition, in calculating the growth direction of ice that freezes while extending like a branch, it is difficult to realistically express the growth direction with an isotropic kernel [32]. Our technique can anisotropically generate or grow ice even when the flow of water changes dynamically due to external forces or VOLUME 9, 2021 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ collisions with objects. Therefore, it is possible to express not only the growth of ice according to the flow of water, but also the phenomenon of rapidly freezing water.

A. PROBLEM STATEMENT
In general, the flow of fluid has a significant influence on determining the shape of the ice during simulation. Various techniques have been proposed to perform freezing simulation such as icicles in the field of computer graphics, but in most cases, the form of freezing is determined only by gravity. This causes patterned freezing, making it difficult to produce realistic freezing in various virtual environments. When modeling the winter environment, ice is created in a curved form not only by gravity but also by various environmental factors such as wind and water flow, and this feature is important when creating virtual environments or games (see Figure 2). Typical examples of external forces affecting freezing surfaces are wind and liquid flow. In particular, supercooled water is more susceptible to liquid motion because it freezes immediately (see Figure 2). However, since the existing method explicitly calculates the force based on water particles, if the speed of the water particles is high or the impact is strong due to the collision with the solid, the water particles become unstable in the process of changing to ice. As a result, ice surfaces are not generated properly or noise is generated severely. In addition, when water particles collide with a solid with large curvature, the surface tension becomes unstable, resulting in scattering of ice particles, resulting in generation of low-quality ice surfaces. Also, when using the isotropic approach to generate a curved ice surface considering liquid motion, the direction cannot be accurately estimated. Therefore, in this paper, we proposed a method to improve the quality of freezing surfaces by solving the aforementioned stability problem.
As shown in Figure 2, it can be seen that the frozen ice near a tree or house has a curved shape due to the physical environment such as wind or water flow. In this paper, we propose a new framework that can stably express this freezing shape. To implement the proposed freezing solver, the following sub-problems must be solved in the input fluid simulation: 1) Integration with reliable water simulation solvers.
Most of the freezing phase changes from water to ice, and in this paper, we propose a freezing solver based on IISPH (Implicit Incompressible SPH) to stably integrate fluid motion into the growth direction of ice. 2) Anisotropic calculation of freezing direction. In the freezing process, the growth direction is calculated based on the anisotropic kernel to accurately and stably deliver the fluid motion to the freezing solver. 3) Improved quality of freezing surfaces by integrating advanced surface tension and viscosity. When the phase of water particles is changed to ice, scattering motion occurs due to impact by collision, and this instability lowers the quality of ice shape generation. In this paper, improved surface tension and viscosity terms are integrated into the IISPH-based freezing solver to alleviate these problems.

II. RELATED WORK
In this section, we explore various studies related to ice, freezing, and snow simulations similar to the proposed method.

A. ICE AND FREEZING SIMULATIONS
In the field of physics-based simulation, studies related to the formation and growth of ice have been continuously published. Kim et al. [3] realistically expressed frost growing on flat or curved surfaces by expressing the growth of ice through crystallization of small and thin water particles using discrete-continuous method (i.e., diffusion limited aggregation, phase field method, and stable fluid solver). Im et al. [1] and Nishino et al. [4] proposed a method to express detailed features such as opaque areas inside ice. Kim et al. [14] proposed a method to model cracks and scratches that appear on ice surfaces when ice is impacted. Kharitonsky and Gonczarowski [12] realistically expressed the growth of icicles, taking into account the physical properties of water droplets, based on thermodynamics. Unlike existing approaches, Kim et al. [2] expressed the growth of icicles by solving the Stefan problem. This method realistically represents the shape of the icicles being merged together by simulating several icicles at the same time, but it required a lot of computation. Gagnon and Paquette [13] proposed a procedural method that allows users to easily control the ice growth simulation, making it possible to model not only the shape of the ice but also the glaze effects according to the user's intention. Compared to the previously mentioned approaches, their method produces icicles quickly, but only expresses a single form of icicles, and it is difficult to express the merging of multiple icicles. Ishikawa et al. [5] simulated the freezing of water droplets as they fall under gravity based on how glaze effects and icicles are generated. Since the aforementioned methods modeled the ice formation under the precondition that the water droplets slowly flow downwards only by gravity, it is difficult to model the freezing shape generated under other conditions. Several studies have attempted to represent the freezing of flowing water [6]- [8], [15], [16]. These methods represent an opaque freezing phenomenon based on the diffusion of air bubbles, but cannot express icicles or glaze effects. Iwasaki et al. [7] proposed a method to express the melting and freezing effects on a particle basis, and expressed ice spikes generated on the ice surface, but most of the results were focused on melting. Wicke et al. [8] simulated various materials based on implicit connectivity, and produced the result of ice formation in flowing water using a heat transfer method. Makkonen [15], [16] proposed an ice growth theory model based on observations, but it is difficult to attempt simulation or visualization because they do not propose an algorithm to discretize and implement the theory.

B. SNOW SIMULATIONS
Many approaches to simulate ice or snow using continuum mechanics have been studied and used to represent elastoplastic snow models [17]. Since Stomakhin et al. [18] first introduced the material point method (MPM) to the snow simulation technique in the field of computer graphics, the MPM method has been expanded in various ways to handle phase shifts [19]- [23]. Recently, Han et al. solved the frictional contact between snow and hair in the MPM framework [24].
In addition to MPM, there are other discretization techniques used to simulate snow. Wong and Fu [25] applied the discrete element method (DEM) to snow particles and spring mechanics to represent snow in an interactive simulation environment. Mukai et al. [9] simulated snow splitting down from the roof using the extended DEM method, and Dagenais et al. [26] simulated snow movement using position based dynamics (PBD) and level-set. Takahashi and Fujishiro [37] modeled the movement of snow using flow based on smoothed particle hydrodynamics (SPH).
Abdelrazek et al. [27] simulated an avalanche using the Bingham viscosity model, and Goswami et al. expressed the snow movement in real time on a GPU [10].

III. PROPOSED FRAMEWORK
Previous studies explicitly calculate the force of water particles, so if the speed of water particles is high or impact is applied due to collision with solids, the process of changing water particles into ice becomes unstable. As a result, ice surfaces may not be generated properly or noise may occur severely [32]. In this paper, we use IISPH, an implicit method that can reliably calculate fluid motion, rather than an explicit solver that directly calculates force based on position.
The framework proposed in this paper consists of a fluid simulation step and a freezing simulation step. In the fluid simulation stage, fluid motion, fluid-solid interaction, surface tension and contact force are calculated according to the IISPH (Implicit incompressible SPH) [28] based solution. We calculate the interaction between the fluid and the solid by creating boundary particles [29] on the vertices, edges, and faces of the solid mesh. The contact force and surface tension were implemented by the method of He et al. [31], and a viscosity-based fluid was modeled to stably express the glaze of ice and the curved freezing shape. The key point of this paper is to reliably simulate the anisotropic formation and growth of ice.
Simulating the supercooled water phenomenon makes it easy to model not only the shape of ice growing downwards by gravity, but also the more complex freezing shape. Since the temperature of supercooled water is less than zero degrees, freezing proceeds by crystallization between molecules, not by temperature change. The fluid in this paper is also assumed to be supercooled water, and the crystallization step for the fluid to change to ice is newly proposed based on the nucleation energy proposed by Im et al. [32]. A list of symbols is available in Table 1.

A. IMPLICIT INCOMPRESSIBLE SPH SIMULATION
The IISPH technique is one of the simulation methods to improve the incompressible problem in particle-based fluids [32]. IISPH can discretize the LHS (left hand side) of the continuity equation through the RHS (right hand side) of the basic SPH technique and the first-order accuracy of the VOLUME 9, 2021 Euler method as follows: Dρ Dt = ρ · v. As a result, the continuity equation can be rewritten as (see Equation 1).
where S (i) is the particles adjacent to the i-th particle, ρ i and m i are the density and mass of the i-th particle, and v ij is the difference in velocity between the i and j-th W is calculated as the third Spline interpolation as follows (see Equation 2).
where h is the length of the kernel function, q is defined as |r| h , and r is: r = x ij . In this paper, the radius of the simulation domain is calculated as kh to be proportional to h, and k is set to 2. W ij is defined as W | r=x ij , i is the slope for x i (i.e. ∂ ∂x i ), and is written as W ij in this paper. Also, x ij is the position difference between the ith and the jth particle: , v ij (t + t) and ρ i (t + t) are not known, but the condition ρ i (t + t) = ρ 0 is satisfied to maintain the incompressibility of the continuity equation, and then an equation for pressure is derived based on this equation.
Newton's second law is used to express explicitly the pressure and non-pressure due to the force acting at the position of the ith particle. In this process, the RHS for F i is calculated using the first-order accuracy of Euler method, and F i = m i Dv i Dt and F i are calculated separately for the pressure and non-pressure parts (see Equation 3).
where F p and F np represent the force for pressure and nonpressure, respectively. Using the above equation, the velocity is newly defined as follows (see Equation 4).
Velocity is divided into two components according to the pressure contribution of force as follows (see Equation 5).
respectively, and the subscripts l and l + 1 represent discretized values at time t and t + t, respectively. IISPH is similar to Chorin's two-stage projection algorithm [38]: 1) Mid-velocity u * is related to viscosity and velocity in the previous time-step, and 2) new velocity is related to pressure and mid-velocity in the new time-step. In this paper, v np i,l is similar to mid-velocity u * in Chorin's projection method, and v p i,l+1 affects the pressure gradient of Navier-Stokes equation. In Chorin's projection method, the pressure equation can be obtained by applying a divergence-free condition to the velocity at a new time-step. However, the pressure equation in IISPH can be obtained by applying the incompressibility to the density at the new time-step t + t. By dividing the density into the pressure part and the non-pressure part by Equation 1, we can calculate the mid-density ρ np i,l from the mid-velocity v np i,l (see Equation 5).
Applying Equation 5 to Equation 6 results in Equation 7.
where ρ i,l is ρ i,l = j m j W ij , and ρ i,l+1 = ρ 0 must be satisfied within the limits of numerical accuracy to maintain incompressibility. Using Equation 5 and v p ij,l+1 , the above equation is made in a simpler form, and the pressure equation can be expressed implicitly (see Equation 8).
The pressure is then calculated using the above equation. Equation 9 can be rewritten by applying the ρ np i of Equation 6 as shown below (see Equations 10 and 11).
where v t is used to simplify LHS, and the force on pressure with conserved momentum is calculated as follows (Equation 12).
Then, F np i,l is calculated as follows (see Equation 13).
whereh ij is calculated as:h ij = (hi+hj) 2 , this value is affected by the radius of the kernel, andρ ij is calculated as:ρ ij = (ρi+ρj) 2 . In Equation 11, the values of the LHS variables are known, and the RHS can be further simplified through a simplified version of Equation 12 (see Equation 15). where W ik , and can be rewritten as , the RHS of Equation 11 can be written as follows (see Equation 16). (16) In the above equation, the contribution to the ith particle is extracted from the underlined term in Equation 16, and can be written as follows (see Equation 17). The underlined term in Equation 16 expresses the contribution of water particle i, and expresses the influence of the i-th particle when the pressure is updated with an iterative numerical solution. (17) Since the above equation must be solved with a numerical solution such as the Jacoby iteration method, the RHS of Equation 11 is rewritten as follows (see Equation 18 and 19). (19) where a ii , d ii , and d ik are defined as follows (see Equation 20).

RHS
To update the pressure with an iterative numerical solution, the RHS of Equation 11 is replaced by Equation 19 (see Equation 21).
where the superscripts r − 1 and r represent the previous and new iteration values, respectively, and ρ np i is defined as follows (see Equation 23).
The above equation must be solved repeatedly until the residual of p r i,l+1 becomes sufficiently small. In this paper, the SOR (Successive over relaxation) method was used, and the convergence rate of the iterative method was improved using the following method (see Equation 24).
where w was determined through repeated experiments, and in this paper, we set it to 0.8. Pressure can be calculated with an iterative numerical solution based on Equation 21, and this value is used to calculate F p i,l+1 in Equation 15. Also, F np i,l is calculated using Equation 13. These two values are used in Equation 3 to calculate the new velocity. The particle's position, x i , is updated using the first-order Euler method. Figure 3 shows the IISPH-based fluid simulation results, and we integrate this fluid solution with the freezing simulation to create a new ice growth result.

B. ANISOTROPIC FREEZING SIMULATION
Our system is divided into a fluid simulation step and a freezing solver, and before explaining our framework in detail, the definition of basic elements used in the system is briefly explained in the table below (see Table 2).

1) VIRTUAL WATER FILM ON SOLID SURFACE
Ice grows in two properties, wet and dry, both of which are influenced by the freezing fraction, which affects the rate at which water turns into ice. The freezing fraction decreases as the humidity increases, and the water freezes into ice as the wet diffuse [16]. The freezing process is affected by various physical environments. In general, freezing occurs at a low temperature, but, like supercooled water, when the temperature is already below zero, freezing occurs by crystallization between molecules. In this case, curved freezing surfaces are expressed as they are greatly affected by wind or fluid movement. To change the fluid into ice, a nucleation process is required, and this nucleation energy is generated when the intermolecular forces, which are affected by temperature, humidity and other factors, become stronger. It is assumed that the fluid treated in this paper is supercooled water, and the crystallization of the fluid into ice was implemented by modeling the intermolecular force based on the nucleation energy.
In the dry growth process, the fluid particles are converted into ice immediately after colliding with the object, and in the wet growth process, the fluid particles flow on the surface of the object, forming a water film as ice accumulates. The level-set technique [33] can be used to form the water film, but this method is generally expensive to calculate and is not suitable for particle-based simulation, making it difficult to represent the water film.
In this paper, a virtual water film is added to the solid particles to efficiently express the water film. When creating a virtual water film, , the virtual mass of fluid particles adjacent to the solid particles is considered. As a result, the amount of mass transferred from the fluid particles is calculated taking into account the freezing fraction, α ice|solid (see Equation 25).
where ice|solid is the amount of virtual mass obtained from fluid particles.  Figure 5 shows calculated as water and solid particles collide. As shown in Figure 5a, increases on the surface in contact with water particles, and this value diffuses around the solid particles as the contact surface widens (see Figure 5b).

2) ANISOTROPIC GROWTH DIRECTION OF ICE
Due to gravity, ice grows faster, mainly vertically, rather than horizontally. For example, icicle-like shapes are created by slow flowing water droplets. As the water droplets flow down the surface of the ice under the influence of gravity, they form an elongated shape. Because water droplets are attached to the surface for a long time due to contact force and surface tension, icicles mainly lengthen in the vertical direction, and grow little by little in the horizontal direction due to the thin water film remaining on the surface. Depending on the environment, if an external force, such as a strong wind, acts on a water drop, it is more affected than gravity, so icicles can grow in different forms. In this paper, the growth direction of ice is modeled anisotropically by reflecting the flow of water during rapid water freezing. When water collides with an object, its velocity changes, and the growth direction of ice is calculated based on the velocity of the fluid particles at this time (see Equation 26). (26) where gv ice|solid is the growth direction of the ith solid and ice particles, and v j and S(i) are the velocity of the fluid particle and its adjacent particles, respectively. The linear transformation G rotates and increases the radial vector r, resulting in W A being anisotropic weighted kernel (see Equation 27).
In order to calculate the anisotropic kernel from fluid particles, the covariance matrix is first calculated (see Equation 28), and principal component analysis (PCA) is performed to calculate the orientation and stretching of the fluid particles. In this paper, PCA was calculated through SVD (Singular value decomposition), and the eigenvectors and eigenvalues calculated using this are as follows (see Equation 29).
where R is the rotation matrix calculated using principal axes, and is the diagonal matrix containing eigenvalues. From this value, the anisotropic matrix G is calculated. As a result, the growth direction of ice is calculated using W A , an anisotropic kernel. Because ice forms on the surface of the object, each time a fluid particle collides with the object, the growth direction vector, gv ice|solid , is calculated.
Unlike previous techniques that generate ice only in a specific direction, the proposed method uses gv ice|solid to grow ice anisotropically according to the flow of water. Figure 6 represents the growth direction of ice. As fluid particles collide with the object, freezing proceeds rapidly. In this figure, gv ice|solid , which represents the growth direction, is represented by a blue arrow.

3) PHASE SHIFT
In this paper, we propose a method to calculate the freezing fraction that affects the phase shift of fluid particles into ice particles based on gv ice|solid . If the fluid particles around the solid particles have the same freezing fraction, the ice will grow equally in the same direction and only glaze effects will be expressed. Since we want to control the freezing effect according to the flow of water, the freezing fraction is calculated as follows (see Equation 31).
where fr is the freezing fraction of the fluid particle, and x isf is: It should be noted that each vector in the above equation is a normal vector(x isf ). Equation 31 shows that the freezing fraction of each fluid particle is determined by the angle between gv ice|solid and x isf . When this angle exceeds 90 degrees, fr becomes negative, and ice does not grow in the region where the fr value is negative. In this paper, fr is scaled as follows to adjust its range to [0, 1] so that ice grows in regions with high humidity and positive fr values: fr = fr+1 2 . Nucleation is required to change a fluid into ice, and the nucleation energy required for this process occurs when the force between molecules affected by temperature and humidity increases. In this paper, we model the force between molecules through nucleation energy. Since we are assuming supercoolded water, the simulation is performed under the assumption that the temperature of the object is already below zero degrees. When the object gets wet, a water film appears on the surface of the object, so the humidity is approximated by the virtual water film described in Section III-B1, and then nucleation energy is calculated. As a result, the anisotropic nucleation energy proposed in this paper is calculated as follows (see Equation 32).
where ce x ice i , x represent nucleation energy for fluid and ice particles, respectively, and i represents a virtual water film. Since the nucleation energy for ice growth actually appears differently depending on the surrounding environment, in this paper, the nucleation energy for the fluid and ice particles is independently calculated. In addition, glaze and icicle can be expressed individually by adjusting nucleation energy using glaze threshold and icicle threshold. VOLUME 9, 2021 When the ce of the fluid particles in contact with the solid particles exceeds the glaze threshold gz th , they are converted into ice particles to express the glaze generated on the object surface. Similarly, when the ce of water particles in contact with ice particles exceeds the icicle threshold ic th , they are converted into ice particles. Figure 7 shows the result of changing fluid particles into ice particles. As shown in the figure, as fluid particles collide with solids, rapid freezing proceeds, and changes into ice particles according to the previously described algorithm.

FIGURE 7.
Phase shift from fluid particles to ice particles (red particle: ice, blue particle: water).

4) IMPROVING THE VISUAL QUALITY OF FREEZING SURFACE BY REDUCING INSTABILITY OF WATER PARTICLES
Fluid motion is important in modeling the shape of ice. The unnatural or unstable fluid motion degrades the quality of the freezing simulations applied based on it. Like Im et al.'s method [32], we express the effect of fluid particles clumping together during freezing simulation by adding air pressure and viscosity as well as surface tension to the fluid motion.

a: SURFACE TENSION AND AIR PRESSURE
Im et al. [32] used the surface tension method proposed by Akinci et al. [35], but there was a problem of numerically unstable convergence. In this paper, we improve the fluid solver by utilizing He et al.'s method [31] to better represent the shape of ice particles clumping together and to ensure that simulations converge reliably (see Equation 34).
where V j is the volume of the jth fluid particle, and x ij is the distance between the particles i and j. In the process of calculating the surface tension F s i , the average of the two surface tension energy densities is calculated and used in the surface tension equation to preserve the momentum (see Equation 35).
where κ is the squared gradient energy coefficient [31], [36], which is set to 0.02 in this paper.
In this paper, we add air pressure force F a i to fluid particles to improve the effect of freezing the fluid as it flows along the surface [31] (see Equation 36).
where p atm is the air pressure of fluid particles, which is set to 5000 in this paper. Our surface tension method includes both forces defined earlier: Figure 8 is a comparison between the surface tension used in Im et al.'s method [32] and the surface tension we used. This paper did not propose a new surface tension technique, but improved the instability of the surface tension used by Im et al. [32]. In this paper, the surface tension is calculated in the iterative solver to calculate the pressure in IISPH, and this force is used in the integration process to update the position of the particle.

b: FLUID VISCOSITY
Im et al. [32] used only the surface tension of Akinci et al. [35] to express the curved ice that appears when a fluid flows and freezes. However, in this way, if the velocity of fluid particles colliding with the object is very slow, it sticks to the surface of the object and flows down, so you can model ice glaze or curved ice, but the faster the fluid particle speed, it becomes difficult to model the freezing shape.
As shown in Figure 9 showing the results of Im et al. [32], when the speed of the particle increases, it is expressed as a splash, but the contact force or surface tension of the fluid are not properly exerted. In this paper, F vf i∈fluid and F vs i∈solid , which are forces related to viscosity, as well as contact and cohesion forces due to surface tension, are added to reduce simulation instability and improve the quality of freezing shape (see Equations 37 and 40).
where β 0 and β 1 mean bulk visibility and sound wave respectively, and in this paper, we set them to 0.1 and 10. Since the above equation is a viscous effect only for fluid particles, it can express a sticky flow of water. Since this viscous effect helps to better express the surface tension or contact force of the fluid (see Figure 10), the freezing effect according to the fluid flow is also well expressed (see Figures 10c and 10d). If we observe the actual phenomenon, when the water touches a cold solid, the heat is absorbed by the solid, and the water sticks to the solid surface and flows down. The fluid viscosity is well represented by F vs i∈solid , but it is not sufficient to model ice glaze effects or curved ice. This is because fluid particles should flow down by sticking to the surface and not scattering after colliding with the solid surface. As shown in Figures 10c and 10d, after the fluid collides with the wall, the effect of the fluid particles being sticky and clumping together by viscosity is well expressed, but the movement as if sticking to the wall and flowing down is not expressed. To solve this problem, in this paper, we calculate F vs i∈solid to express additional viscous effects based on boundary particles (see Equation 40).
where j is the volume of boundary particles proposed by Akinci et al. [29], and it is recommended to refer to their study for more detailed volume calculation. Figure 11 is the comparison of viscous effects near the boundary by F vf and F vs . Compared to Figure 11a, which is the result of applying only the volume of the fluid, Figure 11b is more sticky and flowing down the wall or arm of the armadillo. When a fluid collides with a solid and freezes, a noisy surface is often created when fluid particles are scattered (see Frame 51 in Figure 9), but our method can stably express glaze effects of a solid surface or curved ice by applying viscosity near the boundary.

5) SURFACE RECONSTRUCTION FROM ICE PARTICLES
When reconstructing the ice surface, the SVD value used to calculate the anisotropic direction of ice growth is used. In this paper, the ice surface was reconstructed in detail using the anisotropic property of SVD in the algorithm proposed by Yu et al. [34], and the level-set of ice particles is calculated as follows (see Equation 41).
where D ice i represents the transformation matrix for the ice particle i, whose properties are shown below (see Equation 42).
where k s is the scaling constant, limiting σ n to within a certain range to preserve the stretching with the greatest variation. For a more detailed explanation, it is recommended to refer to the work of Yu et al. [34]. In this paper, the GPU-based Marching Cubes algorithm [11] is used, and the minimum stretching value is set to be larger than half the width of the grid so that the ice surface is well represented.

Algorithm 1 Pseudo-Code of Our Algorithm
For each frame Do // Base water solver Advect 3D water partices using IISPH // Section III-A Compute the interaction of water & solids // Freezing water solver Viscosity force // Section III-B4 Surface tension & air pressure forces // Section III-B4 Virtual water film // Section III-B1 Anisotropic direction for growthing // Section III-B2 Phase transition from water to ice // Section III-B3 Ice surface reconstruction // Section III-B5 End For

IV. IMPLEMENTATION
The experiments in this paper were implemented in the following environment: Intel i7-7700k 4.20GHz CPU, 32GB RAM, NVIDIA GeForce GTX 1080Ti graphics card. The IISPH-based fluid solution was implemented in the GPU environment and used as an underlying water simulation [28]. The Marching Cubes algorithm was used to reconstruct ice surfaces by adding a grid of 200 × 200×200. The boundary particle method proposed by Akinci et al. [29] was used to handle the collision between fluids and solids. Pseudo-code of our foam simulation algorithm is shown in Algorithm 1.
Autodesk 3ds Max was used for high quality rendering for the final result. Since our results mainly deal with freezing surfaces, we did not render opaque regions or bubbles of ice. In this paper, rendering was performed only on ice surfaces, and rendering quality was improved using opacity, bump, reflection, refraction, and displacement mapping.

V. RESULTS
In order to analyze the technique proposed in this paper from various viewpoints, we experimented in various scenarios, and it is explained in a simulation view and a high quality rendering section.
A. SIMULATION VIEW Figure 12 shows the result of rapid freezing when a spherical liquid is dropped onto the Stanford Bunny model. For the simulation of this scene, the time-step was set to 0.006 and approximately 50,000 fluid particles were used. As mentioned earlier, with Im et al.'s method [32], the spherical liquid is scattered when it collides with the object, so that the virtual water film cannot be properly formed on the solid surface, and as a result, ice cannot be expressed (see inset image in Figure 12a). Even in their method, there is a region where ice is expressed, but because the fluid is scattered, it is more like noise than freezing while flowing through the object surface. On the other hand, the proposed method well expressed the glaze effects of ice produced on solid surfaces. Not only that, but also curved ice that freezes along the solid surface was well represented (see Figure 12b).      [32] did not properly express the freezing effect due to the unstable movement of fluid particles near the boundary (see Figure 13a), but the proposed method stably expressed not only the freezing near the boundary but also the thin sheets of ice expressed during rapid freezing (see Figure 13b). Figure 14 shows the result of rapid freezing when a spherical liquid is dropped onto the Armadillo model. As it freezes downwards by gravity, it is affected by the movement of the fluid and freezes in a zigzag shape. At the top of the  Armadillo model, the glaze effects of ice are expressed, and in Frame 72, ice thin like an icicle is being created. In addition, it freezes and spreads directionally according to gv, which is the anisotropic growth direction of ice (see red rectangle in Figure 14). Figure 15 shows that when the liquid is dropped, it freezes rapidly, similar to the results shown earlier. Unlike the unstable freezing modeling by Im et al.'s method [32] (see inset image in Figure 15), the proposed method realistically represents the glaze effects of ice and the ice surface frozen by the flow of fluid.
As seen in the results shown earlier, the proposed technique showed more stable and realistic freezing quality than the state-of-the-art freezing technique, Im et al. [32].

B. HIGH QUALITY RENDERING
In this section, we look at how water particles are frozen by throwing them in various directions. Figure 16 is a scene where water particles fall downward, and glaze effects and curved ice surfaces are well represented on the head and back of the Armadillo model. In this scene, gz th and ic th were set to 0.004 and 0.0001, respectively. Since water was dropped from the upper middle of the Armadillo model, the ice surfaces spread isotropically and freeze after the impact on the back. Figure 17 is a scene in which water particles are thrown in a diagonal direction from the top left, and curved ice surfaces appear well on the head and back of the Armadillo model as in the previous result. Compared to the scene in     Figure 16, the gz th was set to be larger so that glaze effects were less expressed, and an icicle with a relatively anisotropic shape appeared. In particular, the curved ice shape, which is well seen in the supercooled water phenomenon, was clearly shown. This shape is frozen in the red arrow direction, as shown in Figure 17, because it is affected by the throwing direction of the water. Figure 18 shows the result of throwing water particles in various directions. We rendered from various viewpoints to observe whether the ice surfaces were well represented according to the throwing direction of the water particles. Simulation stability is essential to create such a complex scene, and Figure 18 shows the glaze effects and icicles of ice surfaces stably. Figure 19 is a scene where water particles are dropped on the rotating Happy Buddha model, and the ice surfaces are well represented not only on the surface of the 3D model, but also on the ground. In particular, in this scene, we tried to express the fluctuating movement of water by rotating the Buddha model (see FigureFigure 19). Due to this, irregular patterns appeared on the freezing water surfaces, unlike the results shown earlier, and the ice surfaces connected from the surface of the Buddha model to the ground in the form of a bridge were well expressed. Table 3 shows the parameters and scene configuration used to produce the experimental results in this paper. The user can control the freezing surfaces using the corresponding parameters. The growth of the glaze is more constrained as the gz th increases. gz th is set lower than the ic th to show the ice glaze more clearly. As the icicle threshold becomes larger, the icicle grows more slowly.

VI. CONCLUSION AND FUTURE WORK
In this paper, we propose a framework for stable and realistic representation of glaze effects and curved ice surfaces expressed in supercooled water phenomena. The growth direction of ice was calculated anisotropically according to the flow of the fluid, and the stability of freezing simulation was improved by adding surface tension, air pressure, and viscosity to the movement of the fluid to express the glaze effects and curved ice surfaces. Unlike the existing technique that expressed rapid freezing only in specific scenes, the proposed method showed stable and improved freezing results in various scenes.
In this paper, the freezing shape was reliably modeled by integrating particle-based water simulations and a freezing solver. Our method is physically inspired because it uses a physics-based fluid solver, but it cannot be called a complete physics-based simulation technique because pressure and heat transfer were not considered in the process of changing the phase from water to ice. Nevertheless, as shown in the results, it shows stable integration with particle-based fluids. As mentioned in the problem statement, in the previous technique, it is difficult to model the ice shape of a curved surface because only gravity is affected. We think that our technique can be improved by conserving physical momentum as in the Stefan problem, or by using the existing physics-based approach that can accurately solve the boundary problem between phases during phase transition.
In this paper, humidity is calculated when water particles and solid particles come into contact, and a phase shift is carried out due to crystallization-based energy''. However, in the process of particle contact, heat transfer to each material was not considered. Since the algorithm was initially designed on the premise that the main target is supercooled water, the accuracy or results may vary because heat transfer is not considered. We will come up with a supplement to this problem in the future.
In the future, we plan to study how to efficiently express freezing simulation using the flow of surrounding air or wind, not fluid particles. The flow of fluid is an important factor in expressing the detailed curved ice surface, but because of the large amount of calculation, we will study how to efficiently process it. In addition, we will study how to physically express glaze effects caused by water or ice sticking to cold objects based on heat transfer.