# Problem with NDsolve in a simple hydraulic simulation

Posted 3 months ago
638 Views
|
7 Replies
|
3 Total Likes
|
 Dear Community,I try to solve a simple physics problem: two water tanks are connected with a pipe, and a pump pumps water uphill. During the solution however NDSolve gives some strange error message, and I have no idea, what's the problem. In an other application this case runs without any problem. I've attached the notebook, the correct results are visible at the bottom of it. Could you pls. have a look, what this error message means, and what should be corrected?Tx very much in advance, best regards Andras Attachments:
7 Replies
Sort By:
Posted 3 months ago
 Andras,Your problem is with the Friction function. You have a discontinuity at 0 flow and using Return[] is not correct. I would change it to friction[ q_ ] := Module[ { rey = 0. , fl = 0. , ft = 0. } , If[q >= 0.001, rey = ( 4. * q * \[Rho] ) / ( N[ \[Pi] ] * dia * \[Mu] ) ; fl = 64. / rey ; ft = 1. /(( 1.8 * Log10[ 6.9/rey + (ee /(3.7 * dia))^1.11 ] )^2) ; Piecewise[ {{ (fl), ( 0. < rey && rey < 2000. ) } , { ( fl + ( ft - fl ) * ( rey - 2000.) / 2000. ), ( 2000. <= rey && rey < 4000. ) } , { ( ft), ( 4000. <= rey ) }} ] , friction[0.001]] ] or if you want to use Which: friction[ q_ ] := Module[ { rey = 0. , fl = 0. , ft = 0. } , If[q>=0.001, rey = ( 4. * q * \[Rho] ) / ( N[ \[Pi] ] * dia * \[Mu] ) ; fl = 64. / rey ; ft = 1. /(( 1.8 * Log10[ 6.9/rey + (ee /(3.7 * dia))^1.11 ] )^2) ; Which[ ( 0. < rey && rey < 2000. ) , ( fl ) , ( 2000. <= rey && rey < 4000. ) , ( fl + ( ft - fl ) * ( rey - 2000.) / 2000. ) , ( 4000. <= rey ) , ( ft ) ] (* Which[ *), friction[0.001]] ] Either of these works. Your plot is still not matching the goal plots you included, however, you can now debug from here. Note your initial conditions are not correct. For example, Tank 2 seems to start at 2.5 Regards,Neil
Posted 3 months ago
 Dear Neil,Tx for the kind help. You're right, I made a mistake in the definition of the friction function - now it works fine.Thank you very much, best regardsAndras
Posted 3 months ago
 Dear Neil,May I ask for a little more help? I've strongly simplified the model, and now I would like to shut down fluid flow in the 2000 - 2500 s simulation interval ( i.e. q1(t) = 0. and consequently fluid levels h1 and h2 must stay constant ). I try to achieve this with the WhenEvent function in NDSolve, but I get some convergence problems. Do you possibly know how to make it correctly? I've attached the notebook.Tx in advance, best regards Andras Attachments:
Posted 3 months ago
 Andras,What you are asking to do may not make physical sense -- you can't just instantaneously shut down the tanks -- they oscillate because the physics says they should oscillate. You can, however, model whatever you want including a system with discrete behavior (such as a sudden alteration of the system -- closing a valve, etc.). To do this, I would introduce a discrete variable that is either 1 or 0 and affectively "turn off" your system when you hit a condition. For example, odesys = { D[ h1[t] , t ] == ( Minus[ q1[t] ] / a1 )* d[t], D[ h2[t] , t ] == ( q1[t] / a2 ) * d[t] , D[ q1[t] , t ] == pipefactor * ( h1[t] - h2[t] ) + pumphead[ t ] - frictionfactor * Abs[ q1[t] ] * q1[t] , h1[0] == 1.5 , h2[0] == 1.5 , q1[0] == 0., d[0] == 1 } ; Where d[t] is the discrete variable, set it to 1 initially.Now integrate it with sol = NDSolve[ {odesys, WhenEvent[ (h1[t]>1.5), d[t] ->0]} , { h1 , h2 , q1 } , { t , 0. , 3000. } , DiscreteVariables->d[t]] ; When your condition is hit (in this case When h1[t] crosses through the neutral point at 1.5, it shuts a valve and the two tank height derivative equations go to zero.I hope this helps.Regards,Neil
Posted 3 months ago
 In rereading your file, I am now not sure if you wanted to set the valve at time t>2000 instead of when the tanks cross. You can do the same thing based on time: sol = NDSolve[ {odesys, WhenEvent[ 2000. <= t <= 2500., d[t] -> 0, "LocationMethod" -> "LinearInterpolation"]} , { h1 , h2 , q1 } , { t , 0. , 3000. } , DiscreteVariables -> d[t]] ; Note that I had to add the LocationMethod to the WhenEvent for a time based event because the event detection had trouble finding the crossover (although it still gave the right answer with some warnings). I eliminated the warnings by changing the method used to find the event.Regards
 Andras,If you want that then I would use two events -- one turns off the valve and the other turns it on. You get the Plot you want. sol = NDSolve[ { odesys, WhenEvent[ 2000. <= t <= 2500., d[t] -> 0. , "LocationMethod" -> "LinearInterpolation" ], WhenEvent[ t > 2500., d[t] -> 1. , "LocationMethod" -> "LinearInterpolation" ] } , { h1 , h2 , q1 } , { t , 0. , 3000. } , DiscreteVariables -> d[t]] ; Regards