The Direct Problem

For solving an inverse problem of radiation in a participating media, we need to solve the problem of direct radiation. A fundamental quantity in the study of radiation transfer for a participating media is the spectral radiation intensity, ^Iλ

( )

^s,Ω) , where Ω)

is the direction of propagation and s is the path of propagation. It represents the flow of radiation energy per unit area normal to the direction of propagation of the radiation beam, per unit wavelength, per unit solid angle, per unit time. If the energy per unit time is measured in Watts, the wavelength λ is measured in μm, and the solid angle in steradian, sr, then the dimension of the spectral radiation intensity, ^Iλ

( )

^s,Ω) , becomes W/(m². μm. sr) where the area is measured perpendicular to the direction of propagation of the radiation beam [106]. In the study of radiation transfer, the radiation intensity is

the fundamental quantity, which is obtained through solving the Equation of Radiative Transfer (ETR).

In this study the horizontal coordinate x is used to measure linear distances in the y directional discretized cross-sections of the jet, see Fig. 64. The one dimensional radiative transport equation for spectral intensity I_λ passing through an absorbing, emitting and scattering medium is [72]

( )

⁺

( ( )

⁺

( ) ) ( )

⁼

( ) ( )

⁺

( ) ∫

^Φ

(

^′

) (

^′

)

^Ω′ the scattering phase function giving the scattered intensity from direction μ′ to μ; I_λ⁰ is the spectral radiation intensity of the blackbody at medium temperature T and the λ in the subscript refers to the wavelength. At local thermodynamic equilibrium, the emission gain is proportional to Plank’s law, ^I_λ^{0 =2C}₁

(

^λ5

(

^exp

(

^C₂

( )

^λT

)

⁻¹

) )

^{where C}₁

and C₂ are Planck’s radiation constants given by C₁=0.59552197×10⁸Wμm⁴/m²sr and C₂ =14387.69μmK, respectively.

The sum of the spectral scattering and absorption coefficients is called the spectral extinction coefficient k_λ^ext, i.e. k_λ^ext =α_λ +σ_λ and the ratio of the scattering to extinction coefficient, i.e. ω_λ =σ_λ /k_λ^ext is called the spectral single scattering albedo [99].

For infrared radiation the scattering process is often neglected, since the scattering effect is small in the infrared domain [89]. The line of sight is considered orthogonal to the symmetry axes of the jet

(

μ =1

)

. With these assumptions, in each cross-section (z=const.) of the jet, the ERT for instantaneous spectral radiation intensity is simplified and can be written according to the notations in Fig. 64

( )

_x^x,y

(

^{x2 y2}

)

^I

^{( )}

^x,y

(

^{x2 y2}

) (

^I0

(

^{T x2 y2}

) )

The thickness of the jet is dependent on the radial distance from the centreline y, and on the axial distance from the exit z, i.e. l

( )

z,y =2 R

( )

z ²−y² . Since the spectral interval of the camera is large, the detected radiation cannot be validly considered as elementary. The expression of I_λ must be integrated over the spectral sensitivity range of the device, and effective radiation properties over the device sensitivity range must be used [97]. For example the mean absorption coefficient inside the band is

∫

where Δλ is the spectral range of the detector. The weighted integral presented in (50) is also a Planck-averaged quantity. The expression of ERT in (49) has to be rewritten into the band averaged form

( ) (

^{x2 y2}

)

^I

^{( )}

^x,y

(

^{x2 y2}

) (

^I0

(

^{T x2 y2}

) )

After time-averaging, the ERT can be obtained in the following form (dependences assigned between parentheses in (51) are not indicated in the next but conditionally are implied.)

where the overbar denotes the mean value obtained by the time-averaging, and the prime denotes the fluctuating component. The time-averaged equation contains unknown fluctuating correlated terms. The correlation α_Δ′_λI_Δ′_λ is negligible in the case of optically thin fluctuation approximations [94]. If the optical thickness of the radiating medium, based on the macro-scale of turbulence, is small compared to unity, then

=0

′

′ _Δ

Δλ λ

α I . Physically, the correlation α_Δ′_λI_Δ′_λ represents the effect of the local fluctuations in the emission process on the transport of the time-averaged radiation intensity. The last correlation term in (52) can be negligible if it can be assumed that the spectral absorption coefficient is statistically independent of the temperature [94].

Making use of the approximation stated above, the time- and band-averaged equation of radiative transfer is reduced to the form

0λ intensity at x=x₀ is identical to the background radiation intensity I_Δ^bg_λ. Assuming that the background emissivity is close to unity, the background radiation intensity can be expressed with the Planck function. I_Δ^bg_λ =I_Δ⁰_λ

( )

T₀ .

The time-averaged blackbody radiation intensity I_Δ⁰_λ in the medium is a function of the time-averaged temperature and of the fluctuating intensity of temperature. The assumed probability density functions (PDF) are often used in the study of turbulent flow, for

example Gaussian distribution [94]. In this study I have used a Gaussian distribution to simulate the turbulent fluctuation of temperature as well. The Gaussian PDF of the temperature is

where T is the time-averaged temperature and σ_T is the standard deviation of temperature. The time and band-averaged blackbody spectral radiation intensity can be approximately expressed as

This means that the time and band-averaged blackbody radiation intensity is a function of the time-averaged temperature and its standard deviation. The standard deviation σ_T can be calculated from the time series of fictitious temperature images.

If scattering is neglected, the only optical property required in the radiative transfer calculation is the absorption coefficient. The parameterization of the absorption coefficient in terms of physical quantities is applied in terms of the cloud’s liquid water content and the effective droplet radius is defined from the droplet size distribution in Appendix G.

The parameterisation for the optical properties of water clouds in the infrared spectrum have been presented in [89]. The coefficients for the rational functions of optical properties can be determined by performing a least square fitting to points calculated by the exact Mie theory [99]. The accuracy of the parameterization is within 5% compared to the exact Mie calculation. I applied a band averaged absorption coefficient parameterization in this experiment. The selected band limit (10.2-12.5 μm) corresponds to the operating band of the detector used for the measurements. The parameterization is valid for an effective radius in the region of 2-40 μm. The band-averaged extinction coefficient k_Δ^ext_λ is normalized by dividing it with the liquid water content (LWC), ψ_Δλ =k_Δ^ext_λ LWC. The band-averaged value for the specific extinction is defined with respect to the Planck-function

∫

The simple rational function to parameterize the specific extinction coefficient ψ_Δλ from [89] has the form

3 4 2 3 2

1

0 are a /re a /re a /re

a + + + +

Δλ =

ψ . (57)

The coefficients are the following: a₀ =2.9⋅10⁻²m²g⁻¹, a₁ =−3.9657⋅10⁻⁴μm⁻¹m²g⁻¹,

1 2 2 =1.4902μm m g⁻

a , a₃ =−5.0777μm²m²g⁻¹ and a₄ =5.2170μm³m²g⁻¹.

The rational function for the band-averaged single scattering coalbedo,

(

1−ω_Δλ

)

has been formulated as

2 3 2 1

1−ω_Δλ =b0+b /r_e+br_e+br_e , (58)

with coefficients of b₀ =4.9787⋅10⁻¹, mb₁ =7.4581⋅10⁻¹μ , b₂ =−3.7379⋅10⁻³μm⁻¹ and b₃ =9.0555⋅10⁻⁵μm⁻².

The band-averaged absorption coefficient can be determined from the band-averaged single scattering albedo and extinction coefficient, α_Δλ =

(

1−ω_Δλ

)

k_Δ^ext_λ. The expression for the liquid water content can be derived in the form

(

_λ

)

_λ

λ ω ψ

α_Δ − _{Δ Δ}

= / 1

LWC ,

(

^{g m3}

)

^. (59)

The effective droplet radius has to be obtained before the calculation of the LWC. 3.6.2 The Inverse Analysis

In this experiment I approximated the unknown absorption coefficient and temperature distributions in a bell shape function form. By exploiting the axisymmetric round jet assumption, the function of α

( )

δ and T

( )

δ , δ = x²+y² , in the jet cross-sections can be approximated well by a bell-shape function, see Fig. 77. Several functions are suitable for approximating a bell-shape curve. I approximated the time-averaged absorption coefficient profile by the function of a path variable δ as

( )

δ = p0

(

arctan

(

p1

(

δ+r12

) )

−arctan

(

p1

(

δ −r12

) ) )

α , (60)

where δ∈

[ ]

0,R and R=2r₁₂. The time-averaged temperature distribution in the cross-section can be approximated by a similar function

( )

q0 q1

(

arctan

(

q2

(

r12

) )

arctan

(

q2

(

r12

) ) )

T δ = + δ+ − δ− ,q₀ >T_bg, (61)

where T_bg is the background temperature. Because the radiation of the jet can be detected only when there is a temperature difference between the jet and the background, z₀ >T_bg has to be ensured.

The problem defined by (53) with the absorption coefficient and temperature are unknown, but the measured exit intensities are known, is an inverse problem, which can be solved by the minimization of the following objective function

( ) _∑ ( ( ) ( ) )

= Δ −

= ^M

i I i x; Y x

S i

1

2 1 1 p,q

q

p, _λ , (62)

where Y_i

( )

l is the measured exit radiation intensities, M is the total number of measurements, I_Δλi

(

x₁;p,q

)

is the estimated exit radiation intensity at the surface of the jet and p,q are the vectors of unknown parameters. The estimated exit radiation intensities are obtained through solving the direct problem, by using an estimate for the vector of unknown parameters p,q.

Considering the assumptions above, the vectors of unknown parameters in the inverse problem are p=

[

p0;p1

]

and q=

[

q0;q1;q2

]

. Liu et al. in [93] have published that a simultaneous estimation of the temperature and the absorption coefficient from the exit radiation measurement results in an objective function which is not a unimodal function of all unknown parameters, therefore temperature and the absorption coefficient cannot be estimated simultaneously in a single minimization procedure.

In this experiment the following iterative procedure is developed for solving the minimization problem in (62).

(i) First cycling: estimation of an initial absorption coefficient distribution

Initially it is assumed, that the temperature of water droplets in a cross-section can be approximated as constant T_c with a standard deviation of σ_T. T_c is the expected maximum temperature in the cross-section, which must be an input parameter.

Therefore, only the parameters of the absorption coefficient function p are unknown and have to be estimated by the minimization of the objective function

( ) _∑ ( ( ) ( ) )

= Δ −

= ^M

i I i x; ,Tc, T Y x

S i

1

2 1 0

1

0 _λ p σ

p . (63)

The solution is an initial guess to parameters of function α

( )

δ in (60). Superscript refers to the iteration.

(ii) Second cycling: estimation of an initial temperature distribution

A coarse estimation of the temperature distribution function parameter q¹ is derived from solving the following minimization problem with the initialized parameters of

Tbg

q¹₀ = , (64)

( ) ( ( ( ) ) ( (

12

) ) )

The objective function is

( ) _∑ ( ( ( ) ) ( ) )

calculation of q¹, I utilized the correlation between the absorption coefficient and the temperature profiles.

(iii) Third (iterative) cycling: simultaneous approximation of the absorption coefficient and temperature distributions

An iterative procedure follows the preliminary steps. The parameters of the temperature distribution function are collected from the preceding two iteration cycles

bg objective function is the same as in (65), the minimization searches for p^j in the form

( ) _∑ ( ( ( ) ) ( ) )

(iiii) Fourth cycling: refinement of the absorption coefficient and temperature distributions

The iteration procedure for p^j is followed by the minimization of the objective function of

( ) _∑ ( ( ( ) ) ( ) )

= Δ −

= ^M

i

j j

i x; , , Y x

I

S _i

1

2 1

1 α δ

λ q p

q . (68)

By solving this inverse problem the temperature distribution function can be refined by applying the absorption coefficient distribution function ^{α j}

(

^pj,δ

)

. Initial values for q have to move away a little bit from the last q^j, to enforce the solver to begin a searching process. The complex inverse algorithm is completed with a last minimization of the objective function in (65) to obtain the appropriate parameters of the absorption distribution function. The schematic flowchart of the above iterative minimization process can be seen in Fig. 81.

Every minimization means a complete solution of an inverse algorithm. I applied the conjugate gradient method. The computational algorithm for the conjugate gradient inversion method consists of the following main steps [106]:

Step 1. Guess initial values of unknown parameters;

Step 2. Solve the direct problem;

Step 3. Compute the sensitivity coefficients;

Step 4. Compute the gradient direction and the search step size;

Step 5. Compute the new estimate to the unknown parameters;

Step 6. Check the stopping criterion, if it is not satisfied, go to step 2.

For solving the direct problem I used an ODE solver, because the neglected scattering reduced the direct problem to solution of an ordinary differential equation with variable (non-linear) coefficients.

minimize

Fig. 81. Flow chart of the iterative inverse algorithm

3.7 Results and Discussion

I proved the accuracy of the presented inverse algorithm through test data. The exit radiation intensities were generated from known absorption coefficient distribution functions, i.e. ^p_test⁼

[

^0.22;0.0001

]

. The temperature was set to constant in the cross-section. There differences between the original and the inversion procedure resulted parameters were p^j −p_test <4⋅10⁻²⁸, therefore the resulting function and the test curve are identical, see Fig. 82a for intensity and Fig. 82b for the absorption coefficient.

0 0.2 0.4 0.6 0.8 1

8.14 8.16 8.18 8.2 8.22 8.24 8.26 8.28

x/r

Intensity (W/m2)

Test data Approximation

0 0.2 0.4 0.6 0.8 1

0 1 2

x 10^-4

x/r

Absorption coefficient (m-1)

Test data Approximation

Fig. 82. Verification of the numerical accuracy of the inversion algorithms

3.7.1 Inverse Analysis of FEM Generated Data

Selecting a cross-section at z/D=4, the FEM generated temperature distribution with the estimated radial temperature distribution function are shown in Fig. 83.

0 0.2 0.4 0.6 0.8 1

290 300 310 320 330 340 350

y/R

Temperature (K)

Approximation Simulation

Fig. 83. The simulated and estimated temperature distributions at the cross-section of z/D=4

The simulated and the estimated absorption coefficient profiles can be compared in

near the boundary of the jet. The curve of the assumed temperature goes over the simulated data because in this region a larger absorption coefficient means more liquid water content and results in a higher average temperature. The initial temperature was set to the simulated temperature at the centre line of the selected cross-section. The residual norm determined by the minimization algorithm was equal to 1.7⋅10⁻³.

0 0.2 0.4 0.6 0.8 1

0 0.5 1 1.5 2 2.5 3

x/r Absorption coefficient (m-1 )

Approximation Simulation

Fig. 84. The absorption coefficient distributions of the FEM simulation and the estimation of inversion at the cross-section of z/D=4

Sensitivity Analysis

The most important element of the conjugate gradient method is the evaluation of the sensitivity coefficients ∂I ∂p_i ,∂I ∂q_j with i=0,1 and j=0,1,2, which indicates the variation of the exit radiation intensity with respect to the unknown parameters p_i and qj. Since the unknown parameters can assume different orders of magnitude, the sensitivity coefficients with respect to the various parameters may also differ by several orders of magnitude, creating difficulties in their comparison. These difficulties can be alleviated through the analysis of relative sensitivity coefficients defined as

i i

p p

p I

J i ∂

= ∂ , 2i=1, and

j j

q q

q I

J j ∂

= ∂ , j=0,1,2. (69)

Fig. 85 and Fig. 86 show the radial distributions of the relative sensitivities. Comparing the figures, it can be seen, that the system is more sensitive –about one order of magnitude– to q than to p. Figures reveal that the magnitudes of the sensitivity coefficients are much larger at the centreline than close to the boundary. This could be the reason for the largest deviations between simulation and approximation near the boundary. The sensitivity of the intensity field to q₂ is much larger than that to q₀ and q1. This fact could make the estimation of q₀ and q₁ difficult, but fortunately the correlation between the absorption coefficient and temperature distributions help to

approximate their values in an iterative manner that is presented in the inversion algorithm above.

0 0.2 0.4 0.6 0.8 1

-0.18 -0.16 -0.14 -0.12 -0.1 -0.08 -0.06 -0.04 -0.02 0 0.02

y/R

Relative Sensitivity

p0 p1

Fig. 85. The relative sensitivity coefficients for the absorption coefficient function

0 0.2 0.4 0.6 0.8 1

-1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0

y/R

Relative Sensitivity

q0 q1 q2

Fig. 86. The relative sensitivity coefficients for the temperature distribution

Estimation of Liquid Water Content

(

LWC

)

The formula of the liquid water content (LWC) calculation can be derived from (59).

The liquid water content could be strongly influenced by the assumed effective droplet radii. The LWC distribution of the examined cross-section of simulation at various assumed effective droplet radius can be seen in Fig. 87. Error-bars show the ambiguity in the optical property values defined in [26]. It can be seen that there is a good agreement between the simulation and calculation, the largest deviation is at the half width of the jet. These positive/negative deviations can be smoothed by integration.

0 0.2 0.4 0.6 0.8 1 0

10 20 30 40 50 60 70 80 90 100

x/r LWC( g/m3 )

r_e

r_e=5μm

r_e=35μm

FEM generated data with assumption

of r_e=10μm

Fig. 87. The estimated LWC distribution of the examined cross-section at various assumed effective droplet radii and the the LWC distribution in FEM simulation with assumption of r_e =10μm

The integration of the LWC curve gives the total liquid water content in the examined cross-section. The total liquid water content is important in the environmental load calculations. Assuming that the effective droplet radius is r_e =10μm, The total liquid water content in the cross-section is 0.9058g/m. The integration of the liquid water concentration in the cross-section of FEM simulation results in a total liquid water content of 0.8884g/m. The absolute error between the estimated and simulated data is

g/m

0.017 i.e. 1.9%.

3.7.2 Inverse Analysis of IR Thermography Data

The examined cross-section of the jet is discretized by equally spaced points corresponding to the image-pixels. The path-length end points

[

x0,x1

]

can be calculated from the cross-section half with radius r₁₂, according to R=2r₁₂, see Fig. 64. The pixel-size has to be determined from the focal-plane distant of the measurement. In this case the resolution is 0.0017 m/pixel.

The inverse analysis can be executed at several cross-sections, of which I have presented some examples. Selecting the cross-section at z/ D=5.7, the detected and calculated cross-section intensities are compared in Fig. 88. In this figure the box-plot of the cross-section time-series data is shown. The box has lines at the lower quartile, median, and upper quartile values. The whiskers are lines extending from each end of the box to show the extent of the rest of the data. Outliers are data with values beyond the ends of the whiskers. A thin line represents the time-averaged data. The bold line shows the intensity curve obtained by the presented inverse algorithm. The residual norm between the averaged and the approximated data is 1.56⋅10⁻⁴. It can be seen that the approximation curve lies inside the quartile boxes close to the average values through the cross-section.

1 2 3 4 5 6 7 8 91011121314151617181920212223242526272829303132333435363738 8.05

8.1 8.15 8.2 8.25 8.3 8.35

Intensity (W/m2 )

Radial Pixel Number

Approximation Average

Fig. 88. Box-plot and average of time series data completed by the approximation of inverse analysis at cross-section of z/D=5.7

Estimated temperature and absorption coefficient distributions at two distinct cross-sections are plotted in Fig. 89 and Fig. 90, respectively.

0 0.2 0.4 0.6 0.8 1

0 0.05 0.1 0.15 0.2 0.25 0.3

Absorption Coefficient (m-1 )

300 310 320 330 340 350 360 370

y/R

Temperature (K)

α T

Fig. 89. Resulting temperature and absorption coefficient distributions at z/D=4.2

The temperature near the boundary is higher than the background temperature which agrees with the assumption that radiation intensity could be detected only if there is a temperature difference between the inspected medium and the background. It follows from this that the whole boundary layer could not be detected with passive IR techniques.

The liquid water concentration distribution depends on the effective droplet radius as well. The estimated concentration distributions of the examined cross-section with various assumed effective droplet radii can be seen in Fig. 91.

0 0.2 0.4 0.6 0.8 1 0

0.05 0.1 0.15 0.2 0.25 0.3 0.35

Absorption Coefficient (m-1)

300 310 320 330 340 350 360 370

y/R

Temperature (K)

α T

Fig. 90. Resulting temperature and absorption coefficient distributions at z/D=5.7

0 0.2 0.4 0.6 0.8 1

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

y/R LWC (g/m3 )

r_e= 2 μm r_e= 4 μm r_e= 6 μm r_e= 8 μm r_e=10 μm

Fig. 91. The estimated LWC distributions of the examined cross-section z/D=5.7 with various assumed effective droplet radii

In this study the total liquid water content (TLWC) in cross-sections between z/D=5 and z/D=8 proved to be practically identical. By assumption, the effective droplet radius is equal to 10μm (this assumption is based on the typical cloud droplet sizes at the end of the droplet formation cycle [42]), the total liquid water content in the examined jet transitional region is equal to 2.4⋅10⁻²g/m. This value shows that relatively low liquid water content can be detected with the suggested method.

3.8 Summary and Conclusion

In this study I have proved that thermograms of semitransparent free jets, provided by a common, high temperature resolution IR camera, are suitable not only for qualitative, but quantitative analyses as well. The long-IR imaging is a new technique for off-surface monitoring of turbulent phenomenon. The detected long-IR radiation of the

condensed water droplets in the jet represents the turbulent scalar field. I have used the instantaneous images for qualitative analysis, i.e. for visualisation of turbulent free jet flow. I have analysed the time-averaged images to estimate the liquid water content in the jet.

The fictitious jet temperature provided by the IR camera was only slightly higher than the background temperature, therefore in the instantaneous raw images the jet appears as a semitransparent, diffuse cloud. Removal of the high-frequency noise from the image requires special smoothing techniques, since the observed turbulent phenomenon contains a high frequency component as well. Smoothing of this component is not desirable. I have proved that a fuzzy-rule based diffusion filter with an intensity dependent diffusion coefficient could remove the background noise without decreasing the intensity variations in the jet area of the image. I investigated the filtered images

The fictitious jet temperature provided by the IR camera was only slightly higher than the background temperature, therefore in the instantaneous raw images the jet appears as a semitransparent, diffuse cloud. Removal of the high-frequency noise from the image requires special smoothing techniques, since the observed turbulent phenomenon contains a high frequency component as well. Smoothing of this component is not desirable. I have proved that a fuzzy-rule based diffusion filter with an intensity dependent diffusion coefficient could remove the background noise without decreasing the intensity variations in the jet area of the image. I investigated the filtered images

In document Nemlineáris hőterjedés és hőképanalízis (Pldal 94-0)