# Test Case for Transient Heat Conduction

## Description

This test case represents transient heat conduction in an one-dimensional
composite slab. The "slab" is represented by two regions with stationary gas1
and gas2 with different heat conductivity, but otherwise identical properties,
that are separated by an interface (single cell band with r_alpha = 0.5).

The intent is to verify functionality of homogeneous heat transfer model, with
large heat transfer coefficient, which effectively equalises the phase
temperatures T.gas1 = T.gas2, thus mimicking homogeneous model with a single
temperature field.

Without the heat transfer model, the two gas regions are decoupled and the heat
is not properly transferred between gases with different heat conductivities,
giving an incorrect solution.

The homogeneous model requires the partial elimination algorithm in the energy
equation, without which the test will fail.

The theoretical solution was obtained numerically with Mathematica and is based
on the paper:

De Monte, F. (2000).
Transient heat conduction in one-dimensional composite slab.
A natural analytic approach.
International Journal of Heat and Mass Transfer, 43(19), 3607-3619.

The solution is derived for the limit of high Biot number (representing fixed
boundary temperatures), in contrast to the paper, where a more general solution
is given.

## Technical Information

* 1D
* 20 cells (1x20x1)
* 1m x 1m x 0.1m
* Transient heat conduction
* Nu = 1000
* nEnergyCorrectors = 3

Both gases are initially at low temperature T.gas1 = T.gas2 = T0 = 300 K, Gases
are surrounded by high temperature boundaries Tinf = 310 K and begin to heat up
over time. The two gas regions are equal in length L1 = L2 = 0.5m.

Simulation runs to t = 100 s, which is equal to dimensionless
tau = alpha_k1*t/h^2 = 0.4, where alpha_k1 = 0.001 m^2/s is the thermal
diffusivity of gas1 and h = 0.05 is the cell size.

## Post-Processing

* verticalLine of T.gas1 and T.gas2
* residuals of p_rgh

are evaluated via functionObjects at run time.

Three tests are performed:

* Check that T.gas1 = T.gas2
* Validate with reference solution T.gas1.ref.gz
* Validate with analytical solution T.gas1.analytical.gz

For comparison with theoretical solution, maxAbsDev = 0.1 K.
With Tinf = 310 K, this means ~0.03% relative deviation.
This can be reduced with higher nEnergyCorrectors, at longer run time
(currently around 1.3 s).

The results may be visualised with the gnuplot script.

## Keywords

1D, transient heat conduction, composite slab, two-phase, \
stationary enclosed gas,

## Appendix - Mathematica script to produce the theoretical solution

Copy-Paste to a notebook and run with shift-enter

```Mathematica
k1 = 0.001;
k2 = 0.002;
kappa = k2/k1;
alpha1 = k1;
alpha2 = k2;
delta = alpha2/alpha1;
L1 = 0.5;
L2 = 0.5;
gamma = L2/L1;

Eq46 = Tan[beta] == -Sqrt[delta]/kappa*Tan[(gamma/Sqrt[delta])*beta];
pguess = 20;
s = NSolve[Eq46 && 0 < beta < 2*pguess, beta, Reals];
p = Length@s;
BM = Array[bm, {Length@s}];
bm[m_] /; MemberQ[Range[Length@s], m] := s[[m, 1, 2]]

Pi1m[bm_] := Tan[bm]
Pi2cm[bm_] := Tan[(gamma/Sqrt[delta])*bm]

Npm[bm_] := kappa/2*(1 + Pi1m[bm]^2) + gamma/2*(1 + Pi2cm[bm]^2)
cpm[bm_] :=
 kappa/(Npm[bm]*bm)*(Cos[bm] - Cos[(gamma/Sqrt[delta])*bm] +
    Pi1m[bm]*Sin[bm] - Pi2cm[bm]*Sin[(gamma/Sqrt[delta])*bm])
Xc1m[xsi_, bm_] := Sin[bm*xsi] + Pi1m[bm]*Cos[bm*xsi]
Xc2m[xsi_, bm_] :=
 Sin[bm*xsi/Sqrt[delta]] - Pi2cm[bm]*Cos[bm*xsi/Sqrt[delta]]

Theta1[xsi_, tau_] :=
 Sum[cpm[bm[m]]*Xc1m[xsi, bm[m]]*Exp[-bm[m]^2*tau], {m, 1, p}]
Theta2[xsi_, tau_] :=
 Sqrt[delta]/kappa*
  Sum[cpm[bm[m]]*Xc2m[xsi, bm[m]]*Exp[-bm[m]^2*tau], {m, 1, p}]
Theta[xsi_, tau_] :=
 Piecewise[{{Theta1[xsi, tau], -1 <= xsi <= 0}, {Theta2[xsi, tau],
    0 <= xsi <= gamma}}]

(* Dimensionless solutions from tau 0.1 to 1.0, dtau = 0.1: *)
(*
solutions[xsi_] := Table[Theta[xsi, tau], {tau, 0.1, 1.0, 0.1}];
Plot[solutions[xsi],{xsi, -1, gamma}, PlotRange\[Rule]{Automatic,Automatic}]
thetas = Table[solutions[xsi], {xsi, -1, gamma, 0.05}];
Export["thetas.dat", thetas] *)

(*Setup for compositeSlab test case*)
(*Temperature solutions in terms of T(xsi): *)
T1inf = 310;
T2inf = 310;
T0 = 300;
Theta1Dim[xsi_, tau_] := T1inf - Theta1[xsi, tau]*(T1inf - T0)
Theta2Dim[xsi_, tau_] := T2inf - Theta2[xsi, tau]*(T2inf - T0)
ThetaDim[xsi_, tau_] :=
 Piecewise[{{Theta1Dim[xsi, tau], -1 <= xsi <= 0}, {Theta2Dim[xsi,
     tau], 0 <= xsi <= gamma}}]

(* Solution only at time tau = 0.4
(t = 100 s in the test case, dt = 1s, dx = 0.1m *)
dimSolutions[xsi_] := Table[ThetaDim[xsi, tau], {tau, 0.4, 0.4, 0.4}];

(* Definition from the paper:

xsi = x/L1

dxsi = dx/L1, dx = 0.1, L1 = 0.5, dxsi = 0.2

OpenFOAM domain goes from 0 to 1 (L1+L2, or gamma=L2/L1, same value
in this case). We want output at cell centres, to compare these values
with fields in OpenFOAM.
So limits are :
*)
xsi = x/L1;
dx = 0.05;
dxsi = dx/L1;
xsi1 = -1 + dxsi/2;
xsi2 = gamma - dxsi/2;
tgas1 = Table[{xsi*L1 + L1, dimSolutions[xsi]}, {xsi, xsi1, xsi2,dxsi}]
Export["tgas1.dat", tgas1]
```
