[overture] Re: cgins with low Reynolds number flow (Stokes flow)
- From: Yongsheng Lian <yongshenglian@xxxxxxxxx>
- To: freelists <overture@xxxxxxxxxxxxx>
- Date: Mon, 11 Oct 2010 11:18:33 -0400
Hi Bill,
The full implicit version works for low speed flow. As you said, the
dtMax needs to be adjusted a few times to make it work.
Since Overture uses nu and assumes unit density in the simulation, I have
some problems to make it work for the dimensional calculation. I tried to
define mu in the cg file but it did not look it was read properly.
I notice that viscosity mu and density in InsParameters.C. And the normal
force is calculated as follows
fn(Ib1,Ib2,Ib3,0)=(
fluidDensity*ug(Ib1,Ib2,Ib3,pc)*normal(Ib1,Ib2,Ib3,0)
-(nu*((ux(Ib1,Ib2,Ib3,uc)+ux(Ib1,Ib2,Ib3,uc))*normal(Ib1,Ib2,Ib3,0)+
(uy(Ib1,Ib2,Ib3,uc)+ux(Ib1,Ib2,Ib3,vc))*normal(Ib1,Ib2,Ib3,1)) ) );
In the manual you almost always work with the dimensionless equations, I
am not sure how did you factor the density into the final solution. It seems
that the viscous forces should be multiplied by fluidDensity also.
Yongsheng
On Sat, Oct 9, 2010 at 2:05 PM, Bill Henshaw <henshaw@xxxxxxxx> wrote:
> Hi Yongsheng,
>
> Sorry for not calling you back yesterday. I had many visitors.
>
> There is another "full" implicit solver option that might
> work for your very low Reynolds case. It treats the advection
> terms implicitly too.
>
> The attached command file shows how to use it. I used the
> command
>
> cgins cic -g=cice2.order2 -inflow=parabolic -wall=noSlip -nu=1000.
> -tf=1000. -tp=10. -dtMax=.1 -iv=full
>
> which uses the grid from "ogen noplot cicArg -factor=2 -interp=e" (I
> think). Here
> the Reynolds number is about Re=1.e-3
>
> The full implicit solver can take a very large time step but if you take
> too big
> a time step then you may get funny results. Therefore you need to set
> "dtMax" to
> something not too big (e.g. dtmax=.1*(max velocity)). Once the solution
> settles
> down you can probably increase dtMax.
>
>
> Regards,
> Bill.
>
>
> #
> # cgins command file for flow past a cylinder
> #
> # Usage:
> # cgins [-noplot] cic -g=<name> -tf=<> -tp=<> -nu=<> -ts=[pc|im]
> -iv=[viscous|full] -tm=[les] ...
> # -solver=[best|yale|mg] -project=[0|1] ...
> # -psolver=[best|yale|mg] -pc=[ilu|lu] -inflow=[uniform|parabolic]
> -wall=[noSlip|slip] -cyl=[noSlip|slip] ...
> #
> # Options:
> # -project : 1=project initial conditions
> # -solver : implicit solver
> # -iv : implicit variation: viscous=viscous terms implicit. full=all terms
> implicit.
> # -psolver : pressure solver
> # -inflow : inflow boundary condition
> # -wall : boundary condition for the top and bottom walls
> # -cyl : boundary condition for the cylinder
> #
> # Examples:
> # cgins cic -g=cice2.order2
> # cgins cic -g="cice2.order2.ml2.hdf" -solver=mg
> # cgins cic -g=cice4.order2 -inflow=parabolic -wall=noSlip -ad2=1 [ok
> # cgins cic -g=cice4.order2 -inflow=parabolic -wall=slip -ad2=1 [ok
> # cgins cic -g=cice4.order2 -inflow=parabolic -wall=slip -ad2=1
> -ad22=.001 [OK
> # cgins cic -g=cice4.order2 -inflow=parabolic -wall=noSlip -ad2=0 [OK
> # cgins cic -g=cice4.order2 -inflow=parabolic -wall=slip -cyl=slip -ad2=1
> [*BAD* -> wrong number of imp solvers
> # cgins cic -g=cice4.order2 -order=2 -wall=slip -tp=.05 -nu=1.e-5 -ad2=1
> -ad22=2. -ts=im [OK. note nu is "too" small
> # -- 2nd-order + fourth-order diss:
> # cgins cic -g=cice4.order2 -order=2 -wall=slip -tp=.01 -nu=1.e-5 -ad2=0
> -ad4=1 -ad42=2. -project=0 -ts=pc [OK
> # cgins cic -g=cice4.order2 -order=2 -wall=slip -tp=.01 -nu=1.e-5 -ad2=0
> -ad4=1 -ad42=2. -project=0 -ts=im [BAD... why ??
> # -- Full implicit method: (good for low Reynolds number)
> # cgins cic -g=cice2.order2 -inflow=parabolic -wall=noSlip -nu=1000.
> -tf=1000. -tp=10. -dtMax=.1 -iv=full
> # -- Fourth-order
> # cgins cic -g=cice2.order4 -wall=slip -tp=.1 -ad2=0 -nu=.01 [ok
> # cgins cic -g=cice4.order4 -wall=slip -tp=.05 -ad2=0 -nu=.005 [ok
> # cgins cic -g=cice2.order4 -wall=slip -tp=.01 -ad2=0 -nu=.001 -ad4=1
> -ad42=2. -ts=pc [ok
> # cgins cic -g=cice4.order4 -inflow=parabolic -wall=noSlip -ad2=0
> -nu=.005
> #
> # -- LES
> # cgins cic -g=cice2.order2 -ts=pc -tm=les -lesOption=0 -go=halt
> # cgins cic -g=cice2.order2 -ts=pc -tm=les -nu=.05 -lesOption=1
> -lesPar1=.001 -go=halt
> # -- les+implicit : finish me...
> # cgins cic -g=cice2.order2 -ts=im -tm=les -lesOption=0 -nu=.05 -go=halt
> #
> # mpirun -np 2 $cginsp cic -g=cicLongChannele2 -nu=.01 -recomputeDt=2
> -debug=3 [ok]
> # mpirun -np 2 $cginsp cic -g=cice2.order2.ml2 -nu=.01 -recomputeDt=200
> -debug=3 [ok]
> # mpirun -np 2 $cginsp noplot cic -g=cice2.order2.ml2 -nu=.01
> -recomputeDt=200 -debug=3 -go=og [ok
> # srun -N1 -n2 -ppdebug memcheck_all $cginsp cic -g=cice2.order2 -nu=.01
> -recomputeDt=200 -debug=3 -go=halt
> # srun -N1 -n2 -ppdebug $cginsp -noplot cic -g=cicLongChannele2 -nu=.01
> -recomputeDt=200 -debug=3 -tf=.2 -go=go
> # srun -N1 -n4 -ppdebug $cginsp -noplot cic -g=cicLongChannele4 -nu=1.e-3
> -recomputeDt=200 -debug=3 -tf=.2 -go=go
> #
> # mpirun -np 2 $cginsp cic -g=cicLongChannele2 -nu=.01 -recomputeDt=2
> -debug=3 -pc=lu
> # mpirun -np 2 $cginsp noplot cic -g=cicLongChannele2 -nu=.01
> -recomputeDt=2 -debug=3 -solver=mg -go=og
> # mpirun -np 2 $cginsp noplot cic -g=cice2.order2.ml2 -nu=.01
> -recomputeDt=2 -debug=3 -solver=mg -go=og
> # srun -N1 -n2 -ppdebug $cginsp -noplot cic -g=cicLongChannele2 -nu=.01
> -solver=mg -debug=3 -tf=.1 -go=go
> #
> $tFinal=10.; $tPlot=.05; $nu=.1; $cfl=.9; $debug=1; $show=" "; $go="halt";
> $project=1; $dtMax=.1;
> $solver="yale"; $rtoli=1.e-5; $atoli=1.e-6; $idebug=0;
> $psolver="yale"; $rtolp=1.e-5; $atolp=1.e-6; $pdebug=0;
> $recomputeDt=200; $refactorFrequency=1000; $userDefinedOutput=0;
> $grid="cice2.order2.hdf"; $outflowOption="neumann";
> $order = -1; # use default order of accuracy based on the grid
> $ad2=1; $ad21=.5; $ad22=.5; $ad4=0; $ad41=1.; $ad42=1.;
> $pc="ilu";
> $ts="im"; # time stepping
> $iv=""; # implicit variation
> $tm = "#"; # turbulence model
> $lesOption=0; $lesPar1=.01;
> $inflow = "uniform"; $wall="noSlip"; $cyl="noSlip";
> # ----------------------------- get command line arguments
> ---------------------------------------
> GetOptions(
> "g=s"=>\$grid,"tf=f"=>\$tFinal,"tp=f"=>\$tPlot,"nu=f"=>\$nu,"recomputeDt=i"=>\$recomputeDt,\
>
>
> "debug=i"=>\$debug,"pdebug=i"=>\$pdebug,"idebug=i"=>\$idebug,"order=i"=>\$order,"show=s"=>\$show,\
>
>
> "solver=s"=>\$solver,"psolver=s"=>\$psolver,"pc=s"=>\$pc,"project=i"=>\$project,"dtMax=f"=>\$dtMax,\
>
>
> "go=s"=>\$go,"tm=s"=>\$tm,"ts=s"=>\$ts,"lesOption=i"=>\$lesOption,"lesPar1=f"=>\$lesPar1,\
>
>
> "ad2=i"=>\$ad2,"ad21=f"=>\$ad21,"ad22=f"=>\$ad22,"ad4=i"=>\$ad4,"ad41=f"=>\$ad41,"ad42=f"=>\$ad42,\
>
>
> "inflow=s"=>\$inflow,"wall=s"=>\$wall,"cyl=s"=>\$cyl,"outflowOption=s"=>\$outflowOption,\
> "iv=s"=>\$iv );
> #
> if( $solver eq "best" ){ $solver="choose best iterative solver"; }
> if( $solver eq "mg" ){ $solver="multigrid"; }
> if( $psolver eq "best" ){ $psolver="choose best iterative solver"; }
> if( $psolver eq "mg" ){ $psolver="multigrid"; }
> if( $pc eq "ilu" ){ $pc = "incomplete LU preconditioner"; }elsif( $pc eq
> "lu" ){ $pc = "lu preconditioner"; }else{ $pc="#"; }
> if( $tm eq "les" ){ $tm ="LargeEddySimulation"; }
> if( $ts eq "fe" ){ $ts="forward Euler";}
> if( $ts eq "be" ){ $ts="backward Euler"; }
> if( $ts eq "im" ){ $ts="implicit"; }
> if( $ts eq "pc" ){ $ts="adams PC"; }
> if( $ts eq "pc4" ){ $ts="adams PC order 4"; }
> if( $ts eq "mid"){ $ts="midpoint"; }
> if( $ts eq "afs"){ $ts="approximate factorization";}
> if( $iv eq "full" ){ $iv = "useNewImplicitMethod\n use full implicit system
> 1\n implicitFullLinearized"; }else{ $iv="#"; }
> if( $order eq "2" ){ $order = "second order accurate"; }elsif( $order eq
> "4" ){ $order = "fourth order accurate"; }else{ $order="#"; }
> if( $go eq "halt" ){ $go = "break"; }
> if( $go eq "og" ){ $go = "open graphics"; }
> if( $go eq "run" || $go eq "go" ){ $go = "movie mode\n finish"; }
> #
> # specify the overlapping grid to use:
> $grid
> # Specify the equations we solve:
> incompressible Navier Stokes
> $tm
> # Define LES parameters that are accessed by
> getLargeEddySimulationViscosity.bf
> define integer parameter lesOption $lesOption
> define real parameter lesPar1 $lesPar1
> exit
> # Set the order of accuracy
> $order
> #
> # Next specify the file to save the results in.
> # This file can be viewed with Overture/bin/plotStuff.
> show file options
> compressed
> open
> $show
> frequency to flush
> 3
> exit
> # display parameters
> turn off twilight zone
> #
> # choose the time stepping:
> $ts
> #
> # -- choose the implicit variation (implicitViscous is default)
> $iv
> # ---------
> #
> # but integrate the square explicitly:
> choose grids for implicit
> all=implicit
> ## square=explicit
> done
> dtMax $dtMax
> #
> final time $tFinal
> times to plot $tPlot
> plot and always wait
> #* no plotting
> #
> recompute dt every $recomputeDt
> refactor frequency $refactorFrequency
> #
> pde parameters
> nu
> $nu
> # turn on 2nd-order AD here:
> OBPDE:second-order artificial diffusion $ad2
> OBPDE:ad21,ad22 $ad21 $ad22
> # turn on 4th-order AD here:
> OBPDE:fourth-order artificial diffusion $ad4
> OBPDE:ad41,ad42 $ad41,$ad42
> OBPDE:expect inflow at outflow
> if( $outflowOption eq "neumann" ){ $cmd = "use Neumann BC at outflow";
> }else{ $cmd="use extrapolate BC at outflow"; }
> $cmd
> #
> done
> #
> # OBPSF:show variable: vorticity 1
> # OBPSF:show variable: divergence 1
> #
> cfl $cfl
> #
> # -- Here we assume that the grid was made with boundary conditions:
> # 1=inflow, 2=outflow, 5=bodies
> boundary conditions
> if( $wall eq "noSlip" ){ $cmd="all=noSlipWall"; }else{
> $cmd="all=slipWall"; }
> $cmd
> if( $inflow eq "uniform" ){ $cmd = "bcNumber1=inflowWithVelocityGiven,
> uniform(p=1.,u=1.)"; }\
> else{ $cmd = "bcNumber1=inflowWithVelocityGiven ,
> parabolic(d=.5,p=1.,u=1.)"; }
> $cmd
> #* square(0,0)=inflowWithVelocityGiven , uniform(p=1.,u=1.),
> oscillate(a0=.5,a1=.5,t0=.0,omega=.5)
> ## square(0,0)=inflowWithVelocityGiven , parabolic(d=.2,p=1.,u=1.),
> oscillate(a0=1.,t0=.0,omega=.5)
> #* square(0,0)=inflowWithPressureAndTangentialVelocityGiven, uniform(p=1.)
> #* square(1,0)=outflow, pressure(1.*p+1.*p.n=0.)
> #* square(1,0)=outflow, pressure(1.*p+0.*p.n=1.)
> #
> bcNumber2=outflow , pressure(1.*p+1.*p.n=0.)
> if( $cyl eq "noSlip" ){ $cmd="bcNumber5=noSlipWall"; }else{
> $cmd="bcNumber5=slipWall"; }
> $cmd
> done
> #
> pressure solver options
> $ogesSolver=$psolver; $ogesRtol=$rtolp; $ogesAtol=$atolp; $ogesPC=$pc;
> $ogesDebug=$pdebug;
> include ogesOptions.h
> exit
> #
> implicit time step solver options
> $ogesSolver=$solver; $ogesRtol=$rtoli; $ogesAtol=$atoli; $ogesPC=$pc;
> $ogesDebug=$idebug;
> include ogesOptions.h
> exit
> #
> initial conditions
> # read from a show file
> # cylinder.show
> # 9
> uniform flow
> p=1., u=1.
> exit
> allow user defined output $userDefinedOutput
> if( $project eq "1" ){ $cmd="project initial conditions"; }else{ $cmd="#";
> }
> $cmd
> debug $debug
> continue
> #
> $go
>
>
>
>
Other related posts:
