Geophysical Modelling - Math Problem Example

Geophysical modelling - practical V - solutions Problem ment and simplification In the most common case, the conductive heat transfer can be described analytically via the general conduction equation: , or (1) . (2) Here, the temperature () is a function of , , (), and (); () is the thermal conductivity, () is the density, () is the specific heat per unit mass, and () is the rate of internal energy conversion (so-called heat generation) per unit volume. This equation is a direct consequence of the first law of thermodynamics (Pitts & Sissom 1998, p. 30-31). The principal limitations of the general conduction equation (2) are: incompressibility of the heat conduct medium, and absence of convection (Lienhard & Lienhard 2003, p. 54-56). Therefore, we may use this equation for description of the heat transfer through soil as incompressible and non-convective medium (1st assumption). The heat conduction equation (2) can be simplified essentially by following assumptions. Firstly, we may declare the thermal conductivity of the soil layers as constant value (2nd assumption). Indeed, if heat transfer is plane and steady, and if is not large (in our case ), one can make a reasonably accurate approximation using a constant average value of (Lienhard &Lienhard 2003, p. 51). Then, we may declare the density and the specific heat of the soil layers as constants and (3rd assumption). It is necessary to note that such assumption is somewhat groundless, especially for the cases of strictly inhomogeneous soils. However, only after this we may introduce a constant diffusion coefficient () of the soil (Bird et al. 2002, p. 268). Finally, we may declare that the rate of internal energy conversion is negligibly small (4th assumption). In fact, such assumption is declaration of absence of the heat generation or consumption within the soil. Again we note that this assumption can be groundless for the inhomogeneous soils with stone inclusions, fluid- or air-filled interstices with internal convective flows. Moreover, the heat transfer in such porous and composite media is very difficult to analyze (Bird et al. 2002, p. 281-283). Only after all these assumptions we may derive Fourier equation , (3) and to use it in one-dimensional (5th assumption) simplified form: . (4) Equation (4) has well-known equilibrium or steady-state solution (5) for the steady boundary conditions and ; is the thickness of the soil upon the rock background. Solution (5) is easy to derive analytically (Haberman 1983, p. 13-14), so we will use it for checking our numerical model by approximation at . Heat transfer model parameters The simplified problem (4) is stated by following values: m, m2s-1. Boundary conditions are: , (1st case), and (2nd case). Initial values are stated by equation . Model discretization is stated by number of soil layers , their thicknesses m, and the timestep (in seconds) which we can modify (240s, or 550s). Model geometry is shown at the figure 1. In our model zero-level () is located at the surface of the rock background because the soil thickness is rather unstable parameter. Indeed, thickness of the real soil cover is a function and for the small areas only. Therefore, we will use more "stable" rock surface to count out -values of the soil layers. Model dynamics (i.e. heat conduction process) is described by (4) which is transformed in a form of difference equation ; (6) here, denotes ; , ; , . The work equation (6) is derived from (4) by FTCS scheme, when forward differentiation was used for and centered differentiation was used for (Boyce & DiPrima 2001, p. 419f). Figure 1 - Model geometry Solutions 1. Let us transform the FTCS scheme (6) into an explicit form: . (7) For the bottom soil layer () we have (8) because of . For the layer near the soil surface () we have (9) because of . 2. To create a Matlab script for solving equation (4) in the explicit FTCS form (7), we can use both initial script and examples of (Mathews & Fink 1999, p. 526-536). Work model code is in Appendix. For the first time step we obtain distribution (see fig. 2) which shows some heat diffusion down into the top () soil layer. The rest of the soil (i.e. the layers ) is at the initial temperature which equals to temperature of the rock bottom layer. Such situation is correct physically. Figure 2 - Temperature distribution for the first time step 3. After model running for 2 and 24 hours we obtain distributions (blue) and (green; see figure 3) which show gradual temperature increasing from the free soil surface to the rock bottom. In fact, we observe warming-up process which is caused by stable afflux (for boundary conditions stated for the case 1) of heat from the air, and diffusion of this heat through the soil medium. This process is provided by molecular mechanisms of the heat transfer (Bird et al. 2002, p. 265f). Here, we can examine our model by approximation at . Correct simulation must produce linear temperature distribution in a form of analytic equation (5). So, we need to obtain linear temperature increase from near the bottom soil level () to near the surface soil level (). Numerical model was tested by approximation at days. Computed temperature distribution is similar to analytic (see figure 4), so we can conclude that our model is correct dynamically. Difference between numeric solution and analytic approximation (i.e. steady-state solution) can be explained by too rough discretization of the soil in -direction, or too small time interval (24015=3600 iterations in total). Figure 3 - Temperature distributions within the soil after 2 (blue) and 24 (green) hours. Figure 4 - Temperature distributions within the soil after 10 days (green) and at (red). 4. If we will increase the time step from 240s to 550s, we obtain following results (see figure 5). Computations of temperature distributions after 1 (blue) and 2 (green) hours give unstable and oscillating "solutions". This is caused by existence of so-called "stability areas" for differential equations which are solved numerically (Boyce & DiPrima 2001, p. 445f). Such areas of numerical stability are strongly dependent on the time step values. Usually, increasing time steps destroy stability of the numerical processes for ODEs or PDEs solving. Figure 5 - Temperature distributions within the soil after 1 (blue) and 2 (green) hours are unstable numerically at s. 5. For the case of slowly oscillating boundary condition we obtain non-trivial temperature distribution (see figure 6; after 3 days). Computed distribution is almost homogenic. Indeed, the average daily value of is equal to . Little (less than ) inhomogeneity in temperature distribution is caused by oscillating character of the heat transfer from the air to the soil. So, we can observe cyclic heat transfer from the air to the soil and then to the rock bottom (by day), and from the heated soil to the cold air and the rock bottom (by night). Figure 6 - Temperature distribution within the soil after 3 days for the case of oscillating boundary condition . Bibliography Bird, RB, Stewart, WE & Lightfoot, EN 2002, Transport phenomena, 2nd edn, John Wiley & Sons, New York. Boyce, WE & DiPrima, RC 2001, Elementary differential equations and boundary value problems, 7th edn, John Wiley & Sons, New York. Haberman, R 1987, Elementary applied partial differential equations with Fourier series and boundary value problems, 2nd edn, Prentice Hall, Englewood Cliffs, New Jersey. Lienhard IV, JH & Lienhard V, JH 2003, A heat transfer textbook, 3rd edn, Phlogiston Press, Cambridge, Massachusetts. Mathews, JH & Fink, KD 1999, Numerical methods using MATLAB, 3rd edn, Prentice Hall, Upper Saddle River, New Jersey. Pitts, DR & Sissom, LE 1998, Theory and problems of heat transfer, McGraw-Hill Schaum's Outline Series, New York. Appendix - Model code % Solution of the 1D diffusion equation using FTCS scheme clc; clear; close all; % Set model parameters and constants dt=240; N=24*15; J=10; D=1e-7; L=0.1; dz=L/J; DD=D*dt/(dz*dz); % Set up z-grid and t-grid j=[1:J]; n=[1:N]; z(j)=(j-0.5)*dz; js=[2:J-1]; % Define initial conditions Tair=40; % Tair=10; % for the case #2 only Trock=10; Ti(j)=Trock; T(n,j)=NaN; % Calculate 1st timestep T(1,1)=DD*(Ti(2)-2*Ti(1)+Trock)+Ti(1); T(1,js)=DD*(Ti(js+1)-2*Ti(js)+Ti(js-1))+Ti(js); T(1,J)=DD*(Tair-2*Ti(J)+Ti(J-1))+Ti(J); % Calculate subsequent timesteps for n=2:N % t=n*dt; % for the case #2 only % Tair=10+10*sin(2*pi*t/(24*60*60)); % for the case #2 only T(n,1)=DD*(T(n-1,2)-2*T(n-1,1)+Trock)+T(n-1,1); T(n,js)=DD*(T(n-1,js+1)-2*T(n-1,js)+T(n-1,js-1))+T(n-1,js); T(n,J)=DD*(Tair-2*T(n-1,J)+T(n-1,J-1))+T(n-1,J); end % Show results plot(T(24*15,:),z,'g-o'); grid on; xlabel('Temperature (^oC)'); ylabel('Distance from the rock bottom (m)'); axis([0 40 0 0.1]); print -djpeg99 -r200 -zbuffer MyFigure Matlab ver. 7.0.1.24704 (R14) SP1 was used for executing this script on the Windows XP Professional SP2 PC platform. Read More