Passing out additions parameters after ODE solver.
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
Thomas Stone-Wigg
am 23 Jan. 2024
Kommentiert: Thomas Stone-Wigg
am 25 Jan. 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 Kommentare
Torsten
am 23 Jan. 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" ?
Akzeptierte Antwort
Torsten
am 23 Jan. 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 Kommentare
Torsten
am 24 Jan. 2024
Bearbeitet: Torsten
am 24 Jan. 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 ?
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Applications finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!