Programmiererforum

Code: Select all

import numpy as np
from scipy.integrate import quad
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt

# Constants
kb = 1.380649e-23  # Boltzmann constant (J/K)
hbar = 1.0545718e-34  # Reduced Planck's constant (J*s)
theta_D = 280  # Debye temperature (K)
v = 2300  # Average phonon group velocity (m/s) -

# Define the inverse relaxation times
def tau_PD_inv(omega, A):
return A * omega**4

def tau_U_inv(omega, B, T):
return B * omega**2 * np.exp(-theta_D / (3 * T))

def tau_GB_inv(omega):
D = 0.7e-6  # Grain boundary characteristic length (m)
return v / D

def tau_total_inv(omega, A, B, T):
return tau_PD_inv(omega, A) + tau_U_inv(omega, B, T) + tau_GB_inv(omega)

# Debye Callaway model for thermal conductivity
def kappa_lattice(T, A, B):
def integrand(x):
omega = x * kb * T / hbar
tau_total = 1 / tau_total_inv(omega, A, B, T)
return x**4 * np.exp(x) / (np.exp(x) - 1)**2 * tau_total

integral, _ = quad(integrand, 0, theta_D / T)
prefactor = kb**4 * T**3 / (2 * np.pi**2 * v * hbar**3)
return prefactor * integral

# Vectorize the model for curve fitting
def kappa_lattice_vectorized(T, A, B ):
return np.array([kappa_lattice(Ti, A, B ) for Ti in T])

# Example thermal conductivity and temperature data for Ge WQ
T_data = np.array([3.0748730372823156, 3.1962737114707203, 3.4441443065638015, 3.785528026732028,
4.043999986888003, 4.524033859711473, 4.902329259251676, 5.270E+00,
5.902534836759174, 6.721645553370004, 7.690683028568884, 9.031E+00,
1.103E+01, 12.811754075061538, 16.07373842258417, 19.67675447459597,
23.951202854256202, 29.710517518825778, 3.684E+01, 45.030707655028536,
56.92469139605285, 6.929E+01, 85.62796640159988, 105.41803279520694,
129.5096325800293, 159.77670578646953, 200.9314386671018, 240.34130109997898,
291.63127697234717])  # Temperature in K
kappa_data = np.array([ 0.5021456742832973, 0.5907837911587943, 0.7331309897165834, 0.9570228747931655,
1.1281866396926783, 1.4904675008427906, 1.8682223847710366, 2.1856825810901013,
2.8515011626042894, 4.057791636291224, 5.328111104734368, 7.645627510566391,
11.045409412114584, 14.277055520040301, 19.54065663430087, 24.738891215975602,
28.776422188514683, 32.562485211039274, 33.37572711709385, 32.455084248250884,
28.56802589740088, 24.937473590597776, 20.887192910896978, 17.37073097716312,
14.783541083880793, 12.26010669954394, 9.991973525044475, 8.531955488387288,
7.347906481962733])  # Thermal conductivity in W/mK

# Initial guess for parameters A, B, and C
initial_guess = [4.75E-42, 1.054E-17]

# Curve fitting
popt, pcov = curve_fit(kappa_lattice_vectorized, T_data, kappa_data, p0=initial_guess, maxfev=8_000)
#popt, pcov = curve_fit(kappa_lattice_vectorized, T_data, kappa_data, p0=initial_guess, bounds=(0, np.inf))

# Extract fitted parameters
A_fit, B_fit = popt

# Print the fitted parameters
print(f"Fitted parameters:  A = {A_fit}, B = {B_fit}")

# Generate curve fit data
T_fit = np.linspace(min(T_data), max(T_data), 100)
kappa_fit = kappa_lattice_vectorized(T_fit, *popt)

# Plot the original data, fitted data, and extrapolated data
plt.figure(figsize=(8, 6))
plt.plot(T_data, kappa_data, 'o', label='Experimental data')
plt.plot(T_fit, kappa_fit, '-', label='Fitted model')
#plt.plot(T_extrap, kappa_extrap, '--', label='Extrapolated (400-600K)')
plt.xlabel('Temperature (K)')
plt.ylabel('Thermal Conductivity (W/mK)')
plt.legend()
plt.title('Thermal Conductivity Fitting ')
plt.show()
< /code>
import numpy as np
from scipy.integrate import quad
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt

# Constants
kb = 1.380649e-23  # Boltzmann constant (J/K)
hbar = 1.0545718e-34  # Reduced Planck's constant (J*s)
theta_D = 483.0  # Debye temperature (K)
v = 4618.0  # Average phonon group velocity (m/s) -

# Define the inverse relaxation times
def tau_PD_inv(omega, A):
return A * omega**4

def tau_U_inv(omega, B, T):
return B * omega**2 * T * np.exp(-theta_D / (3 * T))

def tau_GB_inv(omega):
D = 1e-6  # Grain boundary characteristic length (m)
return v / D

def tau_new_inv(omega, C):
return C * omega**2

def tau_total_inv(omega, A, B, C, T):
return tau_PD_inv(omega, A) + tau_U_inv(omega, B, T) + tau_GB_inv(omega) + tau_new_inv(omega, C)

# Debye Callaway model for thermal conductivity
def kappa_lattice(T, A, B, C):
def integrand(x):
omega = x * kb * T / hbar
tau_total = 1 / tau_total_inv(omega, A, B, C, T)
return x**4 * np.exp(x) / (np.exp(x) - 1)**2 * tau_total

integral, _ = quad(integrand, 0, theta_D / T)
prefactor = kb**4 * T**3 / (2 * np.pi**2 * v * hbar**3)
return prefactor * integral

# Vectorize the model for curve fitting
def kappa_lattice_vectorized(T, A, B, C):
return np.array([kappa_lattice(Ti, A, B, C) for Ti in T])

# Example thermal conductivity and temperature data for Ge WQ
T_data = np.array([31.21198032, 31.63863233, 33.88415583, 36.02698918, 38.39872854,
39.56437842, 43.09065509, 44.50006634, 46.62080745, 43.20705304,
47.77194375, 51.1374173, 53.33894737, 55.68342411, 57.7839973,
60.44084481, 63.31137245, 66.34735373, 69.16983511, 74.09615366,
77.63998236, 82.66616204, 87.95740211, 93.19748533, 99.59297955,
105.6903923, 113.1334233, 120.6776065, 127.0279836, 134.1693093,
141.3911575, 148.1783285, 155.5512459, 162.631427, 169.7085111,
177.378102, 183.5898785, 190.9367188, 197.6917958, 205.7338436,
212.8065032, 220.0698325, 226.7609477, 234.0571384, 241.1287763,
247.8786187, 255.2722521, 262.3425521, 269.4163916, 276.4872816,
283.5510925, 297.1613034, 304.7919793, 311.8668132, 318.9352276,
326.0927208, 333.0893085, 340.164332, 347.2375209, 354.9542919,
362.0246034, 369.0955138, 375.8247965, 383.5799572])  # Temperature in K
kappa_data = np.array([ 1.085444844, 1.13592882, 1.176429149, 1.217086899,
1.258678341, 1.305290878, 1.357081371, 1.398857871, 1.26084173,
1.436622858, 1.483042548, 1.528031021, 1.561551609, 1.598710119,
1.632418606, 1.665656346, 1.703762245, 1.742520105, 1.774186181,
1.820109169, 1.855169951, 1.893680154, 1.932568555, 1.962172074,
1.999232142, 2.029595098, 2.0490363, 2.067859052, 2.092450212,
2.119053495, 2.131413546, 2.145488003, 2.158842932, 2.169084281,
2.169487956, 2.184424479, 2.196747698, 2.20620033, 2.217744549,
2.223537925, 2.236483003, 2.258522251, 2.273748326, 2.278514432,
2.28270469, 2.288247862, 2.291668986, 2.298648488, 2.302662675,
2.299560106, 2.311410793, 2.339970141, 2.347949378, 2.349474819,
2.354437976, 2.369869851, 2.378039597, 2.384365082, 2.39183344,
2.395266095, 2.399300875, 2.40766296, 2.409888354, 2.442162261])  # Thermal conductivity in W/mK

# Initial guess for parameters A, B, and C
initial_guess = [4.75E-43, 3.6979E-19, 1e-16]

# Curve fitting
popt, pcov = curve_fit(kappa_lattice_vectorized, T_data, kappa_data, p0=initial_guess,maxfev=8_000)
#popt, pcov = curve_fit(kappa_lattice_vectorized, T_data, kappa_data, p0=initial_guess, bounds=(0, np.inf))

# Extract fitted parameters
A_fit, B_fit, C_fit = popt

# Print the fitted parameters
print(f"Fitted parameters:  A = {A_fit}, B = {B_fit}, C = {C_fit}")

# Generate curve fit data
T_fit = np.linspace(min(T_data), max(T_data), 100)
kappa_fit = kappa_lattice_vectorized(T_fit, *popt)

# Plot the original data, fitted data, and extrapolated data
plt.figure(figsize=(8, 6))
plt.plot(T_data, kappa_data, 'o', label='Experimental data')
plt.plot(T_fit, kappa_fit, '-', label='Fitted model')
#plt.plot(T_extrap, kappa_extrap, '--', label='Extrapolated (400-600K)')
plt.xlabel('Temperature (K)')
plt.ylabel('Thermal Conductivity (W/mK)')
plt.legend()
plt.title('Thermal Conductivity Fitting ')
plt.show()