Passing out additions parameters after ODE solver.
4 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
Thomas Stone-Wigg
il 23 Gen 2024
Commentato: Thomas Stone-Wigg
il 25 Gen 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 Commenti
Torsten
il 23 Gen 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" ?
Risposta accettata
Torsten
il 23 Gen 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 Commenti
Torsten
il 24 Gen 2024
Modificato: Torsten
il 24 Gen 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 ?
Più risposte (0)
Vedere anche
Categorie
Scopri di più su Gas Dynamics in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!