Question 13.4: A plane layer of radiating and isotropically scattering medi...

With conduction, radiation, and internal heat generation, the energy equation (Equation 10.1)

\rho c_p\frac{DT}{Dt}=\beta T\frac{DP}{Dt}+\nabla . (k\nabla T-q_r)+\dot{q} +\phi _d              (10.1)

is \nabla .(k\nabla T-q_r)+\dot{q}=0 or, using a summation form for the radiation and conduction fluxes,

\sum\limits_{i=1}^{3}{\frac{\partial q_{r,i}}{\partial x_i} } +\sum\limits_{i=1}^{3}{\frac{\partial q_{e,i}}{\partial x_i}}=\dot{q}                    (13.55)

For the present one-dimensional problem, since q_r = I^{(1)} , where I^{(1)}  is the first moment of the intensity, the energy equation can be put in the dimensionless form

\frac{d \widetilde{I} ^{(1)}}{d\tau _1}  =4N\frac{d^2\vartheta }{d\tau ^2_1} +\frac{\dot{S} }{\tau _D}                      (13.56)

where N =k  (κ + σ_s )/4σ T^3_{w1},\dot{S} =\dot{q} D/\sigma T^4_{w1},\tau _1=(k+\sigma _s)x_1,\tau _D=(k+\sigma _s)D,\widetilde{I} ^{(1)}=I^{(1)}/\sigma T^4_{w1}, and ϑ = T/ T_{w1}. This is the defining equation for the derivative of I^{(1)} . In the pure radiation solution of Example 13.3 without internal heat sources, the derivative of the first moment in Equation 13.43

\frac{dI^{(1)}}{d\tau _1}  =(1-\omega )(4\pi I_b-I^{(0)})                     (13.43)

was set equal to zero, because for these conditions \widetilde{I}^{(1)}=q_r/\sigma T^4_{w1} is constant and {d\widetilde{I}^{(1)}}/{d\tau _1}=0 .This is not the case here with conduction and an internal heat source included. The presence of a second derivative of dimensionless temperature requires two boundary conditions for solving the energy equation: \vartheta (τ = 0) = 1 and \vartheta (τ_1 = \tau _D) = T_{w2}/T_{w1} .

To proceed with the solution, two coupled second-order differential equations are derived for\widetilde{I}^{(0)} and\widetilde{I}^{(1)} (the nondimensional intensity integrated over all solid angles, and the radiative flux). The first is obtained by equating Equation 13.56 and the first-moment differential equation (Equation 13.43) in the P_N method. For the present one-dimensional problem, this yields

\frac{4N}{1-\omega } \frac{d^2\vartheta }{d\tau ^2_1}+\frac{\dot{S} }{\tau _D(1-\omega )}-4\vartheta ^4=-\widetilde{I} ^{(0)}                   (13.57)

where \widetilde{I} ^{(0)}=I^{(0)} /\sigma T^4_{w1} . The second equation is found by substituting the closure Equation 13.38 into the second-moment differential Equation 13.34 to obtain the relation between \widetilde{I} ^{(0)}   and \widetilde{I} ^{(1)} 

I^{(ij)} =\frac{1}{3}\delta _{ij}I^{(0)}          (13.38)

\sum\limits_{i=1}^{3}{\frac{\partial I^{(ij)}}{\partial \tau _i} } =-I^{(j)}             (3 equations: j=1,2,3)            (13.34)

\frac{d\widetilde{I} ^{(11)}}{d\tau _1} =\frac{1}{3}\frac{d\widetilde{I}^{(0)} }{d\tau _1} =-\widetilde{I}^{(1)}                  (13.58)

Now, Equation 13.58 is differentiated with respect to τ_1, and the result is substituted into Equation 13.56 to yield

\frac{d^2 \widetilde{I} ^{(0)}}{d\tau^2 _1}  +12N\frac{d^2\vartheta }{d\tau ^2_1} +\frac{3\dot{S} }{\tau _D}=0                     (13.59)

Equations 13.57 and 13.59 can be solved simultaneously for \widetilde{I} ^{(0)}   and  ϑ.They can be combined into a single fourth-order equation in ϑ by differentiating Equation 13.57 twice with respect to τ_1  and the result is substituted into Equation 13.59 to eliminate the second derivative of  \widetilde{I} ^{(0)} . The resulting fourth-order equation, or the two second-order equations, require two boundary conditions in addition to the known boundary surface temperatures. These are generated from the conditions in Equations 13.40 and 13.41 using Equation 13.50 to eliminate J, which results in

\int_{\Omega =0}^{2\pi }{I_0(\Omega )} Y_l^m(\Omega ) d\Omega =\int_{0}^{2\pi }{\left[\epsilon (\Omega )I_{b,w}+\frac{1}{\pi }\int_{\Omega _i=0}^{2\pi }{\rho (\Omega ,\Omega _i)I(\Omega _i)l_id\Omega _i}  \right]d\Omega  }                         (13.40)

\int_{\Omega =0}^{2\pi }{I_0(\Omega )} Y_l^m(\Omega ) d\Omega =\int_{\Omega =0}^{2\pi }{I_0(\Omega ) l_i^{2n-1}d\Omega } =  \int_{\Omega =0}^{2\pi }{I_0(\Omega )l_id\Omega }=j                (13.41)

J(\tau _1=0)=\sigma T^4_{w1}-\frac{\left(1-\epsilon _w\right) }{\epsilon _w}q_r                  (13.50)

\frac{1}{4} \widetilde{I}^{(0)} (\tau _1=0)=1-\left\lgroup\frac{1}{\epsilon _w}-\frac{1}{2} \right\rgroup \widetilde{I}^{(1)} (\tau _1=0)
\frac{1}{4} \widetilde{I}^{(0)} (\tau _1=\tau _D)=\vartheta ^4_{w2}+\left\lgroup\frac{1}{\epsilon _w}-\frac{1}{2} \right\rgroup \widetilde{I}^{(1)} (\tau _1=\tau _D)                  (13.60)

or

\frac{1}{4} \widetilde{I}^{(0)}_i =\vartheta ^4_{wi}\pm \left\lgroup\frac{1}{\epsilon _w}-\frac{1}{2} \right\rgroup \widetilde{I}^{(1)}_i           (13.61)

where the i subscript denotes walls 1 or 2 and the positive sign applies at wall i = 2. Inserting Equation 13.58 to eliminate \widetilde{I}^{(1)} results in the final boundary relations for \widetilde{I}^{(0)}

\pm \left\lgroup\frac{d\widetilde{I}^{(0)} }{d\tau _1} \right\rgroup_i=-\frac{3}{4\left\lgroup\frac{1}{\epsilon _w} -\frac{1}{2} \right\rgroup }[\widetilde{I}_i^{(0)} -4\vartheta ^4_{wi}]             (13.62)

These can be directly applied as the boundary conditions for Equation 13.59. The problem is completely specified with two second-order differential equations, that is, Equations 13.57 and 13.59 and the four boundary conditions, the specified boundary temperatures for Equation 13.57, and the conditions (i = 1, 2) in Equation 13.62 for Equation 13.59. An iterative numerical solution can be used to obtain ϑ(τ ) and \widetilde{I}^{(0)}_i (\tau _1). Then \widetilde{I}^{(1)} at the boundaries (i = 1, 2) is found from Equation 13.61, which gives the desired radiative fluxes at the boundaries. The total energy transfer also requires the amount of heat conduction. Since the temperature distribution has been determined, this can be found by evaluating −kdT /dx at the boundaries.