A Generalised RBF Finite Difference Approach to Solve Nonlinear Heat Conduction Problems on Unstructured Datasets

Radial Basis Functions have traditionally been used to provide a continuous interpolation of scattered data sets. However, this interpolation also allows for the reconstruction of partial derivatives throughout the solution field, which can then be used to drive the solution of a partial differential equation. Since the interpolation takes place on a scattered dataset with no local connectivity, the solution is essentially meshless. RBF-based methods have been successfully used to solve a wide variety of PDEs in this fashion. Such full-domain RBF methods are highly flexible and can exhibit spectral convergence rates Madych & Nelson (1990). However, in their traditional implementation the fully-populated matrix systems which are produced lead to computational complexities of at least order-N2 with datasets of size N. In addition, they suffer from increasingly poor numerical conditioning as the size of the dataset grows, and also with increasingly flat interpolating functions. This is a consequence of ill-conditioning in the determination of RBF weighting coefficients (as demonstrated in Driscoll & Fornberg (2002)), and is described by Robert Schaback Schaback (1995) as the uncertainty relation; better conditioning is associated with worse accuracy, and worse conditioning is associated with improved accuracy. Many techniques have been developed to reduce the effect of the uncertainty relation in the traditional RBF formulation, such as RBF-specific preconditioners Baxter (2002); Beatson et al. (1999); Brown (2005); Ling & Kansa (2005), or adaptive selection of data centres Ling et al. (2006); Ling & Schaback (2004). However, at present the only reliable methods of controlling numerical ill-conditioning and computational cost as problem size increases are domain decomposition Hernandez Rosales & Power (2007); Wong et al. (1999); Zhang (2007); Zhou et al. (2003), or the use of locally supported basis functions Fasshauer (1999); Schaback (1997); Wendland (1995); Wu (1995). In this work the domain decomposition principle is applied, forming a large number of heavily overlapping systems that cover the solution domain. A small RBF collocation system is formed around each global data centre, with each collocation system used to approximate the governing PDE at its centrepoint, in terms of the solution value at surrounding collocation points. This leads to a sparse global linear system which may be solved using a variety A Generalised RBF Finite Difference Approach to Solve Nonlinear Heat Conduction Problems on Unstructured Datasets 13


Introduction
Radial Basis Functions have traditionally been used to provide a continuous interpolation of scattered data sets.However, this interpolation also allows for the reconstruction of partial derivatives throughout the solution field, which can then be used to drive the solution of a partial differential equation.Since the interpolation takes place on a scattered dataset with no local connectivity, the solution is essentially meshless.RBF-based methods have been successfully used to solve a wide variety of PDEs in this fashion.Such full-domain RBF methods are highly flexible and can exhibit spectral convergence rates Madych & Nelson (1990).However, in their traditional implementation the fully-populated matrix systems which are produced lead to computational complexities of at least order-N 2 with datasets of size N.In addition, they suffer from increasingly poor numerical conditioning as the size of the dataset grows, and also with increasingly flat interpolating functions.This is a consequence of ill-conditioning in the determination of RBF weighting coefficients (as demonstrated in Driscoll & Fornberg (2002)), and is described by Robert Schaback Schaback (1995) as the uncertainty relation; better conditioning is associated with worse accuracy, and worse conditioning is associated with improved accuracy.Many techniques have been developed to reduce the effect of the uncertainty relation in the traditional RBF formulation, such as RBF-specific preconditioners Baxter (2002); Beatson et al. (1999); Brown (2005); Ling & Kansa (2005), or adaptive selection of data centres Ling et al. (2006); Ling & Schaback (2004).However, at present the only reliable methods of controlling numerical ill-conditioning and computational cost as problem size increases are domain decomposition Hernandez Rosales & Power (2007); Wong et al. (1999); Zhang (2007); Zhou et al. (2003), or the use of locally supported basis functions Fasshauer (1999); Schaback (1997); Wendland (1995); Wu (1995).In this work the domain decomposition principle is applied, forming a large number of heavily overlapping systems that cover the solution domain.A small RBF collocation system is formed around each global data centre, with each collocation system used to approximate the governing PDE at its centrepoint, in terms of the solution value at surrounding collocation points.This leads to a sparse global linear system which may be solved using a variety of standard solvers.In this way, the proposed formulation emulates a finite difference method, with the RBF collocation systems replacing the polynomial interpolation functions used in traditional finite difference methods.However, unlike such polynomial functions RBF collocation is well suited to scattered data, and the method may be applied to both structured and unstructured datasets without modification.The method is applied here to solve the nonlinear heat conduction equation.In order to reduce the nonlinearity in the governing equation the Kirchhoff integral transformation is applied, and the transformed equation is solved using a Picard iterative process.The application of the Kirchhoff transform necessitates that the thermal property functions be transformed to Kirchhoff space also.If the thermal properties are a known and integrable function of temperature then the transformation may be performed analytically.Otherwise, an integration-interpolation procedure can be performed using 1D radial basis functions, as described in Stevens & Power (2010).In recent years a number of local RBF collocation techniques have been proposed, and applied a wide variety of problems (for example; Divo & Kassab (2007); Lee et al. (2003); Sarler & Vertnik (2006); Wright & Fornberg (2006)).A more comprehensive review of such methods is given in Stevens et al. (2009).Unlike most local RBF collocation methods that are used in the literature, the technique described here utilises the Hermitian RBF collocation formulation (see section 2 for more details), and allows both the PDE-boundary and PDE-governing operators to be included within in the local collocation systems.This inclusion of the governing PDE within the basis functions is shown in Stevens et al. (2009) to significantly improve the accuracy and stability of solutions obtained for linear transport problems.Additionally, the incorporation of information about the convective velocity field into the basis functions was shown to have a stabilising effect, similar to traditional upwinding methods but without the requirement to alter the stencil configuration based on the local convective field.The standard approach to the solution of linear and nonlinear heat conduction problems is the use of finite difference and finite volume methods with simple polynomial interpolants Bejan (1993); Holman (2002); Kreith & Bohn (2000).Due to the dominance of diffusion in most cases, central differencing techniques are commonly used to compute the heat fluxes.However, limiter methods (such as the unconditionally stable TVD schemes) may be used for nonlinear heat conduction problems where the effective convection term, which results from the non-zero variation of thermal conductivity with temperature, can be expected to approach the magnitude of the diffusive term (see, for example, Shen & Han (2002)).Full-domain RBF methods have also been examined for use with nonlinear heat conduction problems (see Chantasiriwan (2007)), however such methods are restricted to small dataset sizes, due to the computational cost and numerical conditioning experienced by full-domain RBF techniques on large datasets.The present work demonstrates how local RBF collocation may be used as an alternative to traditional finite difference and finite volume methods, for nonlinear heat conduction problems.The described method retains freedom from a volumetric mesh, while allowing solution over unstructured datasets.A central stencil configuration is used in each case, and the solution is stabilised via the inclusion of the governing and boundary PDEs within the local collocation systems ("implicit upwinding"), rather than by adjusting the stencil configuration based on the local solution field ("traditional upwinding").The method is validated using a transient numerical example with a known analytical solution (see section 4), and the ability of the formulation to handle strongly nonlinear problems is demonstrated in the solution of a food freezing problem (see section 5).

RBF method formulation
The Hermitian RBF collocation method operates on a domain which is covered by a series of N scattered data points, including a distribution of points over all domain boundaries.The solution is constructed using N distinct basis functions, which are composed of a partial differential operator applied to a radial basis function Ψ j which is centred on the data point j.The partial differential operator that is applied to each basis function is the PDE boundary operator for points lying on the domain boundary, and the governing PDE for points lying within the domain (see equation 3).A polynomial term is required to complete the underlying vector space.In the present work, the multiquadric radial basis function is used throughout, with m = 1.
The multiquadric RBF is a conditionally positive definite function of order m,whichrequires the addition of a polynomial term of order m − 1, together with a homogeneous constraint condition, in order to obtain an invertible interpolation matrix.For the m = 1c a s eu s e d in this work, the polynomial is simply a constant term.The 'c' term is known as a 'shape parameter', and describes the relative width of the RBF functions about their centres.The tuning of the shape parameter can have a significant effect upon the accuracy of the solution and the conditioning of the numerical system.Consider a typical linear boundary value problem where the operators L [] and B [] are linear partial differential operators on the domain Ω and on the boundary ∂Ω, describing the governing equation and boundary conditions respectively.Data points ξ j are distributed over the boundary and inside the domain, and the solution is constructed from basis functions centred around the ξ j .At data points lying on ∂Ω the boundary operator is applied to the RBF, in order to form the basis function, while the PDE governing operator is applied at those points inside the domain: The RBF formula (3) is then collocated at each of the data points ξ j applying the PDE boundary equation, B [u] = g, at points on the domain boundary, and the PDE governing equation, L [u] = f , at points within the domain.This leads to a symmetric collocation system, as represented by equation ( 4), which can be solved to obtain the RBF weighting coefficients

283
A Generalised RBF Finite Difference Approach to Solve Nonlinear Heat Conduction Problems on Unstructured Datasets

Fig. 1. Local system representation
The basic Hermitian RBF collocation method, as outlined above, is highly accurate and offers excellent convergence rates.However due to the presence of a fully-populated collocation matrix (4) solving even moderately-sized problems becomes computationally expensive.The poor numerical conditioning of this system also places limits upon the size of the dataset that can be solved.To address these issues of computational cost and numerical conditioning, the Hermitian method is applied here in "finite difference mode".In this way a large number of small and highly overlapping collocation systems are formed, leading to a numerical method that is of order-N computational complexity, while maintaining many of the flexibilities present in the global Hermitian RBF collocation method.
The Hermitian RBF method in finite difference mode takes a scattered dataset of internal ("solution") and boundary centres, and forms a local stencil around each non-boundary node (see Figure 1).Additionally, within each stencil a number of PDE centres are distributed.
At the boundary centres the boundary operator is enforced within the Hermitian collocation system, at the solution centres the unknown solution value is collocated, and at the PDE centres the PDE governing operator is enforced.In this way a series of small symmetric collocation systems are formed, which can be solved independently to interpolate the solution over the domain in a piecewise fashion.
A series of local RBF collocation systems are constructed for the interpolation coefficients α (k) ,where Here the u i are the unknown values of the solution field u, at the solution centres contained within the local system k.The shape parameter, c, in the multiquadric RBF ( 1) is kept constant within the local system, but may be allowed to vary from system to system.By defining the shape parameter in terms of the local system support domain size, R, a dimensionless shape parameter may be described (7) which is consistent over varying stencil sizes.
The value of u within the bounds of the local system may be reconstructed as follows where the "reconstruction vector" H (x) (k) is given by By inverting the matrix system (6), i.e. α (k) = A (k) −1 d (k) , it is possible to express the field variable u at any point within the stencil of local system k in terms of the data vector d (k) ,i.e.
Similarly, the PDE governing operator, L, can be applied to the reconstruction vector H (k) in order to reconstruct its value at any location within the valid domain: As before, by expressing the vector α (k) in terms of the data vector d (k) , the value of the PDE governing operator can be given by where W

(k)
L is a stencil weights vector for the differential operator L at local system k.B y reconstructing the PDE operator at the system centrepoint, a relation can be obtained linking the values of u i within the local system Applying the above reconstruction to each local system k,aseriesofN simultaneous equations are produced for u i , i = 1,...,N,w h e r eN is the global number of solution centres.In the resulting global linear system the corresponding boundary conditions of the problem have already been imposed, at the local interpolation, in those stencils containing boundary points.The resulting linear system is sparse, with the number of non-zero entries in each row equal to the number of solution centres in the corresponding collocation system.Therefore, if the local system size is kept constant, the method may be scaled efficiently to very large datasets.The method of solution for time dependent problems requires the creation of a modified PDE-operator via a finite difference approximation of the time derivative.The procedure is illustrated here using a Θ-weighted Crank-Nicolson approach, but is easily extensible to any number of finite difference time advancement schemes.

285
A Generalised RBF Finite Difference Approach to Solve Nonlinear Heat Conduction Problems on Unstructured Datasets www.intechopen.com For a general initial-boundary value problem a finite difference approximation is made to the time-derivative From this approximation, a modified PDE operator is obtained where The time stepping algorithm implies that the original initial boundary value problem reduces at each time step to the solution of a boundary value problem defined by the non-homogeneous partial differential equation ( 16), with the non-homogeneous term given in terms of the solution of the problem at the previous time step.As such the solution procedure at each time step is performed in the same way as for the steady problem, using the modified operator, L, and the corresponding non-homogeneous term, L u n−1 + S 2 (x, t n ).
After solution of the global linear system, and using the current collocation systems and updated data-vectors, a reconstruction of L [u n ] at the solution and PDE centres must be made, ready for the next time step.This is done via the creation of further reconstruction arrays for every x i at which the reconstruction is required, within local system k.
The initial time step is always performed using the Θ = 1 first-order implicit time stepping formulation.Were a value of Θ < 1tobeused,L u 0 would be required in order to calculate the non-homogeneous term L u 0 .This quantity is unknown at the initial configuration, and in order to be calculated would require an interpolation of the initial solution field, using a different interpolation system which does not rely on a non-existent previous time step.At subsequent time steps any value of Θ can be chosen without altering the interpolation system.

Solution procedure
The nonlinear heat conduction equation is solved, with an arbitrary body source: where T is the temperature ρ is the material density C p is the specific heat capacity k is the thermal conductivity For convenience, the density and specific heat capacity terms may be combined into the volumetric heat capacity; c v (T) = ρ (T) C p (T).
The Kirchhoff integral transformation is applied, in order to reduce the degree of nonlinearity in the governing equation ( 19).The Kirchhoff transformation is taken as: By applying the Kirchhoff transformation to the temperature field, equation ( 19) can be rewritten as: The transformation of the nonlinear diffusion term into Kirchhoff space reduces the equation from a strongly nonlinear to a weakly nonlinear form, by removing the multiplication of first derivatives in the diffusive term: In order to complete the transformation it is necessary to obtain a closed-form representation of ψ (T) and its inverse function T (ψ),aswellasforthetwofunctionsk (ψ) and c v (ψ).Ifk (T) and c v (T) are known analytically and are integrable, it is possible to describe these functions exactly.In other cases the functions must be computed numerically .A procedure for this is suggested in Stevens & Power (2010), which can be applied to functions of k and c v that are sampled in pointwise fashion.The surface heat flux q s (x, t) has a corresponding value in Kirchhoff space, which must be used for applied heat-flux boundary conditions: In this case the transformation to Kirchhoff space removes the nonlinearity in the boundary condition; a nonlinear heat-flux condition in temperature-space is reduced to a linear Neumann condition in Kirchhoff-space (see equation ( 23)).
The solution to the nonlinear equation ( 21) requires a process of nonlinear iterations.The solution procedure performs Picard iterations on the Kirchhoff transform variable and functions thereof, in order to iteratively approach the solution at each timestep.An approximation is taken to the Kirchhoff transform variable ψ, written henceforth as ψ.F r o m This equation can be solved directly, using the procedure outlined in section 2. Convergence is examined by comparing the maximum change in ψ over the nonlinear iteration against a user specified convergence parameter ǫ NL .Once convergence has been reached, the L ψ n quantity (see equation ( 16)) can be reconstructed ready for the next time step.This procedure is summarised in Figure 2. The procedure for steady problems is similar, except for the absence of the time-advancement loop and the L ψ n reconstruction.

Validation test case
To validate the implementation of the HRBF-FD method, the unsteady nonlinear heat conduction equation ( 19) is solved in 3D over a domain, using hypothetical values for the thermal properties and source term.Taking with α and λ as arbitrary scalar parameters, the governing equation ( 19) becomes which has an analytical solution given by: At the origin, the value of the thermal conductivity becomes zero, and as such the solution becomes singular at this location.The computational domain is considered to be a cuboid, represented by The analytical solution field ( 27) is applied at the x max , y max and z max boundaries as a Dirichlet boundary condition.The gradient of the analytical temperature field is imposed at the x min , y min and z min boundaries, i.e.
where n i represents the surface normal at the boundary.At the locations where the Neumann boundaries converge, each of the converging temperature gradient conditions is imposed simultaneously, taking advantage of the double collocation property of the RBF Hermitian method (see LaRocca & Power (2007)).The initial condition is obtained by applying the analytical solution ( 27) over the interior of the solution domain, at t = 0. Two separate domains are examined.The first domain takes x min = y min = z min = 0, and x max = y max = z max = 2.As such, the singularity at the origin is included within each of the three Neumann boundaries of the problem.However, the normal temperature gradient at each of these three Neumann surfaces is considered to be zero at the origin, where the singularity lies.This can be considered valid, since the temperature gradient is analytically zero over each of these three surfaces in the limit r → 0. The second domain takes x min = y min = z min = 1, and x max = y max = z max = 3, leading to a non singular solution throughout the domain, and a non-zero temperature gradient over each of the Neumann boundaries.The two domains described above will be henceforth referred to as the singular and the translated domains respectively.The solution domain is discretised using (11 × 11 × 11) uniformly distributed solution and boundary centres.Local systems are formed by connecting each of the solution centres to the solution or boundary centre within a single Cartesian index of itself, leading to a stencil of 27 solution or boundary centres.PDE centres are placed at every Cartesian half-index, leading to 8 PDE centres within each local stencil.A dimensionless shape parameter of value c * = 5 is used throughout.The value of the parameters α and λ are both taken to be 0.5.For the transient case, local systems are reformed after every timestep.For the steady case local systems are reformed after every nonlinear iteration.In both cases, the convergence parameter is taken to be ǫ NL = 10 −6 in the L ∞ norm.For the transient case a timestep of size ∆t = 0.01 is used, with a variety of time advancement schemes.
Figure 3 shows the variation of the L 2 solution error with time, for the singular case, using three different time advancement schemes.In each case the profile of the error appears similar.
The error rises rapidly from the initial condition, reaching its maximum value at around t = 0.6.From here the error decreases towards a minimum value at around t = 5, before rising again towards a steady value.The second-order Θ = 0.5 time advancement scheme offers the best accuracy overall, with the first-order Θ = 1 scheme providing the least accurate solution for most runtimes.The Θ = 0.65 mixed scheme offers accuracy intermediate to the other two schemes.
Figure 4 shows the equivalent error variations for the translated case.Here the accuracy is significantly improved for each of the time advancement schemes, in comparison to the singular case.This is most likely a consequence of the singularity not being present within the solution domain; with the singular case the maximum error is always found at the solution centre closest to the singularity, whereas with the translated case the maximum error location may change as the solution progresses.Once again, the Θ = 0.5 case provides the most accurate solution and the Θ = 1 scheme the least accurate.The error profile appears similar to the singular case, with the main difference being that the peak error is achieved at a much earlier runtime (around t = 0.2).In both the singular and translated cases, the solution is replicated to a high degree of accuracy throughout the time advancement procedure.
When the steady solution is obtained directly, using the steady solution procedure, the L 2 error at the solution centres is 1.49 × 10 −3 for the singular case, and 3.44 × 10 −4 for the translated case.Therefore it appears that approaching the steady solution using any of the transient solution schemes offers a higher degree of accuracy than can be achieved by using the steady solution procedure, when a consistent shape parameter value is used.This is likely a consequence of providing an accurate initial condition to the transient solver.The steady solver begins with an initial guess of T (x) = 0.

Phase change example
To demonstrate the capability of the method to handle rapid changes in thermal properties, the freezing of mashed potato is considered.The functions for heat capacity and thermal conductivity typically vary rapidly during phase-change, which leads to strong nonlinearity in the PDE governing equation.In this case, a piecewise-linear approximation is taken to the thermal properties in order to facilitate their tuning to experimental results.The thermal properties for different foodstuffs may vary significantly, however they all share common features (see Figure 5).As their temperature is reduced they go from an unfrozen "liquid" state, through a transitional state, to a fully frozen state at some temperature several degrees below zero.During the transition zone the thermal conductivity changes significantly, and a large spike is observed in the heat capacity, representing the latent heat of fusion.The rapid change in the magnitude of the heat capacity makes the accurate simulation of freezing processes challenging.Experiments performed at the University of Palermo, Dipartimento di Ricerche Energetiche ed Ambientali, provided data for the freezing of a hemispherical sample of mashed potato.
The experiment was then replicated numerically, adjusting the functions for k and c v using the piecewise linear approximations described above in order to better represent the experimental data.More detail on the experimental setup, the functional parameterisation, and the optimisation procedure are given in Stevens et al. (2011).
To model the freezing process, a 3D hemispherical dataset was created.The dataset is represented in Figure 6, and consists of an unstructured, though fairly regular, distribution of 3380 nodes in total.The base surface of the hemisphere consists of 367 nodes, and at these locations a zero heat-flux boundary condition is applied, representing contact with the insulating material beneath the sample.The upper surface of the hemisphere consists of 1164 nodes, over which a time-varying temperature profile is enforced, as obtained from the (smoothed) experimental results.Additionally, 66 nodes are present along the base edge, where the top surface meets the bottom surface.Over these nodes, both boundary conditions are enforced simultaneously, taking advantage of the double collocation property of the local Hermite collocation method.The local system size varies slightly, however the modal number of boundary and solution centres present in each local system is 14.Additionally, PDE centres are added to each local system.A tetrahedralisation is performed on each local system, using the boundary and solution centres as nodes, with PDE centres placed at the centre of each resulting tetrahedron.The modal number of PDE centres present in each local system is 24.It is important to note that the tetrahedralisation is performed only to provide suitable staggered locations for the PDE centres, and plays no part in the actual solution procedure, which is entirely meshless.Since the tetrahedralisation is local, it may be performed very cheaply.It is also possible to collocate the PDE centres with the solution and boundary centres, however previous research (see Stevens et al. (2009)) indicates that a staggered placement leads to the most accurate results in the majority of cases.The simulation is performed using a second-order Crank-Nicholson implicit time advancement scheme, and a timestep of size 50 seconds.The nonlinear convergence parameter is set to ǫ NL = 10 −5 .The shape parameter is taken as c * = 1.0; significantly lower than in the validation example of section 4. It is typical among RBF methods that cases involving irregular datasets and rapid variations in governing properties will tend to favour  gradual drop in temperature between T max = −1.5 o C and T s = −4.5 o C, occurring between t = 18000 and t = 18500.In contrast, the numerical results predict a near-instantaneous drop in temperature, from T max down to well below T s , at a slightly later time.This "sudden dropoff" behaviour was replicated across a wide range of thermal parameters, and could represent a limitation in the piecewise-linear approximation to the thermal properties.The local Hermitian method was able to produce stable results using a wide range of thermal parameters, and convergence at each timestep was typically relatively fast.The size and intensity of the spike in the function for c p (see Figure 5(b)) is the feature that has most impact upon numerical stability.By increasing the height of the spike sufficiently, it is possible to find configurations where the method is unstable at any shape parameter.This is not unexpected, as an increasingly sharp spike will represent increasingly strong nonlinearities in the governing equation ( 19), within the phase transition zone.Tests were performed using stencil configurations without PDE centres, i.e. without the "implicit upwinding" feature.However, it was not found to be possible to obtain a stable solution for spikes of intensity close to that which was required to match the experimental results.That inclusion of PDE centres provides a stabilising effect has previously been demonstrated for convection-diffusion problems Stevens et al. (2009), and the stabilising effect appears to be present here also.

Discussion
The use of local radial basis function methods in finite difference mode (HRBF-FD) appears to be a viable option for the simulation of nonlinear heat conduction processes, particularly when irregular datasets are required.Traditional polynomial-based finite difference methods are difficult to implement on irregular datasets, and RBF collocation allows a natural generalisation of the principle to irregular data.The inclusion of arbitrary boundary operators within the local collocation systems allows the flexibility to enforce a wide variety of boundary conditions, and the double-collocation property of the Hermitian RBF formulation allows multiple boundary operators to be enforced at a single location where required (such as on converging boundaries).
The inclusion of the governing PDE operator within the local collocation systems is optional, but when present introduces an "implicit upwinding" effect, which stabilises the solution and improves accuracy, at the expense of larger local systems and hence higher computational cost (discussed further in Stevens et al. (2009)).The stabilisation effect is similar to that of stencil-based upwinding, but operates on a centrally defined stencil.Therefore, the HRBF-FD method may be of benefit to problems which may otherwise require upwinding schemes, in particular with unstructured datasets, where the selection of appropriate upwinding stencils may be particularly challenging.
The application of the Kirchhoff transformation greatly simplifies the PDE governing equation and linearises heat-flux boundary conditions, at the cost of requiring thermal property functions to be transformed to Kirchhoff space.Using this Kirchhoff formulation the HRBF-FD method is able to solve a benchmark heat transfer problem to a high degree of accuracy, using both steady and transient solution procedures.Additionally, the method was able to produce stable results for a phase change model involving the freezing of food, in the presence of strongly varying thermal properties.By tuning the thermal properties it was possible to replicate the experimental data to a good degree of accuracy, potentially allowing the calibrated thermal properties to be used in further numerical simulations.

294
Convection and Conduction Heat Transfer www.intechopen.com

286
Convection and Conduction Heat Transfer   www.intechopen.com

Fig. 2 .
Fig. 2. Schematic description of the nonlinear solution procedure this, the current guess-values of k = k ψ and c v = c v ψ can be obtained.Within the nonlinear iteration these quantities are considered constant, leading to the linearised equation:

Fig. 5 .
Fig. 5. Typical variation of thermal properties with temperature (food freezing case)

292
Convection and Conduction Heat Transfer   www.intechopen.com