import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import root
from scipy.integrate import quad
from scipy.interpolate import interp1d
import os
# =============================================================
# 物理常数
# =============================================================
epsilon_0 = 8.854e-14 # 真空介电常数 (F/cm)
q = 1.602e-19 # 电子电荷 (C)
k = 1.380649e-23 # 玻尔兹曼常数 (J/K)
T = 300 # 温度 (K)
beta = q / (k * T) # β = q/(kT): 热电压的倒数 (V^{-1})
phi_t = k * T / q # 热电压 (V)
# =============================================================
# 材料参数
# =============================================================
params = {
'N_D': 1e16, # cm^{-3}, 衬底施主掺杂浓度
't_F': 200e-7, # cm, 铁电层厚度 (200 nm)
'epsilon_F': 200, # 铁电层相对介电常数
'epsilon_i': 3.9, # 绝缘层相对介电常数
'epsilon_Si': 11.7, # Si相对介电常数
'EOT': 4e-7, # cm, 等效氧化层厚度
'A_f': 100e-8, # cm^2, 栅面积
'W': 10e-4, # cm, 沟道宽度
'L': 1e-4, # cm, 沟道长度
'mu': 200, # cm^2/V·s, 载流子迁移率
'n_i': 1.5e10, # cm^{-3}, 本征载流子浓度
}
# =============================================================
# 计算派生参数
# =============================================================
params['C_i'] = params['A_f'] * epsilon_0 * params['epsilon_i'] / params['EOT']
params['L_D'] = np.sqrt(epsilon_0 * params['epsilon_Si'] * phi_t / (q * params['N_D']))
params['psi_b'] = -phi_t * np.log(params['N_D'] / params['n_i'])
# =============================================================
# 读取外部P-E关系数据
# =============================================================
pe_file = r"D:\code\phi\pe.txt"
if not os.path.exists(pe_file):
raise FileNotFoundError(f"P-E关系文件未找到: {pe_file}")
pe_data = np.loadtxt(pe_file)
E_ext = pe_data[:, 0] # MV/m
P_ext = pe_data[:, 1] # C/m²
max_idx = np.argmax(E_ext)
E_up = E_ext[:max_idx+1] # 上扫分支
P_up = P_ext[:max_idx+1]
E_down = E_ext[max_idx:] # 下扫分支
P_down = P_ext[max_idx:]
# 转换为代码使用的单位
E_up = E_up * 1e4 # MV/m -> V/cm
P_up = P_up * 1e-4 # C/m² -> C/cm²
E_down = E_down * 1e4
P_down = P_down * 1e-4
# 创建插值函数
pe_interp_up = interp1d(E_up, P_up, kind='linear', fill_value="extrapolate")
pe_interp_down = interp1d(E_down, P_down, kind='linear', fill_value="extrapolate")
# =============================================================
# 外部极化函数
# =============================================================
def external_polarization(E, direction):
if direction == 'up':
return pe_interp_up(E)
else: # down
return pe_interp_down(E)
# =============================================================
# C-V 特性仿真
# =============================================================
def Q_s(psi_s):
"""
半导体表面电荷密度 (C/cm^2)
"""
# psi_s = np.clip(psi_s, -20.0, 20.0)
#print(psi_s)
term1 = (params['n_i']**2 / params['N_D']**2) * (np.exp(-beta * psi_s) + beta * psi_s - 1)
term2 = np.exp(beta * psi_s) - beta * psi_s - 1
total_term = term1 + term2
if total_term < 0:
total_term = 1e-30
prefactor = (np.sqrt(2) * epsilon_0 * params['epsilon_Si']) / (beta * params['L_D'])
# print(psi_s)
return -np.sign(psi_s) * prefactor * np.sqrt(total_term)
def cv_simulation(V_g, polarization_func, direction='up'):
"""
C-V特性仿真
"""
C_total = np.zeros_like(V_g)
psi_s0, E_f0 = 0.1, 0.1
for i, V_g_val in enumerate(V_g):
def equations(vars):
psi_s, E_f = vars
#psi_s = np.clip(psi_s, -20.0, 20.0)
# E_f = np.clip(E_f, -1e7, 1e7)
P = polarization_func(E_f, direction)
Q_s_val = Q_s(psi_s)
#print(Q_s(psi_s))
V_f = E_f * params['t_F']
V_i = -Q_s_val * params['EOT'] / (epsilon_0 * params['epsilon_i'])
eq1 = P + Q_s_val
eq2 = V_g_val - (psi_s + V_i + V_f)
return [eq1, eq2]
sol = root(equations, [psi_s0, E_f0], method='lm', tol=1e-4)
if sol.success:
psi_s, E_f = sol.x
psi_s0, E_f0 = psi_s, E_f
d_psi = 1e-3
Q_s_plus = Q_s(psi_s + d_psi)
Q_s_minus = Q_s(psi_s - d_psi)
dQ_s_dpsi = (Q_s_plus - Q_s_minus) / (2 * d_psi)
C_si = params['A_f'] * np.abs(dQ_s_dpsi)
C_f = params['A_f'] * epsilon_0 * params['epsilon_F'] / params['t_F']
C_i = params['C_i']
C_total[i] = 1 / (1/C_f + 1/C_i + 1/C_si)
else:
if i > 0:
C_total[i] = C_total[i-1]
else:
C_total[i] = 1e-12
return C_total
# =============================================================
# I-V 特性仿真 - 改进版本 (使用quad积分)
# =============================================================
def solve_surface_potential(V_g, polarization_func, direction):
"""
求解表面势 - 增强鲁棒性
"""
def equations(vars):
psi_s, E_f = vars
psi_s = np.clip(psi_s, -20.0, 20.0)
E_f = np.clip(E_f, -1e7, 1e7)
P = polarization_func(E_f, direction)
Q_s_val = Q_s(psi_s)
#print(Q_s(psi_s))
V_i = -Q_s_val * params['EOT'] / (epsilon_0 * params['epsilon_i'])
V_f = E_f * params['t_F']
eq1 = P + Q_s_val
eq2 = V_g - (psi_s + V_i + V_f)
return [eq1, eq2]
# 尝试不同初始值
initial_guesses = [
[0.5, 1e5],
[-0.5, -1e5],
[params['psi_b'], 0],
[-params['psi_b'], 0],
[0, 0]
]
for initial_guess in initial_guesses:
sol = root(equations, initial_guess, method='lm', tol=1e-4)
if sol.success:
return sol.x[0]
return abs(params['psi_b'])
def pao_sah_integrand(psi, V_ds):
"""
Pao-Sah被积函数 - 改进数值稳定性
"""
psi = np.clip(psi, -20.0, 20.0)
V_ds = np.clip(V_ds, -1.0, 1.0)
beta_psi = beta * psi
beta_Vds = beta * V_ds
# 限制指数范围防止溢出
beta_psi = np.clip(beta_psi, -30, 30)
beta_Vds = np.clip(beta_Vds, -30, 30)
# 计算各项
exp_beta_psi = np.exp(beta_psi)
exp_beta_Vds = np.exp(beta_Vds)
exp_neg_beta_psi = np.exp(-beta_psi)
term1 = exp_beta_psi - beta_psi - 1
term2 = (params['n_i']**2 / params['N_D']**2) * exp_beta_Vds * \
(exp_neg_beta_psi + beta_psi * np.exp(-beta_Vds) - 1)
total_term = term1 + term2
# 处理小值情况
if total_term <= 1e-40:
return 0
prefactor = (2 * params['N_D'] * k * T) / (epsilon_0 * params['epsilon_Si'])
xi = np.sqrt(prefactor * total_term)
if xi < 1e-20:
return 0
exponent = -beta * (psi - V_ds)
exponent = np.clip(exponent, -30, 30)
numerator = (params['n_i']**2 / params['N_D']) * np.exp(exponent)
return numerator / xi
def pao_sah_integral(V_ds, V_g, polarization_func, direction):
"""
使用quad计算Pao-Sah积分
"""
try:
psi_s = solve_surface_potential(V_g, polarization_func, direction)
except:
psi_s = abs(params['psi_b'])
# 确定积分限
a = min(params['psi_b'], psi_s)
b = max(params['psi_b'], psi_s)
# 避免积分区间过小
if np.abs(b - a) < 1e-6:
return 0
# 使用quad进行自适应积分
try:
# 处理被积函数在端点附近的奇异行为
# 在端点附近添加微小偏移
a_adj = a + 1e-5 if a == params['psi_b'] else a
b_adj = b - 1e-5 if b == psi_s else b
# 使用quad进行积分
integral, _ = quad(
lambda psi: pao_sah_integrand(psi, V_ds),
a_adj, b_adj,
points=[params['psi_b'], psi_s], # 指定关键点
epsabs=1e-10, # 绝对误差容限
epsrel=1e-6, # 相对误差容限
limit=1000 # 最大子区间数
)
return integral
except Exception as e:
print(f"积分失败: V_ds={V_ds}, V_g={V_g}, 错误: {e}")
return 0
def iv_simulation(V_g, V_ds, polarization_func, direction='up'):
"""
I-V特性仿真 - 使用quad积分
"""
I_d = np.zeros_like(V_ds)
# 预计算表面势（减少重复计算）
psi_s = solve_surface_potential(V_g, polarization_func, direction)
#print(psi_s)
for i, V_ds_val in enumerate(V_ds):
if np.abs(V_ds_val) < 1e-6:
I_d[i] = 0
continue
try:
integral = pao_sah_integral(V_ds_val, V_g, polarization_func, direction)
if not np.isfinite(integral) or integral < 0:
integral = 0
# 计算电流
I_d[i] = q * params['mu'] * (params['W'] / params['L']) * integral
if not np.isfinite(I_d[i]) or I_d[i] < 0:
I_d[i] = 0
except Exception as e:
print(f"IV求解失败: V_ds={V_ds_val}, 错误: {e}")
I_d[i] = 0
return I_d
# =============================================================
# 主程序
# =============================================================
if __name__ == "__main__":
# 设置扫描电压范围
V_g_cv = np.linspace(-1.6, -0.4, 30) # C-V栅压扫描
V_ds_iv = np.linspace(-0.5, 0, 20) # I-V漏压扫描
V_g_iv = -0.3 # I-V固定栅压
print("开始C-V仿真 (上扫)...")
C_up = cv_simulation(V_g_cv, external_polarization, direction='up')
print("开始C-V仿真 (下扫)...")
C_down = cv_simulation(V_g_cv[::-1], external_polarization, direction='down')
C_down = C_down[::-1]
print("开始I-V仿真 (ON状态)...")
I_d_on = iv_simulation(V_g_iv, V_ds_iv, external_polarization, direction='up')
print("开始I-V仿真 (OFF状态)...")
I_d_off = iv_simulation(V_g_iv, V_ds_iv, external_polarization, direction='down')
# 绘制结果
plt.figure(figsize=(15, 10))
# 1. P-E曲线
plt.subplot(2, 2, 1)
plt.plot(E_up, P_up, 'b-', label='Up Sweep')
plt.plot(E_down, P_down, 'r--', label='Down Sweep')
plt.xlabel('Electric Field (V/cm)')
plt.ylabel('Polarization (C/cm²)')
plt.title('Ferroelectric P-E Curve')
plt.grid(True)
plt.legend()
# 2. C-V曲线
plt.subplot(2, 2, 2)
plt.plot(V_g_cv, C_up * 1e6, 'b-', label='Up Sweep')
plt.plot(V_g_cv, C_down * 1e6, 'r--', label='Down Sweep')
plt.xlabel('Gate Voltage (V)')
plt.ylabel('Capacitance (μF)')
plt.title('MFIS Capacitor C-V Characteristics')
plt.grid(True)
plt.legend()
# 3. I-V曲线 (ON/OFF状态)
plt.subplot(2, 2, 3)
plt.plot(V_ds_iv, I_d_on * 1e6, 'g-', label='ON State')
plt.plot(V_ds_iv, I_d_off * 1e6, 'k--', label='OFF State')
plt.xlabel('Drain Voltage (V)')
plt.ylabel('Drain Current (μA)')
plt.title(f'FeMFET I-V Characteristics (V_g = {V_g_iv}V)')
plt.grid(True)
plt.legend()
# 4. 不同栅压下的I-V曲线
plt.subplot(2, 2, 4)
V_g_list = [-1.0, -0.5, 0, 0.5] # 修改为实际有意义的栅压值
for V_g_val in V_g_list:
I_d = iv_simulation(V_g_val, V_ds_iv, external_polarization, direction='up')
plt.plot(V_ds_iv, I_d * 1e6, label=f'V_g = {V_g_val}V')
plt.xlabel('Drain Voltage (V)')
plt.ylabel('Drain Current (μA)')
plt.title('FeMFET I-V at Different Gate Voltages')
plt.grid(True)
plt.legend()
plt.tight_layout()
# 保存图像
output_dir = r"D:\code\phi"
if not os.path.exists(output_dir):
os.makedirs(output_dir)
plt.savefig(os.path.join(output_dir, "mfisfet_simulation_results.png"), dpi=300)
print(f"结果已保存至: {os.path.join(output_dir, 'mfisfet_simulation_results.png')}")
plt.show()
# 保存数据
np.savetxt(os.path.join(output_dir, "cv_data.txt"),
np.column_stack((V_g_cv, C_up * 1e6, C_down * 1e6)),
header="V_g(V) C_up(μF) C_down(μF)")
np.savetxt(os.path.join(output_dir, "iv_data.txt"),
np.column_stack((V_ds_iv, I_d_on * 1e6, I_d_off * 1e6)),
header="V_ds(V) I_on(μA) I_off(μA)")
print("仿真数据已保存")半导体电容是在哪部分计算的
扫一扫
952
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
