1. c1 = [1]*100 + [0]*100
2. nc1=c1
3. DtildeList=[]
4. DCuList=[]
5. DNiList=[]
6. DCuS=2.03
7. DtildeS=1.45
8. DNiS=2
9.  
10.  
11. dt = 100
12. dx=0.01
13. xlist=np.linspace(-1,1,200)
14.  
15.  
16. def DtildeL(X,DA,DB):
17.     return X*DB+(1-X)*DA
18.  
19.  
20.    
21.  
22. Xlist=np.linspace(0,1,200)
23. for i in range(200):
24.     DCuList.append(10**(-12.92-Xlist[i]*DCuS))
25.     DNiList.append(10**(-13.5-Xlist[i]*DNiS))
26.     DtildeY=DtildeL(Xlist[i], DNiList[i],DCuList[i])
27.     DtildeList.append(DtildeY)
28.  
29. def DtildeFunc(X):
30.     return DtildeL(X,10**(-12.92-(1-X)*DCuS),10**(-13.5-(1-X)*DNiS)) #mm^2/sec
31.    
32. TlistN=np.linspace(973,1273,1080)
33. for k in range(10800):
34.    
35.     for i in range(1, len(c1)): #This loop calculates the numerical solution for the concentration
36.         if i== len(c1)-1:
37.             nc1[i] =c1[i] +(dt * DtildeFunc(c1[i])*10**6)/(dx**2) * (c1[i] - 2 * c1[i] + c1[i-1])+((DtildeFunc(c1[i])*10**6 - DtildeFunc(c1[i-1])*10**6)*(c1[i]-c1[i-1])*(1/(dx**2)))
38.         else:
39.             nc1[i] = c1[i]+(dt * DtildeFunc(c1[i])*10**6)/(dx**2) * (c1[i+1] - 2 * c1[i] + c1[i-1])+((DtildeFunc(c1[i])*10**6 - DtildeFunc(c1[i-1])*10**6)*(c1[i]-c1[i-1])*(1/(dx**2)))
40.     c1=nc1.copy()
41.     if k==1080:
42.         c30=c1
43.     if k==3600:
44.         c100=c1
45.     if k==7200:
46.         c200=c1
47.  
48. plt.figure(figsize=(12,8))
49. plt.plot(xlist,c30, label="Diffusion profile @30h")
50. plt.plot(xlist,c100, label="Diffusion profile @100h")
51. plt.plot(xlist,c200, label="Diffusion profile @200h")
52. plt.plot(xlist,c1, label="Diffusion profile @300h")
53. plt.legend()
54. plt.ylabel('Concentration of Cu',fontsize=12)
55. plt.xlabel('position [mm]', fontsize = 12)
56. plt.title(r'Diffusion profile with use of $\tilde{D}$')
57. plt.grid()
58.  
59. ;