Solve a pde equation with finite differences for Simulink

2 次查看（过去 30 天）
Raúl Rivera
Raúl Rivera 2025-1-30
编辑： Raúl Rivera 2025-1-30
Ran in:
Hello , I want to transform this code that solves a pde equation with the ode solver into finite diferences, because I want to take the code as a matlab function block in simulink so it stands no ode solver(since it is an iterator take much time every time step so never ends simulation ) thats why i want to take it into finite differences .The equations are the following
The inital code is the following with ode solver:
L = 20 ; % Longitud del lecho (m)
eps = 0.4; % Porosidad
u = 0.2; % Velocidad superficial del fluido (m/s)
k_f = 0.02; % Constante de transferencia de masa (1/s)
c0 = 0;
Kf = 4; % Constante de Freundlich
rhop = 1520;
n = 2; % Exponente de Freundlich
% Concentración inicial del fluido (kg/m³)
q0 = 4.320; % Concentración inicial en el sólido (kg/m³)
% Densidad del adsorbente (kg/m³)
tf = 10; % Tiempo final de simulación (horas)
Nt = 100;
t = linspace(0, tf*3600, Nt);
Nz = 100;
z = linspace(0, L,Nz);
dz = z(2) - z(1);
% Initial conditions
ICA = max(ones(1, Nz) * c0, 1e-12); % Evitar valores negativos o cero
ICB = ones(1, Nz) * q0;
IC = [ICA ICB];
options = odeset('RelTol', 1e-6, 'AbsTol', 1e-8, 'InitialStep', 1e-4, 'MaxStep', 100);
[t, y] = ode15s(@fun_pde, t, IC, options, Nz, eps, n, Kf, k_f, u, rhop, dz);
% Define value
cc = y(:, 1:Nz);
qq = y(:, Nz+1:end);
% Recalculate new limit conditions
cc(:, 1) = 0;
cc(:, end) = cc(:, end-1);
% Plotting
cp = cc(:, end) ./ c0;
qp = qq(:, :) ./ q0;
%q_promedio = mean(qq, 2); % Promedio de q en el lecho para cada instante de tiempo
%conversion = 1 - (q_promedio / q0); % Conversión normalizada
figure;
subplot(2, 1, 1);
time = t / 3600; % Convertir a horas
plot(time, 1- qp, 'b', 'LineWidth', 1.5);
xlabel('Tiempo (horas)');
ylabel('Conversion');
title('Curva de conversión durante la desorción');
grid on;
subplot(2, 1, 2);
plot(t / 3600, (cc(:,:)), 'LineWidth', 1.5);
xlabel('Tiempo (horas)');
ylabel('Soluciòn kg/m3');
title('Curva de carga de la solucion durante la desorciòn');
grid on;
% PDE function
function dydt = fun_pde(~, y, Nz, eps, n, Kf, k_f, u, rhop, dz)
dcdt = zeros(Nz, 1);
dqdt = zeros(Nz, 1);
c = y(1:Nz);
q = y(Nz+1:2*Nz);
% Boundary conditions
c(1) = max(c(1), 0); % Asegurar que c(1) sea no negativo
c(end) = c(end-1); % Asegurar que c(1) sea no negativo
% Interior nodes
qstar = zeros(Nz, 1);
dcdz = zeros(Nz, 1);
for i = 2:Nz-1
qstar(i) = Kf .* max(c(i), 1e-12).^(1/n); % Evitar problemas numéricos
dqdt(i) = k_f .* (qstar(i) - q(i));
% if i < Nz
dcdz(i) = (c(i+1) - c(i-1)) / (2 * dz);
%else
% dcdz(i) = (c(i) - c(i-1)) / dz;
%end
dcdt(i) = -u * dcdz(i) - rhop * ((1 - eps) / eps) .* dqdt(i);
end
dydt = [dcdt; dqdt];
end
next is a try to solve with finite diferences but get someting different:
L = 20 ; % Longitud del lecho (m)
eps = 0.4; % Porosidad
u = 0.2; % Velocidad superficial del fluido (m/s)
k_f = 0.02; % Constante de transferencia de masa (1/s)
c0 = 0;
Kf = 4; % Constante de Freundlich
rhop = 1520;
n = 2; % Exponente de Freundlich
% Concentración inicial del fluido (kg/m³)
q0 = 4.320; % Concentración inicial en el sólido (kg/m³)
% Densidad del adsorbente (kg/m³)
tf = 10; % Tiempo final de simulación (horas)
Nz = 100; % Número de nodos espaciales
% Discretización espacial y temporal
z = linspace(0, L, Nz);
t = linspace(0, tf*3600, Nt);
dz = z(2) - z(1);
dt = t(2) - t(1); % Paso temporal
% Condiciones iniciales
c = ones(Nt, Nz) * c0; % Concentración en el fluido
q = ones(Nt, Nz) * q0; % Concentración en el sólido
% Iteración en el tiempo (Diferencias Finitas Explícitas)
for ti = 1:Nt-1
for zi = 2:Nz-1
% Isoterma de Freundlich
qstar = Kf * max(c(ti, zi), 1e-12)^(1/n);
% Transferencia de masa (Desorción)
dqdt = k_f * (qstar - q(ti, zi));
% Gradiente espacial de concentración (Diferencias centradas)
dcdz = (c(ti, zi+1) - c(ti, zi-1)) / (2 * dz);
% Ecuación de balance de masa en el fluido
dcdt = -u * dcdz - rhop * ((1 - eps) / eps) * dqdt;
% Actualizar valores asegurando que sean positivos
c(ti+1, zi) = max(c(ti, zi) + dcdt * dt, 0);
q(ti+1, zi) = max(q(ti, zi) + dqdt * dt, 0);
end
end
% Condiciones de frontera
c(:, 1) = c0; % Entrada con concentración baja
c(:, Nz) = c(:, Nz-1); % Gradiente nulo en la salida
% Cálculo de la conversión normalizada
qp = q(:, :) ./ q0;
% Graficar resultados
figure;
subplot(2, 1, 1);
plot(t / 3600, 1-qp, 'b', 'LineWidth', 1.5);
xlabel('Tiempo (horas)');
ylabel('Conversion');
title('Curva de conversión durante la desorción');
grid on;
subplot(2, 1, 2);
c_salida = c(:, :); % Concentración en la salida del lecho
plot(t / 3600, c_salida, 'r', 'LineWidth', 1.5);
xlabel('Tiempo (horas)');
ylabel('Soluciòn kg/m3');
title('Curva de carga de la solucion durante la desorciòn');
grid on;
I dont know why is wrong
Thanks in advance

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 PDE Solvers 的更多信息

标签

产品


版本

R2024a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by