TMM4100 Materialteknikk
Eulers metode¶
Ficks 2. lov (1) for diffusjon og ligningen for varmeledning (2) i en dimensjon har samme format:
\begin{equation} \frac{\partial C}{\partial t}= D \frac{\partial^2 C}{\partial x^2} \tag{1} \end{equation}\begin{equation} \frac{\partial Q}{\partial t}= \alpha \frac{\partial^2 T}{\partial x^2} \tag{2} \end{equation}Ligningene er vanskelig (eller ikke mulig) å løses analytisk for annet enkle og veldefinerte grensebegingelser og initialbetingelser.
Forward Euler en enkel numerisk metode som er godt egnet for akkurat disse ligningen, og vil alltid gi en stabil løsning så lenge grensebetingelsene på begge sider av en lengde er kjent og tidssteget i iterasjonen er tilstrekkelig lite.
Forward Euler for løsning av 1-dimensjonalt diffusjonsproblem¶
Kort fortalt: ligning (1) kan formuleres med endelige inkrement for både tid og avstand:
\begin{equation} \frac{C_{i}^{n+1} - C_{i}^{n}}{\Delta t} = D \frac{C_{i+1}^{n}-2C_{i}^{n}+C_{i-1}^{n} }{\Delta x^2} \tag{3} \end{equation}Parameteren $F$ er definert som
\begin{equation} F = \frac{D \Delta t}{\Delta x^2} \tag{4} \end{equation}og denne har en viktig rolle med hensyn på stabilitet (må være <= 0.5 for dette tilfellet). Dermed kan vi skrive (3) som
\begin{equation} C_{i}^{n+1} = C_{i}^{n} + F(C_{i+1}^{n}-2C_{i}^{n}+C_{i-1}^{n}) \tag{5} \end{equation}Kode:
In [1]:
import numpy as np
def diffusjonForwardEuler(D,X,nx,t,Ca,Cb,C0,Cinit=[]):
'''
Input:
D : diffusjonskoeffisient,
X : tykkelse
nx : antall segment (oppdeling langs x)
t : tid for beregningen
Ca : konsentrasjon ved x=0
Cb : konsentrasjon ved x=X
C0: : initiell konsentrasjon når denne er uniform
Cinit : initiell konsentrasjonsfordeling. Må ha samme dimensjon som C
Output:
xs : x-koordinater, C : konsentrasjonen langs x
'''
F=0.5 # vilkår for stabilt tidsinkrement er at F<=1/2
dx=X/nx # lengde av segment
dt=(F*dx**2)/D # stabilt tidsinkrement gitt av (3)
nt = int(round(t/float(dt))) # antall tidssteg
ts = np.linspace(0, nt*dt, nt+1) # tider fra 0 til t
xs = np.linspace(0, X, nx+1) # x-koordinater gjennom tykkelsen
if len(Cinit)==0: # betyr at initiell konsentrasjon er uniform
C_n = C0*np.ones(nx+1) # initiell/forrige konsentrasjonsfordeling
else: # bruker heller en ikke-uniform konsentrasjonsfordeling
C_n = Cinit # initiell/forrige konsentrasjonsfordeling som en array blir brukt
C = np.zeros(nx+1) # nåværende/ny konsentrasjonsfordeling
for ti in range(0, nt): # FE iterasjoner:
for xi in range(1, nx): # over x-verdier
C[xi] = C_n[xi] + F*(C_n[xi+1] - 2*C_n[xi] + C_n[xi-1]) # ligning (4)
C[0] = Ca # grensebetingelse, x=0
C[nx] = Cb # grensebetingelse, x=X
C_n,C = C,C_n # switcher mellom forrige og nåværende
return xs, C
Eksempler:
In [5]:
import matplotlib.pyplot as plt
x, C = diffusjonForwardEuler(D=1E-9,X=0.01,nx=200,t=1000,Ca=1,Cb=0,C0=0)
fig,ax=plt.subplots(figsize=(6,4))
ax.plot(x,C,color='red',alpha=0.2)
ax.fill_between(x,C,0,color='red',alpha=0.1)
ax.set_title('Eksempel-1')
ax.set_xlabel('x [m]')
ax.set_xlim(0,x[-1])
ax.set_ylim(0,1)
ax.set_ylabel('Konsentrasjon')
#ax.grid(True)
plt.tight_layout()
plt.show()
In [7]:
x, C = diffusjonForwardEuler(D=1E-9,X=0.01,nx=200,t=1000,Ca=0.8,Cb=0.1,C0=0.4)
fig,ax=plt.subplots(figsize=(6,4))
ax.plot(x,C,color='red',alpha=0.2)
ax.fill_between(x,C,0,color='red',alpha=0.1)
ax.set_title('Eksempel-2')
ax.set_xlabel('x [m]')
ax.set_xlim(0,x[-1])
ax.set_ylim(0,1)
ax.set_ylabel('Konsentrasjon')
#ax.grid(True)
plt.tight_layout()
plt.show()
