{'fraction_non_alphanumeric': 0.028055796359051147, 'fraction_numerical': 0.011821262510836157, 'mean_word_length': 4.875, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Modelling of localised shrinkage pressure in a directionally solidifying domain

A novel approach to modelling shrinkage during solidification has been proposed. It combines Phase-field (PF) simulation with fluid mechanics for a solidifying medium. A mass balance method, based on a viscous creeping flow in two phases (liquid, solid and interface) with an applied penalisation coefficient is proposed. The role of this coefficient is to continually separate the fluid from the solid domain in correspondence to any change due to a phase transformation process. The advantage of the method is that the medium is modelled with no explicit tracking of the interface which significantly improves computational cost and the ease of implementation. This is achieved via the application of a monolithic ‘single-field’ formulation to the system of equations. The method has the potential to locally calculate shrinkage pressure by taking into account multiple parameters, including thermal gradient, permeability, viscosity and morphology.


Introduction
During solidification, a metal melt undergoes a transition from liquid to solid under the influence of temperature and/or composition gradient. This transition generally occurs by the nucleation and growth of dendrites. Heat flow and crystallographic anisotropy dictate the most energetically favourable directions of dendritic growth and are the primary cause for the formation of complex dendritic networks. Deep within those networks, between the solid and liquid phases, an interface exists in which a drastic variation in behaviour and thermo-physical parameters such as density, viscosity and motion are present. In particular, the density variation across the interface causes shrinkage to occur due to mass conservation between both phases. This effect can generate situations where liquid pressure is substantial enough to cause deformation of the already solidified melt.
The deformation of the dendritic networks strongly depends on the coherency state and the inflow and outflow of liquid. The two most common defects formed due to shrinkage during solidification, porosity and hot cracking, are both contributed to the hydrostatic depression of the semi-solid (mushy) zones [1]. This depression is caused by the insufficient liquid feeding in the shrinking mushy zones. Due to the complexity of the mechanisms involved in porosity and hot-tearing formation, most of the developed models are semi-empirical or based on the solidification interval, most popular of which are Niyama criterion [2,3] and RDG criterion [4]. In those models, solidification is generally approximated via Scheil equation and mass flow is taken into consideration via Darcy's law with an empirical determination of the permeability term [5,6]. However, permeability, which is determined from the tortuosity of the dendritic networks and the solidification interval, is very sensitive to the solidification conditions. These semi-empirical models are also highly sensitive to the coherence value of the array that can be in the range from 0.9 to 0.99 volume fraction of solid [7]. There is a large ongoing research IOP Conf. Series: Materials Science and Engineering 861 (2020) 012020 IOP Publishing doi:10.1088/1757-899X/861/1/012020 2 on multiphysics solidification models but most of them mainly focus on natural and forced convection flows during solidification [8,9], with shrinkage effects ignored. The current understanding of the effects of fluid flow and pressure deep within the mushy zone is still not completely understood, especially when predictive methods are attempted to be developed.
The aim of the current work is to demonstrate a method for fluid-pressure calculation specially aimed for domains undergoing phase transition. The formulation coupled with a micro-solidification model such as Phase-field (PF), allows for local calculation of flow and pressure within the dendritic arrays. First, a simple PF model capable of dendritic solidification is built and then, a methodology for fluid dynamics is described. Shrinkage is implemented within a growing array via coupling the PF field variable and the Navier-Stokes equation for local mass conservation.

Phase-field
In phase-field simulations, the sharp interface is replaced by a diffuse interface with a finite thickness, where the phase-field variable φ changes smoothly from φ = 0 in liquid to φ = 0 in solid. Such an interface enabled us to simulate the interface migration without tracking the interface. The significance of this method for modelling dendritic growth was brought to light by Kobayashi [10] and it still remains one of the most powerful and versatile methods for simulating growing structures. Kobayashi's model dealt with a very thin steep interface layer in a dimensionless relationship between the field variable and the physical parameters, latent heat and anisotropy.
The computations presented in this work are based on his model in order to model two impinging solidifying fronts. Kobayashi's formulation offers the complexity of dendritic arrays but short computational times. Pure material solidification with incorporated interface anisotropy is simulated by the following evolution equation: where is the driving force for solidification, -anisotropy and ε ́ = ∂ε/∂θ.
The dimensionless equation describing latent heat release during dendrite growth of a pure material is: ∂T ∂t = c∇ 2 T+K ∂ϕ ∂t (2) where K is a dimensionless latent heat, proportional to the undercooling and c is diffusion constant. All simulation parameters and operating conditions are listed in table 1 [10].

Fluid flow
Fluid flow formulations generally require explicit tracking of the interface or complicated meshing techniques such as level-set [11], IBM [12] or ALE [13,14]. In the current work, however, a monolithic solution for prescribing interface boundaries is developed by borrowing an idea from the area of topology optimisation problems [15][16][17]. The unique idea is to model the entire domain as fluid but enforce a location based damping term in order to simulate a rigid body surface, as shown in figure 1. Numerically, this causes the velocity of the fluid to approach zero at the interface, implicitly enforcing 'no-slip' condition. This term can be included into the steady-state Navier-Stokes equation as follows: where μ is viscosity, ρ -density. The term 'αv' is a 'volume' or 'friction' force and the coefficient α can be treated as a penalisation constant: In this work, the interface is considered as a part of the liquid domain with regards to the coefficient α. However, a more in-depth research can be performed to investigate the relationship between k and μ. Research already done on this area by Guest et al. [17], showed that k is inversely proportional to μ, similarly to Darcy's equation. However, no investigation has been performed on solidifying structures. The shrinkage of the moving interface is expressed in the mass conservation equation via the change of phase field variable φ with time, representing the change in solid fraction in the domain. β is shrinkage coefficient and it is a function of the solid and liquid density ratios. Since only a localised interdendritic flow is modelled, convection is assumed negligible.
As shown in figure 1, the penalisation constant determines the phases in the Navier-Stokes equation. There is no need to decouple the equation into solid and liquid parts, since a solution is obtained only for the fluid regions, whereas the solid remains rigid.

Configuration
The model is set to simulate two solidifying fronts growing opposite to each. A square domain of size L0 = 8×8 was created, with cell size h0 = 0.08, and a time step of t0 = 0.002 was set. As the entire formulation is non-dimensionalised, the time step and the grid size are artificial. The domain and the main boundary conditions are shown in figure 2. Symmetry boundaries are applied to the phase field equation. Heat is extracted from the left and right side of the domain with heat flux of q0 = -2, where the top and bottom boundaries are set to be symmetric. For the Navier-Stokes' equation, fluid is allowed to freely leave or enter at the top and bottom and 'no-slip' condition is imposed on the rest of the boundaries. Two seeds with radial size less than 3 times the grid size are positioned as shown in the diagram below. Phase-field, heat equation and Navier-Stokes formulation are solved via finite element analysis (FEA) with a linear Lagrange discretisation on COMSOL Multiphysics® software. The total simulation time was 11 min 26 s on a university PC (Intel® Core™ i5-8500 CPU @3.00GHz, 16GB RAM).

Figure 2. Schematic of computational domain and boundary conditions
The fluid flow is simulated on the same grid as the PF simulation, hence no interpolation or multi-level approaches are needed. The coefficient α is set as a function of the field variable ϕ where α(ϕ>0.99)=10 4 is considered solid and α(ϕ=0)=0 is liquid. The magnitude of α determines how impervious a material is. If αs is too small, the solid domain allows fluid to enter while if αs is too large, a numerical singularity would ensue. An optimum parameter exists such that solid permeability is optimised. This parameter has been investigate in literature for different physical processes and was found to be dependent on the mesh size as α=1/(0.05hx 2 ), for square grid [17][18][19]. It was found that this value does indeed provide a good first approximation but also that a range of values exists spanning across several magnitudes. Any of these optional values for α, cause an insignificant difference in the solution for pressure.

Results
The four main processes heat transfer, phase-field, fluid flow and pressure are plotted in colour map format in figure 3. Three time steps are chosen, representing the three distinct stages from the process of impinging dendrite fronts. Up to 0.3 s both solidifying fronts have grown freely in the domain. High solidification rates and high velocity can be observed. It is important to note that at this stage, low pressure zones occur both at the tips of the dendrites and deep into the network. This can be contributed to the high solidification rates which give rise to high shrinkage at the dendrite tips. However, these pressure zones are easily accessible by fluid, hence pressure differences are relatively low. At 0.6 s, the two fronts are impinged. At this stage, heat has not been dissipated yet and the fronts have little space to grow, leading to reduced solidification rates. Fluid velocity, caused by shrinkage, also decreases and mostly occurs at the zones close to the boundaries. Pressures deep into the dendritic networks have begun decreasing. The pressure depression in these zones occurs due to the lack of feeding caused by their remoteness from the open boundaries and complexity of the channels. Results are also plotted at 2.6 s in order to represent the final stages of solidification, where the domain is 98% solidified. The channel no longer exists and only a few liquid pockets remain in the now solidified melt. There is barely any flow from the open boundaries as they are cut off by a solid phase. Important to observe is that the pressure depression in the liquid pockets is 5-6 orders of magnitude lower that in the other two cases. Such low pressures can potentially lead to pore formation and hot cracking.
Solid fraction and normalised fluid pressure in the entire domain, both plotted with time, are shown in figure 4. The three stages explained above can clearly be observed. Up to 0.3 s, the domain solidifies rapidly until it reaches 78% solid fraction. Barely any hydrostatic depression occurs due to the small pressure differences between tips and network zones. At 0.5 s however, when most of the dendrite arms impinge, a sudden drop in pressure can be observed. After this a steady-state occurs, where the domain slowly solidifies and pressure slowly recovers. At 2.6 s or 98% solid fraction, a sudden drop in pressure is observed. This, as already explained, is due to the cut-off channels and the formation of liquid pockets inside the solid.
The above described processes have been discussed in literature and proven experimentally, however, it is promising to see that our simulation results can clearly divide them into separate solidification stages. The advantage of the method, once applied to real alloy system, is that the exact time and solid fraction for each stage can be obtained. Furthermore, the exact positions for pore formation can investigated and identified with regards to depreciation pressure within dendritic arrays. (b) (a)