System and Control
AE 315,231
Assignment Number 4
Fluid system
The system required to establish in Matlab
The given parameters
A=5,
Rf=0.5,With
h0=0m
and tend=20s.
Part 1 Equation of motion
The equation of continuity
AV=constant
In other words
Mass per unit volume is constant
∆V
=Constant
∆t
∆V
Q=
∆t
Δx1 = v1Δt
Now, at the lower end of the pipe, the volume of the fluid that will flow into the pipe will be:
V = A1 Δx1 = A1 v1 Δt
It is known that mass (m) = Density (ρ) × Volume (V). So, the mass of the fluid in Δx1 region will
be:
Δm1= Density × Volume
=> Δm1 = ρ1A1v1Δt ——–(Equation 1)
Now, the mass flux has to be calculated at the lower end. Mass flux is simply defined as the mass
of the fluid per unit time passing through any cross-sectional area. For the lower end with cross-
sectional area A1, mass flux will be:
Δm1/Δt = ρ1A1v1 ——–(Equation 2)
Similarly, the mass flux at the upper end will be:
Δm2/Δt = ρ2A2v2 ——–(Equation 3)
Here, v2 is the velocity of the fluid through the upper end of the pipe i.e. through Δx2 , in Δt time
and A2, is the cross-sectional area of the upper end.
In this, the density of the fluid between the lower end of the pipe and the upper end of the pipe
remains the same with time as the flow is steady. So, the mass flux at the lower end of the pipe is
equal to the mass flux at the upper end of the pipe i.e. Equation 2 = Equation 3.
Thus,
ρ1A1v1 = ρ2A2v2 ——–(Equation 4)
This can be written in a more general form as:
ρ A v = constant
The equation proves the law of conservation of mass in fluid dynamics. Also, if the fluid is
incompressible, the density will remain constant for steady flow. So, ρ1 =ρ2.
Thus, Equation 4 can be now written as:
A1 v1 = A2 v2
This equation can be written in general form as:
A v = constant
Now, if R is the volume flow rate, the above equation can be expressed as:
R = A v = constant
Part 2 function of the flow
function [h v] = boundary(h,v, leftboundary, rightboundary,N,Choice)
% Sets boundary conditions (either zero gradient or reflective)
if Choice == 1 % Lax-Friedrichs
if leftboundary == 0 && rightboundary == 0
h = [h(1) h h(N)];
v = [v(1) v v(N)];
elseif leftboundary ==1 && rightboundary == 0
h = [h(1) h h(N)];
v = [-v(3) v v(N)];
elseif leftboundary == 0 && rightboundary == 1
h = [h(1) h h(N)];
v = [v(1) v -v(N-2)];
elseif leftboundary == 1 && rightboundary == 1
h = [h(1) h h(N)];
v = [-v(3) v -v(N-2)];
end
elseif Choice == 4 % Adams Average
if leftboundary == 0 && rightboundary == 0
h = [h(1) h h(N)];
v = [v(1) v v(N)];
elseif leftboundary ==1 && rightboundary == 0
h = [h(1) h h(N)];
v = [-v(3) v v(N)];
elseif leftboundary == 0 && rightboundary == 1
h = [h(1) h h(N)];
v = [v(1) v -v(N-2)];
elseif leftboundary == 1 && rightboundary == 1
h = [h(1) h h(N)];
v = [-v(3) v -v(N-2)];
end
elseif Choice == 2 % MacCormack
if leftboundary == 0 && rightboundary == 0
h(1) = h(2);
h(N+2) = h(N+1);
v(1) = v(2);
v(N+2) = v(N+1);
elseif leftboundary ==1 && rightboundary == 0
h(1) = h(2);
h(N+2) = h(N+1);
v(1) = -v(3);
v(N+2) = v(N+1);
elseif leftboundary == 0 && rightboundary == 1
h(1) = h(2);
h(N+2) = h(N+1);
v(1) = v(2);
v(N+2) = -v(N);
else
h(1) = h(2);
h(N+2) = h(N+1);
v(1) = -v(3);
v(N+2) = -v(N);
end
elseif Choice == 3 % Lax-Wendroff
if leftboundary == 0 && rightboundary == 0
h(1) = h(2);
h(N+2) = h(N+1);
v(1) = v(2);
v(N+2) = v(N+1);
elseif leftboundary ==1 && rightboundary == 0
h(1) = h(2);
h(N+2) = h(N+1);
v(1) = -v(3);
v(N+2) = v(N+1);
elseif leftboundary == 0 && rightboundary == 1
h(1) = h(2);
h(N+2) = h(N+1);
v(1) = v(2);
v(N+2) = -v(N);
else
h(1) = h(2);
h(N+2) = h(N+1);
v(1) = -v(3);
v(N+2) = -v(N);
end
Part 3 Laplace Transform
function[h,v] = redependent(U1,U2,N, leftboundary, rightboundary,Choice)
if Choice == 1
h = U1(2:N+1); % Calculates new h and v values for Lax-F
v = U2(2:N+1)./U1(2:N+1); % (Old ghost values are deleted)
elseif Choice == 4
h = U1(2:N+1); % Calculates new h and v values for Adams Average
v = U2(2:N+1)./U1(2:N+1);
elseif Choice == 2
h = U1; % < Calculates new h and v values for MacCormack
v = U2./U1;
elseif Choice == 3
h = U1; % < Calculates new h and v values for lax-Wendroff
v = U2./U1;
end
[h, v] = boundary(h,v, leftboundary, rightboundary,N,Choice); % Recalculates boundary values
end
%=========== Function to calculate timestep ==============================
function [dt, max1, max2, maxoverall] = timestepinequality(dx,s,v,g,h)
max1 = max(abs(v + sqrt(g*h)));
max2 = max(abs(v - sqrt(g*h)));
maxoverall = max(max1,max2); % Calculates maximum wave speed
dt = s * dx / maxoverall; % Calculates timestep
end
Part4 Simulink function
% example_3 - inversion of a fractional F(s) in symbolic form
clear, close all
syms D1 alfa1 R1 s
% parameters of the network with constant phase element
I=0.25; Rs=0.1; R1=100; D1=1; alfa1=-0.7;
% F(s) in symbolic form
F1=I*(Rs+(R1*D1*s^alfa1/(R1+D1*s^alfa1)))*(1-exp(-4000*s))/s;
% F(s) as a string
F1=char(F1);
% parameters of the fractional control system
k=20.5; a1=3.7343; alfa1=1.15; a2=0.8; alfa2=2.2; a3=0.5; alfa3=0.9;
F2=(k+a1*s^alfa1)/(k+1+a1*s^alfa1+a2*s^alfa2+a3*s^alfa3)/s; % F(s) in symbolic form
F2=char(F2); % F(s) as a string
[t1,ft1]=INVLAP(F1,0.01,1e4,1000,6,39,89);
[t2,ft2]=INVLAP(F2,0.01,5,200);
figure(4)
set(4,'color','white')
subplot(2,1,1)
plot(t1,ft1), grid on, zoom on
xlabel('t [s]'), ylabel('f(t)')
title('response to the input current impulse')
subplot(2,1,2)
plot(t2,ft2), grid on, zoom
xlabel('t [s]'), ylabel('f(t)')
title('step response of a fractional control system')
