<figure><img src="image_1.png" /></figure>

HAL
open science

# 3D fluid-structure simulation of innovative composites for the design and thermal management of electronic devices

Fabien Salmon, Delphine Lacanette, Marie Duquesne, Alexandre Godin

<figure><img src="image_2.png" /></figure>

## To cite this version:

Fabien Salmon, Delphine Lacanette, Marie Duquesne, Alexandre Godin. 3D fluid-structure simulation of innovative composites for the design and thermal management of electronic devices. Energy Conversion and Management, 2023, 280, pp.116824. 10.1016/j.enconman.2023.116824. hal-04828228

HAL Id: hal-04828228

https://hal.science/hal-04828228v1

Submitted on 9 Dec 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

## Abstract

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

# 3D fluid-structure simulation of innovative composites for the design and thermal management of electronic devices

F. Salmon¹*, D. Lacanette¹, M. Duquesne², A. Godin²,³

¹Université de Bordeaux, CNRS, Bordeaux INP, I2M, Bât A11, 351 cours de la Libération, 33400 Talence, France

²Université de La Rochelle, LaSIE UMR CNRS 7356, Avenue Michel Crépeau, CEDEX 1, 17042 La Rochelle, France

³ev Lab, EDF R&D, CNRS, LaSIE, La Rochelle University, Avenue Michel Crépeau, CEDEX 1, 17042 La Rochelle, France

*Corresponding author: Fabien.Salmon@u-bordeaux.fr

Thermal management of electronic devices remains a challenging issue, especially the cooling of electronic components. Current methods mostly involve noisy, cumbersome and underperforming fans. It has recently been shown that phase change materials (PCMs) can perform better and offer more efficient damping of temperature variations in electronic components. Here, we present three novel contributions. First, we consider a triply-periodic minimal surface, called a gyroid, as a porous architecture containing the PCM. Such a topology mathematically maximizes the exchange surface between the copper structure and the PCM. Heat transfer and thus the performance of such a composite are therefore optimized. Second, we set up a complete 3D coupled fluid-structure model of the composite. This model simulates the heat transfer (conduction and convection) inside both the PCM and the porous structure, as well as the mechanical stresses induced on the container by thermal expansion of the PCM. The 3D coupled theoretical model is based on open-source software codes, OpenFOAM, CalculiX and preCICE. With this novel model, deformation of the structure and mechanical stresses are now predictable. It is therefore possible to anticipate whether a configuration will be prone to leakage due to mechanical failure. Thermal contact resistance is also computed, giving information on heat transfer quality at the interface between the composite and the electronic component. Third, the application of this model to the innovative composite shows that reducing the initial filling level of PCM is a better solution than strengthening the container to maximize thermal management.

**Keywords:** Electronic component; Gyroid; Phase change materials and conductive matrix composite; Fluid-structure interaction; Conjugate heat transfer; CFD
## List of symbols

$c_p$ Specific heat at constant pressure (J.kg⁻¹.K⁻¹)

$d$ Structure displacement (m)

$D$ Thermal diffusivity (m².s⁻¹)

$E$ Young's modulus (Pa)

$g$ Gravitational acceleration (m.s⁻²)

$L$ Characteristic length (m)

$m$ Mass (kg)

$p$ Pressure (Pa)

$R_{th}$ Thermal resistance (K.m².W⁻¹)

$s$ Phase change interface position (m)

$T$ Temperature (K)

$u$ Velocity (m.s⁻¹)

$V$ Volume (m³)

## Greek symbols

$\alpha$ Solid fraction

$\beta$ Coefficient of thermal expansion (K⁻¹)

$\delta$ Peak to valley (m)

$\Delta h_f$ Heat of fusion (J.kg⁻¹)

$\varepsilon$ Temperature range where melting occurs (K)

$\varepsilon$ Infinitesimal strain tensor

$\zeta$ Square root of the liquid thermal diffusivity over the solid one

$\eta$ Poisson's ratio

$\kappa$ Penalization parameter (10²⁰ kg.m⁻³.s⁻¹)

$\lambda$ Thermal conductivity (W.m⁻¹.K⁻¹)

$\mu$ Dynamic viscosity (Pa.s)

$\nu$ Cinematic viscosity (m².s⁻¹)

$\xi$ Filling level of PCM

$\rho$ Density (kg.m⁻³)

$\sigma$ Stress tensor (Pa)

$\tau$ Characteristic time (s)

$\chi_T$ Adiabatic compressibility (Pa⁻¹)

## Subscripts

0 Reference configuration

air Air properties

conv Convection

diff Diffusion

l Liquid properties

m Melting

s Solid properties

## Acronyms

CFD Computational Fluid Dynamics
CHT Conjugate Heat Transfer

FSI Fluid-Structure Interaction

OF OpenFOAM

PCM Phase Change Material

51

52

# 1. Introduction

53

54

55

56

57

58

59

60

61

Phase change materials (PCMs) are taking an increasingly prominent place in many fields relating to energy management. With current and future energy issues, applications of these materials are emerging or improving. Thermal energy storage is of great interest for renewable energies [1]. Solangi et al. [2] have proposed to use MXene PCMs to store solar energy. Palomba et al. [3] assessed the performance of a prototype based on PCMs for renewable energy purposes. More broadly, many studies have been conducted on thermal storage using PCMs over the years [4], [5]. PCMs are also widely used in buildings to improve their insulation [6], [7]. The impact of such materials on electricity consumption has been studied by Qureshi et al. [8], among others.

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

Thermal management with PCMs is also of great interest in electronics [9], [10], with cooling of electronic components being a key issue in this area. Currently, most of these components are used in combination with fin heat exchangers that dissipate heat by natural or forced convection and radiation. The coolers are mainly composed of noisy moving parts such as fans, which have many disadvantages in terms of performance, ergonomics and safety. Sahoo et al [11] have shown that use of PCMs is particularly suitable for cooling electronic components with intermittent operation. These observations have been confirmed by numerous studies that have shown these technologies to be useful [12], [13], [14]. PCMs considerably delay the onset of overheating and limit the energy requirements inherent in moving parts [11], [15], [16]. In this regard, Gharbi et al. [17] compared the thermal performances of several PCMs, while a cascade configuration of PCMs to cool devices was studied experimentally and numerically by Huang et al. [18]. The cooling performance of carbon nanotubes associated with paraffin was also investigated experimentally [19]. Ho et al. tried to optimize geometries to improve thermal management [20].

77

78

79

80

81

82

83

84

85

86

87

88

89

90

Most studies in this field are experimental because theoretical predictions are made difficult by the complexity of the physical phenomena involved. However, experimentation can become expensive and some configurations may be out of reach. Nowadays, thanks to improved computational performance, Computational Fluid Dynamics (CFD) is increasingly helpful as a complement to experimentation. In particular, simulation can predict the physical behavior of thermal management based on PCMs. Many authors manage heat transfer in PCMs considering only conduction, without convection [21], [22]. Other works take convection into account, which modifies the temperature field and thus thermal performance [23], [24]. These numerical studies were carried out with commercial codes. To reduce the cost of simulations, open-source software should be used. There are a few investigations with free codes such as OpenFOAM [25], since the set of governing equations cannot be solved directly but must be implemented personally. 2D numerical simulations of PCMs have been conducted with OpenFOAM [26], [27].

91

92

93

94

95

Although the technology is promising, containers can suffer from leakage due to the high mechanical stresses induced by thermal expansion [28], [29]. Such issues are not theoretically predictable at present, however, and there are no studies focusing on this limiting problem. To address it, a fully coupled model should set up between a PCM, which is either solid or liquid,
and solid structures containing the PCM. This model should consider conjugate heat transfers (CHT) and fluid-structure interactions (FSI). There are two ways of solving such a complex problem. First, the monolithic approach [30] considers one set of equations involving every physical phenomenon (fluid and solid). This method is numerically stable but very time-consuming, and is therefore very scarcely used by authors. Second, the partitioned approach [31] considers several sets of equations, which are solved one at a time, often by different software. This method requires an interface coupling between the codes that are used. For configurations involving a fluid and a solid, there is often one code handling the fluid, a second one managing the solid, and a third one coupling them. Numerical stability issues can occur, but this approach is often preferred in the literature since it is less time-consuming.

Critical temperatures and power losses depend on the electronic components. Based on the requirements specified for the component, a PCM and a conductive porous matrix must be selected. The investigation of a composite and the modeling of its thermal behavior could make it possible to optimize thermal performance by selecting the porous matrix topology that maximizes heat transfers. The objective of this paper is therefore to propose and model the thermal and mechanical behavior of new composite materials (architected conductive porous matrix containing phase change materials) that meet precise specifications, such as the highest possible energy density, adapted phase change temperatures, limited volume expansion, properties allowing the intensification of transfers, and the lowest possible thermal resistance between the electronic constituents of the composite materials. Some technical and technological obstacles could be overcome thanks to a better fundamental understanding of the phase transition and transfer phenomena involved in these new heterogeneous materials. Development of innovative solutions for the thermal management of electronic components therefore necessarily requires particular attention to the material scale.

Here, we present coupled simulations of a PCM inside gyroid structures [32]. Gyroids are triply-periodic minimal surfaces, which maximizes the heat exchange surface between PCMs and the surrounding structure. The simulations were managed by a fully open-source software suite composed of OpenFOAM [25], CalculiX [33] and preCICE [34]. For the first time, coupled fluid-solid simulations of PCMs are performed. 3D conjugate heat transfer is considered first between the solid container and melting PCM. Convection of the liquid PCM is taken into account. Then, 3D fluid-structure simulations of the global system are performed. In particular, the fluid pressure caused by thermal expansion and then the induced stresses in the structure are calculated. This global simulation allows us to predict potential failures and leakage, which is a major step forward. In addition, with the simulated solid deformation, contact thermal resistance can be calculated. This theoretical advance allows the assessment of heat transfer losses due to imperfect contact between PCM systems and electronic components, for instance. While most studies focus on the development of new materials increasing heat transfers and exchanges [35], [36], [37], [38], [39], contact thermal resistance plays a significant role in the thermal management of components and is the first obstacle to be overcome for good thermal management. The fluid-structure coupling will thus allow us to design the porous structure correctly to maintain the same contact between the component and the innovative composite throughout the different thermal cycles undergone by the component. The present study therefore challenges the following points:

* creating a 3D open-source coupled FSI and CHT modeling of PCMs and porous architecture

* computing stresses in porous architecture induced by thermal expansion of PCMs to settle mechanical failures

* computing structure displacements to assess the contact thermal resistance
* applying the modeling to gyroid structures which mathematically maximize the exchange surface

* assessing the influence of the initial filling level of PCM and the thickness of the gyroid branches on the structure displacements and mechanical stresses

In the paper, we begin by describing the theoretical modeling of both a PCM and a conductive gyroid solid structure in section 2. The CHT and FSI couplings are described more particularly. The PCM model was implemented in OpenFOAM and must thus be validated. Section 3 focuses on the comparison between our PCM model and the analytical solution of the Stefan problem, as well as experimental data from a specific configuration. Section 4 presents the numerical configuration of the simulations. In particular, a parameter study is performed on the geometry and filling level of the PCM. Section 5 concerns the results of these simulations. Temperatures, velocities, mechanical stresses and structure displacements are studied specifically. The last section discusses the most suitable geometry given the simulated mechanical stresses and displacements.

## 2. Numerical methods

The paper focuses on 3D numerical simulation of the fusion of a Newtonian PCM within a porous architectured structure. First, we present the equations governing the PCM, then those describing the thermo-mechanical behavior of the structure, and finally the interface coupling between both solvers.

### 2.1 PCM modeling

#### 2.1.1 Equations

After melting, the liquid PCM is Newtonian and the incompressible Navier-Stokes equations (1) thus govern its behavior reliably. The variation in density with temperature is considered by the Boussinesq approximation. The influence of this assumption will be addressed in section 3.2. To avoid any velocity in the remaining solid material, a penalization term is added into the momentum equation. The energy equation is based on the equivalent enthalpy approach, where the melting process is considered through an additional source term.

$$ 
\left\{ 
\begin{aligned}
\frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} &= -\frac{\nabla p}{\rho_0} + (1 - \beta(T - T_0))\mathbf{g} + \nu \Delta \mathbf{u} - \frac{\kappa \alpha}{\rho_0} \mathbf{u} \\
(\rho c_p)_{eq} \left(\frac{\partial T}{\partial t} + \mathbf{u} \cdot \nabla T\right) &= \nabla \cdot (\lambda_{eq} \nabla T) - \rho_s \Delta h_f \frac{d\alpha}{dt}
\end{aligned}
\right. 
\quad (1) $$

In equations (1), **u** is the velocity (m.s⁻¹), **p** is the pressure (Pa), **T** is the temperature (K), **T**₀ is a reference temperature (K), **ρ**₀ is the density at the reference temperature (kg.m⁻³), **β** is the coefficient of thermal expansion (K⁻¹), **g** is the gravitational acceleration (m.s⁻²), **ν** is the kinematic viscosity (m².s⁻¹), **α** is the solid fraction, **κ** is a huge penalization parameter (10²⁰ kg.m⁻³.s⁻¹) aiming to avoid velocity in the solid domain, **ρ** is the density which depends on temperature (kg.m⁻³), **c**ₚ is the specific heat capacity (J.kg⁻¹.K⁻¹), **λ** is the thermal conductivity (W.m⁻¹.K⁻¹), **ρ**ₛ is the solid density at melting temperature (kg.m⁻³) and Δ**h**f is the heat of fusion (J.kg⁻¹). The index eq in the energy equation refers to equivalent properties giving

$$ (\rho c_p)_{eq} = (1 - \alpha)\rho_l c_{p,l} + \alpha \rho_s c_{p,s} \quad (2) $$
$$ \lambda_{eq} = (1 - \alpha)\lambda_{l} + \alpha\lambda_{s} \qquad (3) $$

where the $l$ and $s$ indices refer to liquid and solid respectively.

Due to the temporal derivative of the solid fraction in the energy equation, a $C^1$ function is required for $\alpha$ and must be equal to 0 for temperatures greater than the melting temperature, and equal to 1 for lower temperatures, with as sharp a transition as possible around the melting temperature. Equation (4) presents the chosen function, already used in other models [40], [41]:

$$ \alpha = \frac{1}{2} \left( \tanh \left( \frac{T_m - T}{\varepsilon} \right) + 1 \right) \qquad (4) $$

where $T_m$ is the melting temperature of the solid PCM (K) and $\varepsilon$ (K) a small parameter controlling the sharpness of the transition. Here we chose $\varepsilon = 10^{-2}K$, which is a common choice to avoid numerical issues while maintaining a sharp transition [40], [41]. The derivative of the solid fraction in the energy equation becomes:

$$ \frac{d\alpha}{dt} = -\frac{1}{2\varepsilon} \left(1 - \tanh^2\left(\frac{T_m - T}{\varepsilon}\right)\right) \frac{dT}{dt} \qquad (5) $$

To be efficient, the penalization term to avoid velocity inside solid parts must be treated implicitly in the linear resolution, in order to get a huge value on the diagonal of the matrix and thus avoid velocity. The sign of this term also has an impact on numerical stability. In order for the matrix to be strictly diagonally dominant, it is better to add a positive penalization term to the matrix diagonal.

Turbulence is not considered in the resolution. The greatest velocities encountered are about 1 cm.s⁻¹ (section 5) while the structure is also of the order of 1 cm. This leads to a Reynolds number approximately equal to $Re = \frac{uL}{\nu} \sim \frac{10^{-2} \times 10^{-2}}{10^{-6}} \sim 100$, which allows us to assume laminar flows.

### 2.1.2 Numerical schemes

The PCM equations are solved with OpenFOAM where they have been implemented. The discretization of the equations is based on the finite volume method. Resolution of the velocity-pressure coupling relies on the classical PISO (Pressure Implicit with Splitting of Operator) algorithm. The temporal term is discretized with a blend between the Euler and Crank-Nicolson schemes. In the simulations, the weight associated with the Crank-Nicolson scheme is 0.9, which ensures an almost second-order time discretization with enough numerical diffusion to avoid divergence. The advection term in the momentum equation (1) is discretized with the second-order upwind scheme. In the energy equation, the advection term is discretized with a TVD (Total Variation Diminishing) scheme switching between first (upwind) to second order (centered) depending on the velocity gradient. The central difference discretization is applied to the diffusion term. The mesh discretization and the boundary conditions are discussed in section 4.

## 2.2 Porous architecture structure modeling

The structure is governed by the linear elastic thermo-mechanical equations (6). The first equation computes heat conduction while the second describes mechanics.
$$ egin{cases} \frac{\partial T}{\partial t} = D_s \Delta T \\ \rho_s \frac{\partial^2 d}{\partial t^2} = \nabla \cdot \boldsymbol{\sigma} \end{cases} \qquad (6) $$

$D_s = \frac{\lambda_s}{\rho_s c_{p,s}}$ is the thermal diffusivity of the solid, which is assumed to remain constant, $d$ is the structure displacement and $\boldsymbol{\sigma}$ is the stress tensor based on Hooke's law and an additional thermal term $\boldsymbol{\sigma} = \frac{E}{1+\eta}\boldsymbol{\varepsilon} + \frac{E\eta}{(1+\eta)(1-2\eta)}\mathrm{tr}(\boldsymbol{\varepsilon})\boldsymbol{I} - \frac{E\beta}{1-2\eta}(T-T_0)\boldsymbol{I}$. $E$ is Young's modulus of the structure, $\eta$ is Poisson's ratio, $\boldsymbol{I}$ is the second rank unity tensor and $\boldsymbol{\varepsilon} = \frac{1}{2}(\nabla^T d + \nabla d)$.

These equations were solved with CalculiX, which is based on the finite element method. The temporal terms are discretized with the second-order backward Euler scheme. Section 4 presents a mesh convergence on the structure.

## 2.3 Coupling

### 2.3.1 Physical considerations

In order to optimize the calculation time, we first conduct a dimensional analysis to predict which phenomena can be neglected in the simulations. Section 4 presents the properties of the materials involved in the simulations in section 5. We use these properties in the following preliminary calculations. The results presented in section 5 show that the velocity is about 1 cm.s⁻¹. The convection time is $\tau_{conv} = \frac{L}{u} \sim \frac{10^{-2}}{10^{-2}} \sim 1$ s, the thermal diffusion characteristic time is $\tau_{diff,s} = \frac{L^2}{D} \sim \frac{10^{-4}}{10^{-4}} \sim 1$ s in the porous architectured structure and $\tau_{diff,s} = \frac{L^2 \rho c_p}{\lambda} \sim \frac{10^{-4} \times 10^3 \times 10^3}{10^{-1}} \sim 10^3$ s in the PCM, and the mechanical characteristic time is $\tau_{mech} = \frac{L}{\sqrt{\frac{E}{\rho}}} \sim \frac{10^{-2}}{\sqrt{10^{11} \times 10^{-3}}} \sim 10^{-6}$ s. The deformation velocity of the porous architectured structure is thus far higher than the other phenomena velocities in the problem. Then, we can consider quasi-static calculations for the mechanical equations.

Now, let us estimate the impact of mechanical deformations on liquid convection. The calculated liquid pressures in the simulations discussed in section 5 are of the order of 100 Pa.

The ratio of the dynamic pressure over the pressure is $\Pi = \frac{1}{2} \frac{\rho u^2}{p} \sim \frac{10^3 \times 10^{-4}}{10^2} \sim 10^{-3}$. The pressure variations are thus weakly related with the velocity field, and rather linked with density variations due to the increase in temperature. The simulated structure displacements in section 5 are smaller than $10^{-4}$ m. This weak value can globally impact neither heat transfer nor convection in the melted PCM.

Therefore, the impact of solid deformations is negligible on the PCM melting process. The calculations can thus be divided into two simulations and the mechanical and thermal equations can be uncoupled. First, a transient Conjugate Heat Transfer (CHT) simulation is carried out, without solving the mechanical equations. Second, steady Fluid-Structure Interaction (FSI) simulations are performed.

### 2.3.2 Conjugate heat transfer coupling
In the CHT coupling, mesh deformation is neglected. This allows us to perform the simulation in an explicit way, whereas considering mesh motion would have required implicit coupling with several sub-iterations. The boundary conditions at the interface for both software codes (OpenFOAM and CalculiX) are based on a Robin-Robin coupling. This means that the temperatures and heat fluxes at the interface are exchanged between the PCM and the structure. Then, a Robin boundary condition is applied to the interface in both codes. This part of the process was handled with preCICE.

### 2.3.3 Fluid-structure interaction coupling

The volume of the solid PCM is smaller than that of the liquid PCM after melting. In addition, as temperature increases, the liquid PCM expands. If the container volume is smaller than the liquid volume, then pressure can surge. The paper aims to compute the mechanical impact of the thermal expansion of the PCM on the surrounding porous architecture structure. However, the impact of the expansion of the volume of the liquid on pressure is not calculated in the model described in the prior sections since an incompressible fluid has been considered. Density is therefore assumed to be independent of pressure, but does depend on the temperature thanks to the Boussinesq approximation. This assumption is sufficient for the previous CHT simulation, but not for the FSI one. The compressibility of the melted PCM must be considered to estimate the pressure variation induced by the fluid thermal expansion. Solving the compressible Navier-Stokes equations requires an equation of state and there is no such a relation in the literature.

We will thus not be able to compute the liquid pressure during the melting process. However, we are interested in the maximum mechanical stresses in the porous architecture structure. This happens at the end of the simulation, when all the PCM has melted and involves the maximum liquid expansion. When the porous architecture structure is only filled with a liquid at a similar temperature thanks to convection, the pressure can be assumed homogeneous. As noticed previously, the pressure gradients induced by the liquid velocity are indeed very weak, and the temperature heterogeneities in the liquid PCM are very low at the end of the melting process.

The liquid pressure can then be calculated with the following analytical approach. Let us note $\chi_T = -\frac{1}{V} \frac{\partial V}{\partial p} \Big|_T$ the adiabatic compressibility ($Pa^{-1}$). Under the assumption that this property remains constant during compression,

$$p - p_0 = \frac{1}{\chi_T} \ln \frac{V_0}{V} \qquad (7)$$

where $V_0$ is the volume of the liquid phase under pressure $p_0$. With $p_0 = 0 Pa$, $V_0$ would be the volume of the melted PCM without constraint. We assume that the encapsulating structure is initially partially filled with the PCM at the melting temperature such that $m = \xi_0 \rho_l(T_m) V_p$ where $V_p$ is the volume of the encapsulating structure, $\rho_l(T_m)$ is the liquid PCM density at the melting temperature and $\xi_0$ is the initial filling level of liquid PCM. Due to thermal expansion, the volume at temperature $T$ becomes $V' = \frac{m}{\rho_l(T)} = \xi_0 \frac{\rho_l(T_m)}{\rho_l(T)} V_p$. This corresponds to the volume of the melted PCM without compression. So $V' = V_0$. The pressure is therefore given by
$$p = \frac{1}{\chi_T} \ln \xi \frac{\rho_l(T_m) V_p}{\rho_l(T)} \frac{V}{V} \qquad (8)$$

The pressure on the structure walls will expand the initial liquid volume, while the stresses induced in the elastic porous architecture will tend to avoid such a displacement. To calculate the mechanical balance between pressure and solid stresses, an implicit FSI coupling is necessary.

A static simulation based on the temperatures at the final CHT time step is considered. The pressure in the liquid is first calculated with equation 8 based on OpenFOAM. The forces undergone by the solid are exchanged and a mechanical computation is performed by CalculiX. Note that the thermal expansion of the solid structure is also considered in the simulation. The displacement of the solid PCM structure interface is sent to OpenFOAM, which first deforms the mesh and then computes the new pressure value based on the new liquid volume. The process iterates until convergence. This procedure is summarized in Fig. 1.

<figure><img src="image_4.png" /><figcaption>Fig. 1. Resolution procedure of the fluid-structure interaction simulation. The interface corresponds to the surface between the PCM and the porous architecture structure. Blue refers to OpenFOAM, ochre to CalculiX and green to preCICE.</figcaption></figure>

The implicit coupling is based on the IQN-ILS method (Quasi-Newton Inverse Least Squares) to accelerate the convergence [42].

## 3. Validation of the PCM model

The theoretical model described in section 2.1 was implemented in version 7 of OpenFOAM. Here, we present two test cases to validate this implementation, both analytically and experimentally.

### 3.1 Stefan problem

We first compare the model with the analytical solution of the Stefan problem. This consists in a semi-infinite one-dimensional geometry, which is initially considered as solid at temperature $T_0$. On the left boundary, a Dirichlet temperature boundary condition $T_d$ is enforced in order to melt the material. When $x$ tends to infinity, the temperature is supposed to be equal to the initial temperature $T_0$. Fig. 2 displays this configuration.
<figure><img src="image_5.png" /></figure>

<caption>**Fig. 2.** Stefan problem corresponding to a 1D semi-infinite material initially solid at the temperature $T_0$. The left side is then heated, which makes the material melt.</caption>

The phase change interface position $s$ and the temperature field $T$ are known analytically according to the two following relations [43]:

$$ T(x,t) = \begin{cases} T_d + (T_m - T_d) \frac{\operatorname{erf}\left(\frac{x}{2\sqrt{D_l t}}\right)}{\operatorname{erf}(\psi)} & \text{if } x < s(t) \\ T_0 + (T_m - T_0) \frac{\operatorname{erfc}\left(\frac{x}{2\sqrt{D_s t}}\right)}{\operatorname{erfc}(\psi\zeta)} & \text{otherwise} \end{cases} \quad (9) $$

$$ T(x,t) = \begin{cases} T_d + (T_m - T_d) \frac{\operatorname{erf}\left(\frac{x}{2\sqrt{D_l t}}\right)}{\operatorname{erf}(\psi)} & \text{if } x < s(t) \\ T_0 + (T_m - T_0) \frac{\operatorname{erfc}\left(\frac{x}{2\sqrt{D_s t}}\right)}{\operatorname{erfc}(\psi\zeta)} & \text{otherwise} \end{cases} \quad (10) $$

$$ s(t) = 2\psi\sqrt{D_l t} \quad (11) $$

where $\zeta = \sqrt{D_l/D_s}$ and $\psi$ is the solution of the transcendental equation

$$ \frac{\lambda_s}{\lambda_l} \zeta (T_m - T_0) \frac{e^{-\psi^2 \zeta^2}}{\operatorname{erfc}(\psi\zeta)} + (T_m - T_0) \frac{e^{-\psi^2}}{\operatorname{erf}(\psi)} + \psi \Delta h_f \sqrt{\pi} \rho_s \frac{D_l}{\lambda_l} = 0 \quad (12) $$

Let us consider the following properties and parameters of water: $\rho_l = 1000 \text{ kg.m}^{-3}$, $\rho_s = 920 \text{ kg.m}^{-3}$, $\lambda_l = 0.5 \text{ W.m}^{-1}.\text{K}^{-1}$, $\lambda_s = 1.3 \text{ W.m}^{-1}.\text{K}^{-1}$, $c_{p,l} = 4.2 \text{ kJ.kg}^{-1}.\text{K}^{-1}$, $c_{p,s} = 2 \text{ kJ.kg}^{-1}.\text{K}^{-1}$, $\Delta h_f = 333.4 \text{ kJ.kg}^{-1}$, $T_m = 0^{\circ}\text{C}$, $T_0 = -5^{\circ}\text{C}$ and $T_d = 5^{\circ}\text{C}$.

Fig. 3.a displays the phase change interface position against time, while Fig. 3.b shows the relative error with the analytical solution. The analytical solution is compared with the results from the corresponding simulation with OpenFOAM. In the figure, $dx$ refers to the spatial step and $dt$ to the temporal step in the simulation. The two most refined simulations lead to similar accurate results.
366

<figure><img src="image_6.png" /><figcaption>**Fig. 3.** a) Comparison of the phase change interface position displacement calculated analytically and with three OpenFOAM simulations. b) Relative error of the phase change interface position.</figcaption></figure>

367

368

369

370

Fig. 4 presents the temperature field one hour after the start of heating. The theoretical interface position is also presented. The OpenFOAM (OF) solution corresponds to the intermediate refinement level in Fig. 3. The simulated temperature is almost equal to the analytical solution, all along the 1D geometry. The energy equation in the set of equations (1) thus gives consistent results.

371

372

373

374

375

376

<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Comparison of the analytical and the simulated temperature field, after one hour.</figcaption></figure>

377

378

379

## 3.2 2D validation with experimental and numerical results

380

381

382

383

384

385

386

This section concerns the validation of the model, which can be 3D, based on the 2D numerical results presented in [26], which have been successfully compared with the corresponding experiment carried out by the same authors [44]. The test case consists in melting a PCM within a square which is heated on the left side at 308.15 K and maintained at 298.15 K, the initial temperature, on the right side. The bottom and top edges are considered to be adiabatic in the simulation discussed in [26], whereas they actually result from
387

conduction in thin plates in the experiment. In our simulation, adiabatic walls are considered.

388

The same PCM properties given in [26] are considered in the present simulation. Note that

389

there is a small opening at the upper left corner in the experimentation to drain the excess

390

amount of fluid resulting from the volume expansion during melting. Fig. 5 compares the

391

phase change interface at different times from our simulation with the simulation and the

392

experimentation achieved by Faden et al. [26] The mesh is composed of 40,000 cells; the

393

same mesh is used in the reference simulation.

<figure><img src="image_8.png" /></figure>

394

<caption>**Fig. 5.** Phase change interface 1, 2, 3 and 4 hours after the start of heating. The dot points and the crosses stem from [26].</caption>

395

396

Fig. 6 also displays the velocity field in the physical domain after four hours. The same figure is proposed in [26] and is very similar. It is worth noting that our results are close to the numerical simulation managed by Faden et al. [26], and they are also consistent with the experimental results [26], [44].

397

398

399

400
<figure><img src="image_9.png" /><figcaption>Fig. 6. Velocity magnitude in the geometry four hours after the start of heating. The OF simulation is based on a mesh composed of 40,000 cells. The left side is heated to 10 K warmer than the initial temperature.</figcaption></figure>

The reference simulation is based on the compressible Navier-Stokes equations, while ours is based on the Boussinesq approximation. The prior results clearly show that the Boussinesq approximation is reliable for this range of temperature variation and is satisfactory in comparison with compressible equations, which are more time-consuming to solve.

# 4. Numerical setup

## 4.1 Configuration

The geometry of the encapsulating structure is based on triply-periodic minimal surfaces, which optimize the heat transfer exchange surface [32]. Here, we consider a gyroid surface, given by the implicit equation $cos(x) sin(y) + cos(y) sin(z) + cos(z) sin(x) = 0.15$. Giving a thickness to this surface leads to the geometry in Fig. 7. This geometry contains two isolated volumes.
<figure><img src="image_10.png" /></figure>

420

421

422

**Fig. 7.** Thick gyroid volume where the thickness of branches is equal to 2% of the bounding box side length.

423

424

425

426

427

428

Based on this idea, we propose to build a porous architecture containing a volume like in Fig. 7 enclosed by a cube with the same thickness as the thick faces of the gyroid branches. This geometry is then filled with PCM. Fig. 8 displays the general geometry that is investigated in the remaining part of the paper. The side length of the surrounding box is equal to 5 cm.

<figure><img src="image_11.png" /></figure>

429

430

431

**Fig. 8.** Geometry investigated in the paper. The blue parts and the surrounding box are solid, while the remaining part is liquid.

432

433

The porous architecture will be made of copper, while the PCM will be palmitic acid. The melting temperature of palmitic acid is slightly below the usual critical temperature
of electronic components. Copper allows a fast thermal conduction process which makes heat transfer efficient. Tab. 1 itemizes the considered properties for both materials. The palmitic acid properties mainly stem from [45], except for thermal conductivities [46]. The gravitational acceleration is assumed to be equal to $g = 9.81 \text{ m.s}^{-2}$ in the simulations.

<table><caption>**Tab. 1.** Properties of the materials involved in the simulations. The blank cells are not considered or have no physical meaning.</caption><thead><tr><th>Properties</th><th>Solid palmitic acid</th><th>Liquid palmitic acid</th><th>Copper</th></tr></thead><tbody><tr><td>$T_m$ (K)</td><td>335.35</td><td>-</td><td>-</td></tr><tr><td>$\Delta h_f$ (kJ.kg<sup>-1</sup>)</td><td>211.85</td><td>-</td><td>-</td></tr><tr><td>$\rho$ (kg.m<sup>-3</sup>)</td><td>989.6</td><td>-0.7349(T - 273.15) + 907.54</td><td>8950</td></tr><tr><td>$c_p$ (J.g<sup>-1</sup>.K<sup>-1</sup>)</td><td>Fig. 3 [45]</td><td>2.39</td><td>380</td></tr><tr><td>$\lambda$ (W.m<sup>-1</sup>.K<sup>-1</sup>)</td><td>0.22</td><td>0.16</td><td>401</td></tr><tr><td>$\mu$ (Pa.s)</td><td>-</td><td>4.381 x 10<sup>-6</sup>e<sup>2488.8/T</sup></td><td>-</td></tr><tr><td>$\beta$ (K<sup>-1</sup>)</td><td>-</td><td>0.11</td><td>1.7 x 10<sup>-6</sup></td></tr><tr><td>$E$ (GPa)</td><td>-</td><td>-</td><td>130</td></tr><tr><td>$\eta$</td><td>-</td><td>-</td><td>0.33</td></tr></tbody></table>

In sections 5 & 6, several thicknesses for the gyroid volume will be tested to compare the mechanical stresses. The thicknesses will vary from 2 to 10% of the box side length. The filling level of the structure $\xi_0$ will also vary between 98.5% and 100%. Below these values, we can check that the thermal expansion of the liquid does not involve compression under 85°C (equation 8), which will not be exceeded in this paper. Higher temperatures damage electronic components, so there is no interest in exceeding this temperature.

In section 5, the bottom face of the box will be heated with a constant power equal to 30 W which is a common power experienced by electronic components. Even if this power involves temperatures greater than 85°C, the simulations will be stopped at this temperature. The other external faces of the box will undergo a convective heat transfer with a coefficient $h = 5 \text{ W.m}^{-2}.\text{K}^{-1}$ and an outside temperature of 22°C (which is also the initial temperature). The faces will be free to move.

A few studies focus mainly on the measurement of the adiabatic compressibility of palmitic acid. Due to uncertainties on its value, two FSI simulations will be performed for each configuration. They will rely on $\chi_T = 4 \times 10^{-10} \text{Pa}^{-1}$ and $\chi_T = 10^{-9} \text{Pa}^{-1}$ [47], [48], [49], [50].

## 4.2 Mesh convergence

### 4.2.1 Porous architecture structure

Mesh convergence is achieved on the solid part of the geometry (Fig. 8). Only for mesh convergence, the bottom face is fixed and heated to 90°C. The same boundary conditions as in section 4.1 are considered elsewhere. The properties of the copper are summed up in Tab.1. The static results relying on six meshes are compared in Fig. 9. The meshes consist of 185, 320 and 580 thousand first or second order elements respectively. The element sizes are respectively equal to 1 mm, 0.8 mm and 0.65 mm. The comparison is made along a curvilinear line in red located in the middle plane in green in Fig. 9.a. The displacement, von Mises stress and temperature along this curve are presented in Fig. 9. There is almost no difference on the displacement (Fig. 9.b). The von Mises stress also varies weakly with the meshes (Fig. 9.c). The greatest deviations appear close to the bottom face, where the stress is most significant. The simulated temperatures are grouped according to the number of elements rather than the order of the elements (Fig. 9.d). The order does influence the results
473
474
475
476
477

considerably for the three variables. The number of elements also has a weak impact. Since it does not lead to a prohibitive calculation time, we have chosen the second mesh with 320,000 order 1 elements rather than the coarsest one, where there are only 2 tetrahedra in the thickness of the surrounding box, which decreases the stress values close to the bottom face.

<figure><img src="image_12.png" /><figcaption>b</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>c</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>d</figcaption></figure>

<figure><img src="image_15.png" /></figure>

478
479
480
481
482

**Fig. 9.** a) Comparison of the results from six different meshes along the red curve in the green plane. The origin of the curve is on the bottom face (red point). b) Evolution of the structure displacement along the red curve. c) Evolution of the von Mises stress along the red curve. d) Change in the temperature along the red curve (1 k = 1000).

### 4.2.2 PCM

483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498

Here, only the volume associated with the PCM is studied. The material properties in the simulations are gathered in Tab. 1. Like in section 4.2.1, the bottom face is heated to 90°C while zero Neumann boundary conditions are imposed on the other faces. The PCM is initially solid at 22°C. The phase change interface position is plotted along the red line in Fig. 10.a. It should be noted that it was not possible to choose a planar curve considering the winding PCM geometry. Three meshes were compared with 40 thousand, 800 thousand and 6.8 million cells, corresponding to cell sizes of 1 mm, 0.5 mm and 0.25 mm respectively. There are expected differences in the temperature field during the melting process, 33 minutes after the start of heating (Fig. 10.b). Except at this moment, the gaps remain small. Concerning the change phase interface position, the deviations increase with time (Fig. 10.c). They are linked with the observation of the temperature around the delay on the melting process. Since the computation time of the finest mesh is very long (4 days with 192 processors), we consider that the accuracy given by this mesh is not sufficient. In addition, coupled CHT-FSI simulations will be explored, which increases calculation time even more. We thus chose the medium mesh.
a

<figure><img src="image_16.png" /></figure>

b

<figure><img src="image_17.png" /></figure>

c

<figure><img src="image_18.png" /></figure>

500

<caption>**Fig. 10.** a) PCM geometry. The phase change interface position is plotted along the red line inside the geometry. The origin of the curve is on the bottom face (red point). b) Temperature according to time at point (0.02, 0.005, 0.02). This point was chosen close to the bottom face, where the heating process begins early. c) Phase change interface position according to time.</caption>

501

502

503

504

505

506

# 5. Results

507

## 5.1 Conjugate Heat Transfer simulations

508

509

510

511

512

513

514

515

516

517

518

This section focuses on the thermal coupling exclusively, i.e., on the heat transfer and the melting process of the palmitic acid in the porous architectured conductive structure. Five simulations with different gyroid thicknesses (from 2 to 10% of the box side length) are performed. The general behavior of the coupled phenomena is described in specific configurations. Fig. 11 displays the temperature in the middle cross-section of the gyroid box (8% thickness), 17 minutes after the start of heating. The copper temperature, between 76 and 78°C, is far greater than that of the PCM. This stems from the thermal conductivity of copper which is 2,000 times higher than for the PCM (Tab. 1). Therefore, even if only the bottom face is heated, the PCM is actually heated mainly by the whole interface with the copper gyroid. Since the PCM melts at 62.2°C, a significant part of the palmitic acid is still solid as it is quite far from edges.
a

<figure><img src="image_19.png" /></figure>

b

<figure><img src="image_20.png" /></figure>

520

<caption>**Fig. 11.** Temperature field both in the PCM and copper, 17 minutes after the start of heating, in the middle cross-section of the porous architecture. The configuration corresponds to a gyroid thickness of 8%. The white contour corresponds to the limit between PCM (OpenFOAM) and copper (CalculiX). a) The temperature scale is adapted to the PCM temperature. b) The temperature scale is adapted to the solid temperature.</caption>

521

522

523

524

525

526

527

528

529

530

531

532

533

The convection process thus begins only from approximately 15 minutes after the start of heating, when the PCM close to the solid interface is melted. Fig. 12 shows local velocity vectors within the liquid PCM at the same time as in Fig. 11. Due to temperature gradients, the hot liquid moves upward, following the winding structure geometry. The top part of the PCM therefore melts faster than the bottom, despite the fact it is the bottom face that is heated. This process is also observed in Fig. 11.a, in the top left and right corners of the PCM. The upper parts of these areas melt while the lower parts are still solid. The velocities remain smaller than $4 \text{ mm.s}^{-1}$ at this time but can exceed $1 \text{ cm.s}^{-1}$ during the simulation.
<figure><img src="image_21.png" /><figcaption>**Fig. 12.** Velocity vectors in a melting area of the PCM, 17 minutes after the start of heating the structure composed of an 8% thickness gyroid.</figcaption></figure>

Fig. 13 shows the change in the solid PCM inside the gyroid structure composed of 6% thick branches. Following the previous results, a melting process starting from the PCM-copper interface towards the center of the PCM is expected. This behavior is observed in Fig. 13. After 10 minutes (Fig. 13.a), only the PCM close to the edges begins to melt, as in Fig. 10. Note that the temperature in the liquid PCM exceeds 65°C. Seven minutes later (Fig. 13.b), the upper part of the PCM is hotter than the lower part, due to convection. The volume of the solid PCM has decreased drastically. Three minutes later (Fig. 13.c), the solid PCM is shrinking sharply, while the temperature varies weakly in comparison with Fig. 13.b. Only the central part of the PCM remains solid after 20 minutes. Finally, after 23 minutes, only a small part of the PCM is still in the solid state. The temperature will then increase faster without PCM to melt. A few seconds later, the solid PCM disappears.
<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /><figcaption>**Fig. 13.** OpenFOAM temperature field and solid PCM (in black) in a copper gyroid structure consisting of 6% thick branches. a) 10 minutes; b) 17 minutes; c) 20 minutes; d) 23 minutes; after the start of heating.</figcaption></figure>

Fig. 14.a shows the average temperature in the whole PCM according to time, for the five thickness configurations. Before most of the PCM has melted, the lower the thickness, the higher the temperature. During this conduction time, the energy is transferred through the interface. The lower the thicknesses, the greater the exchange interface area. Therefore, there is more conductive heat transfer for copper gyroid structures with thin branches. In addition, for thick gyroid structures, the bounding box is also thicker, so energy coming from the bottom face must diffuse in a thicker copper volume before reaching the PCM. Energy thus attains the PCM faster in thin gyroid structures. Very quickly, however, after several hundreds of seconds (Fig. 14.b), the heating rate becomes lower for thin gyroid structures. Since there is

548

549

550

551

552

553

554

555

556

557

558

559

560
561

more PCM, the melting process plays a more significant role than in thicker configurations.
562
The heating rate thus reduces substantially down to a quasi-constant heating rate of around
563
1°C.min⁻¹. This second phase (melting phase) occurs earlier for thin copper gyroid structures
564
since there is more PCM. Due to the quantity, the melting phase is also longer in thin
565
configurations. Fig. 14.a shows greater damping of temperature variation inside thin copper
566
gyroid structures, thanks to the PCM. Finally, a third phase happens, when the PCM is mostly
567
or totally liquid. The heating rate increases again (Fig. 14.b). Since this last phase begins later
568
for the 2% thick gyroid, the temperature is lower than in the other simulations. At the end, the
569
thinner it is, the colder the temperature.

<figure><img src="image_26.png" /></figure>

<figure><img src="image_27.png" /></figure>

570

<caption>Fig. 14. a) Average temperature in the PCM according to time for the five thickness cases. b) Heating rate (time derivative of temperature) in the PCM according to time.</caption>

571

572

## 5.2 Fluid-structure interaction simulations

573

As described in section 2, we consider a fluid-structure simulation at the end of the melting
574
process of the PCM, when the temperature is almost homogeneous. In the simulation, the final
575
temperature field in the porous architectured structure is considered while the liquid
576
temperature is assumed to be equal to 85°C. The liquid temperature is indirectly involved in
577
equation (8) through the density, and thus has a significant impact on thermal volume
578
expansion. The copper temperature is also involved in the thermal dilatation and generates
579
mechanical stresses in the structure.
580

575

576

577

578

579

580

The deformation of the porous architectured structure (4% thick), completely filled with PCM
581
(ξ₀ = 100%), is presented in Fig. 15 with a scale factor of 30. Fig. 15.b compares a section
582
before and after deformation. The greatest displacements occur at the box sides due to the
583
pressure gradients. The pressure is indeed significant inside the box while the atmospheric
584
pressure is negligible outside. Instead, the pressure is constrained inside gyroid branches on
585
both sides. Note that the side displacement remains relatively small (240 µm).
586

581

582

583

584

585

586
<figure><img src="image_28.png" /><figcaption>**Fig. 15.** Displacement field of the 4% thick gyroid structure initially fully filled with the PCM ($\xi_0$ = 100%), from CalculiX. The deformations presented here correspond to the genuine deformation multiplied by 30. a) Displacement of the 3D copper gyroid structure. b) Comparison of the deformed structure with the initial geometry in the central section of the copper gyroid structure.</figcaption></figure>

Fig. 16 displays the von Mises stresses inside the porous architecture for the same configuration. The greatest stresses do not occur at the box sides but in the thick branches, especially close to the junctures. These parts undergo compressive forces due to the pressure. The stresses reach as much as 570 MPa at several locations. The value of the compressive yield strength of copper can vary according to the material properties. In the paper, we assume that its value is 138 MPa [51]. Based on the von Mises criterion, this configuration thus widely exceeds safe stresses and should lead at least to plastic deformations and perhaps failures.
<figure><img src="image_29.png" /><figcaption>**Fig. 16.** The von Mises stresses from CalculiX in the 4% thick gyroid structure initially fully filled with the melted PCM ($\xi_0$ = 100%).</figcaption></figure>

# 6. Discussion

## 6.1 Stresses involved in the encapsulating solid

The model presented allows us to compute the stresses generated by the thermal expansion of the PCM in the encapsulating structure. Based on these results, it is possible to predict whether a configuration is reliable or not from this point of view. In the paper, we focus on the yield stress and not on failure. Since such a product must undergo several cycles, plastic deformation should not occur for safe and efficient use. Like in the prior section, we follow the von Mises criterion. Fig. 17 displays the maximal von Mises stress in copper for each configuration that was tested. We recall that simulations with different filling levels and thicknesses of gyroid branches were carried out. For each configuration, a range is provided for the maximal von Mises stress, corresponding to the FSI simulations with two values of adiabatic compressibility $\chi_T = 4 \times 10^{-10} Pa^{-1}$ and $\chi_T = 10^{-9} Pa^{-1}$ (§4.1). With a branch thickness of only 2% of the box size length, the stresses are very high and exceed the yield stress for every filling level greater than 99%. With less PCM in the host architecture matrix, its volume expansion does not involve sufficient stresses to deform the copper plastically. This is an important issue since this geometry is the hollowest one, so filled with the greatest amount of PCM. With a thickness of 4% to 8%, the stresses decrease, as expected, but to a level which is not sufficient to change the previous conclusions on 2% thick gyroids. Indeed, from a filling level of 99%, a risk of plastic deformation still exists. With a thickness of 10%, the solid structure can undergo high pressures with von Mises stresses remaining smaller than the yield stress. With this thickness, only a completely filled structure could be deformed plastically. However, such a geometry contains a smaller amount of PCM than for a gyroid composed of 2% thick branches. Considering the maximum filling level of each geometry without jeopardizing the structure, the 2% thick geometry remains the structure that can be filled with the greatest amount of PCM. With a filling level of 98.5%, this geometry can
contain 113 $cm^3$ of PCM while remaining safe according to the simulations. In comparison, the thicker geometry filled to a level of 99.5% can only contain 82 $cm^3$ of PCM. It appears that the filling level is a more interesting parameter since it drastically impacts the stresses involved, while decreasing the quantity of PCM only slightly.

<figure><img src="image_30.png" /><figcaption>**Fig. 17.** von Mises stresses for each configuration. The yield stress is displayed in red. The error bars are associated with the uncertainties on adiabatic compressibility. They correspond to two FSI simulations based on two different values of adiabatic compressibility.</figcaption></figure>

## 6.2 Thermal contact resistance

Excessive heat can severely impact electronic devices. The proposed geometry filled with PCM is used to reduce the temperature. However, to be efficient, the contact between the innovative composite and devices must be efficient, i.e., the contact surface must be maximized. The deformation of the solid structure decreases the exchange area and leads to a loss in heat transfer (Fig. 18). Thermal contact resistance is defined as the ratio between the peak and valley $\delta$ (see Fig. 18) and the thermal conductivity of the medium between the innovative composite and the electronic component (often air) $R_{th} = \frac{\delta}{\lambda_{air}}$. On the example in Fig. 18, $\delta \sim 1 \text{ mm}$, the temperature is approximately 85°C, which yields a thermal conductivity of $\lambda_{air} \sim 0.03 \text{ W.m}^{-1}.\text{K}^{-1}$, and $R_{th} \sim 0.033 \text{ K.m}^2.\text{W}^{-1}$. The thermal contact resistance is very significant in this configuration which is the worst one tested in the paper and the most adverse case from electronic devices. Fig. 19 displays the displacements according to the gyroid thickness for every filling level.
<figure><img src="image_31.png" /><figcaption>**Fig. 18.** Deformed 2% thick gyroid structure completely filled with PCM (scaled up by a factor of 10). The contact with an electronic component is less efficient due to deformation. $\delta$, the so-called peak to valley, is the distance between the component and the structure.</figcaption></figure>

As seen in the previous section, the gyroid with 2% thick branches is the most efficient from the stress point of view. However, with greater displacements, one could expect significant thermal contact resistance. Fig. 19 shows that displacements are indeed far larger in this configuration than with thicker gyroids, except with a filling level of 98.5%. The displacements are of the order of $10^{-4} m$, like in every other configuration with thicknesses of 4 to 10%. With such displacements, the thermal contact resistance is approximately $R_{th} \sim 0.0033 \text{ K.m}^2\text{.W}^{-1}$. Note that this value remains an important issue and must impact the experimental performance. In our model, this feature is not considered and some differences between the theoretical results presented in the paper and experimentation could occur on account of this. Similarly, the contact between the PCM and copper before melting is not taken into account in the model. The heat transfer between every part of the device could therefore be slower than in the simulation, but should modify neither the structure displacement nor the stresses significantly.
<figure><img src="image_32.png" /><figcaption>**Fig. 19.** Maximal displacement for each configuration. The error bars are associated with the uncertainties on adiabatic compressibility. They correspond to two FSI simulations based on two different values of adiabatic compressibility.</figcaption></figure>

# 7. Conclusion

Electronic cooling relies on different techniques. Among them, PCMs mitigate temperature increases. This approach is therefore widely used experimentally and is a subject of huge interest in the community. However, the theoretical modeling of the coupled physical phenomena remains a challenge. We have addressed this very complex problem and proposed a complete 3D model involving CHT and FSI simulations, based on 3 open-source software codes: OpenFOAM, CalculiX and preCICE. The PCM model, which has been implemented in OpenFOAM, takes account of phase changes, conduction and convection. This model has been validated on the analytical solution of the Stefan problem and experimental data of a 2D configuration. The CHT simulations relied on the PCM model (OpenFOAM) and the conduction in the solid architectured structure (CalculiX). The interface coupling between the two codes was managed by preCICE. The FSI simulations aimed to compute the structure displacements and induced stresses. Such information is very useful to prevent mechanical failures and improve the efficiency of heat transfer at the interface with electronic components, as deformation of the structure can worsen the contact.

This modeling was applied to gyroid structures, triply-periodic minimal surfaces, which are expected to behave efficiently. Such shapes do maximize the exchange surface, thus
691 intensifying heat transfer during heating and improving thermal damping in electronic
692 components. We have simulated several configurations with different geometry thicknesses
693 and various filling levels of PCM (from 98.5 to 100%). The thicknesses of the gyroid
694 branches ranged from 2% to 10% of the container box side length. Depending on the filling
695 rate, thermal expansion of liquid PCM led to pressure increases inside the 3D architectured
696 porous conductive matrix, and therefore significant mechanical stresses. For the first time, we
697 have shown that the filling ratio has a notable impact on the stresses that are induced. Only
698 1% short of the maximum volume reduces the mechanical stresses and failure risk
699 considerably. Compared to this low value, increasing the thickness of the container cuts a lot
700 more the PCM volume, and thus thermal performance. To optimize thermal management,
701 decreasing the quantity of PCM is thus preferable to strengthening the container geometry,
702 which scales the volume down too much. As the fluid pressure remains low, the solid
703 deformation is found to be small even for thin gyroids. The simulated maximal displacement
704 with 2% thick branches is of the same order (0.1 mm), as with consolidated gyroids filled
705 with more PCM. The contact thermal resistance is therefore also of the same order of
706 magnitude.

707 With this 3D fluid-structure modeling of the heat and mass transfers in a complex system
708 dedicated to the thermal management of electronic devices, we were able to select the 2%
709 thick gyroid as the most efficient configuration. A future study based on both experimentation
710 and simulation will be carried out with this design of architectured porous structures based on
711 metallic additive manufacturing impregnated with palmitic acid. This innovative composite
712 will be characterized to specify the input parameters of the simulations and validate the
713 theoretical results. Following this study, the material will be integrated into a power supply in
714 order to evaluate its performance and its capacity to manage power peaks more effectively.

715

## Acknowledgments

717 We thank the University of Bordeaux and University of Pau and the Adour Region for
718 providing access to MCIA (Mésocentre de Calcul Intensif Aquitain).

719

## 8. References

720

721

[1] A. Sharma, V. Tyagi, C. Chen et D. Buddhi, «Review on thermal energy storage with phase change materials and applications,» *Renewable and Sustainable Energy Reviews*, vol. 13, n°12, 318-345, 2009.

[2] N. Solangi, N. Mubarak, R. Karri, S. Mazari, A. Jatoi, J. Koduru et M. Dehghani, «MXene-based phase change materials for solar thermal energy storage,» *Energy Conversion and Management*, vol. 273, 116432, 2022.

[3] V. Palomba, V. Brancato et A. Frazzica, «Thermal performance of a latent thermal energy storage for exploitation of renewables and waste heat: An experimental investigation based on an asymmetric plate heat exchanger,» *Energy Conversion and*
*Management*, vol. 200, 112121, 2019.

[4] M. Farid, A. Khudhair, S. Razack et S. Al-Hallaj, «A review on phase change energy storage: Materials and applications,» *Energy Conversion and Management*, vol. 45, n°19-10, 1597-1615, 2004.

[5] R. Jia, K. Sun, R. Li, Y. Zhang, W. Wang, H. Yin, D. Fang, Q. Shi et Z. Tan, «Heat capacities of some sugar alcohols as phase change materials for thermal energy storage applications,» *Journal of Chemical Thermodynamics*, vol. 115, 233-248, 2017.

[6] A. Khudhair et M. Farid, «A review on energy conservation in building applications with thermal storage by latent heat using phase change materials,» *Energy Conversion and Management*, vol. 45, n°12, 263-275, 2004.

[7] E. Mohseni, W. Tang et S. Wang, «Development of thermal energy storage lightweight structural cementitious composites by means of macro-encapsulated PCM,» *Construction and Building Materials*, vol. 225, 182-195, 2019.

[8] W. Qureshi, N.-K. Nair et M. Farid, «Impact of energy storage in buildings on electricity demand side management,» *Energy Conversion and Management*, vol. 52, n°15, 2110-2120, 2011.

[9] H. e. a. Nazir, «Recent developments in phase change materials for energy storage applications: A review,» *International Journal of Heat and Mass Transfer*, vol. 129, 491-523, 2019.

[10] G. Alva, Y. Lin et G. Fang, «An overview of thermal energy storage systems,» *Energy*, vol. 144, 341-378, 2018.

[11] S. K. Sahoo, M. K. Das et P. Rath, «Application of TCE-PCM based heat sinks for cooling of electronic components: A review,» *Renewable and Sustainable Energy Reviews*, vol. 59, 550-582, 2016.

[12] Y. Liu, R. Zheng et J. Li, «High latent heat phase change materials (PCMs) with low melting temperature for thermal management and storage of electronic devices and power batteries: Critical review,» *Renewable and Sustainable Energy Reviews*, vol. 168, 112783, 2022.

[13] C. Kinkelin, S. Lips, U. Soupremanien, V. Remondière, J. Dijon, H. Le Poche, E. Ollier, M. Zegaoui, N. Rolland, P.-A. Rolland, S. Lhostis, B. Descouts, Y. Kaplan et F. Lefèvre, «Theoretical and experimental study of a thermal damper based on a CNT/PCM composite structure for transient electronic cooling,» *Energy Conversion and Management*, vol. 142, 257-271, 2017.

[14] W. Li, F. Wang, W. Cheng, X. Chen et Q. Zhao, «Study of using enhanced heat-transfer flexible phase change material film in thermal management of compact electronic device,» *Energy Conversion and Management*, vol. 210, 112680, 2020.

[15] S. F. Hosseinizadeh, F. L. Tan et S. M. Moosania, «Experimental and numerical studies
on performance of PCM-based heat sink with different configurations of internal fins,»
*Applied Thermal Engineering*, vol. 31, n°117-18, 3827-3838, 2011.

[16] J. Krishna, P. S. Kishore et A. B. Solomon, «Heat pipe with nano enhanced-PCM for electronic cooling application,» *Experimental Thermal and Fluid Science*, vol. 81, 84-92, 2017.

[17] S. Gharbi, S. Harmand et S. Jabrallah, «Experimental comparison between different configurations of PCM based heat sinks for cooling electronic components,» *Applied Thermal Engineering*, vol. 87, 454-462, 2015.

[18] P. Huang, G. Wei, L. Cui, C. Xu et X. Du, «Experimental and numerical optimization of cascaded PCM heat sink by using low melting point alloys,» *Energy Conversion and Management*, vol. 269, 116149, 2022.

[19] A. Farzanehnia, M. Khatibi, M. Sardarabadi et M. Passandideh-Fard, «Experimental investigation of multiwall carbon nanotube/paraffin based heat sink for electronic device thermal management,» *Energy Conversion and Management*, vol. 179, 314-325, 2019.

[20] J. Ho, Y. See, K. Leong et T. Wong, «An experimental investigation of a PCM-based heat sink enhanced with a topology-optimized tree-like structure,» *Energy Conversion and Management*, vol. 245, 114608, 2021.

[21] N. Soares, N. Rosa, J. Costa, A. Lopes, T. Matias, P. Simões et L. Durães, «Validation of different numerical models with benchmark experiments for modelling microencapsulated-PCM-based applications for buildings,» *International Journal of Thermal Sciences*, vol. 159, 106565, 2021.

[22] A. Diani, C. Nonino et L. Rossetto, «Melting of phase change materials inside periodic cellular structures fabricated by additive manufacturing: Experimental results and numerical simulations,» *Applied Thermal Engineering*, vol. 215, 118969, 2022.

[23] Y. Zhuang, T. Chen, J. Chen, J. Li, M. Guan et Y. Chen, «Thermal uniformity performance of a hybrid battery thermal management system using phase change material and cooling plates arrayed in the manner of honeycomb,» *Thermal Science and Engineering Progress*, vol. 26, 101094, 2021.

[24] T. Zhou, Y. Xiao, Y. Liu, J. Lin et H. Huang, «Research on cooling performance of phase change material-filled earth-air heat exchanger,» *Energy Conversion and Management*, vol. 177, 210-223, 2018.

[25] «OpenFOAM-7,» 2022. [En ligne]. Available: https://openfoam.org/.

[26] M. Faden, A. König-Haagen et D. Brüggemann, «An optimum enthalpy approach for melting and solidification with volume change,» *Energies*, vol. 12, 868, 2019.

[27] A. Saraswat, R. Bhattacharjee, A. Verma, M. Das et S. Khandekar, «Investigation of diffusional transport of heat and its enhancement in phase-change thermal energy storage systems,» *Applied Thermal Engineering*, vol. 111, 1611-1621, 2017.
[28] Y. Harmen, Y. Chhiti, F. M’Hamdi Alaoui, F. Bentiss, C. Jama et S. B. M. Duquesne, «Thermal performance of PEG-MWCNTs composites as shape-stabilised phase change materials for thermal energy storage,» *Fullerenes Nanotubes and Carbon Nanostructures*, vol. 29, n°19, 732-738, 2021.

[29] E. Alehosseini et S. Jafari, «Nanoencapsulation of phase change materials (PCMs) and their applications in various fields for energy storage and management,» *Advances in Colloid and Interface Science*, vol. 283, 102226, 2020.

[30] C. Michler, S. Hulshoff, E. van Brummelen et R. de Borst, «A monolithic approach to fluid-structure interaction,» *Computers & Fluids*, vol. 33, n°15-6, 839-848, 2004.

[31] H. Matthies et J. Steindorf, «Partitioned strong coupling algorithms for fluid-structure interaction,» *Computers & Structures*, vol. 81, n°18-11, 805-812, 2003.

[32] S. Blanquer, M. Werner, M. Hannula, S. Sharifi, G. Lajoinie, D. Eglin, J. Hyttinen, A. Poot et D. Grijpma, «Surface curvature in triply-periodic minimal surface architectures as a distinct design parameter in preparing advanced tissue engineering scaffolds,» *Biofabrication*, vol. 9, 025001, 2017.

[33] G. Dhondt, 2022. [En ligne]. Available: http://www.calculix.de.

[34] H.-J. Bungartz, F. Lindner, B. Gatzhammer, M. Mehl, K. Scheufele, A. Shukaev et B. Uekermann, «preCICE - A fully parallel library for multi-physics surface coupling,» *Computers & Fluids*, vol. 141, 250-258, 2016.

[35] W. Alshaer, S. A. Nada, M. A. Rady, C. Le Bot et E. Palomo Del Barrio, «Numerical investigations of using carbon foam/PCM/Nano carbon tubes composites in thermal management of electronic equipment,» *Energy Conversion and Management*, vol. 89, 873-884, 2015.

[36] Z. Rao, S. Wang et G. Zhang, «Simulation and experiment of thermal energy management with phase change material for ageing LiFePO4 power battery,» *Energy Conversion and Management*, vol. 52, 3408-3414, 2011.

[37] W. Li, T. Zhang, B. Li, F. Cui et L. Liu, «Experimental investigation on combined thermal energy storage and thermoelectric system by using foam/PCM composite,» *Energy Conversion and Management*, vol. 243, n°1114429, 2021.

[38] W. Wu, G. Zhang, X. Ke, X. Yang, Z. Wang et C. Liu, «Preparation and thermal conductivity enhancement of composite phase change materials for electronic thermal management,» *Energy Conversion and Management*, vol. 101, 278-284, 2015.

[39] W. Wu, X. Yang, G. Zhang, K. Chen et S. Wang, «Experimental investigation on the thermal performance of heat pipe-assisted phase change material based battery thermal management system,» *Energy Conversion and Management*, vol. 138, 486-492, 2017.

[40] E. Palomo Del Barrio et J.-L. Dauvergne, «New method to characterize phase change materials,» *Advanced Sciences and Technologies*, vol. 74, 243-252, 2010.
[41] F. Jelassi, M. Azaïez et E. Palomo Del Barrio, «A substructuring method for phase change modelling in hybrid media,» *Computers and Fluids*, vol. 88, 81-92, 2013.

[42] F. Salmon et L. Chatellier, «3D fluid-structure interaction simulation of an hydrofoil at low Reynolds number,» *Journal of Fluids and Structures*, vol. 111, 103573, 2022.

[43] D. W. Hahn and M. Necati Ozisik, *Heat Conduction*, Third Edition, Wiley, 2012.

[44] M. Faden, C. Linhardt, S. Höhlein, A. König-Haagen et D. Brüggemann, «Velocity field and phase boundary measurements during melting of n-octadecane in a cubical test cell,» *International Journal of Heat and Mass Transfer*, vol. 135, 104-114, 2019.

[45] M. Duquesne, C. Mailhé, S. Doppiu, J.-L. Dauvergne, S. Santos-Moreno, A. Godin, G. Fleury, F. Rouault et E. Palomo del Barrio, «Characterization of Fatty Acids as Biobased Organic Materials for Latent Heat Storage,» *Materials*, vol. 14, 4707, 2021.

[46] J. Wang, H. Xie, Z. Xin, Y. Li et L. Chen, «Enhancing thermal conductivity of palmitic acid based phase change materials with carbon nanotubes as fillers,» *Solar Energy*, vol. 84, 339-344, 2010.

[47] A. Balcerzak, «Comparison of High-Pressure Behavior of Physicochemical Properties of the Di- and Triacylglycerols Established by Ultrasonic Methods,» *Journal of the American Oil Chemists' Society*, vol. 94, 1261-1268, 2017.

[48] A. J. Rostocki, R. Tarakowski, P. Kielczynski, M. Szalewski, A. Balcerzak et S. Ptasznik, «The Ultrasonic Investigation of Phase Transition in Olive Oil,» *Journal of the American Oil Chemists' Society*, vol. 90, 813-818, 2013.

[49] M. Banu, Siddaramaiah et N. Prasadc, «Radical Scavenging Activity of Tertiary Butyl Hydroquinone and Assessment of Stability of Palm Oil (Elaeis guineensis) by Ultrasonic Studies,» *Chemical and Biochemical Engineering Quarterly*, vol. 30, n°14, 477-487, 2016.

[50] E. H. I. Ndiaye, D. Nasri et J.-L. Daridon, «Speed of Sound, Density, and Derivative Properties of Fatty Acid Methyl and Ethyl Esters under High Pressure: Methyl Caprate and Ethyl Caprate,» *Journal of chemical & engineering data*, vol. 57, n°110, 2667-2676, 2012.

[51] Y. Tang, X. Yang, R. Wang et M. Li, «Enhancement of the mechanical properties of graphene-copper composites with graphene-nickel hybrids,» *Materials Science and Engineering: A*, vol. 599, 247-254, 2014.