# resolution Poiseuille Turbulent avec equation de la Chaleur # # PYL 12/10 # # http://www.lmm.jussieu.fr/~lagree/COURS/ENSTA/HTMLPC/PC5/ # # Define Rstar 100 Define ulog(y) log(0.82*y + pow(1 + 4*pow(0.41,2)*pow(y,2),0.5))*pow(0.41,-1) + (pow(0.41,-2)*pow(y,-1)*(1 - pow(1 + 4*pow(0.41,2)*pow(y,2),0.5)))/2. # periodic box 1 1 GfsSimulation GfsBox GfsGEdge { y=0.5 } { Time {end = 60 } # attention la boite est L0xL0 PhysicalParams { L = 1. } Refine { if (fabs(y) <= 0.05 ) return 7; if (fabs(y) > 0.05 && fabs(y) <= 0.1) return 6; if (fabs(y) > 0.10 && fabs(y) <= 0.2) return 5; if (fabs(y) > 0.2 && fabs(y) < 0.4) return 4; return 3; } Init { istep = 1 } { DY = dy("U"); moinsupvp= ( 0.41*0.41*(y*y) *sqrt(2.)*D2)*sqrt(2.)*D2 } # use backward Euler to avoid Crank-Nicholson oscillations in time SourceViscosity {} { // # longueur de Prandtl double Mu = 1; Mu = 1./Rstar + 0.41*0.41*sqrt(2.)*(y*y)*D2; return Mu; } { beta = 1 tolerance = 1e-4 } VariableTracer T # Add diffusion to tracer T SourceDiffusion T { // # longueur de Prandtl et Prandtl turbulent double Muk = 1; Muk = 1./Rstar/0.7 + 0.41*0.41*sqrt(2.)*(y*y)*D2; return Muk; } Init {} {U = ulog(y*Rstar) T = 1 } Source U 1 # we need this so that acceleration can be balanced by viscous stress AdvectionParams { gc = 1 } AdvectionParams { scheme = none } # AdaptVorticity { istep = 1 } { maxlevel = 7 cmax = 1e-2 } EventStop { istep = 1 } U 1.5e-6 DU ProjectionParams { tolerance = 1e-6 } OutputSimulation { istep = 5 } stdout OutputTime { step= 0.25 } stderr OutputLocation { step = 0.25 } vals.data cuty.dat OutputSimulation { step = 0.25 } SIM/sim-%g.txt { format = text } EventScript { step = 0.25 } { cp SIM/sim-$GfsTime.txt sim.data } } GfsBox { bottom = Boundary { BcDirichlet T 0 BcDirichlet U 0 BcDirichlet V 0 } top = Boundary { BcNeumann T 1 BcNeumann U 0 BcDirichlet V 0 } } 1 1 right # p"< awk '{if(($1<0.01)&&($1>0)){print $0}}' sim.data" u 2:6 w l # set logscale x;p[0.005:]"< awk '{if($1>30){print $0}}' vals.data" u 3:7 w p