Group Abstract

Message Boards

WOLFRAM COMMUNITY

6.2K Views

7 Replies

3 Total Likes

View groups...

Follow this post

Share this post:

GROUPS:

Problem with NDsolve in a simple hydraulic simulation

agilicz

Posted 5 years ago

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: PumpingWater.nb

POSTED BY: agilicz

7 Replies

Sort By:

Neil Singer

Neil Singer, AC Kinetics, Inc.

Posted 5 years 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

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 BY: Neil Singer

Neil Singer

Neil Singer, AC Kinetics, Inc.

Posted 5 years ago

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

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

POSTED BY: Neil Singer

agilicz

Posted 5 years ago

Dear Neil, Tx for the kind reply. Yes indeed, I would like to mimic a valve, which I close @2000s, and reopen @2500s. This should have the effect, that the fluid levels remain constant in this period, and the fluid flow gets 0. After opening (@2500s) fluid levels equilibrate and the flow oscillates and slowly diminish. As visible this approach still does not give the right answer, because there is still fluid flow after 2000s. The correct answer is visible at the bottom of the attached notebook. How should I proceed, if I would like to repeat this closing/opening several times? Can this also be done with the WhenEvent command? Tx for the kind help in advance, best regards Andras Attachments: TwoTanksAndPumpM...nb

POSTED BY: agilicz

Neil Singer

Neil Singer, AC Kinetics, Inc.

Posted 5 years 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

POSTED BY: Neil Singer

Neil Singer

Neil Singer, AC Kinetics, Inc.

Posted 5 years ago

POSTED BY: Neil Singer

agilicz

Posted 5 years 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: TwoTanksAndPump.nb

POSTED BY: agilicz

agilicz

Posted 5 years 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 regards Andras

POSTED BY: agilicz

Reply to this discussion

Reply Preview

Attachments

Remove Add a file to this post

Follow this discussion

or Discard

Feedback