NaN values after certain steps of simulation.

[gauravity30]

Hello,
I am running a simulation over a flat-plate at Re=21000 in cylinder configuration. I am running the simulation at two different mesh resolutions (321x320x20) and (641x640x20). Both the simulations are giving NaN values after 5000 timesteps. I have attached the input.i3d file for your reference. I had modified the BC-Cylinder.f90 as follows to produce the flat-plate geometry, Can you help what is wrong with this ? I tried lowering the time-step thinking maybe the solution is diverging but it didn't help.

do k=nzi,nzf
       do j=nyi,nyf
          ym=yp(j)
          do i=nxi,nxf
             xm=real(i-1,mytype)*dx
             yy = ym-cey
             xx = xm-cex
             if ((yy.le.0.2).and.(yy.ge.-0.2)) then
                    if((xx.le.(1-(0.75*(yy+0.2)))).and.(xx.ge.-(1-(0.75*(yy+0.2))))) then
                        epsi(i,j,k)=remp
                    endif    
             endif
          enddo
       enddo
    enddo

! -*- mode: f90 -*-

!===================
&BasicParam
!===================

! Flow type (1=Lock-exchange, 2=TGV, 3=Channel, 4=Periodic hill, 5=Cylinder, 6=dbg-schemes)
itype = 5

! Domain decomposition
p_row=20 !20             ! Row partition
p_col=10 !10              ! Column partition

! Mesh
nx=641               ! X-direction nodes
ny=640                 ! Y-direction nodes
nz=20                 ! Z-direction nodes
istret = 0            ! y mesh refinement (0:no, 1:center, 2:both sides, 3:bottom)
beta = 0.259065151    ! Refinement parameter (beta)

! Domain
xlx = 55.      ! Lx (Size of the box in x-direction)
yly = 15.25            ! Ly (Size of the box in y-direction)
zlz = 15.25            ! Lz (Size of the box in z-direction)

! Boundary conditions
nclx1 = 2
nclxn = 2
ncly1 = 0
nclyn = 0
nclz1 = 0
nclzn = 0


! Flow parameters
iin = 1               ! Inflow conditions (1: classic, 2: turbinit)
re  = 1333.33 ! nu=1/re (Kinematic Viscosity)
u1  = 7.875              ! u1 (max velocity) (for inflow condition)
u2  = 7.875              ! u2 (min velocity) (for inflow condition)
init_noise  = 0.01   ! Turbulence intensity (1=100%) !! Initial condition
inflow_noise = 0.    ! Turbulence intensity (1=100%) !! Inflow condition

! Time stepping
dt = 0.0001              ! Time step
ifirst = 1            ! First iteration
ilast = 600000     ! Last iteration

! Enable modelling tools
ilesmod=0             ! if 0 then DNS
iscalar=0             ! If iscalar=0 (no scalar), if iscalar=1 (scalar)
iibm=2                ! Flag for immersed boundary method

/End

!====================
&NumOptions
!====================

! Spatial derivatives
ifirstder = 4         ! (1->2nd central, 2->4th central, 3->4th compact, 4-> 6th compact)
isecondder = 5        ! (1->2nd central, 2->4th central, 3->4th compact, 4-> 6th compact, 5->hyperviscous 6th)
ipinter = 3           ! interpolation scheme (1: classic, 2: optimized, 3: optimized agressive)

! Time scheme
itimescheme = 3       ! Time integration scheme (1->Euler,2->AB2, 3->AB3, 4->AB4,5->RK3,6->RK4)

/End

!=================
&InOutParam
!=================

! Basic I/O
irestart = 0           ! Read initial flow field ?
icheckpoint = 10000    ! Frequency for writing backup file
ioutput = 5000        ! Frequency for visualization
ilist    = 1000          ! Frequency for the output to screen
nvisu = 1              ! Size for visualisation collection

/End

!=================
&Statistics
!=================

nstat = 1              ! Size arrays for statistic collection
initstat = 30000       ! Time steps after which statistics are collected 

/End

!########################
! OPTIONAL PARAMETERS
!#######################

!================
&ScalarParam
!================

numscalar = 0          ! How many scalars?
Sc = 0.2               ! Schmidt number

/End

!================
&LESModel
!================

jles = 0              ! LES Model (1: Phys Smag, 2: Phys WALE, 3: Phys dyn. Smag, 4: iSVV, 5: dyn SEV)
smagcst = 0.14        ! Smagorinsky constant
SmagWallDamp = 1      ! 1: Mason and Thomson Damping function, otherwise OFF
walecst = 0.5         ! WALES Model Coefficient
iconserv = 0          ! Formulation SGS divergence (0: non conservative, 1: conservative)

/End


!================
&ibmstuff
!================

cex=27.5        !when simulating a cylinder x coordinate of the center
cey=7.25        !when simulating a cylinder y coordinate of the center
ra=1        !when simulating a cylinder, radius
nraf=10       !level of refinement for iibm==2 to find the surface of the immersed object
nobjmax=2     !number of immersed objects (DO NOT USE ZERO OBJECTS)
iforces=1
nvol=1        !Number of volumes for computing force balance

/End

&ForceCVs

xld(1) = 26.3 !X left for volume control
xrd(1) = 28.7 !X right for volume control
yld(1) = 6.5 !Y bottom for volume control
yud(1) = 8. !Y top for volume control
/End

&CASE
/End

[slaizet]

Hi,
I would suggest to use a stretched mesh in the direction normal to the walls so you can capture properly the boundary layers. I suspect that you have some stabilities issues. You have also not used a proper normalisation for the simulation so this might create some issues.
Best wishes,
Sylvain

[gauravity30]

Hi,
Thank you for the suggestion of using stretched mesh. Can you please explain a bit about the normalization ? Sorry, I don't know what it means.

Thank You,
Gaurav

[slaizet]

We typically normalise all quantities with the density, a reference length and a reference velocity. As an example for the cylinder, we use a constant density equal to 1 (constant for incompressible flow anyway), the diameter D=1 and the free-stream velocity equal to 1 (hence why Re=1/nu).
Best,
Sylvain

[gauravity30]

Okay, I did the normalization for the Reynold's number by dividing the original Reynolds number by velocity and reference length. Similarly, I divided the domain lengths by the the reference length. I will check the normalization again and try the y-stretching.