import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad
from scipy.optimize import fsolve
plt.rcParams['font.sans-serif'] = ['SimHei'] # 解决中文显示问题
plt.rcParams['axes.unicode_minus'] = False # 解决负号显示问题
# ==============================
# 材料参数设定(以钆镓石榴石为例)
# ==============================
Tc = 2.0 # 居里温度 (K)
H_max = 5.0 # 最大磁场强度 (T)
C_m = 10.0 # 磁热容 (J/(kg·K))
mu0 = 4*np.pi*1e-7 # 真空磁导率
# ==============================
# 磁熵变计算函数 (ΔS_M)
# ==============================
def delta_SM(T, H):
"""磁熵变积分公式 ΔS_M = ∫(∂M/∂T)_H dH"""
# 简化假设:M = M_sat * (1 - (T/Tc)^1.5) ,在H足够大时近似
M_sat = 1.0 # 饱和磁化强度 (A/m)
dMdT = -1.5 * M_sat * (T**(0.5))/Tc**1.5 # 导数项
return quad(lambda h: dMdT, 0, H)[0] # 对H积分
# ==============================
# 温度-熵关系模型
# ==============================
def entropy(T, H=0):
"""总熵 = 磁熵 + 晶格熵"""
S_mag = delta_SM(T, H)
S_lattice = C_m * np.log(T) if T > 0 else -np.inf # 处理边界情况
return S_mag + S_lattice
def temperature(S_target, H):
"""求解给定熵和磁场时的温度"""
def f(T):
return entropy(T, H) - S_target
try:
return fsolve(f, x0=Tc)[0]
except:
return np.nan
# ==============================
# 磁制冷循环模拟
# ==============================
# 初始状态
T_initial = 1.5 # 初始温度 (K)
# 等温磁化 (A→B)
T1 = T_initial
S1 = entropy(T1, H=0)
S2 = entropy(T1, H=H_max)
# 绝热退磁 (B→C)
H_steps_rev = np.linspace(H_max, 0, 100)
S_BC = np.full_like(H_steps_rev, S2)
T_BC = [temperature(s, h) for s, h in zip(S_BC, H_steps_rev)]
# 等温退磁 (C→D)
T2 = np.nanmin(T_BC)
S3 = entropy(T2, H=0)
S4 = entropy(T2, H=H_max)
# 绝热磁化 (D→A)
H_steps_mag = np.linspace(0, H_max, 100)
S_DA = np.full_like(H_steps_mag, S3)
T_DA = [temperature(s, h) for s, h in zip(S_DA, H_steps_mag)]
# ==============================
# 绘图
# ==============================
plt.figure(figsize=(12, 6))
# 温熵图 (T-S图)
plt.subplot(1,2,1)
T_range = np.linspace(0.5, 3, 100)
for H in [0, H_max]:
S = [entropy(T, H) for T in T_range]
plt.plot(S, T_range, label=f'H={H} T')
# 绘制循环路径
plt.plot([S1, S2], [T1, T1], 'r--', label='等温磁化') # A→B
plt.plot([entropy(t, h) for t, h in zip(T_BC, H_steps_rev)], T_BC, 'g--', label='绝热退磁') # B→C
plt.plot([S3, S4], [T2, T2], 'b--', label='等温退磁') # C→D
plt.plot([entropy(t, h) for t, h in zip(T_DA, H_steps_mag)], T_DA, 'm--', label='绝热磁化') # D→A
plt.xlabel('熵 S (J/kg·K)')
plt.ylabel('温度 T (K)')
plt.title('温熵图 (T-S Diagram)')
plt.grid(True)
plt.legend()
# 磁制冷循环图
plt.subplot(1,2,2)
H_cycle = np.concatenate([
np.linspace(0, H_max, 100),
np.linspace(H_max, 0, 100),
np.linspace(0, -H_max, 100),
np.linspace(-H_max, 0, 100)
])
T_cycle = np.concatenate([
[T1]*100,
T_BC,
[T2]*100,
T_DA
])
plt.plot(H_cycle, T_cycle, 'r-')
plt.xlabel('磁场强度 H (T)')
plt.ylabel('温度 T (K)')
plt.title('磁制冷循环温度变化')
plt.grid(True)
plt.tight_layout()
plt.show()
帮我详细讲解这个代码中函数的作用
最新发布