Passing out additions parameters after ODE solver.
1 visualización (últimos 30 días)
Mostrar comentarios más antiguos
Thomas Stone-Wigg
el 23 de En. de 2024
Comentada: Thomas Stone-Wigg
el 25 de En. de 2024
I am trying to pass pass out the velocity array at the end of adsorption model. I would like to create a matrix of time and velocity to plot a graph. with the other graphs where their parameters are from the ODE solver. Here is the code with some of the irrelavent parts taken out:
clear,clc
%% MAIN
% Physical Parameters
L = 1; % Column length m
r = 0.25; % Bed radius m
t0 = 0; % Initial Time
tf = 2; % Final time
dt = 0.1; % Time step
t = t0:dt:tf; % Time vector
dz = 0.05; % Mesh size
z = 0:dz:L; % Mesh generation
n = numel(z);
%%%%%%% blended gas parameters %%%%%%%%%%%
%%% calc u input
TPF = [1 300 8e5 2e5 1e5]; % Feed: Velocity (m/s), tempurature (K) and pressure (Bar) [Vfeed Tfeed PH PI PL]
yF = 0.7; % Mole fraction for methane
% y is an array size n*5 of y1 = 1:n, q1 = n+1:2*n,
% q2 = 2*n+1:3*n, T = 3*n+1:4*n, P = 4*n+1:5*n
sol = adsorptionSolver(t,z,yF,TPF);
% sol.x is time steps
purityh = 100*Moley2(n,end) / (sol.y(n,end) + Moley2(n,end));
fprintf('Purity of hydrogen is %f%% after %4.2fs\n', purityh, tf)
figure(1)
subplot(1,2,1)
mesh(sol.x,z,sol.y(1:n,1:end))
xlabel('t')
ylabel('bed length')
zlabel('mole fraction y1')
title('mole fraction of methane')
subplot(1,2,2)
mesh(sol.x,z,Moley2)
xlabel('t')
ylabel('bed length')
zlabel('mole fraction y2')
title('mole fraction of hydrogen')
figure(2)
subplot(1,2,1)
mesh(sol.x,z,sol.y(n+1:2*n,1:end))
xlabel('t')
ylabel('bed length')
zlabel('loading q1 mol/kg')
title('loading of methane')
subplot(1,2,2)
mesh(sol.x,z,sol.y(2*n+1:3*n,1:end))
xlabel('t')
ylabel('bed length')
zlabel('loading q2 mol/kg')
title('loading of hydrogen')
figure(3)
subplot(1,2,1)
mesh(sol.x,z,sol.y(3*n+1:4*n,1:end))
xlabel('t')
ylabel('bed length')
zlabel('Temp')
title('temp of system')
subplot(1,2,2)
mesh(sol.x,z,sol.y(4*n+1:5*n,1:end))
xlabel('time')
ylabel('bed length')
zlabel('Pressure')
title('Pressure of system')
%% Solver function
function out = adsorptionSolver(t,z,yFeed,TPFeed)
n = numel(z); % Size of mesh grid
h = diff(z);
h = h(1);
% start up condidtions
Tw = 300; % Ambient Tempurature K
Pw = 8e5; % Ambient Pressure Pa
y1w = 0.7;
% Start up conditions Conditions
y1_0 = zeros(n,1) + y1w;
q1_0 = zeros(n,1);
q2_0 = zeros(n,1);
T_0 = zeros(n,1) + Tw;
P_0 = zeros(n,1) + Pw;
% y is an array size n*4 of c1 = 1:n, q1 = n+1:2*n, c2 = 2*n+1:3*n, q2 = 3*n+1:4*n
y0 = [y1_0; q1_0; q2_0; T_0; P_0];
% % Solving
out = ode15s(@(t,y) adsorptionModelODE(t,y,h,n,yFeed,TPFeed),t,y0);
end
Ive tried to use another function to separate the two adsorption model fucntion outputs but i havent got it to work and it seems inefficient
function dydt = adsorptionModelODE(t, y, h, n, yiFeed, TPFeed)
[dydt, ~] = adsorptionModel(t, y, h, n, yiFeed, TPFeed);
end
%% Adsorption model
function [dydt, Velocity] = adsorptionModel(t,y,h,n,yiFeed,TPFeed)
% t is used to calculate boundary pressure during a few cycle steps
% rest of code and odes
Velocity=cat(1,uhalf,u);
Velocity([1:2:end,2:2:end])=Velocity;
%%%%%%%%%%%%%%%% Concatenate vectors of time derivatives
dydt = [dy1dt;dq1dt;dq2dt;dTdt;dPdt];
end
Ive seen various forum pages but being a beginner they have confused me. I think persistent is the best way to implement it but ive tried and failed.
Thanks in advance,
Tom
4 comentarios
Torsten
el 23 de En. de 2024
Hopefully this will show what im trying to get at with my orginial question
No. I can only guess: You want to make a surface plot of "purityh" ?
Respuesta aceptada
Torsten
el 23 de En. de 2024
sol = adsorptionSolver(t,z,yF,TPF)
for i = 1:size(sol.y,2)
[~,velocity(:,i)] = adsorptionModel(sol.x(i),sol.y(:,i),z(2)-z(1),numel(z),yF,TPF);
end
Since "velocity" is an array of size 43x61 and not 21x61, I don't know how the 43 should be interpreted.
9 comentarios
Torsten
el 24 de En. de 2024
Editada: Torsten
el 24 de En. de 2024
Sorry, but it means too much effort for me to check.
But it cannot be all you have to do to switch to first-order upwind by defining
Phalf(2:n) = P(1:n-1);
I always build up a code in steps and check whether it produces reasonable results from step to step. So I think it cannot be true that you started with the difficult WENO scheme, can it ?
Más respuestas (0)
Ver también
Categorías
Más información sobre Applications en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!