Flowobs : Stationary Flow of an Incompressible Fluid in a Rectangular Area

Reference

• Neculai Andrei, Nonlinear Optimization Applications Using the GAMS Technology,Springer Optimization and Its Applications, Model Flowobs (8.12) in chapter Heat Transfer and Fluid Dynamics , 2013

Category : GAMS NOA library

Mainfile : flowobs.gms

\$ontext
Stationary flow of an incompressible fluid in a rectangular area in the
presence of an obstacle.

McKinney, D.C., Savitsky, A.G., Basic optimization models for water and
energy management. June 1999 (revision 6, February 2003).
http://www.ce.utexas.edu/prof/mckynney/ce385d/papers/GAMS-Tutorial.pdf
\$offtext

set X /u1*u20/;
set Y /u1*u20/;

* Determination of zone for water movement equation
Set vyside(X,Y);
Set vxside(X,Y);
vxside(X,Y)                  = yes;
vxside(X,Y)\$(ord(X)=1)       = no;
vxside(X,Y)\$(ord(X)=card(X)) = no;
vxside(X,Y)\$(ord(Y)=1)       = no;
vxside(X,Y)\$(ord(Y)=card(Y)) = no;
vxside('u10','u10')            = no;
vxside('u10','u11')            = no;
vyside(X,Y) = vxside(X,Y);

* Parameters
scalar dx  step space in x direction   /1/;
scalar dy  step space in y direction   /1/;
scalar r   density of the fluid        /1000/;

parameter m(X,Y) kinematic viscosity ;
m(X,Y) := 0.5;

variables
obj        objective variable
D(X,Y)     error
P(X,Y)     pressure
Vx(X,Y)    x-direction velocity
Vy(X,Y)    y-direction velocity  ;

* Variable bounds and initialization
Vx.l(X,Y)     =  0.0;
Vy.l(X,Y)     =  0.0;
Vy.l(x,y)\$(vxside(x,y)) = 0.0;

D.lo(X,Y)     = 0.0;
D.up(X,Y)     = 10.0;

P.up(X,Y)     =  20000;
P.lo(X,Y)     = -20000;
P.l(X,Y)      =  0.0;

*Boundary conditions
Vx.fx('u1',Y)   = 0.5;
Vx.fx('u20',Y)  = 0.5;
Vx.fx(X,'u1')   = 0;
Vx.fx(X,'u20')  = 0;
Vy.fx('u1',Y)   = 0;
Vy.fx('u20',Y)  = 0;
Vy.fx(X,'u1')   = 0;
Vy.fx(X,'u20')  = 0;

* Obstacle description
vx.fx('u10','u10') = 0;
vx.fx('u10','u11') = 0;

Equations
For_Vx(X,Y)
For_Vy(X,Y)
Div_Vxy(X,Y)
eobj  objective function ;

For_Vx(X,Y)\$(vxside(X,Y))..
(P(X+1,Y)-P(X,Y))/(r*dx) =e=
m(x,Y)*((Vx(X+1,Y)-2*Vx(X,Y)+Vx(x-1,Y))/(dx*dy)
+        (Vx(X,Y+1)-2*Vx(X,Y)+Vx(X,Y-1))/(dx*dy));

For_Vy(X,Y)\$(vxside(X,Y))..
(P(X,Y+1)-P(X,Y))/(r*dy) =e=
m(X,Y)*((Vy(X+1,Y)-2*Vy(X,Y)+Vy(X-1,Y))/(dx*dy)
+        (Vy(X,Y+1)-2*Vy(X,Y)+Vy(X,Y-1))/(dy*dx));
*--
*For_Vx(X,Y)\$(Vxside(X,Y))..
*    Vx(X,Y)*(Vx(X+1,Y)-Vx(X-1,Y))/(2*dx) +
*    0.25*(Vy(X+1,Y-1)+Vy(X+1,Y)+Vy(X,Y-1)+Vy(X,Y)) *
*         (Vx(X,Y+1)-Vx(X,Y-1))/(2*dy) +
*         (P(X+1,Y)-P(X,Y))/(r*dx)
*    =e=
*    m(X,Y)*((Vx(X+1,Y)-2*Vx(X,Y)+Vx(X-1,Y))/(dx*dx) +
*            (Vx(X,Y+1)-2*Vx(X,Y)+Vx(X,Y-1))/(dy*dy));
*
*For_Vy(X,Y)\$(Vyside(X,Y))..
*    0.25*(Vx(X-1,Y+1)+Vx(X-1,Y)+Vx(X,Y+1)+Vx(X,Y)) *
*         (Vy(X+1,Y)-Vy(X-1,Y))/(2*dy) +
*         Vy(X,Y)*(Vy(X,Y+1)-Vy(X,Y-1))/(2*dy) +
*         (P(X,Y+1)-P(X,Y))/(r*dy)
*    =e=
*    m(X,Y)*((Vy(X+1,Y)-2*Vy(X,Y)+Vy(X-1,Y))/(dx*dx) +
*            (Vy(X,Y+1)-2*Vy(X,Y)+Vy(X,Y-1))/(dy*dy));
*--
Div_Vxy(X,Y)\$((ord(X) > 1) \$ (ord(Y) > 1))..
(Vx(X,Y)-Vx(X-1,Y))/dx + (Vy(X,Y)-Vy(X,Y-1))/dy =e=
D(X,Y);

eobj.. obj =e= SUM((X,Y),(D(X,Y)*D(X,Y)));

Model flowobs /all/;

\$iftheni x%mode%==xbook
\$onecho >bench.opt
solvers conopt knitro minos snopt
\$offecho

flowobs.optfile=1;
option nlp=bench;
\$endif

Solve flowobs using nlp minimizing obj;

\$iftheni x%mode%==xbook
* Put the solution
file res1 /pressure.dat/
put res1;
loop(X, loop(Y, put P.l(x,y):10:4; ); put /;); put /;

file res2 / vx.dat/
put res2;
loop(X, loop(Y, put Vx.l(x,y):9:5; ); put /;); put /;

file res3 / vy.dat/
put res3;
loop(X, loop(Y, put Vy.l(x,y):5:1; ); put /;); put /;

file res4 /err.dat/
put res4
loop(X, loop(Y, put D.l(x,y):5:1; ); put /;); put /;
\$endif

* end flowobs
*************************************************** November 19, 2008