Search All of the Math Forum:

Views expressed in these public forums are not endorsed by NCTM or The Math Forum.

Notice: We are no longer accepting new posts, but the forums will continue to be readable.

Topic: Problem with heat equation integration
Replies: 1   Last Post: Dec 19, 2012 10:36 AM

 Messages: [ Previous | Next ]
 Torsten Hennig Posts: 2,419 Registered: 12/6/04
Re: Problem with heat equation integration
Posted: Dec 19, 2012 10:36 AM

> Hi everyone, (first excuse me for my english)
> I have to integrate numerically heat equation in
> subsurface, that is:
>
> C*dT/dt=d(lambda*dT/dz)/dz + H
>
> where lambda is thermal conductivity, H volumetric
> heat production and C is volumetric heat capacity.
> These magnitudes are z-dependent. Boundary condition
> at the bottom (zN) is:
>
> dT/dz=-qb, where qb is basal heat flow
>
> On the other hand, I have to impose T(z1) varying
> like a sinusoid (annual wave).
>
> So, I have to use theta method to integrate in time.
> First I solve the stationary problem and I use the
> solution like initial condition in time integration.
> Stationary problem is right but the time stepping is
> s wrong, because I obtain same solution for all the
> times. Nevertheless, I put all the code although the
> important part time stepping at the end
>
> clear all; close all; clc;
>
> % =============== Stationary problem
> =====================================
>
> nz_1=99;% Number of nodes - 1
> dt=0;
> temp_med_top=15+273.15; % Mean temperature at
> surface
> w=(2*pi)/(365*24*60*60); % Annual wave pulsation
> ampli=5;% Wave amplitude
> qb=30e-3;% Basal heat flow;
> top=0; % Up boundary
> bot=100; % Low boundary
>
> % At first, define the position of the nodes
>
> z=logspace(0,log10(bot),nz_1+1); %I add 1 more
> associated with the ghost cell
> z=[0,z]; % Add a node at 0m
>
> nz=length(z)-1;% nz=number of nodes
>
> % Calculate node spacing --> Cell longitude
>
> dz=-diff(z);
> nc=length(dz);
> dz_N1=dz(nc); % Cell longitude of the ghost cell
> dz=dz(1:nc-1);
> nc=length(dz);
>
> % Define lambda at the center of the cells. I suppose
> it increases lineary
> % with depth from 1 to 6
>
> lambda=linspace(1,6,nc+1); % 1 more value associated
> with the ghost cell
> lambda_N1=lambda(nc+1);
> lambda=lambda(1:nc);
>
> % Calculate L matrix
>
> dc=0.5*(dz(2:nc)+dz(1:nc-1));
> dl=lambda./dz;
> % Initialize a,b,c vectors
> a=zeros(1,nz);
> b=zeros(1,nz);
> c=zeros(1,nz);
> c(3:nz) =dl(2:nc)./dc; % Upper diagonal
> a(1:nz-2)=dl(1:nc-1)./dc; % Lower diagonal
> b(2:nz-1)=-(a(1:nz-2)+c(3:nz));
> L=spdiags([a' b' c'],-1:1,nz,nz);
>
> % Impose boundary conditions
>
> L(1,1)=-1;
> L(nz,nz-1)=1/dz(nc);
> L(nz,nz)=-1/dz(nc);
>
> % Define C(z) y H(z) at the centers of the cells. I
> suppose that C or H value at
> % the cell i+1 is the value at the cell i plus a
> thousandth of the latter
> % (for example)
>
> C_cel=zeros(nc+1,1); % nc+1 for the ghost cell
> C_cel(1)=1e6;
>
> for i=2:nc+1
> C_cel(i)=C_cel(i-1)+C_cel(i-1)/1000;
> end
>
> C_cel_N1=C_cel(nc+1);
> C_cel=C_cel(1:nc);
>
> H_cel=zeros(nc,1);
>
> for j=2:nc
> H_cel(j)=H_cel(j-1)+H_cel(j-1)/1000;
> end
>
> % Interpolate to the nodes using weighted mean.
> Suppose larger cells have
> % larger weights
>
> C_nod=zeros(nz,1);
> C_nod(1)=C_cel(1); % Asign to 1st node C value of the
> next cell
> C_nod(nz)=(C_cel_N1*dz_N1+C_cel(nc)*dz(nc))/(dz_N1+dz(
> nc)); % Weighted nean between the last cell and the
> ghost cell
>
> for k=2:nz-1
>
>
>
> C_nod(k)=(C_cel(k-1)*dz(k-1)+C_cel(k)*dz(k))/(dz(k-1)
> +dz(k));
> end
>
> H_nod=zeros(nz,1);
>
> for l=2:nz-1
>
>
>
> H_nod(l)=(H_cel(l-1)*dz(l-1)+H_cel(l)*dz(l))/(dz(l-1)
> +dz(l));
> end
>
> % Define C matrix
>
> C=spdiags(1./C_nod,0,nz,nz);
>
> % Define source term
>
> b=-H_nod;
>
> % Impose boundary conditions at b
>
> b(1)=-(temp_med_top+ampli*sin(w*dt));
> dt=0; % Suppose this is the initial state
> b(nz)=qb/(((lambda(nc)*dz(nc)+lambda_N1*dz_N1)/(dz(nc)
> +dz_N1)));
> % The denominator is the interpolation between the
> last cell and the ghost
> % cell
>
> % Multiply b by C
>
> b=C*b;
>
> % Multiply L by C
>
> L=C*L;
>
> % Solving the system
>
> T=L\b;
> Ti=T;
>
> %======================= Time stepping
> ===================================
>
> % Generate T matrix
>
> nt=1000; % number of time steps
> fac=365*24*60*60; % years to seconds factor
> paso=0.00001; % time step (years)
> dt=paso*ones(1,nt-1);% time step vector
> T=zeros(nz,nt);
> T(:,1)=Ti;
> T0=Ti;
> b0=-b;
> b1=-b; % I think b has the opposite sign regarding
> stationary case, but L
> % is exactly the same
> theta=0; % Suppose for example explicit integration
> I=speye(nz);
>
> for it=1:nt-1
>
>
>
>
> b0(1)=(temp_med_top+ampli*sin(w*(it-1)*dt(it)*fac))*C
> (1,1);
>
>
>
> b1(1)=(temp_med_top+ampli*sin(w*it*dt(it)*fac))*C(1,1
> );
> A=I-theta*dt(it)*L;
>
>
>
> p=(I+(1-theta)*dt(it)*L)*T0+dt(it)*(theta*b1+(1-theta
> )*b0);
>
> T1=A\p;
>
> T0=T1;
> T(:,it+1)=T1;
>
> end
>
> space=linspace(0,-bot,nz);
> time=linspace(0,paso*nt,nt);
> figure
> contourf(time,space,T)
>
>
> Thank you!

Use MATLAB's "pdepe" for your heat equation -
the risk of making coding errors is much lower.

Best wishes
Torsten.

Date Subject Author
12/19/12 Alex
12/19/12 Torsten Hennig