%% FeFET 建模与仿真 - MATLAB 实现
% 本代码实现了铁电场效应晶体管(FeFET)的C-V和I-V特性仿真
% 基于文件2.pdf和3.pdf中的物理模型
%% 主程序
function FeFET_Simulation()
%% 物理常数
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); % 热电压的倒数 (V^{-1})
%% 材料参数
params = struct(...
'N_D', 1e16, ... % cm^{-3}, 衬底施主掺杂浓度
't_F', 200e-7, ... % cm, 铁电层厚度 (200 nm)
'epsilon_F', 200, ... % 铁电层相对介电常数
'epsilon_i', 3.9, ... % 绝缘层相对介电常数
'epsilon_Si', 11.7, ... % Si相对介电常数
'P_s', 10e-6, ... % C/cm^2, 饱和极化
'P_r', 8.5e-6, ... % C/cm^2, 剩余极化
'E_c', 50e3, ... % V/cm, 矫顽电场
'delta', 75e3, ... % V/cm, 铁电材料参数
'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}, 本征载流子浓度
);
%% 计算派生参数
% 热电压
phi_t = k * T / q;
% 德拜长度
params.L_D = sqrt(epsilon_0 * params.epsilon_Si * phi_t / (q * params.N_D));
% 费米势
params.psi_b = -phi_t * log(params.N_D / params.n_i);
% 绝缘层电容
params.C_i = params.A_f * epsilon_0 * params.epsilon_i / params.EOT;
%% 读取外部P-E关系数据
pe_file = 'D:\code\phi\pe.txt';
if ~exist(pe_file, 'file')
error('P-E关系文件未找到: %s', pe_file);
end
% 加载数据
pe_data = load(pe_file);
E_ext = pe_data(:, 1); % MV/m
P_ext = pe_data(:, 2); % C/m²
% 将数据拆分为上扫和下扫分支
[~, max_idx] = max(E_ext);
E_up = E_ext(1:max_idx); % 上扫分支
P_up = P_ext(1:max_idx);
E_down = E_ext(max_idx:end); % 下扫分支
P_down = P_ext(max_idx:end);
% 转换为代码使用的单位 (V/cm, C/cm²)
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 = griddedInterpolant(E_up, P_up, 'linear', 'linear');
pe_interp_down = griddedInterpolant(E_down, P_down, 'linear', 'linear');
%% 设置扫描电压范围
V_g_cv = linspace(-5, 5, 30); % C-V栅压扫描
V_ds_iv = linspace(0, 5, 20); % I-V漏压扫描
V_g_iv = 3.0; % I-V固定栅压
V_ds_transfer = 0.1; % 转移特性固定漏压
fprintf('开始C-V仿真 (上扫)...\n');
C_up = cv_simulation(V_g_cv, 'up', params, epsilon_0, beta, pe_interp_up, pe_interp_down);
fprintf('开始C-V仿真 (下扫)...\n');
C_down = cv_simulation(fliplr(V_g_cv), 'down', params, epsilon_0, beta, pe_interp_up, pe_interp_down);
C_down = fliplr(C_down); % 反转结果以匹配电压顺序
fprintf('开始I-V仿真 (ON状态)...\n');
I_d_on = iv_simulation(V_g_iv, V_ds_iv, 'up', params, epsilon_0, beta, k, T, q, pe_interp_up, pe_interp_down);
fprintf('开始I-V仿真 (OFF状态)...\n');
I_d_off = iv_simulation(V_g_iv, V_ds_iv, 'down', params, epsilon_0, beta, k, T, q, pe_interp_up, pe_interp_down);
fprintf('开始转移特性仿真...\n');
V_g_transfer = linspace(0, 5, 20);
I_d_transfer = zeros(size(V_g_transfer));
for i = 1:length(V_g_transfer)
I_d_transfer(i) = iv_simulation(V_g_transfer(i), [V_ds_transfer], 'up', params, epsilon_0, beta, k, T, q, pe_interp_up, pe_interp_down);
end
%% 绘制结果
figure('Position', [100, 100, 1200, 900]);
% 1. P-E曲线
subplot(2, 3, 1);
plot(E_up, P_up, 'b-', 'LineWidth', 1.5); hold on;
plot(E_down, P_down, 'r--', 'LineWidth', 1.5);
xlabel('Electric Field (V/cm)');
ylabel('Polarization (C/cm^2)');
title('Ferroelectric P-E Curve');
grid on;
legend('Up Sweep', 'Down Sweep', 'Location', 'best');
% 2. C-V曲线
subplot(2, 3, 2);
plot(V_g_cv, C_up * 1e6, 'b-', 'LineWidth', 1.5); hold on;
plot(V_g_cv, C_down * 1e6, 'r--', 'LineWidth', 1.5);
xlabel('Gate Voltage (V)');
ylabel('Capacitance (\muF)');
title('MFIS Capacitor C-V Characteristics');
grid on;
legend('Up Sweep', 'Down Sweep', 'Location', 'best');
% 3. I-V曲线 (ON/OFF状态)
subplot(2, 3, 3);
plot(V_ds_iv, I_d_on * 1e6, 'g-', 'LineWidth', 1.5); hold on;
plot(V_ds_iv, I_d_off * 1e6, 'k--', 'LineWidth', 1.5);
xlabel('Drain Voltage (V)');
ylabel('Drain Current (\muA)');
title(['FeMFET I-V Characteristics (V_g = ', num2str(V_g_iv), 'V)']);
grid on;
legend('ON State', 'OFF State', 'Location', 'best');
% 4. 不同栅压下的I-V曲线
subplot(2, 3, 4);
V_g_list = [1.0, 2.0, 3.0, 4.0];
colors = lines(length(V_g_list));
for i = 1:length(V_g_list)
I_d = iv_simulation(V_g_list(i), V_ds_iv, 'up', params, epsilon_0, beta, k, T, q, pe_interp_up, pe_interp_down);
plot(V_ds_iv, I_d * 1e6, 'Color', colors(i,:), 'LineWidth', 1.5, ...
'DisplayName', ['V_g = ', num2str(V_g_list(i)), 'V']);
hold on;
end
xlabel('Drain Voltage (V)');
ylabel('Drain Current (\muA)');
title('FeMFET I-V at Different Gate Voltages');
grid on;
legend('show');
% 5. 转移特性曲线
subplot(2, 3, 5);
semilogy(V_g_transfer, abs(I_d_transfer) * 1e6, 'mo-', 'LineWidth', 1.5, 'MarkerFaceColor', 'm');
xlabel('Gate Voltage (V)');
ylabel('Drain Current (\muA)');
title(['Transfer Characteristics (V_ds = ', num2str(V_ds_transfer), 'V)']);
grid on;
% 6. 内存窗口分析
subplot(2, 3, 6);
memory_window = abs(C_up - C_down) * 1e6;
plot(V_g_cv, memory_window, 'b-', 'LineWidth', 1.5);
xlabel('Gate Voltage (V)');
ylabel('Memory Window (\muF)');
title('C-V Memory Window');
grid on;
% 保存图像
output_dir = 'D:\code\phi';
if ~exist(output_dir, 'dir')
mkdir(output_dir);
end
saveas(gcf, fullfile(output_dir, 'mfisfet_simulation_results.png'));
fprintf('结果已保存至: %s\n', fullfile(output_dir, 'mfisfet_simulation_results.png'));
% 保存数据
writematrix([V_g_cv', C_up'*1e6, C_down'*1e6], ...
fullfile(output_dir, 'cv_data.txt'), 'Delimiter', 'tab');
writematrix([V_ds_iv', I_d_on'*1e6, I_d_off'*1e6], ...
fullfile(output_dir, 'iv_data.txt'), 'Delimiter', 'tab');
fprintf('仿真数据已保存\n');
end % 主函数结束
%% 外部极化函数
function P = external_polarization(E, direction, pe_interp_up, pe_interp_down)
if strcmp(direction, 'up')
P = pe_interp_up(E);
else % 'down'
P = pe_interp_down(E);
end
end
%% 半导体表面电荷密度
function Q_s_val = Q_s(psi_s, params, epsilon_0, beta)
term = (params.n_i^2 / params.N_D^2) * (exp(-beta * psi_s) + beta * psi_s - 1) + ...
(exp(beta * psi_s) - beta * psi_s - 1);
Q_s_val = -sign(psi_s) * (sqrt(2) * epsilon_0 * params.epsilon_i) / (beta * params.L_D) * sqrt(term);
end
%% C-V 特性仿真
function C_total = cv_simulation(V_g, direction, params, epsilon_0, beta, pe_interp_up, pe_interp_down)
C_total = zeros(size(V_g));
options = optimoptions('fsolve', 'Display', 'off', 'Algorithm', 'trust-region-dogleg', 'FunctionTolerance', 1e-6);
% 初始猜测
psi_s0 = 0.1;
E_f0 = 0.1;
for i = 1:length(V_g)
% 定义方程组
function F = equations(x)
psi_s = x(1);
E_f = x(2);
% 计算极化
P = external_polarization(E_f, direction, pe_interp_up, pe_interp_down);
% 计算表面电荷
Q_s_val = Q_s(psi_s, params, epsilon_0, beta);
% 方程1: P = -Q_s
eq1 = P + Q_s_val;
% 方程2: V_G = psi_s + V_I + V_F
V_i = -Q_s_val * params.EOT / (epsilon_0 * params.epsilon_i);
V_f = E_f * params.t_F;
eq2 = V_g(i) - (psi_s + V_i + V_f);
F = [eq1; eq2];
end
% 使用fsolve求解
try
[sol, ~, exitflag] = fsolve(@equations, [psi_s0; E_f0], options);
if exitflag > 0
psi_s = sol(1);
E_f = sol(2);
% 更新初始猜测
psi_s0 = psi_s;
E_f0 = E_f;
% 计算半导体电容 (公式4.5)
term = (params.n_i^2 / params.N_D^2) * (-exp(-beta * psi_s) + 1) + ...
(exp(beta * psi_s) - 1);
% 数值保护
if term < 0
term = 1e-6;
end
C_si = (params.A_f * epsilon_0 * params.epsilon_i) / (sqrt(2) * params.L_D) * sqrt(term);
% 铁电层电容 (公式4.6)
C_f = params.A_f * epsilon_0 * params.epsilon_F / params.t_F;
% 绝缘层电容 (公式4.7)
C_i = params.C_i;
% 总电容 (公式4.8)
C_total(i) = 1 / (1/C_f + 1/C_i + 1/C_si);
else
if i > 1
C_total(i) = C_total(i-1);
else
C_total(i) = 1e-12;
end
warning('CV求解在V_g=%f时失败', V_g(i));
end
catch
if i > 1
C_total(i) = C_total(i-1);
else
C_total(i) = 1e-12;
end
warning('CV求解在V_g=%f时出错', V_g(i));
end
end
end
%% 求解表面势
function psi_s = solve_surface_potential(V_g, V, direction, params, epsilon_0, beta, pe_interp_up, pe_interp_down)
options = optimoptions('fsolve', 'Display', 'off', 'Algorithm', 'trust-region-dogleg', 'FunctionTolerance', 1e-6);
% 初始猜测
psi_s0 = abs(params.psi_b) + 0.5;
E_f0 = 0.1;
% 定义方程组
function F = equations(x)
psi_s = x(1);
E_f = x(2);
% 计算极化
P = external_polarization(E_f, direction, pe_interp_up, pe_interp_down);
% 计算表面电荷
Q_s_val = Q_s(psi_s, params, epsilon_0, beta);
% 方程1: P = -Q_s
eq1 = P + Q_s_val;
% 方程2: V_G = psi_s + V_I + V_F
V_i = -Q_s_val * params.EOT / (epsilon_0 * params.epsilon_i);
V_f = E_f * params.t_F;
eq2 = V_g - (psi_s + V_i + V_f);
F = [eq1; eq2];
end
% 使用fsolve求解
[sol, ~, exitflag] = fsolve(@equations, [psi_s0; E_f0], options);
if exitflag > 0
psi_s = sol(1);
else
psi_s = abs(params.psi_b); % 使用默认值
warning('表面势求解在V_g=%f, V=%f时失败', V_g, V);
end
end
%% Pao-Sah 被积函数
function f = pao_sah_integrand(psi, V, V_g, direction, params, epsilon_0, beta, k, T, pe_interp_up, pe_interp_down)
% 计算公式(21)中的内部表达式
term1 = exp(beta * psi) - beta * psi - 1;
term2 = (params.n_i^2 / params.N_D^2) * exp(beta * V) * ...
(exp(-beta * psi) + beta * psi * exp(-beta * V) - 1);
% 添加数值保护
total_term = term1 + term2;
if total_term <= 0
total_term = 1e-10;
end
% 计算电场 ξ(ψ,V)
prefactor = (2 * params.N_D * k * T) / (epsilon_0 * params.epsilon_Si);
xi = sqrt(prefactor * total_term);
% 公式(20)被积函数
numerator = (params.n_i^2 / params.N_D) * exp(-beta * (psi - V));
% 避免除零
if xi < 1e-10
f = 0;
else
f = numerator / xi;
end
end
%% I-V 特性仿真
function I_d = iv_simulation(V_g, V_ds, direction, params, epsilon_0, beta, k, T, q, pe_interp_up, pe_interp_down)
I_d = zeros(size(V_ds));
for i = 1:length(V_ds)
% 离散化沟道电势
V_points = linspace(0, V_ds(i), 10);
integral_val = 0;
for j = 1:length(V_points)
V_val = V_points(j);
% 求解该沟道电势下的表面势
psi_s = solve_surface_potential(V_g, V_val, direction, params, epsilon_0, beta, pe_interp_up, pe_interp_down);
% 只在表面势大于费米势时积分
if psi_s > params.psi_b
% 离散化表面势范围
psi_points = linspace(params.psi_b, psi_s, 20);
% 计算被积函数值
integrand_vals = zeros(size(psi_points));
for m = 1:length(psi_points) % 使用m替代k避免冲突
integrand_vals(m) = pao_sah_integrand(psi_points(m), V_val, V_g, direction, ...
params, epsilon_0, beta, k, T, pe_interp_up, pe_interp_down);
end
% 使用梯形法进行内层积分
inner_integral = trapz(psi_points, integrand_vals);
% 计算dV (梯形法)
if j == 1
dV = V_points(2) - V_points(1);
elseif j == length(V_points)
dV = V_points(end) - V_points(end-1);
else
dV = (V_points(j+1) - V_points(j-1)) / 2;
end
integral_val = integral_val + inner_integral * dV;
end
end
% 计算漏极电流
I_d(i) = q * params.mu * (params.W / params.L) * integral_val;
end
end
% 运行主程序
FeFET_Simulation();错误: 文件: FeFET.m 行: 211 列: 9
函数定义放置或嵌套错误。
最新发布