How to Implement Constant Heat Flux with Neumann Boundary Conditions?

[Pato4184]

Hello, I hope you are all doing well.

I am currently conducting simulations of flow in a periodic channel using Incompact3D.

I am attempting to replicate the results presented in the paper by Flageul et al. (https://www.sciencedirect.com/science/article/pii/S0142727X15000910).

Specifically, I want to impose a constant heat flux on both walls using Neumann boundary conditions. Below, I have included my input file for reference.

Additionally, I should mention that in the file transeq.f90, immediately after the line !X PENCILS, where the following statement appears:

dphi1(:,:,:,1) = ta1(:,:,:) + tc1(:,:,:),

I have added a source term as follows:

dphi1(:,:,:,1) = ta1(:,:,:) + tc1(:,:,:) + ux1(:,:,:).

I would like to know if the way I have configured the Neumann boundary conditions is correct, as the results I am obtaining are incorrect, and I honestly do not know what might be wrong.

Thank you in advance for your time and help.

#########################################

! -- mode: f90 --

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

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

! Domain decomposition
p_row=0 ! Row partition
p_col=0 ! Column partition

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

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

! Flow parameters
re = 3420 ! nu=1/re (Kinematic Viscosity)
cpg = F ! if cpg=T, then re is friction Reynolds number

! Time stepping
dt = 0.005 ! Time step
ifirst = 1 ! First iteration
ilast = 50000 ! Last iteration

! Enable modelling tools
ilesmod=0 ! if 0 then DNS
numscalar = 1 ! How many scalars? (Set to zero to disable scalars)
iibm=0 ! Flag for immersed boundary method

! Boundary and initial conditions
iin = 1 ! Inflow conditions (1: classic, 2: turbinit, 4: SEM)
u1 = 1. ! u1 (max velocity) (for inflow condition)
u2 = 1. ! u2 (min velocity) (for inflow condition)
init_noise = 0.125 ! Turbulence intensity (1=100%) !! Initial condition
inflow_noise = 0.0 ! Turbulence intensity (1=100%) !! Inflow condition
idir_stream = 1 ! Index of the streamwise direction (1=X, 3=Z, 2 not available)

nclx1 = 0
nclxn = 0
ncly1 = 2
nclyn = 2
nclz1 = 0
nclzn = 0

/End

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

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

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

! Dissipation control
nu0nu = 4.0 ! Ratio between hyperviscosity/viscosity at nu
cnu = 0.44 ! Ratio between hypervisvosity at k_m=2/3pi and k_c= pi

/End

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

! Basic I/O
irestart = 0 ! Read initial flow field ?
icheckpoint = 25000 ! Frequency for writing backup file
ioutput = 1000 ! Frequency for visualization
ilist = 25 ! Frequency for writing to screen
nvisu = 1 ! Size for visualisation collection
nprobes = 4 ! Number of probes

validation_restart = F

/End

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

wrotation = 0.12 ! rotation speed to trigger turbulence
spinup_time = 5000 ! number of time steps with a rotation to trigger turbulence
nstat = 1 ! Size arrays for statistic collection
initstat = 300000 ! Time steps after which statistics are collected

/End

!=================
&ProbesParam
!=================

flag_all_digits = T ! By default, only 6 digits, 16 when True
flag_extra_probes = T ! By default, gradients are not monitored
xyzprobes(1,1) = 0.1 ! X position of probe 1
xyzprobes(2,1) = 0.1 ! Y position of probe 1
xyzprobes(3,1) = 0.1 ! Z position of probe 1
xyzprobes(1,2) = 0.1 ! X position of probe 2
xyzprobes(2,2) = 0.2 ! Y position of probe 2
xyzprobes(3,2) = 0.1 ! Z position of probe 2
xyzprobes(1,3) = 0.1 ! X position of probe 3
xyzprobes(2,3) = 0.3 ! Y position of probe 3
xyzprobes(3,3) = 0.1 ! Z position of probe 3
xyzprobes(1,4) = 0.1 ! X position of probe 4
xyzprobes(2,4) = 0.4 ! Y position of probe 4
xyzprobes(3,4) = 0.1 ! Z position of probe 4

/End

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

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

! Boundary Conditions --> alpha_sc T +- beta_sc dT/dy = g_sc

! Bottom
alpha_sc(:,1) = 0
beta_sc(:,1) = 1
g_sc(:,1) = 2428.2 !Re * Pr

! Top
alpha_sc(:,2) = 0
beta_sc(:,2) = 1
g_sc(:,2) = -2428.2 !Re * Pr

sc(1) = 0.71 ! Schmidt number (or Prandt Number)

nclxS1 = 0
nclxSn = 0
nclyS1 = 2
nclySn = 2
nclzS1 = 0
nclzSn = 0

/End

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

iles = 0 ! LES Model (1: Phys Smag, 2: Phys WALE, 3: Phys dyn. Smag, 4: iSVV, 5: dyn SEV)
smagcst = 0.14 ! Smagorinsky constant
SmagWallDamp = 0 ! 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

&CASE
/End

###########################################################

[mathrack]

Hello,

Have you tried to run a laminar channel flow in 2D with an imposed heat flux to check your input file ? In your input file, the value selected for the imposed flux is probably incorrect. See https://github.com/mathrack/incompact3d/blob/master/neuma/schemes.f90#L1149. Check the detailed implementation of the boundary condition for implicit time schemes in the current code. I think the code is imposing the outward wall-normal derivative.

If you are not sure to understand what the code is doing, you can try to deactivate the convective term and simply look at the diffusion of the scalar.

Please note that the source code and the input file for the cited paper is available at https://github.com/mathrack/incompact3d and https://repo.ijs.si/CFLAG/incompact3d.

As a side remark, when imposing the flux at the top and bottom walls, one might observe a slow temporal drift of the temperature. To overcome this slow temporal drift, one can enforce a constant value for the domain-averaged temperature. If the equilibrium between the imposed flux and the source term is good, the correction will be small.

Cheers,
Mathrack