Feladat¶
Oldja meg az id??f??gg?? h??vezet??s egyenlet??t n??gysz??gletes s??klemezre a k??vetkez?? hat??rfelt??telekkel:
- baloldalt adott a h??m??rs??klet, az ??l ment??n line??ris a h??m??rs??kletprofil
- alul adott konstans a h????ram
- a m??sik k??t ??l ment??n t??k??letes a h??szigetel??s.Hővezetés egyenlete¶
A feladat elv??gz??s??hez el??sz??r meghat??rozom azt a parci??lis differenci??l egyenletet, melynek megold??sak??nt mekaphat?? a h??terjed??s k??t dimenzi??ban.
Az ??ltal??nos h??vezet??si egyenlet egy dimenzi??ban a k??vetkez??:
\begin{equation*} \frac{\partial u}{\partial t} = \alpha \frac{\partial^{2} u}{\partial x^{2}}, \end{equation*} ahol $\alpha=k/ (\rho c_p) $ a h??vezet??si egy??tthat??, $??$ a s??r??s??g, $c_p$ a h??kapacit??s, $k$ pedig a termikus vezet??k??pess??g.
Ehhez nagyon hasonl??an meg tudjuk adni a h??vezet??si egyenletet k??t dimenzi??ban is, ahol megjelenik m??g egy t??rkoodin??ta:
\begin{equation*} \frac{\partial u}{\partial t} = \alpha (\frac{\partial^{2} u}{\partial x^{2}} + \frac{\partial^{2} u}{\partial y^{2}}), \end{equation*}
A differenciál egyenlet numerikus megoldása¶
A h??vezet??st le??r?? parci??lis differenci??legyenletet megold?? egyik legn??pszer??bb m??dszer az ??gynevezett Forward-Time Central-Space (FTCS). Az id?? v??lzoz??ban az Euler-f??le el??re l??ptetett differenciah??nyados l??p a deriv??lt hely??be, m??g a t??r v??ltoz??k eset??ben egy k??zponti differenciah??nyadossal dolgozunk a deriv??lt helyett.
Az egydimenzi??s h??vezet??si egyenlet a k??vetkez?? alakot veszi fel diszkretiz??lva:
\begin{equation*} \frac{u_i^{n+1} - u_i^{n} }{\Delta t} = \frac{\alpha}{\Delta x^{2}}(u_{i+1}^{n} - 2u_i^{n}+u_{i-1}^{n} ) \end{equation*}
A fenti egyenletben, az $u$ f??ggv??ny fenti indexe jel??li az id??l??p??st, az als?? index pedig az adott r??cspontban vett ??rt??ket.
Az egy dimenzi??s egyelethez hasonl??an fel lehet ??rni k??t dimenzi??ra is az egyenletet:
\begin{equation*} \frac{u_{i,j}^{n+1} - u_{i,j}^{n} }{\Delta t} = \frac{\alpha}{\Delta x^{2}} (u_{i+1,j}^{n} - 2u_{i,j}^{n}+u_{i-1,j}^{n}) + \frac{\alpha}{\Delta y^{2}} (u_{i,j+1}^{n} - 2u_{i,j}^{n}+u_{i,j-1}^{n}), \end{equation*}
ahol $i$ az $x$ beli $j$ pedig az $y$ beli ir??nyt jel??li.
\begin{equation*} u_{i,j}^{n+1} = u_{i,j}^{n} + \frac{\Delta t \alpha}{\Delta x^{2}} (u_{i+1,j}^{n} - 2u_{i,j}^{n}+u_{i-1,j}^{n}) + \frac{ \Delta t \alpha}{\Delta y^{2}} (u_{i,j+1}^{n} - 2u_{i,j}^{n}+u_{i,j-1}^{n}), \end{equation*}
Bevezetve:
\begin{equation*} r =\frac{\Delta t \alpha}{\Delta x^{2}} = \frac{\Delta t \alpha}{\Delta y^{2}} \end{equation*}
\begin{equation*} u_{i,j}^{n+1} = u_{i,j}^{n} + r (u_{i+1,j}^{n} - 2u_{i,j}^{n}+u_{i-1,j}^{n} + u_{i,j+1}^{n} - 2u_{i,j}^{n}+u_{i,j-1}^{n}), \end{equation*}
Megkapjuk a legegyszer??bb alakot amit m??r nagyon egyszer??en leht leprogramozni:
\begin{equation*} u_{i,j}^{n+1} = u_{i,j}^{n} + r (u_{i+1,j}^{n} +u_{i-1,j}^{n} + u_{i,j+1}^{n} +u_{i,j-1}^{n} - 4u_{i,j}^{n}), \end{equation*}
A stabilit??s felt??tele:
\begin{equation*} \Delta t \le \frac{\Delta x^{2} \Delta y^{2}}{2 \alpha (\Delta x^{2}+ \Delta y^{2})} \end{equation*}
illetve $\Delta x^{2} =\Delta y^{2}$ eset??n:
\begin{equation*} \Delta t \le \frac{\Delta x^{2} }{4 \alpha} \end{equation*}
#3D heat equation
import numpy as np
from numpy import pi,exp,sqrt
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d
import matplotlib.animation as animation
from matplotlib.animation import FuncAnimation
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
import numpy as np
%matplotlib inline
#parameters
dt=0.1
dx=0.1
L=50
B=50
c=329
spec_heat=23
density=104
Lx=48
Ly=48
dx=Lx/Nx
dy=Ly/Ny
alpha=c/(spec_heat*density)
i=1
Gr=np.eye(1)*300
Gr2=np.eye(1)*500*(i)+j
def assert_heaters(M, Gr,Gr2):
M[0:50,0:1] = np.where(Gr2 > 0, Gr2, M[0:50,0:1])
M[49:50,0:50] = np.where(Gr > 0, Gr, M[49:50,0:50])
avgT=[]
M=np.zeros([L,B])
assert_heaters(M, Gr,Gr2)
Stabilitás feltétele:¶
if dt<= 1/(2*alpha*((1/dx**2)+(1/dy**2))):
print('Stabil')
else:
print('Nem stabil')
T = np.arange(0,20,dt)
MM = []
M2=np.zeros([L,B])
for i in range(len(T)):
Gr2=np.eye(1)*200*(i)
for j in range(1,L-1):
for i in range(1,B-1):
M[i,j] = M[i,j] + (M[i-1,j] + M[i+1,j] + M[i,j-1] + M[i,j+1] -4*M[i,j])/(dx**2+dy**2)
MM.append(M.copy())
avgT.append(np.mean(M))
len(MM)
plt.imshow(MM[190],cmap=plt.cm.plasma) plt.colorbar() plt.show()
for q in range(len(MM)): plt.imshow(MM[q],cmap=plt.cm.plasma) plt.colorbar()
plt.savefig('gif/'+str(q)+'.png')
plt.show()gif előállítása¶
import imageio
import glob
import re
def atof(text):
try:
retval = float(text)
except ValueError:
retval = text
return retval
def natural_keys(text):
return [ atof(c) for c in re.split(r'[+-]?([0-9]+(?:[.][0-9]*)?|[.][0-9]+)', text) ]
images = []
filenames=glob.glob("gif/*.png")
filenames.sort(key=natural_keys)
for filename in filenames:
images.append(imageio.imread(str(filename)))
imageio.mimsave('gif/heat_gif.gif', images)
Eredmények¶
A szimul??ci?? futtat??sa sor??n minden l??p??sben elk??sz??tett heat-mapokat kimentettem png form??tumban egy mapp??ba ??s a kimentett heatmapokat ell??ttam egy h??m??rs??kleti ??rt??keket tartalmaz?? sz??nsk??l??val. Ezeket a k??l??nb??z?? id??tartamokhoz tartoz?? k??peket egy gif-be rendeztem egy python script seg??ts??g??vel. Az interakt??v vizualiz??ci??t ??s a t??bbi egy??b ??br??t szint??n a jupyter notebook-ban k??sz??tettem el.
Képkockák a szimuláciból¶
plt.imshow(MM[30],cmap=plt.cm.plasma)
plt.colorbar()
plt.show()
plt.imshow(MM[180],cmap=plt.cm.plasma)
plt.colorbar()
plt.show()
Négyzetlap időbeli hővezetése:¶
print('Maxim??lis homerseklet: '+str(np.round(max(avgT),3))+'C')
print('Minim??lis homerseklet: '+str(np.round(min(avgT),3))+'C')
from IPython.display import HTML
HTML('<img src="gif/heat_gif.gif">')
Átlaghőmérséklet időbeli változása:¶
plt.plot(avgT)
plt.xlabel('Iterations',size=15)
plt.ylabel('Average Temperature',size=15)
plt.grid(True)
x=np.linspace(0,50,len(M[49]))
b=plt.plot(x,M)
plt.grid(True)
plt.xlabel('x [dx]')
plt.ylabel('Temp [T]')
