0
$\begingroup$

Thank you for looking at this problem. Is this the correct approach to include neumann boundary conditions? With this solution temperature is not correct, and there´s no diffusion. The model seems unbalanced. However if the area is included in km2 instead of m2, temperature is ok and there is diffusion...

I´m trying to aplly the following 1D heat diffusion equation to a lake:

$$ \frac{∂T}{∂t}= \frac{1}{A(z)}\ \frac{∂}{∂z}\ ((A(z)K(z)\frac{∂T}{∂z}) $$

Where: k(z)= Eddy diffusion coefficient, 1.43x10-3 (m2/s), A(z)= lake area at depth z, m2

The surface boundary condition is the following:

$$ ρcK(z,t)\frac{∂T}{∂z}= -Q $$

Where: ρ= water density, 1000 kg/m3 c= water specific heat, 4186 J/kgºC Q= radiation, w/m2 (or J.s-1.m-2). The lake bottom condition is the following:

$$ \frac{∂T}{∂z}= 0 $$

I´m considering the Crank-Nicolson scheme with centered differences in space and time and I have used Thomas Algorithm to solve the system of equations.

The equation discretization is included in the picture below. But to simplify the analysis I´m including here the equations for each coefficient (f = k(z); n is time and i is space):

$$ α_{i}=\frac{1}{A_{i,n+1}+A_{i,n}} $$ $$ s=\frac{∆t}{2(∆z)^{2}} $$

$$a_{i} = −α_{i}s(f_{i−1,n+1}+ f_{i,n+1})$$

$$b_{i} = 1+α_{i}s(f_{i−1,n+1}+ 2f_{i,n+1}+f_{i+1,n+1})$$

$$c_{i} = −α_{i}s(f_{i,n+1}+ f_{i+1,n+1})$$

$$d_{i,n} = α_{i}s(f_{i,n}+ f_{i-1,n})T_{i−1,n}+α_{i}s(f_{i+1,n}+ f_{i,n})T_{i+1,n}+[1−α_{i}s(f_{i−1,n}+ 2f_{i,n}+f_{i+1,n})]T_{i,n}$$

for the surface boundary:

$$b_{0}=1+α_{0}s(f_{−1,n+1}+ 2f_{0,n+1}+f_{1,n+1})$$

$$c_{0}=−α_{0}s(f_{-1,n+1}+ 2f_{0,n+1}+f_{1,n+1})$$

$$d_{0,n}=[1-α_{0}s(f_{-1,n}+ f_{0,n}+f_{1,n}]T_{0,n}+α_{0}s(f_{-1,n}+2 f_{0,n}+f_{1,n})T_{1,n}+α_{0}s(f_{−1,n}+ 2f_{0,n})\frac{2Q_{n}∆z}{ρcK_{0,n}}+α_{0}s(f_{-1,n+1}+ f_{0,n+1})\frac{2Q_{n+1}∆z}{ρcK_{0,n+1}}$$

for the bottom boundary:

$$a_{NZ} = −α_{NZ}s(f_{NZ−1,n+1}+ 2f_{NZ,n+1}+f_{NZ+1,n+1})$$

$$b_{NZ} = 1+α_{NZ}s(f_{NZ−1,n+1}+ 2f_{NZ,n+1}+f_{NZ+1,n+1})$$

$$d_{NZ,n} = α_{NZ}sT_{NZ−1,n}(f_{NZ-1,n}+ 2f_{NZ,n}+f_{NZ+1,n})+[1−α_{NZ}s(f_{NZ−1,n}+ 2f_{NZ,n}+f_{NZ+1,n})]T_{NZ,n}$$

I have also include a simplification of the Python code. I have removed all the loops and I have considered a constante value for eddy diffusivity K(z). It´s a daily simulation (N=10) for a lake with 6 layers. Each layer is 2m deep (z=2m).

import numpy as np
np.set_printoptions(threshold=np.inf)
import pylab as plt



#------------------------------------------------------------------------------------------------------------------------------------
# Tri Diagonal Matrix Algorithm (Thomas algorithm) solver (source:https://gist.github.com/cbellei/8ab3ab8551b8dfc8b081c518ccd9ada9)
#-------------------------------------------------------------------------------------------------------------------------------------

def TDMAsolver(a, b, c, d):

    nf = len(d)
    ac, bc, cc, dc = map(np.array, (a, b, c, d))
    for it in xrange(1, nf):
        mc = ac[it-1]/bc[it-1]
        bc[it] = bc[it] - mc*cc[it-1] 
        dc[it] = dc[it] - mc*dc[it-1]

    xc = bc
    xc[-1] = dc[-1]/bc[-1]

    for il in xrange(nf-2, -1, -1):
        xc[il] = (dc[il]-cc[il]*xc[il+1])/bc[il]

    del bc, cc, dc

    return xc

#-------------------------------------------------------------------------------------------------------------------------------------
#Crank Nicolson scheme
#-------------------------------------------------------------------------------------------------------------------------------------
#Simulation time period (10 days)
N=10

Nlayers=6

#Space, m
DZ=2.0

#Radiation, J.s-1.m-2 (or W/m2)
Q=300.0

#Eddy diffusion coefficient at 25 C (m2/s)
K=1.43*10**-3

#Water temperature array
TempLake=np.zeros((N+1,Nlayers))

# Initial condition, C
TempLake[0]=[25,25,25,25,25,25]

#Lake surface area, m2
A1=23992000.0 #m2
A2=23992000.0 #m2

alfa=1.0/(A1+A2)

#time interval, seconds, s
DT=86400
#Water density, kg/m3
DH2O=1000.0
#Water specific Heat j/kg C
Cp=4186

s=DT/(2*DZ**2.0)

#Initializing coefficients np.arrays
a=np.zeros(Nlayers-1)
b=np.zeros(Nlayers)
c=np.zeros(Nlayers-1)
d=np.zeros(Nlayers)

for i in xrange(1,N+1):  # 1,2,3,4,5,6,7,8,9,10        
    #Surface
    #b0, c0, d0
    b[0]=1+alfa*s*(K+2.0*K+K)
    c[0]=-alfa*s*(K+2.0*K+K)
    d[0]=((1.0-alfa*s*(K+2.0*K+K))*TempLake[i-1][0])+(alfa*s*(K+2.0*K+K)*TempLake[i-1][1])+alfa*s*(K+K)*((2*Q*DZ)/(Cp*DH2O*K))+alfa*s*(K+K)*((2.0*Q*DZ)/(DH2O*Cp*K)) # 0,1,2,3,4,5,6,7,8,9  

    #Bottom
    #an, bn, dn
    a[-1]=-alfa*s*(K+2.0*K+K)
    b[-1]=1+alfa*s*(K+2.0*K+K)
    d[-1]=alfa*s*TempLake[i-1][-2]*(K+2.0*K+K)+(1-alfa*s*(K+2*K+K))*TempLake[i-1][-1] # 0,1,2,3,4,5,6,7,8,9 

    #ai, bi, ci, di
    a[0]=-alfa*s*(K+K)
    a[1]=-alfa*s*(K+K)
    a[2]=-alfa*s*(K+K)
    a[3]=-alfa*s*(K+K)

    b[1]=1.0+alfa*s*(K+2.0*K+K)
    b[2]=1.0+alfa*s*(K+2.0*K+K)  
    b[3]=1.0+alfa*s*(K+2.0*K+K)
    b[4]=1.0+alfa*s*(K+2.0*K+K)

    c[1]=-alfa*s*(K+K)
    c[2]=-alfa*s*(K+K)
    c[3]=-alfa*s*(K+K)
    c[4]=-alfa*s*(K+K)

    d[1]=alfa*s*(K+K)*TempLake[i-1][0]+alfa*s*(K+K)*TempLake[i-1][2]+(1-alfa*s*(K+2.0*K+K))*TempLake[i-1][1]
    d[2]=alfa*s*(K+K)*TempLake[i-1][1]+alfa*s*(K+K)*TempLake[i-1][3]+(1-alfa*s*(K+2.0*K+K))*TempLake[i-1][2]
    d[3]=alfa*s*(K+K)*TempLake[i-1][2]+alfa*s*(K+K)*TempLake[i-1][4]+(1-alfa*s*(K+2.0*K+K))*TempLake[i-1][3]
    d[4]=alfa*s*(K+K)*TempLake[i-1][3]+alfa*s*(K+K)*TempLake[i-1][5]+(1-alfa*s*(K+2.0*K+K))*TempLake[i-1][4]

    TempTDMA=TDMAsolver(a, b, c, d)

    TempLake[i]=TempTDMA

print TempLake

Z=TempLake 

X,Y=np.meshgrid(range(Z.shape[0]+1),range(Z.shape[1]+1)) 

im = plt.pcolormesh(X,Y,Z.transpose(), cmap='jet') 
ax = plt.gca()
ax.set_ylim(ax.get_ylim()[::-1])
plt.colorbar(im, orientation='horizontal') 

plt.show()

enter image description here

Results:

TEST 1 : simulation code included above

TEST 2 : simulation equal to Test 1 but lake area is included in km2 instead of m2

enter image description here

enter image description here

$\endgroup$
1
$\begingroup$

The choice of the boundary conditions is dependent on the physics you want to model. If you're doing heat diffusion through a lake, and you want constant temp on top you'll have a different boundary condition than for constant heat flux out of the top (which appears to be what you want). I'm assuming that you want the bottom of the lake to have no heat flux (as that's what you wrote). as far as I can tell from your break down you have no mistakes that jump out at me, but I didn't check too closely. When you say the temperature is not correct, what do you mean? Thanks.

$\endgroup$
  • $\begingroup$ Thank you for your help. If the boundary conditions are ok, can the problem be related with the value of alpha=1/(Ai,+ Ai) and/or, s=deltaT/2(deltaz)^2. I have included the results I have. As you can see Test 1 results show only a very small increase in lake surface water. However, when lake area is considerable reduced (included in km2 instead of m2), results seem to be ok. In a daily simulation (delta T=86400 seconds) the value of s is very large, 5400...can this be the problem?. Thank you $\endgroup$ – M. A. Feb 4 at 15:56
  • $\begingroup$ Having looked at it again, one thing that jumps out to me is that the equation you describe is heat conduction. Are you sure you discretized this PDE correctly? Also, what are your plots showing? Please label the axes. $\endgroup$ – EMP Feb 4 at 17:04
  • $\begingroup$ I think so, I have checked this n times, but I´m not 100% sure. I have included a picture with the discretization... $\endgroup$ – M. A. Feb 4 at 17:17
  • $\begingroup$ I found the problem. The code wasn't considering that f(z)=A(z)*K(z)...Thank you so much, for helping me to exclude the definition of the boundary conditions $\endgroup$ – M. A. Feb 5 at 11:54
  • $\begingroup$ Happy to help and to hear it works! $\endgroup$ – EMP Feb 5 at 16:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.