Simulation der Wasserhöhe in einem Tank mit Zu- und Abfluss

WICHTIG: Der Betrieb von goMatlab.de wird privat finanziert fortgesetzt. - Mehr Infos...

Partner:

Forum

Option

[Erweitert]

• Diese Seite per Mail weiterempfehlen

Gehe zu:

rebbalmat0887

Forum-Anfänger


	Beiträge: 33

	Anmeldedatum: 08.10.15

	Wohnort: ---

	Version: ---

Verfasst am: 08.10.2015, 14:47 Titel: Simulation der Wasserhöhe in einem Tank mit Zu- und Abfluss

Hallo,

ich möchte den Wasserstand in einem Tank simulieren, der einen konstanten
Zufluss (durch Punmpe) und ein Ausflussrohr besitzt. Da der ausfließende
Wasservolumenstrom von der Spiegelhöhe anhängt, muss es ein Gleichgewicht
zwischen Zu- und Ablauf geben, bei dem letztendlich der Wasserspiegel konstant
bleibt.
Verifizieren will ich das ganze mit Messwerte, die mit dem Plot dann
übereinstimmen sollten, was aber nicht der Fall ist. Deswegen möchte ich fragen,
ob ich die DGL evtl. falsch implementiert habe.

Ich möchte über der Zeit simulieren und habe folgende DGL aufgestellt:

Code:

dh/dt = 1/A * [Qzu - Qab] = 1/Atank2 * [Qin - Aout * sqrt(2*g*h)]

mit:
Atank2 - Grundfläche des Tanks
Qin - Zufließender Volumenstrom (const: ~ 0,00010789 m³/s)
Aout - Querschnittsfläche des Ausflussrohres
g - Erdbeschleunigung
h - Spiegelhöhe

Funktion ohne Link?

Zur Realisierung in Matlab:

Code:

%% Simulation Water Tank Model
close all
clear all
clc

%Initialize Parameter

k = 0.000063091; % [m³/s*GPM]
QinGPM = 1.71; % 1.7 GPM
Qin = QinGPM * k; % Umrechnen in m³/s
Atank1 = ((8+1/2+1/8)*2.54)*((10+1/4+1/8+1/16)*2.54); % [cm²]
Atank2 = Atank1 / 10^4; % Umrechnen in m²
g = 9.81; % [m/²]
dout = 1.5; % Rohrdurchmesser für Ablauf in [mm]
Aout = (pi*(dout/1000)^2)/4; % Rohrfläche für Ablauf in [m²]

% Solveing the Differential Equation in Symbolic Manner:
fcn = dsolve('Dh = 1/Atank2 *(Qin - sqrt(2*g*h))','h(0) = hstart','t');

% Solveing the Differential Equation Numerically:
tstart = 0;
tend = 1.0 * 10^6;
tspan = [tstart tend];
hstart = 0; % Startspiegelhöhe in [m]

[t,h] = ode45(@dglWassertank,tspan,hstart);

hgrenz = (2 * Qin^2) / (pi^2 * (dout/1000)^4 * g);
display(['Water Balanced Hight: ' num2str(hgrenz) ' m']);

%Messwerte aus Excel einlesen
dateiname = uigetfile({'*.xlsx';'*.xls'}); % Dialog zur Auswahl des Filname
data=xlsread(dateiname,1,'a1:d841');

t_s = data(:,1);
t_min = data(:,2);
h_inch = data(:,3);
h_cm = data(:,4);

% Plotting the Results:

figure()
subplot(2,1,1);
plot(t,h);
grid on
xlabel('Time [s]');
ylabel('Height of Water [m]');
title('Watertank Simulation ODE','FontSize',12);

subplot(2,1,2);
plot(t_s,h_cm);
grid on
xlabel('Time [s]');
ylabel('Height of Water [m]');
title('Measurement','FontSize',12);

Funktion ohne Link?

Die Funktion "dglWassertank" habe ich extra als Funktion gespeichert.

Code:

function dh_dt = dglWassertank(t,h)

k = 0.000063091; % [m³/s*GPM]
QinGPM = 1.71; % 1.7 GPM
Qin = QinGPM * k; % Umrechnen in m³/s
Atank1 = ((8+1/2+1/8)*2.54)*((10+1/4+1/8+1/16)*2.54); % [cm³]
Atank2 = Atank1 / 10^4; % Umrechnen in m²
g = 9.81; % [m/²]
dout = 1.5; % Rohrdurchmesser für Ablauf in [mm]
Aout = (pi*(dout/1000)^2)/4; % Rohrfläche für Ablauf in [m²]

dh_dt = 1/Atank2 *(Qin - 2*Aout*sqrt((2*g*h)));

end

Funktion ohne Link?

Ich hoffe, es kann mir jemand helfen. Vielen Dank im Voraus für Eure Mühen!

messwerte.xls

Beschreibung:

1. Spalte: t [s]
2. Spalte: t [min]
3. Spalte: h [inch]
4. Spalte: h [cm]

Download

Dateiname:

messwerte.xls

Dateigröße:

147 KB

Heruntergeladen:

395 mal

Winkow

Moderator


	Beiträge: 3.842

	Anmeldedatum: 04.11.11

	Wohnort: Dresden

	Version: R2014a 2015a

Verfasst am: 08.10.2015, 15:29 Titel:

mir fällt so erstmal nicht viel auf. sicher das deine ausgangswerte richtig sind? 1.5mm ausgangsdurchmesser ist recht wenig. sicher das das nicht in cm ist oder so?
woher kommt die 2 vor Aout? ich hätte da ehr mit 0.6 gerechnet oder so.
warum plottest du die höhe in cm und schreibst dann meter ran?
_________________

richtig Fragen

rebbalmat0887

Themenstarter

Forum-Anfänger


	Beiträge: 33

	Anmeldedatum: 08.10.15

	Wohnort: ---

	Version: ---

Verfasst am: 08.10.2015, 15:56 Titel:

Hallo,

vielen Dank für deinen Einwand. Ich hab mich geirrt, es sind 15,3mm. Der Faktor 2
vor der der Fläche Aout kommt daher, weil ich zwei Ausflussrohre selben Durch-
messers habe. Einen Korrekturfaktor hatte ich erstmal noch garnicht berücktsichtigt.
Ja, das mit der Achsenbeschriftung hatte ich nach dem Abschicken bemerkt und ge-
ändert.
Es sieht auf jeden Fall schon mal besser aus, aber immer noch nicht perfekt. Aber
ich kann erstmal die Fehlerquelle 'falsche ode45-Implementierung' ausschließen.
Das ist schonmal viel Wert. Dankeschön!
Als Korrekturfakor ist 0,6 okay?

Winkow

Moderator


	Beiträge: 3.842

	Anmeldedatum: 04.11.11

	Wohnort: Dresden

	Version: R2014a 2015a

Verfasst am: 08.10.2015, 16:13 Titel:

Zitat:

Als Korrekturfakor ist 0,6 okay?

hab ich auch nur aus wiki übernommen.
https://de.wikipedia.org/wiki/Ausflussgeschwindigkeit der hängt denke ich von einigen faktoren ab. bei 0.1 kommst du zb ehr auf deine werte Smile

kann man vielleicht mal genauer untersuchen was jetzt für deinen fall zutrifft. kann aber auch sein das etwas falsch is das wir beide hier übersehen Smile

_________________

richtig Fragen

rebbalmat0887

Themenstarter

Forum-Anfänger


	Beiträge: 33

	Anmeldedatum: 08.10.15

	Wohnort: ---

	Version: ---

Verfasst am: 09.10.2015, 11:56 Titel:

Hallo,

danke erstmal für deine Hilfe. Ich habe meine Simulation noch um den Verlustterm

Code:

dpv = 0.5*roh*v^2*[Zeta_einström + Zeta_ausström + (lambda*lrohr / drohr)]

Funktion ohne Link?

ergänzt. Mit den anderen kleineren Korrekturen sieht die Simulation jetzt erstmal
ganz gut aus. ABER: mein simuliertes System ist zu schnell, verglichen mit den Mess-
werten. Ich habe mal probiert die Grundfläche des Tanks zu vergrößern. Von 580cm^2
auf 2500 cm^2 (ca das Vierfache).
Gibt es einen Aspekt, den ich möglicherweise vergessen habe?

rebbalmat0887

Themenstarter

Forum-Anfänger


	Beiträge: 33

	Anmeldedatum: 08.10.15

	Wohnort: ---

	Version: ---

Verfasst am: 09.10.2015, 11:58 Titel:

Sorry,

habe vergessen das aktuelle Script anzugeben.

m-File:

Code:

%% Simulation Water Tank Model
close all
clear all
clc

%Initialize Parameter
k = 0.000063091; % [m³/s*GPM]
QinGPM = 1.71; % 1.7 GPM
Qin = QinGPM * k; % Umrechnen in m³/s
Atank1 = ((8+1/2+1/8)*2.54)*((10+1/4+1/8+1/16)*2.54); % [cm³]
Atank2 = Atank1 / 10^4; % Umrechnen in m²
g = 9.81; % [m/²]
dout = 15.3; % Rohrdurchmesser für Ablauf in [mm]
Aout = (pi*(dout/1000)^2)/4; % Rohrfläche für Ablauf in [m²]
kappa = 0.6; % Korrekturfaktor
lambda = 0.02;
zeta = 2;
lrohr = (22*2.54)*10;

% Solveing the Differential Equation Numerically:
tstart = 0;
tend = 4500;
tspan = [tstart tend];
hstart = 0; % Startspiegelhöhe in [m]

[t,h] = ode45(@dglWassertank,tspan,hstart);

%*(1+lambda*lrohr/dout+zeta)

hgrenz = (2 * Qin^2 *(1+lambda*lrohr/dout+zeta)) / (pi^2 * (dout/1000)^4 * g * kappa^2);
display(['Water Balanced Hight: ' num2str(hgrenz*100) ' cm']);

%Messwerte aus Excel einlesen
dateiname = uigetfile({'*.xlsx';'*.xls'}); % Dialog zur Auswahl des Filname
data=xlsread(dateiname,1,'a1:d841');
t_s = data(:,1);
t_min = data(:,2);
h_inch = data(:,3);
h_cm = data(:,4);

% Plotting the Results:
% (1) Beide Graphen in einem Plot
figure(1)
plot(t,h*100);
hold on
grid on
plot(t_s,h_cm);
axis([0 2500 0 5]);

Funktion ohne Link?

ode45-function:

Code:

function dh_dt = dglWassertank(t,h)

k = 0.000063091; % [m^3/s*GPM]
QinGPM = 1.71; % 1.7 GPM
Qin = QinGPM * k; % Umrechnen in m^3/s
Atank1 = ((8+1/2+1/8)*2.54)*((10+1/4+1/8+1/16)*2.54); % [cm³]
Atank2 = Atank1 / 10^4; % Umrechnen in m^2
g = 9.81; % [m/s^2]
dout = 15.3; % Rohrdurchmesser für Ablauf in [mm]
Aout = (pi*(dout/1000)^2)/4; % Rohrfläche für Ablauf in [m²]
kappa = 0.6; % Korrekturfaktor
lambda = 0.02;
zeta = 2;
lrohr = (22*2.54)*10;

dh_dt = 1/Atank2 *(Qin - kappa * 2 * Aout * sqrt((2*g*h)/(1+lambda*lrohr/dout+zeta)));

end

Funktion ohne Link?

rebbalmat0887

Themenstarter

Forum-Anfänger


	Beiträge: 33

	Anmeldedatum: 08.10.15

	Wohnort: ---

	Version: ---

Verfasst am: 12.10.2015, 10:06 Titel:

Ok,

also danke für deine Hilfe. Ich denke, dass das Problem nicht Matlab-basierend ist.
Ich habe aber leider auch keine Ahnung, wo ich noch suchen soll. Trotzdem vielen
Dank für die Hilfe.

Grüße

rebbalmat0887

Themenstarter

Forum-Anfänger


	Beiträge: 33

	Anmeldedatum: 08.10.15

	Wohnort: ---

	Version: ---

Verfasst am: 13.10.2015, 10:54 Titel:

Hallo,

ich habe doch noch eine Idee, woran es liegen könnte. Die Rohrreibungszahl 'Lambda'
ist nicht konstant. Sie muss iterativ berechnet werden, da sie von der aktuellen Höhe
im Wassertank abhängt. Ich werde hier die folgen Schritte posten und hoffe, ihr könnt
mir bei der Implementierung helfen. Ich weiß nämlich nicht, wo genau die Iteration
hinkommt. In die externe Funktion oder ins Hauptskript?

Code:

Iteration:

nue = 1.139 * 10^-6; % kinematische Viskosität Wasser bei 15°C
lambdaStart = 0.02;
lambdaNeu = 0;
d = 0.0153; % Rohrdurchmesser [m]
l = 0.558; % Rohrlänge [m]
g = 9.81; % Gravitationskonstante [m/s²]
zeta = 1; % Einlaufverlust in Rohr [1]

while abs(lambdaStart - lambdaNeu) > 0.0001;

v = sqrt( (2 * g * h) / (1 + zeta + lambdaStart*l/d) ); % Geschwindigkeit
% am Rohraustritt
% berechnen [m/s]
Re = (v * d) / nue; % Berechnen der Reynoldszahl [1]
lambdaNeu = 0.316 / nthroot(Re,4); % neues Lambda berechnen [1]

end

Funktion ohne Link?

Ich denke, ein großes Problem ist, dass ich die DGL mit ode45 löse und daher
keine feste Schrittweite habe. Die bräuchte ich aber um mir zu einem bestimmten
Zeitschritt iterativ erst 'Lambda' und dann die Wasserhöhe zu berechnen.

Könnt ihr mir dabei helfen? Dankeschön!

rebbalmat0887

Themenstarter

Forum-Anfänger


	Beiträge: 33

	Anmeldedatum: 08.10.15

	Wohnort: ---

	Version: ---

Verfasst am: 13.10.2015, 15:35 Titel:

Ok,

also die Iteration zur Berechnung von 'Lambda' habe ich hinbekommen.

Code:

clear all
close all
clc

nue = 1.139 * 10^-6; % kinematische Viskosität Wasser bei 15°C
d = 0.0153; % Rohrdurchmesser [m]
l = 0.558; % Rohrlänge [m]
g = 9.81; % Gravitationskonstante [m/s²]
zeta = 1; % Einlaufverlust in Rohr [1]
h = 0.001;
lambda = 0.02;
diff = lambda;

while diff > 0.0001

lambdaAlt = lambda;
v = sqrt( (2 * g * h) / (1 + zeta + lambda*l/d) ); % Geschwindigkeit
% am Rohraustritt
% berechnen [m/s]
Re = (v * d) / nue; % Berechnen der Reynoldszahl [1]
lambda = 0.316 / nthroot(Re,4); % neues Lambda berechnen [1]

diff = abs(lambda - lambdaAlt);
end
display(['Lambda: ' num2str(lambda)]);

Funktion ohne Link?

Und es ist auch möglich, sich bei ode45 statt nur Anfangs- und Endwert einen
Vektor mit fester Schrittweite vorzugeben (z.Bsp: 0.02s).
Jetzt brauche ich aber eine verschachtelte Schleife, die mir zu jedem Zeitschritt
eine Iteration zur Berechnung von 'Lambda' und h(t,lambda) ausführt.
Kann mir da echt keiner helfen?

Grüße


Einstellungen und Berechtigungen
Beiträge der letzten Zeit anzeigen:	Du kannst Beiträge in dieses Forum schreiben. Du kannst auf Beiträge in diesem Forum antworten. Du kannst deine Beiträge in diesem Forum nicht bearbeiten. Du kannst deine Beiträge in diesem Forum nicht löschen. Du kannst an Umfragen in diesem Forum nicht mitmachen. Du kannst Dateien in diesem Forum posten Du kannst Dateien in diesem Forum herunterladen

Impressum | Nutzungsbedingungen | Datenschutz | FAQ |

RSS

Hosted by:

Copyright © 2007 - 2025 goMatlab.de | Dies ist keine offizielle Website der Firma The Mathworks

MATLAB, Simulink, Stateflow, Handle Graphics, Real-Time Workshop, SimBiology, SimHydraulics, SimEvents, and xPC TargetBox are registered trademarks and The MathWorks, the L-shaped membrane logo, and Embedded MATLAB are trademarks of The MathWorks, Inc.