TITLE								        {Natural.pde }
   'Convezione naturale in fenditura sottile'

SELECT
errlim=5e-3

COORDINATES
cartesian2 	{ problema in 2D}

VARIABLES
   vx     	vy	p        Tp
DEFINITIONS

Lx=5e-3         	Ly= 5*Lx		

visc=1.8e-5		dens=1.2	{ Aria }
cond=0.025		rcp=1e3
beta=1/300

v=vector( vx, vy) 	vm=magnitude( v)
div_v= dx( vx)+ dy( vy)	
unit_x=vector(1,0)     unit_y=vector(0,1)
nx=normal( unit_x)     ny=normal( unit_y)

g=9.8		{ Gravita'}	

Ta=300		Tb=290		{ Temperature}
Tm=(Ta+Tb)/2

Fy= dens*beta*g*(Tp-Tm)		{ forza dovuta alla Gravita'}

natp=visc* [ nx* del2( vx)+ny* del2( vy)]+ny*Fy 	{ Condiz. al contorno per P}

L=globalmax(del2( p))

h=Lx/2	
Re= dens*globalmax(vm)*h/visc		{ Reynolds}
Gr=dens^2*g*beta*abs(Tb-Ta)*h^3/visc^2	{ Grashof}

qa=Eval(-cond*dx(Tp),-Lx/2,y)	{ flusso termico alla parete calda}
qb=Eval(-cond*dx(Tp),Lx/2,y)	{ flusso termico alla parete fredda}
 				
van=visc/dens/(h)*(Gr/12*(x/h-(x/h)^3)) { soluzione analitica}

EQUATIONS							
   vx:	dx( p)- visc* del2( vx)= 0
   vy:	dy( p) - Fy - visc* del2( vy)= 0		
   p:	del2( p)- dy(Fy)- 1e4* div(v)=  0		{ Sostituisce l'eq. di continuita'}
   Tp:	-cond*del2( Tp)+ rcp*[ vx*dx( Tp)+ vy*dy( Tp)]=0

CONSTRAINTS    { Integral constraints }
VOL_INTEGRAL(p)=0		{ definisce lo zero per P}

BOUNDARIES
region 1
start (-Lx/2,-Ly/2)
 value( vx)=0
 value( vy)=0
natural( Tp)= 0
natural( p)=natp
 line to (Lx/2,-Ly/2)	{ parete inferiore, adiabatica}
value( vx)=0
 value( vy)=0
value( Tp)=Ta
natural( p)=natp
 line to (Lx/2,Ly/2)		{ parete di destra, a T alta}
value( vx)=0
value( vy)=0
 natural( Tp)= 0
natural( p)=natp
line to (-Lx/2,Ly/2) 	{ parete superiore, adiabatica}
value( vx)=0
value( vy)=0
value( Tp)=Tb
natural( p)=natp
line to close		{ parete di sinistra, a T bassa}

MONITORS         { mostra i progressi }
surface(div(v)) report( Re)
surface(p) report( L)
surface(Tp) report( Gr^.5)
vector(v) zoom(-Lx/2,-Ly/2,Lx,Lx)
vector(v) zoom(-Lx/2,0,Lx,Lx)  report( Re)
vector(v) zoom(-Lx/2,Ly/2-Lx,Lx,Lx)  report( Gr^.5)
elevation(vy,van)  from (-Lx,0) to (Lx,0)
elevation(vy,van)  from (-Lx/2,-Ly/2+Lx) to (Lx/2,-Ly/2+Lx)
elevation(vy,van)  from (-Lx/2,-Ly/2+Lx/2) to (Lx/2,-Ly/2+Lx/2)
elevation(qa,qb,-cond*(Ta-Tb)/Lx) from (-Lx/2,-Ly/2) to (-Lx/2,Ly/2)

PLOTS
contour( Tp) painted report( Re)
contour(vm)  zoom(-Lx/2,-Ly/2,Lx,Ly)
vector(v) zoom(-Lx/2,-Ly/2,Lx,Lx)  report( Gr^.5)
vector(v) zoom(-Lx/2,0,Lx,Lx)  report( Re)
vector(v) zoom(-Lx/2,Ly/2-Lx,Lx,Lx)  report( Gr^.5)
vector(v) zoom(-Lx/2,-Ly/2,Lx,Ly)  report( Re)
elevation(Tp) from (-Lx/2,0) to (Lx/2,0)
elevation(vy,van) from (-Lx/2,0) to (Lx/2,0)
elevation(vy,van)  from (-Lx/2,-Ly/2+Lx) to (Lx/2,-Ly/2+Lx)
elevation(vy,van)  from (-Lx/2,-Ly/2+Lx/2) to (Lx/2,-Ly/2+Lx/2)
elevation(vx) from (-Lx/2,0) to (Lx/2,0)
surface(div(v))
surface(p) report( L)
elevation(qa,qb,-cond*(Ta-Tb)/Lx)  from (-Lx/2,-Ly/2) to (-Lx/2,Ly/2)
surface(vx)export format "#x#y#b#1" file="vx.txt"	{ esporta il file di velocita' lungo x}
surface(vy) export format "#x#y#b#1" file="vy.txt"	{ esporta il file di velocita' lungo y}
surface(Tp) export format "#x#y#b#1" file="T.txt"	{ esporta il file di temperatura}

END