DAE Solver for Height vs Time plot with variable

Question

0 votes

Hi,

I have a Height vs Time plot for a tank which has 1 inlet and 2 outlets. I want to plot Height v Time with multiple lines with varying alpha values. I have created the code using a DAE solver to plot the line for 1 alpha value but am unsure how to add multiple lines to show different values. The alpha values will be 0.01:0.01:0.2.

Any help will be much appreciated.

function dxdt = G_20_tankf4(t,x)
% G_20_tankf4: Tank level model
%Derivative function solution for Method 4
% Unpack the state vector into local variables
h = x(1);
F1 = x(2);
F2 = x(3);
P2 = x(4);
F3 = x(5);
% Tank model parameters and chonstants
A = 8.0;             %M^2
CV1 = 2;             %(m^3/min)/(kPa^0.5)
CV2 = 3*h;           %(m^3/min)/(kPa^0.5)
P0 = 100;            %kPa
P1 = P1function(t);  %kPa
P3 = 100;            %kPa
rho = 1000;          %kg/m^3
g = 9.81;            %m/s^2
alpha = 0.01;
% Evaluate the derivative
dhdt = 1/A*(F1-F2-F3);
% Calculate the residulas of the algebriac equations
r1 = F1 - CV1*sqrt(P1-P2);
r2 = F2 - CV2*sqrt(P2-P3);
r3 = P2 - (P0 + rho*g*h/1000);
r4 = F3 - alpha*g*h;
% Pack the resuts into dx/dt
dxdt = [dhdt r1 r2 r3 r4]';
% G_20_tanksim4: Method 4
% Set extra parameters
tf = 72;         %h
h0 = 2;          %m
%Work out conistent initial condition for F1, F2, F3, P2 and CV2
CV1 =  2;            % (m^3/h)/kPa^0.5
P0 = 100;            % kPa
P3 = 100;            % kPa
rho = 1000;          % kg/m^3
g = 9.81;            % m/s^2
P1 = P1function(0);  % kPa
CV2 = 3*h0;          % (m^3/h)/kPa^0.5
alpha = 0.01;
P20 = P0 + rho*g*h0/1000;
F10 = CV1*sqrt(P1-P20);
F20 = CV2*sqrt(P20-P3);
F30 = alpha*g*h0;
% Construct initial state vector
x0 = [h0, F10, F20, P20, F30];
% Set mass matrix (one differential variable followed by three algebraics)
M = [1 0 0 0 0; 0 0 0 0 0; 0 0 0 0 0; 0 0 0 0 0;0 0 0 0 0];
% Indicate DAE problem by specifying singular mass matrix and supplying M
% (and also specify the tolerance as before)
opts = odeset('Reltol',1e-5,'MassSingular','yes','Mass',M);
% Solve the DAE
[t,x] = ode15s(@G_20_tankf4,[0 tf],x0,opts);
% Extract the level from x
h = x(:,1);