3.2. DISCRETIZATION OF THE LAYER DOMAIN

The heat transfer phenomena is a phenomenon that can be influenced by many different factors. This program does not take into account all the possible situations. Several assumptions have been considered in order to design a simple and practical model.

In the following section, we calculate the temperature profile in the layers in order to demonstrate that, for simple geometries, the confinement mass which is taken for the heat transfer expressions can be treated as an ‘infinite’ surface of thickness ‘d’ (i.e. infinite in two directions -y, z directions - and of layer thickness ‘d’ in the third -x direction). Mathematically expressed, we will establish that inside the bulk the temperature gradients in the x-direction are an order of magnitude higher than the gradients in the y,z-directions at any time and in any layer:

		(3.13)
		(3.14)
The above facts are evident because for a step dy or dz we have:
		(3.15)
		(3.16)
for any layer ‘U’.

Using the generalized heat balance eq. 3.10 over one layer element in the confinement wall, we can relate the heat transfer in each layer. The scheme of the grid-point distribution applied for calculating the temperature profile in each layer is presented in Fig. 3.2.


Fig.3.2: Scheme of the multilayer confinement. The simulation of the whole multilayer confinement reduces to the analysis of a single layer and can then be extended from layer to layer. The grid-point distribution is chosen with variable step lengths in the heat transfer direction as well as in the time direction.

The function of the heat balance (eq. 3.10) can be singular at the interface of the different layers and at the beginning of the heat transfer process (times around 0). Therefore the grid-point distribution must be chosen with variable step lengths (see Fig. 3.2). The generation of adaptive meshes allow the achievement of a desired resolution in localized regions and decreases by orders of magnitude the calculation time. Grid points are added in regions of high gradients to generate a denser mesh in that region and substracted from regions where the solution is decaying or flattening out.

Fig. 3.2 illustrates the discretization of the governing partial differential eq. 3.10 describing the heat transfer processes occurring inside the layer. The layers can be divided into a series of N mesh planes each having IxJ elements. Applying the FEA method, the position of the mesh planes is moved along the time-axis allowing the calculation of the temperature profile at each location for every x, t grid points of each layer.

Discretization: Let us choose ‘n’ for satisfying following conditions:

		(3.19)
		(3.20)
where e.g. describes with arbitrary units the length of a desired thickness inside one layer. After computing a series expansion of the above equations with respect to the variable ‘k’, we obtain:
		(3.21)
		(3.22)
Taking the natural logarithm, we can solve the inequality 3.22 for ‘n’ and round the obtained value to the nearest integers towards minus infinity:
		(3.23)

The temperature profile may be expressed by substituting finite-difference approximations. We can compute the Taylor series expansion of the second derivates, with respect to the variable x, up to the order three. The series data structure represents an expression as a truncated series in one indeterminate node, expanded about the particular point (i ,j) (see Fig. 3.2). The detailed presentation of the Finite Element Analysis rests upon the scope of this help.



Fig.3.3: Slow cook-off experiments of the rocket motor (A) and simulation (B). The predicted temperature of ignition was 124°C. It is in good agreement with the slow cook-off experiments (126°C).

For simplified geometries, the heat balance can be expressed by the following equations:

		(3.17)
or for cylindrical coordinates
		(3.18)

where J is a geometry factor which is dependent on the type of recipient: J=0 for the infinite plate, J=1 for the infinite cylinder and J=2 for the sphere. The above equation has been extended by the consideration of heat produced by the decomposition reaction Qr which rate is derived from the kinetics of the reaction. The heat balance equation can be now solved from r = 0 (centre of cylinder) to r = R (surface of the cylinder) with AKTS-Thermal Safety Software. The temperature profiles have to be considered for all layers. In each layer the initial temperatures at t = 0 have to be introduced. At the centre and if a layer is perfectly isolated on its left or right side, the boundary condition (see I, Fig. 3.2) is derived from the symmetrical properties of the temperature profile at the wall surface. The other boundary conditions (II, Fig. 3.2) are derived from comparison of the heat transfers at the interface between the different layers.

Consideration of the decomposition kinetics and application of the boundary conditions enable the calculation of the heat transfer, the temperature distribution and the reaction progress in a larger amount of material as encountered in the cook-off experiments. The slow cook-off simulation and experiments of the rocket motor are presented in figure 3.3. The ignition temperature of the slow cook-off was 126°C (Fig. 3.3A). The simulation of the cook-off behaviour and the predicted ignition temperature of 124°C (Fig. 3.3B) are in good agreement with the experiments. In the simulation the following parameters were used: TRocket initial = 40°C during 4 hours followed by a heating rate of 3.3°C/h; rocket motor diameter = 125 mm, thickness of the boost propellant layer = 31.5 mm, thickness of the sustain propellant layer = 31 mm; thermal diffusivity of boost and sustain propellant 'l/(r*cp)' = 0.02 cm²/s. Knowledge of the decomposition kinetics, thermal diffusivity and specific heat combined with FEA can be used to determine the ignition temperature very precisely. More generally, applications of FEA and accurate kinetic description enable the determination of the effect of scale, geometry, heat transfer (insulation), thermal conductivity and ambient temperature on the heat accumulation conditions. In fact, the assumption that it is safe to handle an energetic material at any temperature below the first appearance of an exothermic signal on the DSC curve can be false and even dangerous. The highest safe temperature for handling any energetic material depends on several factors such as its size, shape, and prior thermal history. Therefore, safe storage or transport conditions with tailored safety margins can be defined using numerical simulation.