Complexity matters: highly-accurate numerical models of coupled radiative-conductive heat transfer in a laser flash experiment

Thermal diffusivity measurements of samples transmitting thermal radiation require adjustments to the data treatment procedures in laser flash analysis. Conventionally, an unconstrained diathermic model is used. Current results show that the alternative coupled radiative-conductive models produce substantially different results -- for instance, at high temperatures in oxide ceramics. However, care must be taken to ensure accurate implementations of each constituent computational technique. The latter are presented in this work.


Introduction
High-temperature measurements of thermal properties using different experimental techniques such as the laser flash analysis [1] and the guarded hot plate method [2] can be challenging for many reasons -including, for instance, the stability of data acquisition and detector performance [3]. When conducting tests on materials transmitting thermal radiation (hereinafter referred to as the semi-transparent materials), e.g. metal oxides [4,5], the optical properties of the sample material can significantly influence the measurement accuracy. Interestingly, even oxide nuclear fuel exhibits a degree of transparency to thermal radiation at high temperatures, potentially influencing the measurement procedure [6]. Other applications include thermal barrier coatings for the aerospace industry [7]. The non-vanishing interest in accurately measuring thermal properties of semi-transparent materials has instigated the development of mathematical methods aimed at quantifying radiative transfer and its coupling with heat conduction. Although some authors focussed on delivering quick estimates based on non-coupled heat conduction and radiative transfer [8,9,10,11], significant effort has been undertaken to address the coupled problem, primarily based on the works [12,13,14,15] with a recent development reported by Braiek et al. [16]. These models concern radiative transfer in either non-scattering or weakly-scattering media. An approximate solution to the radiative transfer equation (RTE) using the exponential kernel technique and the two-flux method has been given for the latter case. These two methods had been used extensively in the past as follows from the introduction to [17] and could have hardly yielded realistic results for a scattering phase function with strong anisotropy. Moreover, the heating term was deemed small compared to the ambient temperature, which facilitated the solution of the initial problem. This is typically never satisfied under experimental conditions. To treat the radiative part, the three-flux method has also been considered in [18] and an early attempt to use the discrete ordinates method (DOM) -first introduced by Chandrasekhar [19] was reported by da Silva et al. [20]. Currently, DOM is often applied to this kind of problems [21,22,23], although it is still not clear if this has indeed increased the reliability of laser flash analysis on semi-transparent samples. Even using the DOM formalism, some authors still defer to non-scattering transfer blaming difficulties in the estimation of some coefficients [24]. Zmywaczyk and Koniorczyk [25], Lacroix et al. [26] have used the DOM in its rather conventional form with reference to Fiveland [27,28]. On the other hand, simplified non-coupled models are still the most popular choice for experimental data treatment in the majority of cases [29]. It is thus inconclusive if solving a coupled conductive-radiative problem is advantageous compared to using less demanding methods, and whether anisotropic scattering has any measurable effect on the thermal properties determined from a laser flash experiment. This is partially due to numerical heat transfer still being a developing area with many caveats still not addressed sufficiently: particularly, for the spatial and angular discretisation in DOM [30,31]. This paper is aimed at delivering reliable and fast numerical algorithms for one-dimensional coupled conductive-radiative heat transfer with application to the laser flash analysis. The algorithm and procedures outlined in this work are part of the PULsE (Processing Unit for Laser Flash Experiments) software, which is an opensource, cross-platform Java code freely distributed under the Apache 2.0 license [32].  Early models used in laser flash measurements of semi-transparent samples considered radiation and conduction as non-coupled phenomena since this greatly simplifies the mathematical formulation of the problem. Tischler et al. [8] considered an exponential decay of radiation intensity in a solid partially transparent to the laser pulse. McMasters et al. [9] applied the optically thick approximation and introduced an additional source term in the heat equation. These models are useful to gain a crude estimate of thermal diffusivity e.g. in porous samples and semi-conductors with an intermediate band gap. Rather than considering laser penetration in solids -a complex problem associated with the diffusion of charge carriers, their re-combination and thermalisation by phonon emission [33] -it is much easier to manually restrict the laser absorption depth by applying a graphite coating. Blumm et al. [10] proposed the diathermic model specifically to deal with this case; an analytical solution was later developed by Mehling et al. [11]. A variation of this model is currently being used in software packaged with some commercial instruments.
The diathermic model is based on the following propositions: (a) A cylindrically shaped sample is completely transparent to thermal radiation; (b) The front (laser-facing) and rear (detector-facing) sides of the sample are coated by a thin grey absorber; (c) The coatings are in perfect thermal contact with the bulk material; (d) The side surface is free from any coating.
Consequently, the monochromatic laser radiation is largely absorbed at the front face of the sample (y = 0), causing immediate heating. A portion of thermal radiation causes the rear face (y = 1) to start heating precisely at the same time (ahead of thermal conduction). The remainder energy dissipates in the ambient. It is thus sufficient to consider three radiative heat fluxes. The first two correspond to heat dissipation within the furnace chamber [3]. The third flux acts to thermalise the parallel boundaries by radiative transfer only [34]: where the emissivities of both faces are assumed to be equal (ε 1 = ε 2 = ε). Let η = ε/(2 − ε), so that 0 < η ≤ 1. Since nonlinear heat losses can be neglected [Appendix A], the boundary problem is written as: where eqs. (2a) and (2d) and the corresponding notations are the same as in [3]. The standard non-dimensional variables are used, also defined in the same reference.

A finite-difference solution
Let the superscript i and the subscript j = 0, ..., N − 1 denote the time step and the coordinate index respectively. The boundary conditions [eqs. (2b) and (2c)] are expressed in finite differences as follows: The usual Taylor expansion is written down in the h-vicinity of ξ = ξ 0 and ξ = ξ N−1 , thus defining the virtual nodes ξ = ξ −1 and ξ = ξ N needed to evaluate the boundary derivatives. After some elementary algebra, an O(h 2 + ∆t 2 ) accurate scheme is readily obtained: with the heat equation also given in finite differences: where a j = c j = 1. A fully implicit scheme shown previously to work well in most cases [3] corresponds to: b j = 2 + h 2 /∆t, R j = −h 2 /∆tθ i j . Equations (4) and (5) are reduced to the following linear matrix equation: Since z N−1 = 0, the matrix A does not have the required tridiagonal form. This problem can be solved by applying the bordering method [35]. Consider the following equations equivalent to Eq. (6a): where A = A N−1,N−1 is the border minor of the matrix A, the vectors θ T = θ 0 θ 1 · · · θ N−2 and F = r 0 0 · · · 0 R N−2 . The vector A solution to Eq. (7a) is sought in the form θ i+1 = w I + θ N−1 w II , where w I and w II are the solutions of the linear matrix equations A w I = R and A w II = −u . Since A is a Jacobi matrix, the two latter equations can be solved using the standard tridiagonal matrix algorithm. The following relations hold: The sweep algorithm coefficients can then be easily calculated: An example calculation using the diathermic model with the finite-difference scheme described in this section is shown in fig. 1. The calculation used a default grid density N = 30 and a time step ∆t = t F h 2 , where t F = 0.5.

The general form of the coupled conductive-radiative heat transfer problem
The following is the equation of radiative transfer in a plane-parallel geometry with an axially symmetric radiation field for a grey participating (i.e., emitting, absorbing, and scattering) medium compliant with the Kirchhoff's law [36]: where s is the path travelled by radiation; n is the refractive index of the medium; ψ, χ and ε are respectively the linear absorption coefficient, the scattering coefficient and the emissivity -all averaged over the radiation spectrum; Φ(µ, µ ) is the phase function of scattering, such that Φ(µ, µ )dµ /2 = 1; µ = cos Θ is the cosine of the angle between the light propagation direction and the outward normal to an elementary illuminated surface.
It is convenient to express eq. (10) in terms of the optical thickness τ = ψ ds = ψ dy/ cos Θ, which then allows separating the positive (0 < µ ≤ 1) and negative (−1 ≤ µ < 0) streams. After introducing the albedo for single scattering ω 0 = χ/(ψ + χ), the RTE e.g. for I + takes the form: where the source function is defined as S(τ, µ) = (1 − ω 0 )J + 0.5ω 0 µ IΦ(µ, µ )dµ . A matching equation may be written for I − , thus the RTE may bs solved separately for streams propagating in the positive and negative hemisphere originating at either τ = 0 or at τ = τ 0 := lψ. The complexity of the problem is determined by the source function S(τ, µ), which in some cases, e.g. at ω 0 = 0, allows an analytical solution. Once a solution has been obtained, the net radiative heat flux F(τ) can be calculated using an expression for a radiative field with axial symmetry [19]: Conduction and radiation both contribute to the heat flow, which becomes In the isotropic case dF/dy = τ 0 ×dF/dτ and the dimensionless heat equation may be written as: where q = F/(n 2 σ 0 T 3 0 ) is the dimensionless radiative flux; in addition, the Planck number is introduced: N P = λ / 4σ 0 n 2 T 3 0 l . The re-normalisation of the heat flux simply leads to substituting the emission function J(τ) [eq. (11)] with which is also dimensionless. The boundary radiative fluxes q(0) and q(1) are inferred from the boundary intensities I + (0) and I − (τ 0 ) determined through the conditions of diffuse emission and reflection [e.g. [37]]: where G 0 and G τ 0 is the incident irradiation reaching the respective boundary.

4.
A closer look at the radiation problem 4.1. Useful special cases 4.1.1. Exact solution at ω 0 = 0 In the absence of scattering, the source function s(τ, µ) is simply equal to the emission function j(τ). This then simplifies the equation, which is solved in terms of the exponential integrals E n (t) [see e.g. [36]]. The latter are defined as E n (t) = 1 0 e −t/µ µ n−2 dµ, n > 0. This leads to the following expression for the radiative flux [38]: where i + (0) and i − (0) are the boundary intensities.
Consequently, the radiation fluxes at the boundaries are: Combining eqs. (16) and (18) allows to evaluate i + (0) and i − (τ 0 ) from a set of two linear equations (see e.g. [12]): The heat source term in eq. (14a) can either be calculated using a discrete approximation or exactly using the analytic expression below derived using the properties of the exponential integral (the reader is referred to [36]): A quadrature scheme needs to be used in order to calculate the integrals in eqs. (17) and (20) -this is given in Appendix B.

The two-flux approximation
For a weakly-anisotropic phase function Φ(µ, µ ), when τ 0 is not very large, the two-flux approximation originally introduced by Schuster [39], Schwartzschild and Gesell [40] has been shown to yield sufficiently accurate results [41]. In current notations, this approximation considers I + and I − averaged over the positive and negative hemispheres correspondingly. The governing equations are then [42]: where u is an integral scattering parameter of the model. The phase function can be expanded in a series of Legendre polynomials [19]. In the linear-anisotropic approximation the series is truncated after the second term, which in a axially-symmetric radiation field gives rise to Φ(µ, µ ) = 1 + gµ µ . This corresponds to u = 0.5 − 0.25g.

The general case of strong anisotropic scattering in a nonlinear grey participating medium
The true multi-modal [43] form of the scattering function Φ(µ, µ ) can be derived from the Lorenz-Mie theory (see e.g. [44]). In most practical applications, it is more convenient to use an approximation, which still captures the strongly anisotropic scattering behaviour. This is commonly done using the single-parameter Henyey-Greenstein phase function [45]. Other specialised functions have been discussed in [46,47,48,49].
The phase function of interest is thus: If the integral over Φ(µ, µ ) in eq. (10) cannot be simplified, as in case of Φ(µ, µ ) given by eq. (22), the solution to the RTE becomes quite involved. Some effort in solving the RTE and, indeed, the coupled problem has been undertaken by many authors [26,50,25]. Generally, the discrete ordinates method (DOM) is used for this purpose. Henceforth, the paper is focussed on the numerical implementation of DOM.
Recall the general form of the source function: The idea behind DOM is to evaluate the integral on the right-hand side using a quadrature rule. A discrete set of nodes is introduced: µ m , m = 0, ..., M − 1, with an equal number of negative and positive nodes; each node is assigned a certain weight w m . The discrete form of eq. (23) is: The discrete RTE [eq. (12)] is given by: with the boundary conditions of diffuse emission and reflection [eq. (16)]: The net radiative flux [eq. (17)]:

The solution to the heat problem
It is convenient to first select an appropriate numerical scheme for solving the heat problem outlined in section 3 before launching a full-scale analysis of the radiative transfer problem (section 4.2). For this reason, the analytical solution obtained in section 4.1.1 is used to calculate the heat fluxes q and their derivatives dq/dτ (for details of the calculation method the reader is referred to Appendix B). The current section includes a comparison of various finite-difference scheme for solving the heat problem.

Explicit scheme
The problem can be solved using an explicit finite difference scheme with an embedded fixed-point iteration algorithm. The finite differences are written on a rect- The discretised heat equation serves to calculate the reduced temperature θ j at j = N − 2, N − 3, ..., 0. Let θ j be the temperature value at the previous timestep. The explicit scheme for the heat equation is then: where the second-order differential operator is defined as The equations arising from the boundary conditions are solved iteratively: where Ξ is the pulse function, k is the iteration number.

Implicit schemes
Both the fully-implicit and semi-implicit scheme follow the same solution logic [3]. The discrete heat equation is written as follows [51]: where 0 < σ ≤ 1 is the weight of the scheme (σ = 1.0 for the fully-implicit scheme); φ j is some finite-difference representation of the term −dq/dt.
After some elementary algebra, the following expressions are derived, completing the difference scheme: The scheme is solved iteratively until converged values of k+1 θ 0 and k+1 θ N−1 are obtained (usually a few iterations are required). It is at least O(h 4 + ∆t 2 ) accurate if [51]: where the superscript indicates averaging over two consequent time steps andq = dq/dτ.

Verification and benchmarking
The general method (section 3) and the finite-difference schemes are verified against the reference solutions reported in [13] for τ 0 = 0.1 and τ 0 = 100 at N P = 0.8612. The calculated time-temperature profiles are shown in Figure 2 where a good agreement between the linearised analytical case and the exact numerical solution is observed at δ T m = 0.4. Implicit schemes produce more accurate results compared to the explicit scheme; particularly, the fourth-order accurate semi-implicit scheme described in section 5.2 performs well even for coarse grids with N = 10. This is especially important in the light of a high demand on computational resources expected when solving the inverse coupled radiative-conductive problem. It is also evident that only a numerical scheme is applicable to solving the heat problem eq. (14) at anywhere near realistic δ T m /T 0 values [note the difference between fig. 2 (a) and (c)].

Spatial discretisation and integration
Commonly, the RTE [eq. (25)] is integrated using a diamond-differencing scheme (see e.g. [52]), also known as the central-difference scheme, which for a one-dimensional  problem is exactly the same as the implicit trapezoidal rule -a second-order accurate and A-stable method. All alternative conventional methods are based on the finitevolume methodology [37,53] and include: the first-order step scheme, the secondorder exponential, hybrid and CLAM schemes. Advances in spatial discretisation schemes for RTE, mainly based on NVD and TVD for multi-dimensional radiative transfer, have been reviewed in [31] -however, with no significant progress reported for high-order spatial differencing schemes. More recently, Maginot et al. [54] have used a stiffly-accurate single diagonally implicit Runge-Kutta (SDIRK) method reported originally by Alexander [55]. Notable implementations of SDIRK are included in [56,57]. Despite their advantages, SDIRK methods only allow a stage-order of one [58]. Higher stage-order is useful since this strongly improves accuracy when applied to stiff problems and increases the error-estimate quality [59]. Stage-order two may be achieved with the first-stage explicit SDIRK (ESDIRK). Alternative to ESDIRK is the Rosenbrock method [58], which might be more efficient for some problems [60].

Explicit Runge-Kutta with an adaptive uniform grid
A given ODE can have varying stiffness depending on the parameter values. For the RTE, stiffness is mainly determined by τ 0 . At τ 0 < 1 the problem can be effectively treated as non-stiff. When stiffness is not an issue, explicit embedded Runge-Kutta schemes can be used. If high accuracy is desired, a good fourth-order scheme such as the Dormand-Prince (DP54) [61] scheme with a fifth-order error control and an extended region of absolute stability can be used. Practice shows that for the current use of the RTE, error tolerance can be high, thus a lower order embedded method might be sufficient. In this case, a third-order Bogacki-Shampine (BS32) [62] scheme with second-order error control and good stability can be used. Both schemes are FSAL (first same as last), which saves computational time, and their implementation follows the same pattern described below.
Firstly, let h l denote the signed grid step, which is positive when approaching the right boundary and negative otherwise. The following notations are used: the intensities at each stage n = 1, ..., s are denoted as i (n) with the corresponding coordinate t where i (n) m = i ml +h l ∑ n−1 n =1 a nn f (n ) m are the outward intensities at the node m and stage n. Depending on whether the RTE is solved left-to-right or right-to-left, the angular index m for the outward intensities will run through the indices of either positive or negative nodes (cosines). For the sum over outward intensities, the latter are expressed in the same way using the solution i (n−1) m at the stage n − 1. Inward intensities i m l+c n are not known a priori, which is why the RTE is solved iteratively; this will be described in more detail later in the text. For now these intensities are assumed to be known.
Once the derivative f (n) m becomes known, it is then used to calculate the next stage approximation i (n) m , and so on. This process repeats for all n = 1, ..., s. As soon as all derivatives have been calculated, the intensities i ml+1 may be evaluated using the respective expression. Error control is achieved by evaluating the vector est. the components of which are given by: where m runs through the indices of outward intensities. Absolute and relative tolerances are introduced according to Hairer et al. [63] so that the error threshold is defined via: Thus, max m (est m ) is compared at each subsequent integration step l ±1 against e l±1 -if the former is greater than the latter, integration stops immediately, triggering a grid re-construction with a different segmentation: where [u] indicates the value at current iteration.
As mentioned above, to solve the RTE, one must calculate the intensities corresponding to both the negative and positive µ m . However, when using the method above to solve either of the Cauchy problems, only half of the intensities is readily calculated while the other half is assumed to be known. To solve the RTE for all µ m , an iterative solution is required. Here two techniques are considered [64]: the fixed-point iterations and the successive over-relaxation. In both cases, the intensities at iteration [u + 1] are expressed as: where the relaxation parameter ω R = 1 for fixed-point iterations and 1 < ω R < 2 in the successive over-relaxation technique. The second term on the right-hand side is the solution of the ODEs times the relaxation parameter. For instance, at ω R = 1.7 and for pure isotropic scattering at ω 0 = 1, convergence is reached two times faster than for fixed-point iterations.
The stopping criterion for the iterative procedure regards the relative change to the boundary fluxes q 0 and q N at the left and right boundaries correspondingly: where e it is a relative error tolerance (typically, e it 10 −4 ).

TR-BDF2 with an adaptive stretching grid
For moderately-and highly-stiff problems, e.g. at τ 0 > 10, the use of a uniform grid requires a very small step size h l to make the scheme stable, thus greatly increasing the computational cost of an explicit method. Hence, an adaptive step-size control should be used instead, which achieves true flexibility in the A/L-stable, stiffly-accurate methods, for instance, the TR-BDF2 scheme [65]. The latter can be regarded as a major improvement over the original diamond-differencing scheme for plane-parallel radiative transfer problems, since it includes the same trapezoidal rule (diamond-differencing) at the second stage and uses second-order backward-differencing at the third stage, resulting in stiff accuracy. Furthermore, it provides an asymptotically correct error estimate and allows dense output. TR-BDF2 can be regarded as an ESDIRK scheme [59].
The explicit first stage is calculated in the same way as in section 6.1, noting that TR-BDF2 is also FSAL. The second and third stages are implicit by definition. However, because the ODEs in the DOM are linear, the corresponding intensities can easily be found explicitly from the solution of the following linear set. For instance, the left-to-right sweep at the second stage: where γ = 2 − √ 2 and d = γ/2 [65]. Clearly, this reduces to a linear matrix equation Ai ml+γ , which is solved by matrix inversion. Due to the A matrix usually being low-dimensional (the dimension is equal to a half of the total number of quadrature points), a fast matrix inversion routine has been implemented for the typical quadrature sets. For higher-order quadratures, a matrix inversion tool based on either QR, LU or Cholesky decomposition of the Apache Commons Mathematics Library is used. Since the method is ESDIRK, the final third stage uses the same matrix inverse A −1 . The linear set for the third (and final) stage is (left-to-right sweep): where w = √ 2/4 [65] (this should not be confused with the quadrature weights w m ). The correct error estimate [65] valid for both stiff and non-stiff problems is then simply: Est = A −1 est, where est is given by eq. (34) and [65]: The same general scheme for error control [eq. (35)] is used.
To take advantage of the stability properties of TR-BDF2, an adaptive grid is constructed using stretching functions [66,67]. Since rapid variation of intensities is mainly expected when approaching the boundaries, it is sufficient to maintain a small step in their vicinity. The stiff solver can then use an arbitrary large step in the remainder domain. For this purpose, the grid step h l is defined via a hyperbolic tangent function: where N

[u]
g is the number of segments in a uniform grid and S g is the stretching factor. Figure 3 shows an example grid generated using the above algorithm. When the error becomes higher than the threshold given by eq. (35), the grid is reconstructed by increasing the number of grid points in the same manner as described in section 6.1. The first iteration always starts from a uniform grid with a default of N [u] g = 8 segments. The parameter S g normally does not change during the reconstruction. Finally, the same iterative procedure described in section 6.1 is adopted to obtain convergence.

Interpolation
In each case, knowledge of the dimensionless temperature θ is required at intermediate integration steps t (n) used then to calculate the reduced radiance j(t (n) ). Since the temperature is defined discretely on a different external grid of the heat equation, an interpolation procedure is required to calculate the temperature θ l at the integrator nodes. In this case, the dimensionless temperature is interpolated using natural cubic splines implemented in the Apache Commons Mathematics Library.
Both the explicit [eq. (33)] and implicit [eq. (38)] methods contain summation over the unknown inward intensities i ml+c n . Since all intensities are calculated at the internal grid points l and because l + c n is not a grid point, an interpolation procedure is required here as well to calculate i values obtained at the previous iteration. Additionally for the implicit method, the outward intensities at the intermediate points t l + γh l are not known either, and hence the same procedure needs to be used for their calculation. Because in Runge-Kutta methods both the intensities and their derivatives are calculated, a cheap and convenient method for this interpolation is the globally C 1 Hermite interpolation described in detail in [68]. The Hermite interpolant satisfying the function and derivative values at end points of the segment t ∈ [a, b] is: This allows effective interpolation of both inward and outward intensities at any intermediate point 0 < t < τ 0 .

Angular discretisation
The quadrature choice is central to the DOM as it defines both the overall accuracy of the method and the stability requirements for the spatial integration technique. Chandrasekhar [19] originally considered the Gauss-Legendre and Lobatto (Radau) quadratures for angular discretisation. In modern calculations, the level-symmetric quadratures by Lathrop and Carlson [69] are often used [37]. These and other similar quadratures have been reviewed in [70,71,72,73]. More recently, an extensive review [74] of different quadratures has shown that for problems generating a continuous intensity field, the Gauss-Chebyshev quadrature LC11 derived by Lebedev [75] offers the highest precision. Since in many cases, particularly for the one-dimensional radiative transfer with diffuse emission and reflection conditions, the intensities are discontinuous at µ = 0 (see e.g. [76]), standard quadratures which do not specifically treat the discontinuity would give inaccurate results. The level-symmetric quadratures were designed to cover both the non-continuous and discontinuous case and are applicable to a wide range of problems. However, high-order quadratures (such as S 10 , S 12 etc.) yield negative weights. Although quadratures such as S 8 give sufficiently accurate results in many cases, an alternative should be considered for higher-order calculations. A composite Gaussian quadrature has been considered for Fresnel boundary conditions in [77] where the angular interval was divided in three segments. A similar procedure can be performed for the diffuse emission and reflection boundaries.
Consider the M cosine nodes and weights of a Gauss-Legendre quadrature on [0, 1]: µ m and w m . The goal is to construct a composite quadrature that will work despite the intensities being discontinuous at µ = 0. The 2M cosine nodes of this composite quadrature are then: with the same weights w m = w m . By construction, the composite Gaussian quadrature given by eq. (43) is applicable to discontinuous functions at µ = 0. An example G 16 M ordinate set proposed in this work is given in table 1 (note this quadrature is symmetric).

Verification and benchmarking
To verify the solvers and the discrete ordinate sets, two model cases were considered: (a) a non-scattering grey medium with diffusely emitting and reflecting walls (ε = 0.85,  The deviation is decreased even more when a low error tolerance is selected (rtol = 10 −4 , atol = 10 −5 , e it = 10 −6 ). For comparison with the two-flux model, an artificial quadrature containing two equal-weight symmetric points is examined. An exact match between the approximate analytical model and the discrete ordinates method is shown in fig. 6, thus confirming the reliability of the numeric procedure. Additionally, the performance of different schemes and quadratures was tested for a grey medium with a strong anisotropic scattering (ε = 0.85, ω 0 = 0.4, g = 0.8). The results of different computational methods for the net fluxes shown in fig. 7 show good mutual agreement both in the stiff and non-stiff cases.
Finally, the relative performance of different schemes was assessed in table 2. Here the TR-BDF2 scheme in the high-tolerance mode using the G 8 M quadrature was used as reference, corresponding to the respective 1.00 table entry. Increasing problem stiffness in the high-tolerance mode only marginally increases the computational cost for TR-BDF2. Other schemes do not perform so well in terms of performance, particularly the DP5 at τ 0 = 100 is 50 times slower than the reference. BS23 performs better but still fails to deliver a reasonable computation time for stiff problems. For the G 16 M quadrature there was no fast matrix inversion implemented and hence the TR-BDF2 algorithm relied on a generic decomposition algorithm for the latter. This justifies the considerably more expensive calculations. Problems requiring only a small ordinates set (e.g. S 4 ) show a ≈ 1.5 increase in performance compared to the reference. Same performance for TR-BDF2 and DP5 is achieved at low error tolerance levels for a non-stiff (τ 0 = 0.1) problem, whereas BS23 requires a finer step size, which almost triples the overall cost. The numbers change dramatically even for moderately-stiff problems (τ 0 = 10.0), with the DP5 outperforming the BS23 scheme -as expected, since DP5 is a fourth-order method. On the other hand, both require more resources to achieve the same error tolerance compared to the TR-BDF2 due to the adaptive grid employed for the latter. For the G 16 M quadrature there is an expected drop in performance -and vice versa for the S 4 ordinate set.
With these results in mind, the default settings for calculation are chosen as TR-BDF2 and a G 8 M ordinate set in the high-tolerance mode.

Cross-verification
The goal is to verify the complete solution to the conductive-radiative problem described in sections 5 and 6. Synthetic model parameters used in the tests are listed   in table 3. These correspond to a case of non-scattering grey medium; the latter is especially helpful since it allows an exact solution to the RTE (section 4.1.1), examples of which have previously been shown in fig. 2. The resulting time-temperature profiles generated by solving the boundary problem [eq. (14)] with the radiative fluxes calculated using the discrete ordinates method were compared to the same profiles calculated using the analytical solution to the RTE. No deviation between the two calculation methods is observed (Figure 8), thus indicating a correct implementation of all solvers.

Experimental validation
Experiments conducted with the use of a laser flash analyser (LFA) produce raw data in the form of time-temperature profiles with varying level of noise [3]. Experimental validation requires solving the inverse problem of heat transfer, which boils down to finding a set of parameters (e.g. table 3) corresponding to an optimal solution of the heat problem. A solution is deemed optimal if the objective function (such as the sum of squared residuals) reaches a global minimum in the parameter space. Fortunately, the corresponding optimisation procedure has already been previously implemented and extensively tested in [3]. Nevertheless, some modifications to the procedure are required both for the diathermic model (section 2) and the coupled conductiveradiative problem (section 3). Firstly, the original linear-interpolation procedure has been replaced by spline interpolation. Secondly, the basic procedure in [3] involved only unconstrained optimisation. In case of an ill-posed problem or a tendency of the computational method to fail outside a certain region in the parameter space, the unconstrained optimisation procedure will not behave well. Figure 9 shows two almost identical time-temperature profiles obtained with two very different parameter sets. This is a classical example of an ill-posed problem [78]. To eliminate non-physical solutions, the parameter space should be bounded. The corresponding linear constraints are listed in table 4. The complete solution of the optimisation problem with linear constraints based on the active-set method has been discussed in [79] and the general method of solving ill-posed problems is known as the Tikhonov regularisation. A very simple alternative is considered in this work mainly for demonstration purposes. A one-to-one mapping Y i ∈ R → X i ∈ [a, b] is introduced for each parameter y i in table 4 using hyperbolic functions. This ensures that at each time the parameter x i only takes 'reasonable' values. The optimisation procedure is then effectively the same, except that the search vector is formed of X i rather than Y i . It should also be noted that imposing these constrains is only possible if the thermal properties of the sample (specific heat and density) are known in each experiment -otherwise there is no way of telling whether the parameter value is sensible or not. As a direct consequence, this means that even the diathermic model, which does not require neither the specific heat nor the density values for calculation, will not guarantee physically reasonable results if the thermal properties are unknown and an unconstrained optimisation is used instead.
Finally, a set of experimental data acquired for a synthetic alumina sample (l = 1.181 mm) measured in a laser flash apparatus at high temperatures has been provided for validating the computational procedure. Measurements were conducted using the Kvant instrument at the Moscow Engineering Physics Institute, previously briefly de- scribed in [3]. Specific heat and thermal expansion data have been taken from [80,81]. Density at room temperature was measured using the hydrostatic method. Example time-temperature profiles are shown in fig. 10 along with the solutions to the inverse problem using three different models. A sharp temperature peak at the start of experiment is especially pronounced at the highest ambient temperatures. Only the complete calculation with the Henyey-Greenstein phase function is capable of reproducing this behaviour, although the model deviates from the experiment slightly at the start. Possibly this is due to some residual coating on the side surface of the sample which may have created an easy path for thermal diffusion. Another point to be aware of is the fact that the sample holder used in these experiments covered a significant area of the sample. The holder effectively consisted of two washers pressed against both sides of the sample while typically a three-point contact scheme is used in modern instruments. Thus, laser radiation was non-uniformly absorbed at the front surface, covering approximately 75 %. At each test temperature, thermal diffusivity ( fig. 11) was averaged over three measurements. Results show that a complete calculation produces systematically different values compared to the diathermic model with a maximum deviation of over 10%. The high error margins are due to the optimisation procedure finding different minima depending on the starting conditions. The tendency of the optimiser to slip into a local minimum is due to the objective function being acute and multi-modal, which commonly occurs in multi-variate optimisation; moreover, even though the set of parameters can be sufficiently different, the minima are not. This highlights the necessity   Diathermic Non-scattering Full Figure 11: Thermal diffusivity of synthetic alumina determined from a single set of experimental data using three different models. Note the uncertainties associated with the different starting parameters for the full model.
of introducing additional constraints -relying on e.g. the optical properties.

Conclusions
The numerical method described in this paper combines: (a) a stiffness-aware solver, its error control scheme and an adaptive stretching grid -specifically tailored to solving the initial value problems arising from the discretised radiative transfer equation; (b) a composite Gaussian quadrature designed to treat discontinuous intensities typical to the one-dimensional radiative transfer and (c) a fourth-order semi-implicit finite-difference scheme for numerically solving the heat problem. This combination is applied to enhance the data analysis in laser flash experiments where the material under study scatters thermal radiation anisotropically, such as when conducting measurements on transparent alumina at high temperatures. The calculation procedure reproduces the initial rapid variation of temperature typical to the strongly-scattering medium while still observing physically-reasonable values of secondary model parameters (i.e., of the optical thickness, Planck number, emissivity, scattering albedo and of the anisotropic factor). The estimate quality is benchmarked against a standard diathermic model, where the maximum deviation is observed at high temperatures and pronounced scattering anisotropy. The optimisation procedure has been modified to implement constrained search using a one-to-one mapping of the search variables. This allowed imposing realistic parameter constraints. A further refinement of the search procedure is recommended to correctly address the ill-posed problems often occurring in multi-variate optimisation. The algorithms have been successfully implemented in the PULsE software, with the latest version being immediately available for use.
where Q is the heat absorbed by the thin surface layer and ε is the sample's flat surface emissivity. These equations are then transformed to the dimensionless form: where δ T m = 4Q(πd 2 C p ρl) −1 is the maximum heating of the rear surface in the absence of heat sinks and Bi := 4σ 0 εT 3 0 l/λ is the Biot number, and θ = (T − T 0 )/δ T m . It can be easily seen that if θ δ T m /T 0 is small, the heat loss term becomes simply Bi· θ y , which corresponds to the classical case. When δ T m /T 0 1, using only the first term of the Taylor expansion might not be appropriate; especially at the front surface (y = 0, see Fig. A.12), since θ y=0 θ y=1 at Fo = 0 − 0.15. However, the overall magnitude of the heat sink term is proportional to T 0 /4δ T m . Hence, the significance of this term may be low when the expression in the brackets may be nonlinear.
The finite-difference calculations proceed as follows. The domain is divided into a uniform grid by introducing the coordinate step size h = 1/(N − 1), where N is the number of individual coordinate points on the grid, and the discrete time step τ = τ F h 2 , τ F ∈ R. The grid is used to discretise θ (y, Fo), which becomes θ (ξ j , Fo m ) = θ m j , j = 0, ..., N − 1, m = 0, ..., m 0 , called the grid function. Let Lφ (ξ α ) = (φ α+1 − φ α−1 ) /2h. Then, the finite-difference analog of Eqs. A.2 is: where the time index is implicit. Consider using a Taylor expansion on the grid at j = 0 and j = N − 1 and introducing virtual nodes j = −1 and j = N, thus transforming eq. (A.3) using contraction mapping: φ = ζ (φ ). For a fully-implicit scheme the first coefficients α 1 and β 1 from the tridiagonal matrix equation θ j = α j+1 + θ j+1 β j+1 and the solution at the j = N boundary are calculated at each iteration k + 1 until the scheme converges to a given precision (usually within a few iterations): The solution is shown in fig. A.13 where the heating curves have been normalized. Curves are plotted at different values of ι := δ T max /T 0 , all else being equal. With increasing the r factor, the normalized maximum shifts towards shorter times while the temperature decreases due to heat losses (in this case, Bi = 1.0) becomes more pronounced. For δ T max /T 0 < 5 × 10 −2 (in most practical cases), this effect is so small that the nonlinear behaviour of the heat losses in eq. (A.2) may be completely neglected. Therefore, some care must be taken only when conducting measurements at cryogenic temperatures and at a high laser power applied to poor thermal conductors. Otherwise, keeping nonlinear terms in the boundary conditions is redundant and a simpler (linearised) model of the heat problem may be used instead.

Appendix B. Numerical evaluation of some integrals
The integrand function E 1 (t) is discontinuous at t = 0, which complicates the evaluation of radiative flux derivatives dq/dτ using the standard Newton-Cotes formulae. The latter require significant computational resources, which is inappropriate when the flux derivatives need to be calculated frequently.
The general problem consists in evaluating integrals of the form: The exponential integrals E n (t) are pre-calculated using the midpoint rule with a very large number of integration points by filling a look-up table of typically N tab = 10, 000 − 20, 000 entries, depending on the cutoff value (t exp c = 9.2 − 21.0), which ensures a precision of at least 10 −5 . This table is filled only once at the program start and used later in future calls to the solver. An acceptable accuracy when using a Newton-Cotes formula (e.g. the Simpson's rule) can be achieved at n q = 256 [see table B.5] for integrals of order n ≥ 2 when the integrand is well-defined at zero. Since the exponential integrals rapidly decrease with τ and the emission function j(t) is bounded, the integrand becomes very small where the exponential integrals are near-zero. The integration bounds are calculated as [max{a, (t c − α)/β }, b] at β < 0 and [a, min{b, (t c − α)/β }] at β > 0. This ensures that for large τ 0 , the integration excludes terms smaller in amplitude than a certain threshold defined by the cutoff t c . Additionally, since f (t) is discretised differently to what is used in the quadrature scheme, a natural cubic spline interpolation implemented in the The Apache Commons Mathematics Library is introduced to calculate the function values.
A more effective quadrature has been introduced by Chandrasekhar [19]. It is first noticed that eq. (B.1) may be written as: .
The moments M l are defined as: α+β a x l E n (x)dx.

(B.3)
These can be integrated by parts if the recurrent expression for E n (x) is utilised [19]. After the moments have been calculated, the next step is to calculate the x j ( j = 1, .., m) roots of the monic polynomial x m + ∑ m−1 l=0 c l x l where the coefficients c l form the solution of a linear set:   Table B.6: Comparison of end precision ∆ and computational effort T 10,000 (measured for 10, 000 consecutive calls to the respective integration method) for different quadrature formulae for calculating the integral I 1 = τ 0 τ j[θ (t)]E 1 (α + βt)dt at τ 0 = 3.0, β = 1, τ = −α = 0.5 using the same test temperature profile as in Table B