Finite element formulation to study thermal stresses in nanoencapsulated phase change materials for energy storage

Abstract Nanoencapsulated phase change materials (nePCMs) – which are composed of a core with a phase change material and of a shell that envelopes the core – are currently under research for heat storage applications. Mechanically, one problem encountered in the synthesis of nePCMs is the failure of the shell due to thermal stresses during heating/cooling cycles. Thus, a compromise between shell and core volumes must be found to guarantee both mechanical reliability and heat storage capacity. At present, this compromise is commonly achieved by trial and error experiments or by using simple analytical solutions. On this ground, the current work presents a thermodynamically consistent and three-dimensional finite element (FE) formulation considering both solid and liquid phases to study thermal stresses in nePCMs. Despite the fact that there are several phase change FE formulations in the literature, the main novelty of the present work is its monolithic coupling – no staggered approaches are required – between thermal and mechanical fields. Then, the FE formulation is implemented in a computational code and it is validated against one-dimensional analytical solutions. Finally, the FE model is used to perform a thermal stress analysis for different nePCM geometries and materials to predict their mechanical failure by using Rankine’s criterion.


Introduction
One of the major concerns that society faces currently for its development is producing and supplying energy. In fact, evolution of mankind has been closely related to a progressive increase in energy consumption through history [1]. Therefore, research in energy production appears to be crucial for society. Concerning the production of energy, two different paths seem to arise: searching and exploiting new sources of energy or optimizing the existing facilities of energy production processes to gain in efficiency. In connection with this last alternative, a considerable amount of research in thermal energy storage is being carried out [2][3][4][5]. More precisely, in this field, heat storage systems based on phase change materials are continuously attracting attention, see [6][7][8][9] for more details. These materials change from one state of matter to another one by releasing or absorbing energy and, consequently, they act as regulators: allow storing energy temporarily and freeing it when necessary.
A main application of phase change materials can be found at concentrated solar power plants [10], where they are used together with heat transfer fluids for storing energy. A way for improving the thermal efficiency of these plants consists of adding nanoencapsulated phase change materials (nePCMs) to the heat transfer fluid or to the thermal storage fluid. This mixture, commonly known as nanofluid [11], enables not only to improve the efficiency of heat transfer [12] but also to store energy to overcome the mismatch between supply and demand of energy [5]. Nevertheless and despite the fact that nePCMs have a direct impact in the thermal efficiency and heat storage, their synthesis becomes a difficult task.
From a mechanical point of view and due to the thermal stresses which appear in heating/ cooling cycles [13], one of the major problems to synthesize nePCMs arises in determining the thickness of the shell which confines the phase change material (core) given that a compromise between mechanical reliability and heat storage must be achieved. Both mechanical and thermal capabilities can be measured by the encapsulation ratio g, which is defined as the ratio between the volume of the nePCM core and that of the whole nePCM (core þ shell): g % 1 implies high thermal efficiency but low mechanical reliability, g ( 1 produces high reliability and low thermal efficiency.
Furthermore, increasing the size of the nanoparticle as a way of enhancing its heat storage capacity is discarded given that the colloidal stability of the nanofluid is not guaranteed as early as a threshold value of the nanoparticle radius is overcome [14,15].
Owing to the complexity that this problem entails, different scientific and technical communities are involved in its study. Therefore, together with experimentation, numerical simulations appear to be suitable to gain in understanding while trying to reduce the number and the cost of experiments to be conducted.
Despite the fact that there are several numerical models in the literature concerning different aspects of nanoparticles, their scope of study is rarely devoted to describe the thermomechanical behavior of the nePCMs. For instance, the thermal behavior of the shell is accurately described in Ref. [16], but it does not consider the influence of the thermal stresses on the shell.
Regarding phase change without mechanical interactions, a great variety of numerical schemes are available in the literature; for instance, [17,18] use the finite difference method and [19,20] the finite element (FE) method. According to Refs. [21][22][23][24], materials exhibit two different behaviors when changing their state of matter from solid to liquid or vice versa, see the schematic enthalpy variation for both phase change cases shown in Figure 1: Pure substances present a sharp change in their value of enthalpy H, see Figure 1 (left), which represents H versus temperature T and the two matter states: solid and liquid. Alloys present a smoother variation of H, see Figure 1 (right), since both phases co-exist at the same time when the temperature T 2 ½T s , T l , where T s and T l denote solidus and liquidus temperature [25], respectively. The transition zone is commonly referred as mushy zone.  Numerically, pure substances result more problematic than alloys given that the latent heat released/absorbed leads to a discontinuity in enthalpy. In the framework of the FE, a direct element integration in the presence of jump discontinuities produces errors, which can be solved by regularization techniques [20].
In addition and according to Refs. [21][22][23][24], there are basically two families of numerical schemes to numerically solve phase change: Tracking domain schemes, for which the phase change interface is continuously tracked. Fixed domain schemes, for which the phase change is calculated after the calculation of temperature distributions.
On the one hand, the first scheme is accurate for pure substances but not suitable for alloys. Besides, this method often requires mesh adaptivity or geometric transformations to determine the phase change interface. On the other hand, the second scheme is suitable for both pure and alloy substances and it is easier to implement than tracking methods [19].
Finally, a thermomechanical FE formulation with phase change is reported in Ref. [26]. However, this work uses a staggered approach: first a thermal analysis is performed to obtain the temperature distributions and then a mechanical analysis is conducted. Therefore, the computational time increases and the accuracy and robustness decreases.
In this context, the current work presents a three-dimensional and thermodynamically consistent formulation applied to thermo-elastic phase change pure substances. For this purpose, linear momentum and energy balances are stated and the constitutive equations are obtained from a thermodynamic potential, specifically, from the Helmholtz's potential. Then, the governing equations are discretized in the context of the FE method [27], which is more robust than the finite difference method. In particular, a monolithic (no staggered approach is required) and displacement-based formulation by using eight-noded elements with four degrees of freedom per node is considered. With regard to phase change, a fixed domain scheme is adopted and three implicit numerical schemesequivalent heat capacity, heat source and enthalpywith regularization techniques are implemented and tested by using one-dimensional analytical solutions extended by the authors of the present work.
Finally, the numerical tool developed in the present work is applied to study phase change in nePCMs in order to determine their temperature distribution and asses their mechanical strength. In particular, two nePCM geometries (spherical and cylindrical) and two pair of core@shell materials (Sn@SnO 2 and Al@Al 2 O 3 ) are simulated and the Rankine's criterion is used to predict the mechanical failure of the nePCM shell.
The current work assumes linear elasticity for the solid phase given that, from an experimental point of view, the plastic behavior of the shell should be avoided. For the liquid phase and since the core volume is reduced: i) advection terms are neglected in a first and good approximation as was also adopted in Ref. [26] for modeling welding processes, and ii) the liquid behaves like a liquid at rest, as assumed in Ref. [28]. Constant material properties are considered in each state of matter (solid and liquid). Experimentally, material properties exhibit temperature-dependency, but the lack of data and the dispersion in the measurements reported in literature make the constancy assumption a reasonable modeling choice.

Theoretical formulation
Mathematically, the thermomechanical phase change problem is expressed by a set of two coupled differential equations, called governing equations, which are composed of balance and constitutive equations and of boundary conditions. The mathematical notation used in the present work is presented in Table 1.

Balance equations
Consider a body of domain X, boundary C and its outward normal n containing solid and liquid phases. In order to model the current thermomechanical phase change problem, three balance equations must be considered: linear and angular momentum balances and energy balance.

Mechanical balances
Linear and angular momentum balances for both solid and liquid phases may be expressed as: where q, € u , r , f denote mass density, acceleration, Cauchy stress tensor and body force vector, respectively. Besides, the stress tensor is directly related to the traction vector t by the Cauchy relation: t ¼ r Á n: Finally, the angular momentum balance is automatically satisfied by the symmetry of the Cauchy stress tensor, as expressed in the right equation of (1).

Energy balance
For the sake of convenience, the energy balance is expressed in terms of enthalpy H, which is defined as [23]: where c, L denote heat capacity and latent heat, respectively; T ref , T m are reference temperature at which enthalpy is calculated and melting temperature, respectively; and hðT À T m Þ is the Heaviside step function, which reads: Finally, the energy balance may be expressed as: where q and Q denote heat flux and heat source/sink, respectively.

Constitutive equations
In this section, constitutive equations are obtained by consistent thermodynamic approaches based on equilibrium and non-equilibrium theories, see [29][30][31] for more details.

Thermomechanical constitution
The material constitution for the solid phase is calculated from the Helmholtz energy potential F , which is obtained by combining the first and second law of thermodynamics, by assuming that only reversible processes are considered, by applying a Legendre transformation to exchange the entropy S by T, and by assuming a natural state F ðT ¼ T ref , e ¼ 0Þ ¼ 0 for which the body is undeformed and at a reference temperature T ref : where e ¼ r s u denotes small strain tensor, u the displacement vector with Cartesian components u ¼ ðu, v, wÞ, r s the symmetric part of the displacement gradient, and hot is the abbreviation for high-order terms. The three first terms in the Taylor expansion of (5) vanish since the natural state is zero and there are neither initial stresses nor initial variation of entropy, respectively. Furthermore, Biot coupling [32] is not considered in the current work: a one way thermoelastic coupling is assumed. Finally, C and b denote fourth-order elastic and second-order thermoelastic tensors, respectively, which are explicitly expressed as: where I , I sy denote second-and symmetric part fourth-order identity tensors, respectively [33], and the Lam e parameters are expressed as: where E, and a denote Young's modulus, Poisson's ratio and thermal expansion coefficient, respectively. Finally, the constitutive equation for both solid r s and liquid r l phases is obtained by a standard equilibrium thermodynamics approach [34] to obtain: where it is assumed that the liquid phase change material inside the shell behaves like a liquid at rest (hydrostatics) and then the deviatoric part of stresses in the liquid is not present, as indicated in Ref. [28].

Heat conduction
From a phenomenological point of view, heat flux and its driving forcethe gradient of temperatureare related in a first and good approximation by [34]: where j ¼ jI denotes the isotropic thermal conductivity tensor.

Boundary conditions
The boundary conditions are composed of Dirichlet (also known as first-type) or Neumann (second-type) expressions: where u , T , t and q are the prescribed displacements, temperature, traction vector and thermal flux, respectively.

Outline of numerical phase change schemes
This section briefly describes the three different numerical phase change schemes used in the current work, namely: equivalent heat capacity hc, heat source hs and enthalpy e schemes.

Equivalent heat capacity scheme
In this scheme, the rate of enthalpy is calculated by directly applying the chain rule to (2): where d denotes the Dirac delta function. Introducing (11) in (4), the energy balance becomes: From a numerical point of view and according to Refs. [23,24], a numerical regularization is performed and c(T) reads: where c s and c l denote heat capacity for solid and liquid phases, respectively, is the regularization parameter, which ensures the correct integration of the d function, and T s ¼ T m À and T l ¼ T m þ represent temperatures for solid and liquid phases, respectively.

Heat source scheme
This scheme directly performs the derivative of (2) with respect to time: Now, by applying a backward first-order finite difference with time step Dt to the second term on the right-hand side of (14) and introducing it into (4), the energy balance becomes: where h nþ1 and h n denote the regularized Heaviside step function at current time n þ 1 and at previous time n, respectively. This regularization form at the current time (obviously analogous for h n ) may be expressed as [35]:

Enthalpy scheme
In this scheme, the rate of enthalpy is directly discretized by using a backward first-order finite difference scheme. Then, the energy balance of (4) becomes: where H nþ1 and H n denote the regularized enthalpy at current and previous time, respectively. The regularized enthalpy at the current time (similar for previous time) may be expressed as [23,36]:

Weak forms
Since the strong forms are second-order differential functions of the degrees of freedom u and T, these forms are multiplied in the whole domain by arbitrary test (also called weight) functions du and dT in order to obtain an amenable displacement-based FE formulation. Then, the divergence theorem is applied to the gradient term of both strong forms and the Neumann boundary conditions of (10) are enforced to calculate the weak forms, which are first-order differential equations of the degrees of freedom.
Finally, the mechanical weak form becomes: ð The three thermal weak formsone for each phase change schemeread: ð

Discretization
In order to obtain numerical solutions in the framework of the FE method, the continuum domain X is discretized by n three-dimensional eight-noded brick elements of domain X e and boundary C e . For this purpose, an isoparametric interpolation by using standard shape functions of Lagrange-type N is adopted to interpolate the global coordinates, the test functions and the four degrees of freedom: where the Einstein summation convention is used; a j a denotes the nodal values at a generic node a for each degree of freedom j ¼ ðu, v, w, TÞ; and B s and B denote the discretized form of the symmetric gradient of displacements and gradient of temperature, respectively.

Residuals
Despite the linearity of the problem, a residual-based formulation is adopted in the present work for the sake of completeness. For it, by introducing (21) in (19), the mechanical residual reads: where the constitutive equation of r depends on the phase, namely, solid r s or liquid r l , mathematically: Likewise, by introducing (21) in (20), the thermal residuals for each phase change scheme become: where, as commented, the indexes hc, hs and e refer to the phase change schemes. The discretized form of the heat flux of (9) becomes q ¼ Àj B j a T j :

Assembled tangent matrix
This section presents the final assembled and monolithic matrices at generic nodes a, b for the schemes k ¼ fhc, hsg : and for the case e: where K, C and M denote stiffness, capacity and mass matrices, respectively, and they are explicitly calculated in the Appendix A. In addition, the coefficients c 2 and c 3 are scalar quantities, which result from linearizing the Newmark relations, see [37]. Finally, the numerical formulation is implemented into the research code FEAP [38], which belongs to the University of California at Berkeley (USA). This software holds dummy routines, called user elements, that permit to introduce new modular elements as that of the present work.

Validations
This section presents several comparisons between analytical and numerical solutions in order to check the correct implementation of the numerical formulation. For this purpose, available closed solutions in the literature, which solve phase change problems, are extended by the authors of the current work by including the mechanical field, see Appendix B. Figure 2 shows the geometry and boundary conditions of the numerical model used for the validations. A fixed-free rod of length L y (L y ) L x , L z ) at an initial temperature T i is considered and a time-dependent temperature T 0 is prescribed at the free-end. Since T 0 > T m > T i , the phase change interface will move progressively toward the fixed-end.
For the validations, material properties are those of tin (Sn), which are obtained from [39][40][41][42][43] and summarized in Table 2, and the values T 0 ¼ 573:15 (K) and T i ¼ 303:15 (K) are considered. In addition, Table 2 reports the material properties of aluminum (Al), obtained from [39] and used in Section 6 for further analyses.
Experimentally, materials exhibit temperature-dependency. However, the lack of available data characterizing the temperature-dependency over the desired temperature range, the considerable dispersion of the temperature-dependent values reported in literature and the complexity of measuring some of the temperature-dependent properties are the main reasons to consider constant properties in each phase. Nevertheless, the inclusion of temperature-dependent material properties in the numerical formulation would be straightforward in residual-based FE formulations, like the one developed in the present work. Figure 3 compares analytical (solid, dashed and dotted lines) and numerical solutions (solid circles) for temperature distributions (left column) and axial displacements (right column) along the one-dimensional geometry, for the three different phase change schemesheat capacity (top row), heat source (middle row) and enthalpy (bottom row)and for three different times: 1, 5, 10 (s). For this comparison, the regularization parameter is ¼ 1:25: As observed in Figure 3, analytical and numerical solutions are in good agreement with each other for both temperature and axial displacement and for the three phase changes. In particular, the maximum relative error between analytical and numerical results for temperature and axial displacement for each numerical scheme is reported in Table 3.
In conclusion, any of the phase change schemes can be used to solve thermomechanical phase change problems in pure substances.

Analyses of thermal stresses in nePCMs
In this section, the previously formulated and validated numerical tool is applied to simulate four different scenarios in order to determine the temperature distribution on the nePCM shell and to assess the mechanical reliability and energy density of the nePCMs. For this purpose, two geometries and two pairs (core and shell) of materials are considered. Concerning geometry, spherical and cylindrical nePCM configurations are contemplated, as shown in Figure 4. In both geometries, the diameter of the core is d À 2 th, with th the shell thickness. The height of the cylinder is chosen in such a way to ensure that the total volume Analytical t = 5s Analytical t = 10s Numerical  (core þ shell) of both geometries of nePCMs is the same in order to be able to perform comparative analyses between them. Regarding material properties, two pair of core@shell materials are considered: tin@tin-oxide (Sn@SnO 2 ) and aluminum@alumina (Al@Al 2 O 3 ). Core material properties are reported in Table 2, while shell properties are given in Table 4. Tin oxide properties are obtained from [44][45][46][47][48] and alumina ones from [39,40,49]. Notice that r t denotes the tensile strength.
From a FE point of view, structured meshes of 3584 (sphere) and 3840 (cylinder) eight-noded elements are used. With regard to boundary conditions, the nanoparticle is mechanically fixed at its center and a linearly increasing temperature is prescribed at the outer surface of the shell. The initial temperature of the nanoparticle at t ¼ 0 (s) is T i ¼ 303.15 (K) and the prescribed temperature is linearly increased with time steps Dt¼ 20 (ns) until the final time t ¼ 0.5 (ms). At this time, T 0 ¼ 573:15 and T 0 ¼ 1050:15 (K) are reached for Sn@SnO 2 and Al@Al 2 O 3 nePCMs, respectively. The phase change enthalpy scheme with a regularization parameter of ¼ 1:25 is applied over the present section.

Temperature and Rankine's equivalent stress distributions
The first simulation is aimed to obtain temperature and maximum equivalent stress distributions on the nePCM shell. Shells are normally composed of oxides, which possess a mechanical behavior similar to that of ceramics. Despite the fact that the most adequate failure criterion for   ceramic materials is not clear in the literature [50], the Rankine's equivalent stress is adopted in the present work given that a shear-insensitive criterion is more adequate than a shear-sensitive one to describe the fracture behavior of ceramics [51]. Figure 5 shows contour plots of temperature and Rankine's equivalent stress distributions on half of the nePCM shell for each of the four scenarios of study. According to the experimental work reported in Ref. [13], the diameter d and shell thickness th used in the current section for all the scenarios are: d ¼ 103 and th ¼ 9.78 (nm).
In Figure 5, firstly, it can be observed that all the nanoparticle shells are at uniform temperature and very low gradients of temperature are appreciated since steady state is reached immediately in nanosolids due to their reduced physical size. Concretely, when the prescribed temperature increases its value, the transient temperature distribution disappears quickly and a new equilibrium state (with low gradients of temperature and consequently negligible heat fluxes) is reached for the new boundary condition.
In the second place, thermal stresses appear due to the difference in the thermal expansion coefficients of the core and shell materials, as was experimentally confirmed in [52]. Notice that this result has also been verified numerically.
Thirdly and with regard to the mechanical reliability of the nePCMs, the maximum numerical values of equivalent stresses are compared with their respective tensile strengths r t given in Table 4. From these comparisons, it can be concluded that: Spherical and cylindrical Sn@SnO 2 nePCMs are mechanically reliable during the heating process. In particular, an extra validation has been performed to reproduce the conditions reported in [13] and it is verified that the present numerical tool agrees with the experimental study in that article on the mechanical strength of spherical Sn@SnO 2 nePCMs. Spherical and cylindrical Al@Al 2 O 3 nePCMs are expected to fail.
Finally and for the sake of completeness, Figure 6 shows the time evolution of Rankine's equivalent stress at a point at the outer surface of the shell for each scenario of study. Several conclusions are obtained from these curves: The equivalent stress increases linearly with temperature until the melting temperature is reached and after that, stress decreases. Consequently, the maximum stress developed in the shell occurs just before melting starts. The trend in the time evolution curves is the same regardless of any material property or geometry but the amplitude of the equivalent stresses depends on both these parameters. Stresses in Al@Al 2 O 3 nePCMs are higher than those in Sn@SnO 2 nePCMs due to their difference in core T m , see Table 2. For the shell thickness th ¼ 9.78 (nm), stresses in cylindrical geometries are higher than those predicted in spherical ones. However, this is not always the case for different values of shell thickness, as shown in Section 6.2.

Analysis of the shell thickness
The present analysis shows the influence of shell thickness on energy density (E d ) and maximum Rankine's stress developed at the nePCM shell. The energy density measure used in the present analysis is defined as follows: where V core and V total denote core and total (core þ shell) volume of the nePCM, respectively. Figure 7 shows both the variation of the energy density of a single nePCM and of the maximum Rankine's equivalent stress for three different values of the shell thickness th, namely: {2, 5, 9.78} (nm). Spherical and cylindrical geometries with the same total volume are considered for comparison purposes.
In the first place, in Figure 7, it is observed that both energy density and Rankine's stress decrease with the progressive increase of shell thickness th. That decrease in energy density is caused by the reduction in the volume of the available phase change material (core) and, consequently, the energy efficiency of the nePCM is reduced. However, the increase in the shell thickness improves the mechanical reliability of the nePCM, which, as a result, diminishes the thermal stresses developed at the shell.
Secondly, and from a geometrical standpoint, it is observed that: The energy density of the spherical nePCMs is higher than that of the cylindrical ones because, for equal total volume of both geometries, the volume of core material inside the nePCM is larger in the spherical geometry. The maximum Rankine's stress is slightly higher for spherical nePCMs until a threshold value with increasing shell thickness is overcome and, from that point forward, cylindrical nePCMs are the ones undergoing higher thermal stresses for the same shell thickness. Thirdly, and regarding material properties, it is observed that they exert a direct influence on both energy density and Rankine's stress: Al@Al 2 O 3 nePCMs possess an energy density which is nearly twice the value of that of Sn@SnO 2 nePCMs. The reason of this disparity lies in the difference between the values of latent heat L and mass density q of the core materials, see properties in Table 2. With regard to mechanical reliability, comparing the maximum values of stress in Figure 7 with the r t given in Table 4, it may be concluded that: whilst Sn@SnO 2 nePCMs do not fail under thermal stresses, Al@Al 2 O 3 nePCMs are expected to do it.
Finally, spherical Al@Al 2 O 3 nePCMs posses the best energy performance. However, in terms of mechanical strength, Sn@SnO 2 nePCMs are the only resisting the thermal stresses developed under the previously reported conditions. Since the maximum value of stress is geometry-dependent for a given shell thickness (see Figure 7), a compromise between mechanical strength and energy density has to be achieved for each desired application.
In conclusion, mechanical capability of nePCMs highly depends on: i) the difference between the thermal expansion coefficient of the core-shell, ii) the shell thickness and its tensile strength and iii) the melting temperature necessary to reach the liquid state. In turn, energy capability of nePCMs highly depends on: i) the latent heat and mass density of the core and ii) the core volume of the nePCM.

Conclusions
A three-dimensional finite element formulation has been developed to numerically study thermomechanical phase change problems for pure substances. For this purpose, governing equations for mechanical and thermal fields are stated and discretized within the FE context and three different phase change schemes are considered and compared. The numerical formulation is implemented in a research code, which is validated by comparing numerical results against closed solutions extended by the authors of the present work. From these validations, it is concluded that the three phase change schemes are suitable to deal with phase change phenomena on pure substances. This numerical tool is used to simulate nePCMs in four scenarios of study: two different geometries (spherical and cylindrical) and two core@shell pairs of materials (Sn@SnO 2 and Al@Al 2 O 3 ) are considered. For each scenario, three analyses are performed: i) temperature and maximum Rankine's stress distributions on the nePCM shell, ii) time evolution of Rankine's stress and iii) study of the influence of the shell thickness on stress and energy density. From these analyses, it is concluded that the choice of the nePCM geometry and material pair must respond to a compromise between energy density and mechanical strength, which must be thoroughly examined for each desired application.
Despite the uncertainty associated to the values of some material properties, numerical simulations provide a good estimation of the stresses developed in nePCMs during thermal processes. Hence, this framework appears to be a powerful tool, complementary to experiments, to determine the thickness needed for the nanoparticle shells.

Appendix A. Tangent matrices
According to Ref. [27], tangent matrices are calculated from the residuals of (22) and (24) by solving: where the indexes i, j refer to the degrees of freedom and a, b to two generic nodes.
Applying (A1) to (22), the mechanical matrices for the solid phase become: and, for the liquid phase: Now, applying (A1) to (24), the thermal matrices for the heat capacity hc scheme read: for the heat source hs scheme: dX e , C TT ab, nþ1 ¼ À and, finally, for the enthalpy e scheme:

Appendix B. Analytical solution
This appendix presents an analytical solution for a thermomechanical phase change problem applied to a onedimensional half-space domain. The analytical solution for the thermal field considering phase change is reported in Refs. [53,54]. The authors of the current work have extended that solution by including the mechanical field. For this purpose, it is assumed that the body is not subjected to any traction and, consequently, the axial displacement v of the solid phase may be calculated as: where L y denotes length of the body, as shown in Figure 2. Finally, the expression of the axial displacement reads: a ðT m À T i Þ ffiffiffi p p erfcðnÞ 2 e À y 2 4b s t À e À where t, b s ¼ j s =ðq s c s Þ, erfc and n denote time, thermal diffusivity of the solid phase, the complementary error function and a dimensionless coefficient reported in Refs. [53,54], respectively.
Under restrictive assumptions, an analytical solution for a one-dimensional fluid can be obtained. More precisely, an analytical solution for the case of a non-viscous fluid at rest is provided according to the constitutive law in Eq. (8) (right). In this particular case, pressure in liquid phase can be computed as: where K denotes bulk modulus. By considering (B2), an explicit expression of pressure can be found:

Disclosure statement
No potential conflict of interest was reported by the authors.