Reference
Category : GAMS NOA library
Mainfile : flowobs.gms
$onText
Stationary flow of an incompressible fluid in a rectangular area in the
presence of an obstacle.
Adapted from:
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/;
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