Including Solar Radiation in Sensitivity and Uncertainty Analysis for Conjugate Heat Transfer Problems

With the advancement of computer-aided engineering (CAE) and the development of coupling tools that enable the exchange of data from various codes specialized in their respective fields, the possibility of acquiring more accurate results rises, because some boundary conditions can be calculated. Rauch et al. (2008) stated various other reasons for making use of co-simulation. Hand in hand with performing multi-physics coupling the required knowledge of the underlying theories, model limitations, parameters, and numerical constraints increases. This can lead to involving specialists of various fields but it can also entice someone to make improper choices for parameters of models an engineer is not sufficiently acquainted with. Either way, the number of possible parameter settings usually rises with the number of models employed. To serve as an example, Rauch et al. (2008a) were involved with an underbody simulation of a car. The exhaust system transported the exhaust gas out of the engine and heat was released along the piping system. A commercial thermal radiation solver was used to calculate wall temperatures including conduction and convection. Convective heat transfer coefficients and fluid temperatures at the wall were imported from a commercial computational fluid dynamics (CFD) solver which, in turn, received the wall temperatures. Other parameters in the radiation solver to be set were conductivity, emissivity, and settings for the view factor calculation.


Introduction
With the advancement of computer-aided engineering (CAE) and the development of coupling tools that enable the exchange of data from various codes specialized in their respective fields, the possibility of acquiring more accurate results rises, because some boundary conditions can be calculated.Rauch et al. (2008) stated various other reasons for making use of co-simulation.Hand in hand with performing multi-physics coupling the required knowledge of the underlying theories, model limitations, parameters, and numerical constraints increases.This can lead to involving specialists of various fields but it can also entice someone to make improper choices for parameters of models an engineer is not sufficiently acquainted with.Either way, the number of possible parameter settings usually rises with the number of models employed.To serve as an example, Rauch et al. (2008a) were involved with an underbody simulation of a car.The exhaust system transported the exhaust gas out of the engine and heat was released along the piping system.A commercial thermal radiation solver was used to calculate wall temperatures including conduction and convection.Convective heat transfer coefficients and fluid temperatures at the wall were imported from a commercial computational fluid dynamics (CFD) solver which, in turn, received the wall temperatures.Other parameters in the radiation solver to be set were conductivity, emissivity, and settings for the view factor calculation.
When confronted with deviations from measurements, the simulation engineer is left to choose those parameter changes that will improve the results based on a sound physical footing.This involves either expert knowledge, or a trial and error approach.Sensitivity analysis provides a systematic way for determining the most influential parameters.Global sensitivity analysis uses a statistical approach in understanding the total parameter space of the system.Saltelli et al. (2008) provide an overview on global methods which are suited for optimization.Local methods provide local information for the response of the system under investigation (Saltelli et al., 2008a).Roy and Oberkampf (2011) give an overview of sources of uncertainties that influence the outcome of a simulation.They identified uncertainties that stem from model inputs, such as model parameters, geometry, boundary and initial conditions, and external influences such as forces.They do occur in the manufacturing process, in natural resources, and might change during the life-cycle.Another source of uncertainties is of a numerical nature.To this category belong discretizations in space and time, computational round-off errors, iterative convergence errors, or programming errors.Finally, the model itself might be based on uncertain footings such as assumptions or model simplifications.
In this work first-order local sensitivity analysis is being performed on a conjugate heat transfer case.This method derives the first-order derivatives of the model equations and calculates those derivatives with the results of a simulation.According to Saltelli and Annoni (2010), this method behaves similarly to one-at-a-time (OAT) screening in a highdimensional parameter space.This method is more convenient to the engineer as it deals with results of the actual simulation without having to perform many simulations with parameters that are not of immediate interest.Starting from the simulation results small perturbations are excited in order to discover the surroundings of the solutions for all parameters of interest.In this way, the outcome of a perturbation can be ascribed to a particular parameter.Contrary to stochastic methods, non-influential parameters cannot falsely be identified as important.Furthermore, a model failure can be ascribed to one specific parameter.Finally, it helps the engineer to understand the behavior of the system in an easy to follow way.The draw-back of this method is its limited validity for extrapolation in non-linear systems.
To the author's knowledge, there are only a few papers devoted to sensitivity analysis for thermal cavity radiation including conduction and convection.Blackwell et al. (1998) investigated a control volume finite element system with conjugate heat transfer, but only simplified radiation to the environment without the use of view factors and reflections were considered.Taylor et al. (1993) gave a treatment on the impact of view factors on the outcome of the energy balances for cavity radiation.They were able to show by an academic example, that the accuracy of some form factors might influence the whole system, significantly.They also developed sensitivity balances for the emissivity and temperature for diffuse-gray radiation when a temperature is set as boundary condition.They further elaborated on this topic by expanding the sensitivities to changes in area and heat flux boundaries (Taylor et al., 1994(Taylor et al., , 1995)).They found that enforcing reciprocity and closure brings great improvements to the energy balances.Taylor & Luck (1995) augmented those findings with the study of various methods in reciprocity and closure enforcement for view factors.They found that the weighted least squares correction for view factors gave the most promising result in producing meaningful configuration factors out of approximate view factors.Korycki (2006) derived local sensitivity factors of the first order in the context of the finite element method (FEM).He derived expressions for conjugate heat transfer including the effect of changing view factors due to shape optimizations and heat transfer in participating media.He also offered an adjoint formulation for sensitivity analysis.Bhatia and Livne (2008) developed sensitivity equations for weak coupling.They employed two sets of meshes.On the coarser mesh, external boundaries were calculated using FEM whereas the finer mesh was used for radiation in a cavity.Both meshes were coupled by a connectivity matrix.Since thermal stresses and deformations were investigated, deformed view factors were approximated by a first order Taylor series expansion, which showed some descent results.In another paper (Bhatia & Livne, 2009) they also showed a way to put all temperature dependent parameters into the main diagonal of a matrix with consequent Taylor series expansion to get the inverse of a matrix.Furthermore, transient behavior was introduced.
This article is a continuation of previously published papers (Rauch & Almbauer, 2010a;Rauch, 2011b) that dealt with steady-state uncertainty estimations of some parameters for conjugate heat transfer with cavity radiation.Later, the transient (Rauch, 2011a) case was investigated.The present author also treated view factors in conjugate heat transfer analysis for steady-state (Rauch & Almbauer, 2010) and transient (Rauch, 2011) cases.
In this treatment, this approach of steady considerations of sensitivities based on equations of an element-centered finite difference method (FDM) is extended for a solar heat source term.Therefore, the equation of first order local discrete sensitivities for wall temperature with respect to solar heat flux is being derived.The governing equation is first discretized followed by a first order derivation.Alternatively, the continuous approach would first do the derivation and than discretize (Lange & Anderson, 2010).Two formulations, namely with and without reciprocity relation for view factors, are being considered.The following section is devoted to introducing a simple case resembling a simplified flat plate solar collector.In section four a sensitivity and uncertainty analysis is performed.The numerical stability of the solution matrix for the sensitivity coefficients is investigated.The effect of partial sensitivity analysis on uncertainty coefficients is shown.The difference of incorporating reciprocity for uncertainty coefficients is demonstrated.A comparison to uncertainty factors of other parameters is made.This work aims to contribute to the enhancement of understanding the simulated model at hand by showing where input parameters have a large influence on the results.

Theory
In their work, Siegel & Howell (2002) show how to formulate the net-radiation method, which is based upon first principle thermodynamics in non-participating media.

Basic equations for conjugate heat transfer
To help understand the following treatment the geometry of two thermal nodes is shown in Fig. 1.A thermal node k consists of a center and a virtual extension.The node has an area A k where it can exchange energy with the surroundings by convection and radiation.The dashed line in Fig. 1 marks the border to the thermal back node, if applicable.Through that border exchange from the surface node happens by conduction, only.The area A kt marks the area to the neighboring node t, where heat is transferred from the surface node k to t by conduction.
First, the basic equations for conjugate heat transfer are introduced, starting with the relation for a given boundary heat flux b as a function of radiosity q o , view factor F sk between nodes s and k, and area A for an element k or s: The running index s spans all N thermal nodes of the geometry.This also includes the backside of a face, where applicable.The summation over all surface nodes sets thermal radiation apart from convective and conductive heat transfer.Where for the latter two heat transfer modes only neighboring elements for low order methods have an impact, in thermal cavity radiation all entities have to be included resulting in a dense solution matrix.The boundary heat flux in Eq. ( 1) consists of other heat transfer modes, source, and storage terms: T f,k represents the fluid temperature for node k, T k the wall temperature of node k, h k the convective heat transfer coefficient, k kt the thermal conductivity, which in this study is presumed to be independent of temperature, l kt the distance between two thermal node centers, q solar,k the solar heat flux, N Cond the number of thermally conducting neighbors, and T t the wall temperature of the neighboring node t.The solar irradiance is attenuated by the atmosphere and strikes at the wall.There it can be reflected and in case of transparent material transmitted.The solar heat flux q solar , with dots for heat fluxes being omitted in this treatment, is the actual heat absorbed by the wall.The radiation term of Eq. ( 1) is placed on the left-hand side (LHS) of Eq. ( 2) because the thermal radiation formulation does not follow thermodynamic sign convention.
With wall temperature and boundary heat fluxes as unknowns a second relationship is needed with emissivity  k : Finally, the fourth relationship renders the set of equations solvable with the Stephan-Boltzmann constant : Local sensitivity coefficients of the first order are obtained by deriving the first order differentiation with respect to the solar heat flux.The order of the derivations follows the order of the four basic equations.

Sensitivity of solar heat flux
Starting with Eq. ( 1), total derivation yields: The underlying assumption is that a change in the solar heat flux is affecting the radiosity but not the geometry including view factors.Also, no temperature dependant properties are considered.The partial derivative of the boundary heat flux with respect to radiosity with  as the Kronecker delta function yields: At this point it is noteworthy no mention that in this treatment a derivation of a particular radiosity with respect to another radiosity is only considered when they are identical.In other words, a perturbation is kept local.The transportation of the perturbation information to the entire system is done by solving for the radiosity derivative with respect to the solar flux.Substitution of Eq. ( 6) into Eq.( 5) results in: The Kronecker delta function helps in the implementation and understanding of the equation as the radiosity q o,k in Eq. ( 1) had been drawn into the summation term.Next, the total derivative for Eq. ( 2) is formulated: The simplification of a local perturbation states: ,, ,, 1: , 0: The derivation of the boundary flow with respect to the solar heat flux is straightforward: Likewise, the derivation with respect to the wall temperature T k gives: And for the neighboring wall temperature T t : Substituting Eq. ( 9) to (12) into Eq.( 8) results in: ,, , 11 To bring the radiation term together with convection and conduction Eq. ( 13) is set equal to Eq. ( 7): This equation has two unknown gradients, namely radiosity and wall temperature with respect to solar heat flux.
Since the emissivity is not a function of temperature in this work, Eq. ( 3) is differentiated in the same way as in Eq. ( 5) which gives: Substitution into Eq.( 5) yields: Equation ( 4) is the last relation for total derivation: ,, Equation ( 19) is set equal with Eq. ( 16):  20) provides an expression for the temperature gradient when reformulated: The same relationship will hold for a neighboring node.Therefore, index k in Eq. ( 21) is replaced with t: These two equations can be substituted into Eq.( 14): When rearranging the terms with the unknown radiosity gradients to the left-hand side and the other terms to the right-hand side (RHS) the final relation is obtained:

Reciprocity
The equations of the previous chapter do not assume any interaction between view factors.But for physically meaningful view factors the reciprocity principle needs to be satisfied along with the closure relation.The investigation of closure is beyond the scope of this work.The interested reader may be referred to the work of Taylor & Luck (1995).The reciprocity relation is formulated as follows: It states, that view factors can be transformed into each other by taking into account the size of the nodes.In this work this principle will be applied in reformulating the original equations.A difference of results for sensitivity formulations with and without reciprocity assumption will, apart from numerical issues, reflect the quality of the view factors.
The derivation for the solar flux with respect to radiosity with the reciprocity relation is identical with the one in the previous chapter.Consequently, it suffices to substitute Eq. ( 25) into Eq.( 21) and Eq. ( 24).

Uncertainty analysis
The temperature or radiosity derivatives of the previous sections are called sensitivity coefficients.Once these are obtained they can be further processed to yield uncertainty coefficients by the following equation: U c is the combined standard uncertainty and u is the standard uncertainty for the parameters q solar , T f , and conductance C which is defined as: The terms with the square root are the dimensionless uncertainty factors (UF) and reflect the changeability of the system to a parameter.The standard uncertainties can be supplied as difference for each parameter in form of a scalar or a probability density function.Provided the standard uncertainties for all parameters are the same, the uncertainty factors can be directly compared.Ideally, the standard uncertainty should be drawn into the square root but often a detailed knowledge of standard deviations, biases, or errors for each thermal node is not available.In this paper no distinction between errors, bias, or deviations is being made which would be beyond the scope of this work.A problem with Eq. ( 28) is when parameters become zero, because then the UF will lose some information.This can happen, when the effect of solar heat flux is investigated at nodes where no solar flux actually strikes.In order to circumvent this issue the following relationship is suggested: A further advantage with this formulation is in its usability.The temperature information is already at hand and by multiplying with the standard uncertainties yields a range of temperatures or ranges of temperatures when using a probability density function.

Workflow
When performing an uncertainty analysis the following work-flow can be applied as described in Fig. 2.
The sensitivity and uncertainty analysis used here is performed a posteriori.First a conjugate heat transfer (CHT) simulation has to be run, be it a structural mechanics, a CFD, or a specialized heat transfer code.Then the radiosity gradients of Eq. ( 24) or Eq. ( 27) have to be calculated.This is followed by solving Eq. ( 21) or Eq. ( 26) for the temperature gradients.When uncertainty factors are desired, these temperature gradients need to be summed.This requires a lot of radiosity gradients.Wherefore, it is recommended to employ LU factorization for the LHS of the first set of equations, followed by matrix vector multiplication of the RHS of those relations.When using Eq. ( 28) or Eq. ( 30) a temperature range can be obtained.

Example
For demonstration purposes a simplified test case has been set up as shown in Fig. 3. Some basic data has been obtained from Rodríguez-Hidalgo et. al (2011).The case represents a flat plate solar collector.It consists of a glass panel with a transmittance of 0.76 and a reflectance of 0.08.The absorber sheet is a crucial element.It consists of a thin copper sheet with a special surface coating with high solar absorbance but low emissivity.The heat absorbed is passed on the water pipe made of copper.The frame is made of aluminium.In Fig. 3 a small gap between the absorber sheet and the frame should signify, that those parts are actually insulated.The glass panel also is insulated from the front frame.Normally, between the absorber sheet, water pipe, and back frame there is a thick insulation layer.But because for ease of modelling the system in a finite difference solver and the steady-state calculation the insulation chamber is modelled as a vacuum.Basic parameters are to be found in tables 1a and 1b.As the connection between the absorber sheet and the water pipe would be a line the conduction is modelled by conduction bridges.The conductance as defined in Eq. ( 29) is estimated to be 77.1308W/K.This value is based on the smaller nodal area between sheet and pipe element, the temperature-independent thermal conductivity of copper and the mean distance between the two nodes.The environment is modelled as a desert in Arizona using environmental data from the commercial solver RadThermIR v.9.1.The solar irradiance is 950 W/m² and the zenith angle is 42°.The wind model of McAdam is used with the velocities in tables 1a and 1b which results in a heat transfer coefficient of 14.06 W/m²/K for the exterior.The air chamber is modelled as a single air node and the resulting temperature is found in tables 1a and 1b.Table 2 shows the results of the thermal radiation solver with lowest view factor settings.
The four front nodes are situated at the very centre of the respective part.The total thermal node count amounts to 3160.

Analysis
The results from the previous section form the basis for the following analysis.

Numerical stability
The infinity norm M  of the matrix M is one easily to be obtained indicator for stability.It is defined as follows with max as maximum operator and a ij as a matrix entry of row i and column j: Another method for stability analysis is by using eigenvalues.The Gerschgorin circle theorem estimates eigenvalues  and thus stability of a solver because it gives the maximum possible value of  by the Gerschgorin radius r.The eigenvalue  assures stability if it is less than or equal to 1.This is guaranteed when the radius r is less than or equal 0.5.In that case it can be expected, that the set of equations can be handled by an iterative solver without preconditioning. 1 Figure 6 shows the Gerschgorin radii r and the main diagonal entries for the nodes pertaining to the front side of the absorber sheet when only radiation and convection terms are considered.It reveals that although the radii are greater than 0.5 the diagonal entries are even greater.Thus, the absorber sheet behaves numerically stable because of its diagonal dominance.
Figure 7 shows the absorber sheet when all heat transfer modes are taken into account.The first thing to note, are the high values for both radii and diagonal entries as compared to Fig. 6.This supports the findings of t h e i n f i n i t y n o r m i n F i g .5, that the introduction of conduction has a destabilizing effect.The second finding is that now all radii are greater than their corresponding diagonal entry.This supports, along with the need of multiple solutions for uncertainty factors as mentioned in section 2.5, the use of a direct solver.

Partial uncertainty analysis
Seeing the stability issue when including all heat transfer terms there would be another incentive for omitting a heat transfer mode.For example, when considering conduction all neighboring nodes need to be processed.This cannot be felt while generating the solution matrix.But one needs to remember that the RHS needs to be recalculated for every sensitivity coefficient.In the present example this amounts to 3160 recalculations.This leads to the longer calculation time when estimating the uncertainty factor of, for example, conductance.In fact, the double summation of the last term in Eq. ( 30) represents together with a loop for the RHS calculations an N³ time complexity.The changes in the order of magnitude emphasize the need to include all heat transfer modes into the calculations.

Reciprocity analysis
In section 3.2 the reciprocity relation was introduced together with the formulation of the sensitivity coefficient for solar flux.The difference between the formulations with and without reciprocity in Fig. 9 show that the quality of the view factors behave quite satisfactorily with regards to reciprocity.The differences are not a matter of right or wrong with respect to uncertainties.Which formulation to use is a matter of whether the radiation solver employed uses the reciprocity relation.The solver used in this work does and for this reason all figures use the reciprocity formulation.

Results
The effect of a single perturbation at node 373 shows Fig. 10.At the node itself the impact is the largest and fades quickly with distance.But the information of the excitation is passed on to the whole absorber sheet.

Conclusion
Two formulations for solar heat flux uncertainties were derived.Partial uncertainties should not be used.Uncertainty factors can aid in understanding at what nodes which parameters need more attention.The combined standard uncertainty predicted a change in the solar heat flux of the test case very well.But this is just the beginning of the investigations.Questions of numerical stability, memory and time complexity need to be addressed.Also because of the non-linearity the range of predictability needs to be investigated.Very small values can cause high sensitivity coefficients.A draw-back of this formulation is that it can only be applied where radiosities have been calculated.This hampers the use of this method when systems are modeled without radiation in order to speed-up calculations or when contact between two surfaces is modeled where radiation cannot occur.

Acknowledgment
The author would like to acknowledge the financial support of the "COMET K2 -Competence Centres for Excellent Technologies Programme" of the Austrian Federal Ministry for Transport, Innovation and Technology (BMVIT), the Austrian Federal Ministry of Economy, Family and Youth (BMWFJ), the Austrian Research Promotion Agency (FFG), the Province of Styria and the Styrian Business Promotion Agency (SFG).The author would furthermore like to express his thanks to his supporting scientific partner Graz University of Technology.Vol. 46, No. 3, pp. 578-590, doi:10.2514/1.26236

Fig. 4 .
Fig. 4. Surface mesh of solar panel Figure 4 shows the mesh of the panel with the glass panel and absorber sheet removed.The red lines indicate insulation to prevent conduction.Care is taken to use well shaped quadrilaterals.The outer dimensions are 2 x 0.19 x 0.082 m.The small width is caused because only one pass of the water serpentine is considered.The length of the quadrilaterals is about 0.04 m and the diameter of the pipe is 0.01 m.The water enters the pipe from left with an inlet temperature of 308 K and a flow rate of 0.3x10 -4 m³s -1 .The water is modelled with 50 fluid nodes along the pipe allowing for localized fluid temperatures.

Fig. 5 .
Fig. 5. Infinity norm of the solar panelFigure5shows the infinity norm of the solution matrix for sensitivity coefficients for formulations with reciprocity of the LHS of Eq. (27).The set of bars denote whether only the radiation term (Radiation), the radiation and convection terms (Convection), radiation and conduction terms (Conduction), or all three terms (All) for the sensitivity coefficients had been used.It can be seen that for Radiation the infinity norm is quite low, because many

Fig. 6 .
Fig. 6.Gerschgorin radii of the absorber sheet with convection and radiation term

Figure 8
Figure 8 shows uncertainty coefficients for the solar heat flux with the reciprocity relation, when incorporating various heat transfer mode.The back frame and the glass panel can handle the increase in solar flux quite well when radiation is the only heat transfer mode.But the absorber sheet and the water pipe s h o w u n p h y s i c a l b e havior.Introducing convection helps those nodes in passing the heat on to the water.

Fig. 14 .
Fig. 14.Non-normalized UF for the fluid temperature Finally, the obtained UFs from Fig. 13 shall be used to estimate the temperature range.Solar irradiation data that is available on a daily or monthly data for various locations can hold significant uncertainties as Gueymard & Thevenard (2009) has shown for the location Colorado, USA.To see whether the wall temperatures can be estimated when the solar fluxes change, a simulation is run with a solar irradiance of 1050 W/m².The obtained solar fluxes for the four nodes along with the UFs of Fig. 13 are put into Eq.(30).The positive combined uncertainties calculated for solar flux only is added to the temperatures in table 2 and expressed as percentage difference to the simulated temperatures with the case of solar irradiance of 1050 W/m².

Fig. 15 .
Fig. 15.Difference predicted to simulated temperatures The predicted temperatures with UFs in Fig. 15 are in good agreement with the simulation.

Table 1b .
Boundary conditions continued