MODELLING PHASE CHANGE IN A 3D THERMAL TRANSIENT ANALYSIS

A 3D thermal transient analysis of a gap profiling technique which utilises phase change material (plasticine) is conducted in ANSYS. Phase change is modelled by assigning enthalpy of fusion over a wide temperature range based on Differential Scanning Calorimetry (DSC) results. Temperature dependent convection is approximated using Nusselt number correlations. A parametric study is conducted on the thermal contact conductance value between the profiling device (polymer) and adjacent (metal) surfaces. Initial temperatures are established using a liner extrapolation based on experimental data. Results yield good correlation with experimental data.


INTRODUCTION
Manufacture and assembly of structures is a complex process wherein the final stages of assembly gaps may be found between mating surfaces of components that are to be mechanically-fastened [1]. This fact is inherent in all manufacturing and assembly processes as the control over such processes are always finite [2], however the variations are reported to be larger in components manufactured using composites as opposed to traditional engineering alloys [3]. The gaps arise due to an individual component's manufacturing tolerances, which when placed in their respective final positions present significant deviations from their nominal dimensions due to the summation of the individual variations. These gap heights typically range from 0 to 2.5 mm.
The Gap Sachet ( Figure 1) is a device created for the purpose of profiling gaps arising between mating faces of components and is intended for use in conjunction with one of the many 3D surface digitizing technologies for digital mapping [4]. The device works by injecting a moulding material, plasticine (in its liquid state), into a fabricated thin plastic film (nylon) retainer which is placed into the gap to be profiled. The injection is achieved by a heated motorised-extruder unit, having fixed mass flow rate and injection temperature. The plasticine within the Gap Sachet is allowed to cool and upon hardening it is removed providing a sectional profile of the gap.
This paper aims to model in ANSYS v11.0, the various transient thermal processes involved during the cooling of the Gap Sachet, including the phase change of plasticine and the effect of thermal contact conductance between a polymer-metal interface. The thermal effects taking place during the injection process is approximated and applied as a time dependent thermal load to describe the initial conditions. The results obtained are compared with experimental data.

PROBLEM SPECIFICATION & BOUNDARY CONDITIONS
The gap to be profiled may take on various shapes and sizes and are typically located between components that are required to be fastened; consisting of a wide overhead cover (Top Plate) which attaches onto a smaller structural support component (Bottom Plate) in order to maintain its shape. For the purpose of this paper, the gap in between the two components is assumed to have constant rectangular cross section and a height of 2.5 mm. Figure 2 shows the front and top views of a filled Gap Sachet which has been placed into a typical idealised gap model, indicating the thermal processes involved during cooling. Four main bodies are established; (1) Plasticine (2) Plastic Film (3) Top Plate and (4) Bottom Plate, whose material properties are listed in Table 1.
The ideal placement of the Gap Sachet for gap profiling is along the centre line of the Top and Bottom Plates thus having geometrical symmetry along the z and x-axes, however it is assumed to have thermal symmetry only along the z-axis (axes detailed in Figure 8). Therefore, only half the system is to be modelled.
The computational model of the Gap Sachet is assumed to encompass the entire gap within its designed volume (30 x 7.5 x 2.5 mm, length x width x height) and is defined as two components. A rectangular hollow block (thin film plastic, wall thickness 0.0508 mm) which surrounds a solid block (plasticine). The idealised Top and Bottom Plates have dimensions 50 x 25 x 6 mm and 30 x 15 x 6 mm respectively. Four measurements are made (T 1 , T 2 , Top & Bot) as indicated in Figure 2 whose positions are detailed in Table 2.  (3) Radiation. The rate of conduction is defined by the material properties for each respective body. Convection effects are computed by approximation of the fluid film coefficient with respect to surface and bulk fluid temperature, assuming natural convection. All faces are assumed to be insulated limiting the effects of convection being applied on the exposed face (henceforth referred to as the conv 1 face).
Also, the Top and Bottom Plates are also considered to be insulated on all sides except the region within which the Gap Sachet is to be inserted (Figures 2 and 3). For the experiment, due to the size of the insulation, combined with its overall dimensions a channel of depth 15mm exists perpendicular to the conv 1 face. However, as the insulating material does not produce any heat flux, it is assumed to play no role in the rate of transfer of heat from the exposed face. Further, the injection hole is not modelled in the simulation to reduce the computational costs of the analysis. The effects of radiation are assumed to be negligible.
As the plasticine is to be injected in its fluidic state into the Gap Sachet using the motorised-extruder, during cooling the plasticine mould solidifies taking the shape of the gap. Modelling of this phase change (liquid to solid) is the main aim of the computational analysis through application of an enthalpy based model which describes the energy released   during cooling (latent heat of fusion). The model's accuracy is of interest in this paper and its application to plasticine will be studied. It is well known that for surfaces in contact, perfect conduction is never achieved between the two bodies. A temperature discontinuity exists at the mating interface due to the microasperities on the surface of both components. The thermal contact conductance is a result of contact only being made at a select number of points, as opposed to the entire surface area of the bodies in contact [5].
One method of computationally obtaining the value of thermal contact conductance is based on curve-fitting of the experimental temperature readings (vs. time) at specific locations within a system with the computational results. This is achieved by varying the value of thermal contact conductance across a wide range until a good match is obtained [6]. While it is clear that experimental measurements of thermal contact conductance yield suitable results for specified parameters, the aforementioned indirect method is commonly used in conjunction with FEA packages where only the thermal loads are known [7]. A parametric study is conducted for the values of thermal contact conductance between the range 200 to 800 [W/ K] in steps of 200.
The motorised-extruder has fixed initial temperature (T Injection = 85 [ o C]) and constant mass flow rate of 1.326 [g/s], taking approximately 1.7 [s] to fill the Gap Sachet with plasticine. The flow of plasticine from the extruder into the Gap Sachet is not modelled. However, the temperature profile generated as a result of the flow into the Gap Sachet is approximated from experimental data and applied as initial conditions. The profile is derived by linearly extrapolating temperature values at 2.5 mm intervals between thermocouple readings T 1 and T 2 ( Figure 4) and applying them within quadrants in the computational model, assumed constant along the cross-sectional profile (x-axis). As the duration to fill the Gap Sachet is small in comparison to the overall length of the simulation, it is assumed that the effect of convection within the liquid plasticine is small and that conduction is the dominant heat transfer mode throughout the temperature range of plasticine.
Initial body temperatures are averaged and applied for the computational analyses. Temperature results are extracted from points T 1 , T 2 , Top & Bot and compared with experimental data. Using the Galerkin Weighted Residual Method [11] to integrate around the volume of an element (e), while taking into account the boundary conditions (spatial temporal dependence of temperature, convection and heat flow) [12] and substitution of the shape function of the elements [13], for FEA analysis in ANSYS [10], Eqn. (1) may be written in the following matrix form: Where the subscript e signifies the matrix defining the respective element, [ ( )  Various methods exist to numerically estimate the value of h f for the purpose of computational simulation. One common method is the application of the appropriate Nusselt number correlations to simulate the aid in cooling achieved through the exposed surface due to convection. The procedure involves establishing fluid properties for a range of working film temperatures and determining the corresponding convective parameters, allowing the computation of h f [16][17] [18][19] [20].
The fluid (air) properties are established within a working range of T ∞ < T < 85 [ o C] using correlations developed for dry air at one atmosphere by Ref [21]. The Nusselt number (Nu) expresses the overall heat transfer phenomenon and is defined as the ratio of the convective heat flow to the heat transferred by conduction [17]: Where A is the area [m 2 ] and P is the perimeter [m]. The Nusselt number [17] is described by two dimensionless parameters of the form: Where Gr is the Grashof number, the ratio of the buoyancy forces to viscous forces, Pr is the Prandtl number, the ratio of momentum and thermal diffusivities and m and n are exponents obtained from experiments. The Grashof number [16] is defined as: Where g is the acceleration due to gravity [m/s 2 ], β [1/ o K] is the thermal expansion The product of Gr and Pr defines the dimensionless constant Ra, the Rayleigh number [16]: (9) The Rayleigh number indicates the dominant heat transfer mechanism, whose value when below the order of 10 3 indicates dominant heat transfer through conduction, while an increasing value signifies the takeover of heat transfer by convection [25]. Empirical correlations between Nu and Ra for flat plates [16] are defined in the following form:

PHASE CHANGE
Plasticine is to be injected into the Gap Sachet in its liquid state and allowed to cool to a hardened solid state. Phase change transition initiates when the amplitude of the crystal lattice particles oscillate at a force having value larger than that of the crystal binding energy, thus breaking its bonds and transforming into liquid phase (melting). The removal of energy results in the solidification of the material (crystallization) [27]. The process of phase change during cooling exudes latent heat energy and is taken into account in ANSYS by employing the enthalpy model which defines the enthalpy as a function of temperature. The process is broken down into three phases, (1) solid phase (2) liquid phase and (3) mushy phase. Enthalpies are calculated for different temperature points, based on the latent heat energy dissipated, identifying the three stages detailed [28][29].

Determining Phase Change Enthalpies
Eqn. (2) details the three necessary properties required to be input into the governing equation for computation of a 3D thermal analysis, namely the density (ρ), specific heat capacity (c p ) and thermal conductivity (k). is the element specific heat matrix [J/kg-K] which is defined as: [10] (12) Where {N} is the element shape function, defined as N(x,y,z). The (ρc) terms in the equation above define the enthalpy (H, [J/m 3 ]) through: (13) Thus for the enthalpy method, Eqn (13) is used to describe the relative value of enthalpy in the three states of plasticine (solid (H s ), mushy (H m ) and liquid (H l )) with respect to temperature [30]: Where T s [ o C] is the solidus temperature, T 0 [ o C] is the lower limit reference temperature (taken to be 0 [ºC]), T l [ o C] is the liquidus temperature, T + [ o C] is the upper limit reference temperature (taken to be 90 [ o C]), and c p(m) is the average of the solid and liquid specific heats. Eqn. (14) is applied within the temperature range (T < T s ), Eqn. (15) between (T s ≤ T ≤ T l ) and Eqn. (16) between (T > T l ). Eqns. 14 to 16 are integrated [31][32] to obtain: Where c p (*) is defined as: (20) Eqns 17 to 19 can be used to define the enthalpy [J/m 3 ] vs. temperature [ o C] curve for the plasticine, after defining the limits of temperature.

Phase Change properties of Plasticine
Plasticine is a non-linear, viscoelastic, strain-rate softening/hardening material [33], whose exact composition is unknown, but is known to be composed primarily of calcium carbonate, paraffin wax and long-chain aliphatic acids. Due to the presence of paraffin wax, Plasticine is found to exhibit a phase change phenomenon which corresponds to the crystallization point of paraffin wax, allowing it flow in a liquid state beyond this temperature [34].
Due to the rheology of plasticine, it is commonly used as an analogue material to simulate deformation of geological structures [35], material flow in friction stir welding [36] [37] and material forming processes [38][39] [40]. The non-linear property of plasticine derives from structural changes in the material (through re-orientation of filler chains) and is reported to show softening beyond 200 [ o K], similar to the glass transition phenomenon observed in glass [41]. The physical and thermal properties of plasticine is known to vary across the different brands, and also between the different colours it is manufactured in [35].
It is understood that in crystalline substances melting and solidification take place at the same temperature (T s = T l = T m ), typically assigned to be within 1 [ o C] before the melting peak, while for amorphous substances the phase transition temperature is not located at one point, but takes place across a range of temperatures. (T s < T m < T l ) [42][43] [31]. As such, a wide mushy zone is assigned, assuming T 0 = 0  Table  3, and the resulting Enthalpy vs. Temperature chart is shown in Figure 7.

THERMAL CONTACT CONDUCTANCE
When two nominally flat surfaces make contact, due to the microscopic irregularities on both surfaces contact is only made at a select number of discrete points. As a result, a temperature discontinuity exists at the mating interface where heat transfers not only through conduction at contact points, but also through convection and conduction of the fluid in the interstitial gaps, and radiation [5] Where q [W/m 2 ] is the heat flux (Q) per unit area (A) and ∆T [K] is the temperature drop at the mating interfaces. The value of h c is dependent on parameters such as the surface roughness, hardness, interfacial pressure, thermal and physical material properties, and thickness [6].
Numerous methods exist to obtain experimentally and numerically, the value of thermal contact conductance between numerous interfaces among a for a wide range of the aforementioned parameters [45][47] [48] [49][50] [51]. Results from [52] for tests conducted Modelling Phase Change in a 3D Thermal Transient Analysis It is assumed the value of thermal contact conductance are the same on both sides of the Gap Sachet's mating faces and that the faces are always in contact, from the point the Gap Sachet has been filled (static model). Further it is assumed that pure conduction from the plasticine into the plastic film takes place, and that the resistance due to heat flow is minimal as the plasticine is injected in its liquid state, thus being able to fill the space inside the Gap Sachet confidently, leading to perfect contact.

ESTIMATION OF INITIAL THERMAL LOADING
The initial temperature of plasticine is approximated based on experimental measurements at points T 1 and T 2 , after the injection is complete. The intermediate temperatures are derived using a linear extrapolation of temperatures between the aforementioned points, using: (22) Where T [ o C] is the temperature to be estimated at a point x [mm] along the length of the Gap Sachet, C is the y-intercept set as T 1 and m is the gradient obtained from: (23) A total of 12 intervals (quadrants) are defined from the ranges of 0 [mm] to 30 [mm] at intervals of 2.5 mm. Eqn 22 is evaluated from x = 2.5 [mm] to x =27.5 [mm] and the results obtained (T [ o C]) are used to define the temperatures within the respective quadrants. (Figure 4).

ANSYS TRANSIENT ANALYSIS
The simulation aims to portray in ANSYS, the transient thermal effects of the Gap Sachet during cooling. With the incorporation of Enthalpy (material non-linearity) and thermal contact conductance (contact non-linearity), the transient analysis becomes a non-linear one. The Newton-Raphson method of iteration is employed using the sparse solver to obtain converged results at every sub-step [43]. A summary of the time stepping regime is detailed in Table 4. = + T mx C Using Eqns. 22 and 23, the estimated thermal loads are applied within the plasticine volume as a ramped load for Load Step #1, and is then removed in Load Step #2 & #3, simulating the thermal loads transferred from the plasticine during filling. The iteration convergence tolerance limit (TOLER) was set to the default number of 0.001, to be achieved by the 12 th iteration. Minimum sub-step sizes of 0.1 [s] were found to produce converged and stable results [53] [54].
The SmartSize meshing tool is used to mesh the four bodies created, using mesh setting of 5. For the purpose of the analyses SOLID70 elements (tetrahedral) were used. The incorporation of thermal contact conductance within the computational model involves the application of Target and Contact elements, defined through the Thermal Contact Wizard. TARGE170 and CONTA174 elements were used to define the 3D target and contact bodies (volumes). Figure 8 shows the 3D meshed model, comprising a total of 44,439 elements and 9,227 nodes respectively. Figure 9 show the results of the experimental measurements and computational analyses at points Top, Bottom, T 1 and T 2 for a range of thermal contact conductance values for T injection = 85 [ o C].

RESULTS & DISCUSSION
The Top, Bot and T 1 results show similarities with experimental results in terms of temperature profiles, with the exception of T 2 . As the injection hole has not been modelled on the conv 1 face, the plasticine within the Gap Sachet is not directly exposed to convection in the computational model, thus under predicting the drop in temperature after filling of the Gap Sachet has been completed at 1.7 seconds. Due to the low thermal conductivity of plasticine, this idealization does not manifest into large errors along the length of the Gap Sachet.
The initiation of phase change is seen at 51 [ o C], characterised by a sharp change in slope of the temperature profile of T 1 and T 2 . Good agreement is seen for the curves of thermal contact conductance 400 [W/m-K] with experimental data indicating the true thermal contact conductance to be within this range. However due to the simplification of the Gap Sachet model the computational model maintains a larger volume of plasticine in comparison to the experiment and is therefore expected to over predict results. Further, the enthalpy model  [55] that enthalpies of melting/fusion vary slightly in value and initiation and end temperatures, introducing further errors within the model. Figure 10 shows the Nodal DOF temperature plots at different times for the thermal contact conductance value of 400 [W/m-K], where between t = 0 [s] to 1.7 [s], the linearly ramped thermal load can be seen to maintain the temperature distribution within the defined quadrants during the filling stage in the first Load Step. The subsequent images between t = 5 [s] to 50 [s] shows the cooling of the Gap Sachet in second and third Load Step. The conv 1 face is located on the right hand side of the figures. Figure 11 shows the Nodal Flux Vector Sum at 1.7 seconds, where it can be seen the heat flux has initiated during the loading. At 10 seconds, the heat flux is seen on both top and bottom plates, with peaks located at corners of the Gap Sachet model. This idealization of the Gap Sachet model effectively increases the contact area between the top and bottom faces of the Gap Sachet, therefore the results of the heat flux are simply presented for a qualitative understanding of the flux field.

CONCLUSION
The paper aimed to simulate the transient thermal processes involved during the cooling of the Gap Sachet when placed between two aluminium components. Phase change of plasticine was considered by applying the energy released during solidification over a wide mushy zone. Temperature dependent natural convection on one face was approximated and simulated through Nusselt number correlations, based on surface temperature. A parametric study was conducted on the value of Thermal Contact Conductance between the Gap Sachet and the adjacent contacting faces. Temperature vs. time profiles were shown to show good corroboration with experimental results, except for the exposed face, wherein the injection hole was not modelled, disallowing direct contact of the plasticine within the Gap Sachet with atmosphere. However, due to the low thermal conductivity of components of the Gap Sachet, the inaccuracy due to this simplification is observed to be minimal.