Oplossen differentiaalvergelijking

[Igor Batoukhtine]

Zie onderstaande differentiaalvergelijking. Deze wil ik dusdanig oplossen zodat ik vm = ... kan uitdrukken
 
[sharedmedia=core:attachments:21425]
 

[Professor Puntje]

De differentiaalvergelijking zoals die er nu staat roept vooral veel vragen op. Het is handig de wiskundige vorm van die vergelijking duidelijk te maken. Dan heb je ook kans dat iemand met een wiskundige oplossing komt, zonder dat hij zich eerst in het vakgebied van de baggerpompen hoeft te verdiepen.

 

Van belang is het volgende:

 

- Wat zijn in de vergelijking constanten? Geef die aan als (bijvoorbeeld): a, b, c, etc. 

 

- Wat zijn in de vergelijking de onafhankelijke variabelen? Kennelijk in elk geval t, maar zijn er nog meer?

- Wat zijn in de vergelijking de afhankelijke variabelen? Kennelijk ρ_m en v_m ? Schrijf die dan als y en x. Maar waarvoor staan die streepjes erboven dan?

[Marko]

 
Zie onderstaande differentiaalvergelijking. Deze wil ik dusdanig oplossen zodat ik vm = ... kan uitdrukken
 

 
Dit gaat een lastige uitdaging worden. De dichtheden zijn kennelijk afhankelijk van de tijd, maar het precieze verband is niet gegeven, dus een algemeen bruikbare uitdrukking zal er op deze manier niet uitrollen. Het is niet eens duidelijk of er überhaupt een vergelijking is die dat verband geeft. Dit los van de vraag of die vergelijking tot een integreerbare DV zou leiden.
 
Waarom los je dit niet numeriek op? Je hebt Matlab en Simulink dus dat is een piece of cake. 
 
Symbolen gaan hernoemen is overigens de verkeerde weg. Of er a staat of f maakt voor de uitwerking natuurlijk geen klap uit, maar de hele wereld duidt een frictiefactor aan met f, en niemand met a. 

[Igor Batoukhtine]

Alles is constant behalve pm en vm. De streep erboven bij vm geeft aan dat het om het gemiddelde van vm over de gehele pijpleiding gaat, terwijl vm zonder streep lokaal zal zijn (gok ik). Is er daarnaast geen beginsituatie nodig (vm(0) = ?). De snelheid op tijdstip 0 is in ieder geval 0, en de dichtheid pm op tijdstip 0 is gelijk aan pf = 1025 kg/m². Deze loopt langzaam op tot 1200 kg/m². De implementatie heb ik in simulink geprobeerd (zie onderstaande afbeelding).
 

Diff verg simulink 1244 keer bekeken

 
Als initial condition (dus vm(0)) heb ik dan 0 ingevoerd. Is bovenstaande juist of mis ik wat? bij t --> oneindig gaat vm(t) namelijk naar - oneindig.. en dit is raar, want op een gegeven moment is de pm wel constant (na 1250 seconden). Zie ook onderstaande foutmelding:
 
An error occurred while running the simulation and the simulation was terminated
 
Derivative of state '1' in block 'Simulatie_van_Wijk/Berekenen geleverde druk pomp en snelheid (6.10 en 7.1 v Wijk)/Berekenen snelheid ( 7.1 v Wijk)/Integrator' at time 188.89237727852208 is not finite. The simulation will be stopped. There may be a singularity in the solution.  If not, try reducing the step size (either by reducing the fixed step size or by tightening the error tolerances)

[Marko]

In je vergelijking staan een paar onduidelijkheden.
 
Er staat een paar keer een gemiddelde dichtheid, afhankelijk van t, en een keer een dichtheid zonder streep, maar ook zonder t. 
Als je ervan uitgaat dat die dichtheid steeds constant is, dan kun je zoals je nu doet steeds dezelfde waarde gebruiken.
 
Alleen snap ik dan je opmerking "want op een gegeven moment is de pm wel constant" niet.
 
Sowieso zou ik je dringend aanraden om dichtheid weer te geven als ρ en p te bewaren voor druk (nog liever een hoofdletter om elke verwarring te voorkomen). In Simulink kun je de Griekse letter ρ niet ingeven, maar niets weerhoudt je ervan dit in te typen als rho.
 
Dichtheden worden normaal in kg/m³ gegeven. Soms zie je het in kg/m als het om lineaire dichtheid gaat, bijvoorbeeld bij vezels. Een dichtheid in massa per oppervlak komt ook wel voor, maar dan voor vellen papier of lappen stof. Ga eens na wat je in deze vergelijking nodig hebt. 
 
Ga ook na wat je nu eigenlijk voor L moet invullen, of anders gezegd: voor welke situatie deze DV geldt. 
 
Met andere woorden: zorg dat je eerst helder krijgt wat deze vergelijking beschrijft en wat de oorsprong/achtergrond van iedere term in deze vergelijking is. Dan zul je denk ik ook gaan begrijpen waarom je snelheid in deze situatie naar min oneindig gaat.

[Igor Batoukhtine]

Het is verholpen, heb in Simulink nu een werkende regeling (soortvan) waarbij de vloeistofsnelheid wordt geregeld adv opgegeven vergelijking. De 1e term aan de rechterzijde van de vergelijking is gelijk aan de geleverde druk van de pomp, de 2e term is gelijk aan het drukverlies als gevolg van de verticale hoogte die moet worden overbrugd en de 3e term is het drukverlies als gevolg van wrijving in de leiding.

Het punt dat de snelheid naar min oneindig ging was helaas een foute instelling. De dichtheid is idd in kg/m3, wss een typfout hierboven.

Nou heb ik wel voor de rho,m in de laatste term (die blijkbaar constant is daar) de wisselende rho,m(t) gebruikt, omdat de dichtheid niet constant is dus. Dus ik gebruik overal de term rho,m(t) ipv de constant rho,m.

Als het goed is is alles nu ook wat meer helder(?)

[Marko]

Het was me niet duidelijk waarom je dichtheid veranderde, omdat ik die alleen als input in het model zag, en aannam dat deze daardoor vast lag, maar ik heb dat gewoon niet goed gelezen; jouw input was dat die dichtheid langzaam toenam, en dat kan natuurlijk.
 
Je omschrijving van de termen is juist. De tweede term slaat inderdaad op de druk ten gevolge van de vloeistofkolom. Is die hoogte daadwerkelijk 5000 m? 
 
 
Overigens wordt er met wat je hier schetst nog helemaal niets geregeld. Wat je nu in Simulink hebt staan is het dynamische model van de vloeistofstroming bij een pomp die een vaste druk levert.

[Professor Puntje]

Heb je het verschil tussen de gemiddelde (gestreepte) en de lokale waarden van v_m(t) en ρ_m(t) verwaarloosd?
 
Beschouw je ρ_m = ρ_m(t) als een ons reeds gegeven functie van de tijd?

[Igor Batoukhtine]

 
 
Je omschrijving van de termen is juist. De tweede term slaat inderdaad op de druk ten gevolge van de vloeistofkolom. Is die hoogte daadwerkelijk 5000 m? 

In het proefschrift wordt 5000 m gebruikt, maar in mijn geval gaat het om +- 200 meter. 
 

 
 
Overigens wordt er met wat je hier schetst nog helemaal niets geregeld. Wat je nu in Simulink hebt staan is het dynamische model van de vloeistofstroming bij een pomp die een vaste druk levert.

Hieronder heb ik mijn totale regeling zoals in Simulink.
Eerst een overzicht van het dynamische gedrag van de vloeistofstroming.

Diff verg simulink 1237 keer bekeken

Dit heb ik vervolgens in een subsysteem geplaatst:

Berekening druk 1237 keer bekeken

Wat je hierboven ziet is het volgende: De druk P_e wordt geregeld door de PID regelaar (voor het gemak de waarde Y), waarbij de druk varieert tussen 0<P_e<P_e,fmax. Dit wordt gerealiseerd door de volgende vergelijking:
P_e = min(1,max(0,Y))*(rho,m/rho,f)*P_e,fmax.Hierdoor zal, als de PID regelaar een waarde < 1 uitstuurt, de druk worden geregeld (eigenlijk het toerental gevarieerd) volgens de PID regelaar, als de waarde van de PID regelaar >1 is dan zal de maximale beschikbare druk worden geleverd (maximaal toerental). Zie ook onderstaande afbeelding incl uitleg:

drukregeling 1237 keer bekeken

 
Dit heb ik dan ook weer in een subsysteem geplaatst, waardoor onderstaande afbeelding de regeling toont:

Regeling vloeistofsnelheid 1237 keer bekeken

 
De inertie van de pomp toont het opstartgedrag van de pomp. Dit kan wel of niet worden meegenomen in het model, waarbij de druk van de pomp wordt opgebouwd volgens volgende vergelijking:
 

Inertie pomp 1237 keer bekeken

 
Hierbij is t_e de opstarttijd van de pomp, waarbij de pomp bij maximaal toerental (en waterdraaien) zijn druk bereikt. 
 
Het resultaat van de regeling is in de volgende afbeeldingen te zien
voor 50 seconden looptijd:

Resultaat regeling PID 1237 keer bekeken

Voor 1300 seconden looptijd:

Resultaat regeling PID1 1237 keer bekeken

 
Wat ik alleen niet begrijp is dat de variatie van mijn D in de PID-regelaar geen verschil in uitsturing behaald. Eigenlijk voldoet een PI regelaar al, terwijl in het proefschrift wel degelijk een verschil is te zien tussen verschillende regelaars (overigens krijg ik dat oscillerende gedrag ook niet in mijn model):

PID regelaar proefschrift 1237 keer bekeken

 

 
Heb je het verschil tussen de gemiddelde (gestreepte) en de lokale waarden van v_m(t) en ρ_m(t) verwaarloosd?
 
Beschouw je ρ_m = ρ_m(t) als een ons reeds gegeven functie van de tijd?

 

Verschil inderdaad verwaarloosd en dus ρ_m = ρ_m(t)
 
Graag jullie mening

[Professor Puntje]

Van de gebruikte simulatie en van pompsystemen weet ik te weinig af. Ik ben vooral geïnteresseerd in de wiskundige openingsvraag van dit topic of de gegeven differentiaalvergelijking analytisch oplosbaar is. Om daarmee aan de slag te kunnen wil ik nog graag weten of ρ_m als een gegeven functie mag worden beschouwd. Met andere woorden: neem je al van tevoren aan dat ρ_m als functie van t een gegeven verloop heeft, en wil je vervolgens het verloop van v_m weten die bij dat verloop van ρ_m behoort? 

[Igor Batoukhtine]

Voor het gebruik van de simulatie mag daar vanuit worden gegaan inderdaad (in de realiteit zal het natuurlijk anders lopen en is het van te voren niet bekend hoe het verloop van de dichtheid is). Maar het is nu alleen van belang om te kijken hoe de vloeistofsnelheid reageert op de verandering van de dichtheid. Het verloop dat ik voor ρ_m gebruik is als volgt:
 
ρ_m varieert van 1025 kg/m³ tot 1200 kg/m² in 1250 seconden (lineair)

[Marko]

Je bent in ieder geval heel aardig op weg.
 
Hou in de gaten dat de dichtheid-met-streepje voor jouw situatie per definitie ongelijk is aan de dichtheid-zonder-steepje.
 
De term zonder streepje gaat immers over de pompdruk, en die is het gevolg van de dichtheid ter plekke en op dat moment. De andere termen hebben betrekking op de vloeistof in de buis, en dus draait het om de gemiddelde dichtheid over de hele buis.
 
Maar jij laat de dichtheid van de vloeistof die de pomp ingaat variëren, en de verblijftijd van de vloeistof is relatief lang (dus de dichtheid van pak 'm beet een minuut eerder heeft nog steeds invloed op de pomp beneden). 
 
Je zult dus iets moeten verzinnen voor dat verschil.

[Igor Batoukhtine]

Hmm, stel nou dat ik het volgende doe:
 

1. Ik neem een sampletijd van 1 seconde;
2. Ik neem de tijd die het mengsel erover doet om van de pomp naar einde leiding te komen (L/v_m= 200/4 = 50 s); 
3. Dit betekent dat er 50 samples worden genomen, waarover het gemiddelde van de dichtheid kan worden bepaald : Som rho(t) over 50 seconden / 50.
4. De bepaling voor het volgende gemiddelde zal zijn van rho(t) 1<t<51 etc.

Deze waarde vul ik dan in bij de rho met streepje? De rho zonder streepje houd ik dan op de rho die nu direct varieert.

[Professor Puntje]

Ik zit nog met een ander probleem. De differentiaalvergelijking heeft een partiële afgeleide van v_m (streep) naar t. Betekent dat dat die vergelijking enkel weergeeft hoe v_m (streep) varieert voor constante ρ_m?
 
Edit: Mogelijk wordt er bedoeld dat je enkel naar de tijd moet differentiëren in niet naar de plaats langs de pijpleiding? Hoewel dat gemiddelde over de pijpleiding daar weer onafhankelijk van zou moeten zijn...

[Igor Batoukhtine]

In het proefschrift staat het volgende:
 

Uitleg diff verg 1291 keer bekeken

 
Er wordt dus inderdaad uitgegaan van een constante ρ_m. Maar waarom dat is is mij nog onduidelijk, waar ik uiteindelijk naar toe wil (denk ik) is een variërende ρ_m in combinatie met een variërende ρ_m (met streep)
 
@Marko
 
Heb nu de volgende code voor het bepalen van het gemiddelde over een tijdsperiode:
 
t = 0:200;                      % Tijdsimulatie
rho = 1025;                   % Dichtheid water
inc = 0;
log = zeros(1,200);       % Array voor het loggen van de dichtheden
L = 200;                        % Leidinglengte
v = 4;                             % Leidingsnelheid;
ii = 0;                             % telvariabele voor het berekenen van gemiddelde rho
avgrho = zeros(1,150); % Definiëren van avgrho
 
for tijd = t
    inc = inc + 1;
    log(inc) = rho; % Opslaan van de dichtheden in een array
    rho = rho + 1;  % Lineaire verhoging van de dichtheid (1025 - 1225)
    reistijd = L/v;
 
    if (tijd >= reistijd)
        ii = ii + 1;
        avgrho(ii) = mean(log(inc-reistijd):log(inc));
    end
end