Question 24.4: Although it served our purposes for illustrating the shootin...
Although it served our purposes for illustrating the shooting method, Eq. (24.6) was not a completely realistic model for a heated rod. For one thing, such a rod would lose heat by mechanisms such as radiation that are nonlinear.
0=\frac{d^2 T}{d x^2}+h^{\prime}\left(T_{\infty}-T\right) (24.6)
Suppose that the following nonlinear ODE is used to simulate the temperature of the heated rod:
0=\frac{d^2 T}{d x^2}+h^{\prime}\left(T_{\infty}-T\right)+\sigma^{\prime \prime}\left(T_{\infty}^4-T^4\right)
where σ′ = a bulk heat-transfer parameter reflecting the relative impacts of radiation and conduction = 2.7 × 10^{−9} K^{−3} m^{−2}. This equation can serve to illustrate how the shooting method is used to solve a two-point nonlinear boundary-value problem. The remaining problem conditions are as specified in Example 24.2: L = 10 m, h^′ = 0.05 m^{−2}, T_∞ = 200 K, T (0) = 300 K, \text{ and } T (10) = 400 K.
Just as with the linear ODE, the nonlinear second-order equation is first expressed as two first-order ODEs:
\frac{dT}{dx}= z
\frac{dz}{dx} = − 0.05(200 − T ) − 2.7 × 10^{−9}(1.6 × 10^9 − T^4 )
An M-file can be developed to compute the right-hand sides of these equations:
function dy = dydxn(x,y)
dy = [y(2);−0.05*(200 − y(1)) − 2.7e − 9*(1.6e9 − y(1)^4)];
Next, we can build a function to hold the residual that we will try to drive to zero as
function r = res(za)
[x,y] = ode45(@dydxn,[0 10],[300 za]);
r = y(length(x),1) − 400;
Notice how we use the ode45 function to solve the two ODEs to generate the temperature at the rod’s end: y(length(x),1). We can then find the root with the fzero function:
>> fzero(@res,− 50)
ans =
− 41.7434
Thus, we see that if we set the initial trajectory z(0) = − 41.7434, the residual function will be driven to zero and the temperature boundary condition T (10) = 400 at the end of the rod should be satisfied. This can be verified by generating the entire solution and plotting the temperatures versus x:
>> [x,y] = ode45(@ dydxn,[0 10],[300 fzero(@ res,− 50)]);
>> plot(x,y(:,1))
The result is shown in Fig. 24.7 along with the original linear case from Example 24.2.
As expected, the nonlinear case is depressed lower than the linear model due to the additional heat lost to the surrounding gas by radiation.
