1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
| #!/usr/env python3
# -*- coding: utf-8 -*-
import argparse
parser = argparse.ArgumentParser()
parser.add_argument('--data', action='store_true')
parser.add_argument('--einstein', action='store_true')
parser.add_argument('--debye', action='store_true')
import imgen
args = imgen.setup(parser)
import matplotlib.pyplot as plt
import numpy as np
import scipy.integrate
from scipy.optimize import curve_fit
if args.data:
plotE, plotD = False, False
elif args.einstein and not args.debye:
plotE, plotD = True, False
elif not args.einstein and args.debye:
plotE, plotD = False, True
else:
plotE, plotD = True, True
Tspace = np.linspace(0.1, 400, 300)
def modEinstein(T, Te, Cmax):
w = Te/T
return Cmax * w**2 / (4*(np.sinh(w/2)**2))
def modDebye(T, Td, Cmax):
f = lambda x: x**4 * np.exp(x) / (np.exp(x) - 1)**2
i = np.array([])
for t in T:
i = np.append(i, scipy.integrate.quad(f, 0, Td/t)[0])
return Cmax * 9*(T/Td)**3 * i
# Capacité thermiques massiques (J/kg/K)
# https://ws680.nist.gov/publication/get_pdf.cfm?pub_id=913059
Texp = np.array([ 4, 6, 8, 10, 12, 14, 16, 18, 20, 30, 40, 50, 60, 70, 80, 90, 100, 120, 140, 160, 180, 200, 220, 240, 260, 280, 300, 400])
Cexp_Cu_mass = np.array([0.09942, 0.2303, 0.4639, 0.8558, 1.470, 2.375, 3.640, 5.327, 7.491, 26.40, 57.63, 95.84, 135.2, 171.8, 203.8, 230.9, 253.5, 287.6, 311.6, 329.4, 343.4, 355.0, 364.7, 372.6, 378.6, 382.5, 384.0, 397.9])
Cexp_In_mass = np.array([ 0.9463, 3.505, 8.387, 15.28, 23.56, 32.62, 42.06, 51.61, 61.08, 104.5, 138.6, 163.3, 179.9, 190.5, 196.9, 200.7, 203.1, 206.0, 209.1, 213.5, 225.0, 219.0, 225.0, 230.6, 234.5, 235.7, 233.1, 233.0])
# Capacités thermiques molaires (J/mol/K)
Cexp_Cu = 63.546e-3 * Cexp_Cu_mass
Cexp_In = 114.818e-3 * Cexp_In_mass
# Incertitudes sur Texp
Terr = np.ones_like(Texp) # +/- 1K (source)
for data in [('Cu', Cexp_Cu), ('In', Cexp_In)]:
# Get data and error bars
material = data[0]
Cexp = data[1]
Cerr = Cexp * 0.05 # Incertitude sur Cexp: +/- 5% (source)
### Einstein
popt, pcov = curve_fit(modEinstein, Texp, Cexp, bounds=[(0, 0), (np.inf, np.inf)])
Te = popt[0]
# print(np.sqrt(np.diag(pcov))) # Écarts types de Tmax, Cmax
Ceinst = modEinstein(Tspace, *popt)
### Debye
popt, pcov = curve_fit(modDebye, Texp, Cexp, bounds=[(0), (np.inf)])
Td = popt[0]
Cdebye = modDebye(Tspace, *popt)
# Plot data
plt.errorbar(Texp, Cexp, label=material, linestyle='None', marker='+')
col = plt.gca().get_lines()[-1].get_color() # Récupère la couleur
# Plot models
if plotE:
plt.plot(Tspace, Ceinst, color=col, label=r'$T_E={:3.0f}K$'.format(Te), linestyle="dotted")
if plotD:
plt.plot(Tspace, Cdebye, color=col, label=r'$T_D={:3.0f}K$'.format(Td))
plt.axhline(y=3*8.314, color='k', linestyle='dashed', label=r'$3 k_B N_A$')
plt.xlabel(r'$T$ / K')
plt.ylabel(r'$c_P$ / J/mol/K')
plt.legend()
# Save
imgen.done(__file__)
|