TITLE 'Moto attorno ad un cilindro'     { cilindro.pde }
COORDINATES cartesian2  { problema in 2D}
VARIABLES
 vx      vy       p     Tp
SELECT
errlim=0.5e-2

DEFINITIONS
  {geometria}
Lx=10.0	 Ly=1.0		a=0.1	{Raggio}
  {propr fisiche}
visc=1e-3	  dens=1e3	cond=0.7	Cp=4200 {acqua }
rcp=dens*Cp	 Pr=cp*visc/cond

T_in=400	T_cil=300

cost=visc/dens/2/a
 k=staged	(1,		1,		1,		1)
 v_in =staged	(1e-2*cost, 	1e-1*cost, 	1e-0*cost, 	1e1*cost)

v= vector( vx, vy)
vm= magnitude( v)	vmax=globalmax( vm)
 unit_x= vector(1,0)  unit_y= vector(0,1)
nx= normal( unit_x)	ny=normal( unit_y)
tx=tangential( unit_x)	ty=tangential( unit_y)

Re= dens* vmax*2*a/visc

vxdvx= vx*dx(vx)+ vy*dy(vx)
vxdvy= vx*dx(vy)+ vy*dy(vy)
natp= visc*( nx*del2(vx)+ ny*del2(vy))-k*dens*(nx*vxdvx+  ny*vxdvy)

dist=(x^2+y^2)^.5
costeta=x/dist	sinteta=y/dist
vr=vx*costeta+vy*sinteta
vteta=-vx*sinteta+vy*costeta
tau=-visc*(dx(vteta)*costeta+dy(vteta)*sinteta+vteta/dist)  { Sforzo tangenziale sul cilindro }

fp=2*line_integral(-p*costeta, 'cylinder')  { forza dovuta alla pressione}
fa=2*line_integral(tau*sinteta, 'cylinder')  { forza dovuta all'attrito}
f=(abs(fa)+abs(fp))/(1/2*dens*vmax^2*2*a)  {fattore d'attrito}
fteor=9.689*re^(-.78)*(1+.147*re^.82)  {correlazione per moto attorno a cilindri}

flux=-cond*(dx(Tp)*costeta+dy(Tp)*sinteta)  { flusso termico uscente }

calore=2*line_integral(flux, 'cylinder') { portata termica uscente }
h=calore/(T_cil-T_in) { coeff di scambio termico}
Nu=h*2*a/cond
Nuteor=0.898*(Re*Pr)^.33 {correlazione per moto attorno a cilindri}

EQUATIONS
vx:  k*dens* vxdvx+ dx( p)- visc* del2( vx)= 0
vy:  k*dens* vxdvy+ dy( p)- visc* del2( vy)= 0
p:  del2( p)+ k*dens*(dx(vxdvx)+ dy(vxdvy))- 1e5*visc/Lx^2*div(v)=0 { Sostituisce l'eq. di continuita'}
Tp: -cond*del2( Tp)+ k*rcp*( vx*dx( Tp)+ vy*dy( Tp))=0

CONSTRAINTS
!VOL_INTEGRAL(p)=0		{ definisce lo zero per P}
BOUNDARIES
 region 'domain'  start (-Lx/3,Ly)
value( vx)=v_in		natural( vy)=0		natural(p)= 0		value(Tp)=T_in { In }
  line to (-Lx/3,0)
natural( vx)=0   	value( vy)=0   		natural(p)= natp    	natural(Tp)=0
  line to (Lx,0)
natural( vx)=0   	 natural( vy)=0  	VALUE(p)=0    	natural(Tp)=0 { Out }
line to (Lx,Ly)  		
value( vx)=0  		value( vy)=0   		natural(p)= natp	natural(Tp)=0
     line to close

start 'cylinder' (-a,0)	
 value( vx)=0  		value( vy)=0   		natural(p)=natp 	value(Tp)=T_cil
arc( center= 0,0) to (a,0)

MONITORS
PLOTS
 grid(x,y) zoom(-Ly/4,0, Ly, Ly)
contour( vm) painted report(re)
contour( Tp) painted
elevation( p) on 'cylinder'
elevation(-p*costeta,tau*sinteta) on 'cylinder' as 'sforzi sul cilindro' report(re)
vector( v) norm zoom(-Ly/4,0, Ly, Ly) report(re)
contour( P) painted zoom(-Ly/4,0, Ly, Ly) report(re)
contour( Tp) painted zoom(-Ly/4,0, Ly, Ly) report(re)
elevation( flux) on 'cylinder' as 'flusso termico'  report(re)
elevation( vx,v_in) from(Lx,0) to (Lx,Ly) as 'velocità all uscita'
history(fp,fa) versus (re) as 'forze sul cilindro'
history(f,24/Re,fteor) versus (re) loglog as 'fattore di attrito'
history(Nu,Nuteor) versus (re) loglog as 'Nusselt'
END