Example : Electronic component cooling

This test case has been proposed by Annabelle Le-Hyaric and Michel Fouquembergh formerly at AIRBUS.

We consider a 2D model representative of the neighboring of an electronic component submitted to a cooling air flow. It is described by four geometrical domains in \(\mathbb{R}^2\) named \(\Omega_i,i=1,2,3,4\), see figure. We suppose the velocity \(\vec{v}\) is known in each domain --- for instance in \(\Omega_4\) it is the solution of previous Navier-Stokes computations. --- The temperature \(T\) of the domain \(\Omega = \cup_{i=1}^4 \Omega_i\) is then solution of heat transfer equation :

\[\begin{equation} \label{eq:1} \rho C_i \Big( \frac{\partial T}{\partial t} + \vec{v} \cdot \nabla T \Big) - \nabla \cdot \left( k_i \nabla T \right) = Q_i,\quad i=1,2,3,4 \end{equation}\]

where \(t\) is the time and in each sub-domain \(\Omega_i\), \(\rho C_i\) is the volumic thermal capacity, \(k_i\) is thermal conductivity and \(Q_i\) is a volumic heat dissipated.

eads geometry

One should notice that the convection term in heat transfer equation may lead to spatial oscillations which can be overcome by Petrov-Galerkin type or continuous interior penalty stabilization techniques.

Integrated circuits (ICs) (domains \(\Omega_1\) and \(\Omega_2\) ) are respectively soldered on PCB at \(\mathbf{x1}=(e_{Pcb}, h_1)\) and \(\mathbf{x_2}=(e_{Pcb}, h_2)\). They are considered as rectangles with width \(e_{IC}\) and height \(h_{IC}\). The printed circuit board (PCB) is a rectangle \(\Omega_3\) of width \(e_{PCB}\) and height \(h_{PCB}\). The air(Air) is flowing along the PCB in domain \(\Omega_4\). Speed in the air channel \(\Omega_4\) is supposed to have a parabolic profile function of \(x\) coordinate. Its expression is simplified as follows (notice that \(\vec{v}\) is normally solution of Navier-Stokes equations; the resulting temperature and velocity fields should be quite different from that simplified model), we have for all \(0 \leq y \leq h_{PCB}\)

\[\begin{equation} \label{eq:2} \begin{array}[c]{rl} e_{Pcb} + e_{Ic} \leq x \leq e_{Pcb} + e_a, & \displaystyle \vec{v} = \frac{3}{2 (e_a-e_{Ic})}\ D\ \Big( 1 - \Big( \frac{x - ( \frac{e_a+e_{Ic}}{2} + e_{Pcb})}{\frac{(e_a-e_{Ic})}{2}} \Big)^2 \Big)\ f(t)\ \vec{y}\\ e_{Pcb} \leq x \leq e_{Pcb} + e_{Ic}, & \vec{v} = 0 \end{array} \end{equation}\]

where \(f\) is a function of time modelling the starting of the PCB ventilation, i. e.

\[\begin{equation} \label{eq:13} f(t) = 1-\exp(-\frac{t}{3}), \end{equation}\]

\(D\) is the air flow rate, see table~\vref{tab:1} and \(\vec{y}=(0, 1)^T\) is the unit vector along the \(y\) axis. A quick verification shows that

\[\begin{equation} \label{eq:18} \int_{\Gamma_4 \cap \Omega_4} \vec{v} \cdot \vec{n} = \int_{\Gamma_4 \cap \Omega_4} v_y = D \end{equation}\]

The medium velocity \(\vec{v}_i = \vec{0}, i=1,2,3\) in the solid domains \(\Omega_i, i=1,2,3\).

ICs dissipate heat, we have respectively

\[\begin{equation} \label{eq:12} \begin{array}[c]{rl} Q_1 &= Q \big( 1 - \exp(-t) \big)\quad \text{ in } \Omega_1\\ Q_2 &= Q \big( 1 - \exp(-t) \big)\quad \text{ in } \Omega_2\\ \end{array} \end{equation}\]

where \(Q\) is defined in table~\vref{tab:1}.

We shall denote \(\vec{n}_{|{\Omega_i}} = \vec{n}_i\) denotes the unit outward normal to \(\Omega_i\) and \(\vec{n}_{|{\Omega_j}} = \vec{n}_j\) denotes the unit outward normal to \(\Omega_j\).

1. Boundary conditions

We set

  • on \(\Gamma_3 \cap \Omega_3\), a zero flux (Neumann-like) condition

\[\begin{equation} \label{eq:10} -k_3\ \nabla T \cdot \vec{n}_3\ =\ 0; \end{equation}\]
  • on \(\Gamma_3 \cap \Omega_4\), a zero flux (Robin-like) condition

\[ \begin{equation} \label{eq:7} (-k_4\ \nabla T + \rho C_4 T \mathbf{v} ) \cdot \vec{n}_4\ =\ 0; \end{equation}\]
  • on \(\Gamma_4, (0 \leq x \leq e_{Pcb}+e_a, y=0)\) the temperature is set (Dirichlet condition)

\[\begin{equation} \label{eq:11} T\ = T_0; \end{equation}\]
  • between \(\Gamma_1\) and \(\Gamma_2\), periodic conditions

\[\begin{equation} \label{eq:9} \begin{array}{rl} T_{\big|\mathbf{x} = 0} &= T_{\big|\mathbf{x} = e_{Pcb} + e_{a}} \\ -k_3\ \nabla T \cdot \vec{n}_{3_{\big|\mathbf{x} = 0}} &= k_4\ \nabla T \cdot \vec{n}_4{_{\big|\mathbf{x} = e_{Pcb} + e_{a}}}; \end{array} \end{equation}\]
  • at interfaces between the ICs and PCB, there is a thermal contact conductance:

\[\begin{equation} \label{eq:5} \begin{array}{rl} -k_1\ \nabla T \cdot \vec{n}_1 - k_3\ \nabla T \cdot \vec{n}_3 = r_{13}\big( T_{\partial \Omega_1} - T_{\partial \Omega_3} \big)\\ -k_2\ \nabla T \cdot \vec{n}_2 - k_3\ \nabla T \cdot \vec{n}_3 = r_{23}\big( T_{\partial \Omega_2} - T_{\partial \Omega_3} \big);\\ \end{array} \end{equation}\]
  • on other internal boundaries, the coontinuity of the heat flux and temperature, on \(\Gamma_{ij} = \Omega_i \cap \Omega_j \neq \emptyset\)

\[\begin{equation} \label{eq:6} \begin{array}{rl} T_i &= T_j \\ k_i\ \nabla T \cdot \vec{n}_i &= -k_j\ \nabla T \cdot \vec{n}_j. \end{array} \end{equation}\]

2. Initial condition

At \(t=0s\), we set \(T = T_0\).

2.1. Inputs

The table displays the various fixed and variables parameters of this test-case.

Table 1. Table of fixed and variable parameters

Name

Description

Nominal Value

Range

Units

Parameters

\(t\)

time

\([0, 1500\)]

\(s\)

\(Q\)

heat source

\(10^{6}\)

\([0, 10^{6}\)]

\(W \cdot m^{-3}\)

IC Parameters

\(k_1 = k_2 = k_{IC}\)

thermal conductivity

\(2\)

\([0. 2,150 \)]

\(W \cdot m^{-1} \cdot K^{-1}\)

\(r_{13} = r_{23} = r\)

thermal conductance

\(100\)

\([10^{-1},10^{2} \)]

\(W \cdot m^{-2} \cdot K^{-1}\)

\(\rho C_{IC}\)

heat capacity

\(1. 4 \cdot 10^{6}\)

\(J \cdot m^{-3} \cdot K^{-1}\)

\(e_{IC}\)

thickness

\(2\cdot 10^{-3}\)

\(m\)

\(h_{IC} = L_{IC}\)

height

\(2\cdot 10^{-2}\)

\(m\)

\(h_{1}\)

height

\(2\cdot 10^{-2}\)

\(m\)

\(h_{2}\)

height

\(7\cdot 10^{-2}\)

\(m\)

PCB Parameters

\(k_3 = k_{Pcb}\)

thermal conductivity

\(0. 2\)

\(W \cdot m^{-1} \cdot K^{-1}\)

\(\rho C_{3}\)

heat capacity

\(2 \cdot 10^{6}\)

\(J \cdot m^{-3} \cdot K^{-1}\)

\(e_{Pcb}\)

thickness

\(2\cdot 10^{-3}\)

\(m\)

\(h_{Pcb}\)

height

\(13\cdot 10^{-2}\)

\(m\)

Air Parameters

\(T_0\)

Inflow temperature

\(300\)

\(K\)

\(D\)

Inflow rate

\(7\cdot 10^{-3}\)

\([5 \cdot 10^{-4} ,10^{-2}\)]

\(m^2 \cdot s^{-1}\)

\(k_4 \)

thermal conductivity

\(3 \cdot 10^{-2}\)

\(W \cdot m^{-1} \cdot K^{-1}\)

\(\rho C_{4}\)

heat capacity

\(1100\)

\(J \cdot m^{-3} \cdot K^{-1}\)

\(e_{a}\)

thickness

\(4\cdot 10^{-3}\)

\([2. 5 \cdot 10^{-3} , 5 \cdot 10^{-2}\)]

\(m\)

2.2. Outputs

The outputs are (i) the mean temperature \(s_1(\mu)\) of the hottest IC

\[\begin{equation} \label{eq:3} s_1(\mu) = \frac{1}{e_{IC} h_{IC}} \int_{\Omega_2} T \end{equation}\]

and (ii) mean temperature \(s_2(\mu)\) of the air at the outlet

\[ \begin{equation} \label{eq:4} s_2(\mu) = \frac{1}{e_{a}} \int_{\Omega_4 \cap \Gamma_3 } T \end{equation}\]

both depends on the solution of \eqref{eq:1} and are dependent on the parameter set \(\mu\).

We need to monitor \(s_1(\mu)\) and \(s_2(\mu)\) because \(s_1(\mu)\) is the hottest part of the model and the IC can’t have a temperature above \(340K\). \(s_2(\mu)\) is the outlet of the air and in an industrial system we can have others components behind this outlet. So the temperature of the air doesn’t have to be high to not interfere the proper functioning of these.

Testcases

We build some example tests to compute for this problem.

1. Test 1

This test is the base of all of the other tests, we bold the changes in the next tests. Only these parameters varies in these tests :

Flow rate

\(7 \cdot 10^{-3}\)

Characteristic length

\(5 \cdot 10^{-4}\)

Triangle family

P1

Temporal scheme

BDF 1, \(\Delta t = 1s\)

Stabilisation

GLS

As the periodicity helps to dissipate the heat in both directions of the PCB, the results are attended to be higher than the true results.

The Feel++ toolboxes are configured with two complementary files, the .json file configure the specific materials related to the toolbox for the geometry. The .cfg file configure the generic options, for instance the mesh file, the output directory or the solver. The .json file is common to all examples, except the stationary problem.

// -*- mode: javascript -*-
{
    "Name"      : "Eads",
    "ShortName" : "Eads",

    "Models" : {
	"use-model-name" : 1,
        "fluid" : { "equations" : "Navier-Stokes" }
    },

    "Materials" : {
        "Pcb" : {
	    "markers" : "Pcb",
            "physics" : "heat",
            "rho"     : "1.",
            "k"       : "0.2",
            "Cp"      : "2e6",
	    "mu"      : "1e12"
        },
	"IC1" : {
	    "markers" : "IC1",
	    "physics" : "heat",
	    "name"    : "IC1",
	    "rho"     : "1.",
	    "k"       : "2",
	    "Cp"      : "1.4e6",
	    "mu"      : "1e12"
	},
	"IC2" : {
	    "markers" : "IC2",
	    "physics" : "heat",
	    "rho"     : "1.",
	    "k"       : "2",
	    "Cp"      : "1.4e6",
	    "mu"      : "1e12"
	},
	"Air" : {
	    "markers" : "Air",
	    "physics" : ["heat","fluid"],
	    "rho"     : "1.",
	    "k"       : "3e-2",
	    "Cp"      : "1100",
	    "mu"      : "1.8e-5"
	}
    },

    "BoundaryConditions" : {
	"velocity" : {
	    "Dirichlet" : {
		"Air/Right"           : { "expr" : "{0,0}" },
		"Air/Walls"           : { "expr" : "{0,0}" },
		"Air/Input"           : { "expr" : "{0,(2e-3+2e-3<x)*(x<2e-3+4e-3)*7e-3*(3/2/(4e-3-2e-3))*(1-((x-(4e-3+2e-3)/2-2e-3)/((4e-3-2e-3)/2))^2)*(1-exp(-t/3))}:x:t" }
	    }
	},
	"fluid" : {
	    "outlet" : {
		"Air/Output"          : { "expr" : "0." }
	    }
	},
        "temperature" : {
            "Dirichlet" : {
                "Air/Input"           : { "expr" : "300" },
		"Pcb/Input"           : { "expr" : "300" }
            },
            "Neumann" : {
                "Air/Output"          : { "expr" : "0." },
                "Pcb/Output"          : { "expr" : "0." },
		"Air/Right"           : { "expr" : "0." },
		"Pcb/Left"            : { "expr" : "0." }
            },
	    "VolumicForces" : {
		"IC1"                 : { "expr" : "1e6*(1-exp(-t)):t" },
		"IC2"                 : { "expr" : "1e6*(1-exp(-t)):t" }
	    }
	}
    },

    "InitialConditions" : {
	"temperature" : {
	    "" : {
		"" : { "expr" : "300" }
	    }
	}
    },

    "PostProcess" : {
	"use-model-name" : 1,
	"heat-fluid" : {
	    "Exports" : {
		"fields" : [
		    "fluid.velocity","fluid.pressure",
		    "heat.temperature","fluid.pid"
		]
	    }
	},
	"fluid" : {}, "heat" : {}
    }
}
directory=toolboxes/heatfluid/opus/test1

case.dimension=2

[heat-fluid]
filename=$cfgdir/eads_normal.json
mesh.filename=$cfgdir/eads.geo
gmsh.hsize=5e-4

[heat-fluid.heat]
initial-solution.temperature=300
stabilization-gls=1
bdf.order=1

# gmsh.partition=1

[heat-fluid.fluid]
bdf.order=1
# solver=Newton
pc-type=lu

# gmsh.partition=1

[ts]
time-step=1
time-final=1500
# restart.at-last-save=true