3
$\begingroup$

So I have a PDE that I use to describe how material flows through a volume(2D or 3D). $$\frac{\partial C}{\partial t} + \vec{u} \cdot \nabla C = (D' + D )\nabla^2C$$ Now using finite differences I get to the following: $$C^{i+1}_{j,k} =C^i_{j,k}+\Delta t \biggr[ (D+D')\biggr[\frac{C^i_{j+1,k}-2C^i_{j,k}+C^i_{j-1,k}}{(\Delta x)^2}+\frac{C^i_{j,k+1}-2C^i_{j,k}+C^i_{j,k-1}}{(\Delta y)^2}\biggr]-u_x \frac{C^i_{j,k} -C^i_{j-1,k}}{\Delta x} - u_y \frac{C^i_{j,k} -C^i_{j,k-1}}{\Delta y}\biggr]$$ where D is diffusion coefficient and superscript is temporal, subscript is spatial. (Just imagine a grid). We are given $u_x$ and $u_y$ through some already validated code. Now, we use a boundary condition where $C=0$ at the edges. I know this is a silly boundary condition but it's irrelevant because the problem I'm about to show you persists even if we use small times.(where the concentration doesn't reach the edge)

enter image description here

[Look in the comments...It's telling me I can't post 2 links... Where do I have the second?]

here is the code. It's relatively simple.

import numpy as np 
import matplotlib.pyplot as plt
import time
%matplotlib inline
from scipy import ndimage
import scipy.io as sio
a =sio.loadmat('C:\\Users\\Tsila elbag\\Documents\\IPythonNotes\\New folder\\Cavity10k_1.mat')
b = a['VelField']
N =50
M =300
ux = b[0][0][0]
uy = b[0][0][1]
#ux = np.zeros((50,50))
#uy = np.ones((50,50))
dx = 2/N
dy =2/N
D= .001
Dprime =0
dt = .0001
C = np.zeros((N,N))
#C[np.floor((.25+1)/dx):np.floor((.75+1)/dx),np.floor((.25+1)/dy)\
#  :np.floor((.75+1)/dy)]=1
C[np.floor((.25+1)/dy)\
  :np.floor((.75+1)/dy),.1/dx:np.floor((.25)/dx)]=1
f, axarr = plt.subplots(2, sharex=True)
axarr[0].contourf(C)
Chist= C
Ctemp = C
for T in range(M):
    for j in range(N):
        for k in range(N):
            if k==49  or j ==49 or k ==0 or j ==0:
                Ctemp[j,k]=0
            else:

                Ctemp[j,k]= C[j,k]+dt*((D+Dprime)*(((C[j+1,k]-2*C[j,k]+\
                C[j-1,k])/(dx**2))+((C[j,k+1]-2*C[j,k]+C[j,k-1])/(dy**2)))-\
                ux[j,k]*(C[j,k]-C[j-1,k])/(dx)-uy[j,k]*(C[j,k]-C[j,k-1])/dy)
    #Chist.append(list(C))
    C = Ctemp    
axarr[0].set_title('Original') axarr[2].contourf(C)
axarr[2].set_title('{} seconds later'.format(dt*M))

So the question is why do I have these waves and if anyone knows effective boundary conditions I would be grateful. Thanks in advance.

$\endgroup$
  • $\begingroup$ I can't add the picture until I'm 10 reputation. If you up-vote twice. I will put the image. $\endgroup$ – Aviv Moshe Albeg Aug 25 '16 at 19:47
  • 1
    $\begingroup$ With a picture of the waves, it's possible to give a physical explanation to it -- but odds are good it's a numerical problem. We will have to wait and see what the problem actually is. $\endgroup$ – tpg2114 Aug 25 '16 at 20:41
  • 1
    $\begingroup$ I have 2 possible explanations. (a) You're writing dy=2/N, but bot 2 and N are integers, so dy might be 0 in that case. Try something like dy=2.0/N. The same goes for dx. (b) The problem might also be with your input velocity field ux & uy. Try (a), if that doesn't work, take a look at the velocity field. $\endgroup$ – Hayk Hakobyan Aug 25 '16 at 23:09
  • 4
    $\begingroup$ I'm voting to close this question as off-topic because it is about debugging specific code. $\endgroup$ – ACuriousMind Aug 26 '16 at 0:02
  • 3
    $\begingroup$ Are you using a value of $\frac{u\Delta t}{\Delta x}>1$? $\endgroup$ – Chester Miller Aug 26 '16 at 3:20

Your Answer

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

Browse other questions tagged or ask your own question.