Implementing Stress-Free Boundary Conditions in the 2D Rayleigh-benard Convection Script

[Carla]

Hello,

I’m an undergraduate student, and I am currently starting to work with Dedalus. At the moment I am trying to modify the 2D Rayleigh-Benard convection example script to implement stress-free (free-surface) boundary conditions. In the end, the aim is to verify the value of the critical Rayleigh number (657.5, taken from Turcotte and Schubert’s “Geodynamics”) above which the fluid layer becomes unstable against convection. The problem is now that the outcome is not as expected: When substituting in Rayleigh number values above the critical value and plotting the maximum Reynolds number versus the time, the Reynolds number still decreases exponentially (instead of increasing).

I have tried approaching the problem in two ways: First by defining the components in x and z directions, and then using them and their derivatives to write the boundary conditions as equations in the “Problem” section. The other approach was trying to adapt the code implementing the stress-free BC in the internally heated ball example script, such that it works for the cartesian system in the Rayleigh-Benard scenario.

I have attached the modified version of the script.

If anyone can spot a mistake or help in any other way that would be greatly appreciated!

Thanks,

Carla

[Evan Henry Anders]

Hi Carla,

I think all of the different ways you've tried to implement stress-free boundaries are equivalent and correct (nicely done!) I think the problem here is the value you're taking for the critical Rayleigh number. Per e.g., equation 2.3 of this paper with your aspect ratio Gamma=4, I think the critical rayleigh number for your setup should be ~761. Try using a rayleigh number greater than that and see what happens?

Also, it's probably worth running your code for longer than 800 iterations (try 2000? 10000?); sometimes when you're running a simulation very close to critical, it can take a while for the weak growth mode to overtake the decay of the initial noise conditions.

Best,

Evan

--
You received this message because you are subscribed to the Google Groups "Dedalus Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dedalus-user...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/CAOb829kpV7_56O3rDmbLKY2B1LKO4VyWbNRz-aTaEVJnjrSnRw%40mail.gmail.com.

[Carla]

Hi Evan,

That makes sense, I'll try again and see what happens. Thanks a lot for checking the code and your explanations!

Cheers,

Carla

[Louis-Alexandre Couston]

Hi Carla, Evan,

I want to implement also stress free BCs for a Poiseuille flow and I am curious about the effect of tau correction appearing in the stress free BC (if any).

Approach 1 of Carla doesn't have tau terms in the BC but Approach 2 does (see below). Does the tau term in the latter becomes 0?

Should we have tau terms in the stress free BC? (I'm thinking we shouldn't have tau terms--these are for the bulk right--but could be wrong!)

Thanks,

Louis

--------------------------------------------------------------------------------------------------

grad_u = d3.grad(u) + ez * lift(tau_u1)  # First-order reduction

# For stress-free boundary conditions:

# Approach 1: Defining the components of u, and the z-derivative of u_x:

u_x = u @ ex

u_z = u @ ez

dzu_x = d3.Differentiate(u_x, coords['z'])

# Approach 2: Defining the strain rate matrix, extracting the horizontal shear stress at the boundaries:

strain_rate = grad_u + d3.trans(grad_u)

shear_stress_0 = (strain_rate(z=0) @ ez) @ ex

shear_stress_Lz = (strain_rate(z=Lz) @ ez) @ ex

[Daniel Lecoanet]

Hi Louis,

If your simulation is well-resolved, the tau terms should be small. In practice, I haven’t seen much difference between including the tau terms in the BCs vs not. It is not necessary to include tau terms in the BCs, but it might be more convenient.

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/8d789b1d-bca6-444b-92b0-a0ed47232ec0n%40googlegroups.com.

[Louis-Alexandre Couston]

Hi Daniel,

Ok. Thanks for your quick answer!

Best,

Louis

You received this message because you are subscribed to a topic in the Google Groups "Dedalus Users" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/dedalus-users/7jBhNsRj-XU/unsubscribe.
To unsubscribe from this group and all its topics, send an email to dedalus-user...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/223597EC-1DC5-419D-A57C-7DA4AEE4AA43%40northwestern.edu.

[Carla]

Hi again,

thanks again to Evan, your reply solved my problem and helped me understand. I am now working with a different script, the internally heated full sphere ('Internally heated convection (ball IVP)'), and I'm stuck again: First, I tried to verify the onset of convection again under different boundary conditions: A fixed temperature on the sphere, combined with either stress-free or no-slip boundaries. I've tried implementing the fixed temperature and no-slip boundary like this:

# No-slip condition:

problem.add_equation("u(r=1) = 0")

# Fixed temperature (instead of constant heat flux):

problem.add_equation("T(r=1) = 0")

When plotting the Reynolds number (calculated the same way as in the Rayleigh-Benard script) however, it again does not behave as expected, only this time it increases exponentially at Rayleigh values below the minimum critical values for both the no-slip and stress-free case. The critical values used stem from Chandrasekhar's book. Again, if anyone can spot a mistake and let me know that would be appreciated greatly!

Lastly, I also wanted to determine the Nusselt number of the system. For this, I have tried implementing the Nusselt number like this:

flow.add_property(1/T(r=0), name = 'Nusselt')

(Note: I've used the equation T_0(r=0)/T(r=0) for the Nusselt number in an internally heated sphere, where T_0 is the conduction temperature profile 1-r^2, and T is the overall temperature profile denoted T in the example script).

and then later, in the while-loop

Nus = flow.max('Nusselt')

I only included the 'max' expression because otherwise it would give me an error ('Name Nusselt not defined'), however I figured that since the temperature field is evaluated at a single point, the flow property contains only one value which is consequently also the maximum value.

When computing some numerical values, they seemed to make sense, still I wanted to ask if someone could let me know if this is an appropriate way to get the Nusselt number.

Apart from these changes I left the script pretty much as it was except for adding some things used for plotting. I have attached the script at the bottom, I couldn't figure out how to attach it as a document here.

Thanks a lot in advance, and sorry for pestering again with my beginner questions,

Carla

Script:

import matplotlib

matplotlib.use('Agg')
import matplotlib.pyplot as plt
from textwrap import wrap
import sys
import numpy as np
import dedalus.public as d3
import logging

logger = logging.getLogger(__name__)

fig = plt.figure()
ax = fig.add_subplot(111)

# Allow restarting via command line
restart = (len(sys.argv) > 1 and sys.argv[1] == '--restart')

# Parameters
Nphi, Ntheta, Nr = 128, 64, 96
Rayleigh = 2e03
Prandtl = 1
dealias = 3 / 2
stop_sim_time = 20 + 20 * restart
timestepper = d3.SBDF2
max_timestep = 0.05
dtype = np.float64
mesh = None
t = 0

# Bases
coords = d3.SphericalCoordinates('phi', 'theta', 'r')
dist = d3.Distributor(coords, dtype=dtype, mesh=mesh)
basis = d3.BallBasis(coords, shape=(Nphi, Ntheta, Nr), radius=1, dealias=dealias, dtype=dtype)
S2_basis = basis.S2_basis()

# Fields
u = dist.VectorField(coords, name='u', bases=basis)
p = dist.Field(name='p', bases=basis)
T = dist.Field(name='T', bases=basis)
tau_p = dist.Field(name='tau_p')
tau_u = dist.VectorField(coords, name='tau u', bases=S2_basis)
tau_T = dist.Field(name='tau T', bases=S2_basis)

# Substitutions
phi, theta, r = dist.local_grids(basis)
r_vec = dist.VectorField(coords, bases=basis.radial_basis)
r_vec['g'][2] = r
T_source = 6
kappa = (Rayleigh * Prandtl) ** (-1 / 2)
nu = (Rayleigh / Prandtl) ** (-1 / 2)
lift = lambda A: d3.Lift(A, basis, -1)
strain_rate = d3.grad(u) + d3.trans(d3.grad(u))
shear_stress = d3.angular(d3.radial(strain_rate(r=1), index=1))

# Problem
problem = d3.IVP([p, u, T, tau_p, tau_u, tau_T], namespace=locals())
problem.add_equation("div(u) + tau_p = 0")
problem.add_equation("dt(u) - nu*lap(u) + grad(p) - r_vec*T + lift(tau_u) = - cross(curl(u),u)")
problem.add_equation("dt(T) - kappa*lap(T) + lift(tau_T) = - u@grad(T) + kappa*T_source")
problem.add_equation("shear_stress = 0") # Stress free
# problem.add_equation("u(r = 1) = 0") # No-slip boundary conditions
problem.add_equation("radial(u(r=1)) = 0") # No penetration
# problem.add_equation("radial(grad(T)(r=1)) = -2") # constant heat flux
problem.add_equation("T(r=1) = 0") # constant temperature at the boundary
problem.add_equation("integ(p) = 0") # Pressure gauge

# Solver
solver = problem.build_solver(timestepper)
solver.stop_sim_time = stop_sim_time

# Initial conditions
if not restart:
T.fill_random('g', seed=42, distribution='normal', scale=0.01) # Random noise
T.low_pass_filter(scales=0.5)
T['g'] += 1 - r ** 2 # Add equilibrium state
file_handler_mode = 'overwrite'
initial_timestep = max_timestep
else:
write, initial_timestep = solver.load_state('checkpoints/checkpoints_s20.h5')
initial_timestep = 2e-2
file_handler_mode = 'append'

# Analysis
slices = solver.evaluator.add_file_handler('slices', sim_dt=0.1, max_writes=10, mode=file_handler_mode)
slices.add_task(T(phi=0), scales=dealias, name='T(phi=0)')
slices.add_task(T(phi=np.pi), scales=dealias, name='T(phi=pi)')
slices.add_task(T(phi=3 / 2 * np.pi), scales=dealias, name='T(phi=3/2*pi)')
slices.add_task(T(r=1), scales=dealias, name='T(r=1)')
checkpoints = solver.evaluator.add_file_handler('checkpoints', sim_dt=1, max_writes=1, mode=file_handler_mode)
checkpoints.add_tasks(solver.state)

# CFL
CFL = d3.CFL(solver, initial_timestep, cadence=10, safety=0.5, threshold=0.1, max_dt=max_timestep)
CFL.add_velocity(u)

# Flow properties
flow = d3.GlobalFlowProperty(solver, cadence=10)
flow.add_property(u @ u, name='u2')
flow.add_property(1 / T(r=0), name='Nusselt')
flow.add_property(np.sqrt(u @ u) / nu, name='Re')

# For plotting:
u_list = list()
Re_list = list()
times_list = list()

# Main loop
try:
logger.info('Starting main loop')
while solver.proceed:
timestep = CFL.compute_timestep()
solver.step(timestep)
max_Re = flow.max('Re')
Re_list.append(max_Re)
t += timestep
times_list.append(t)
#Nus = flow.max('Nusselt')
if (solver.iteration - 1) % 10 == 0:
logger.info("Iteration=%i, Time=%e, dt=%e, max(u)=%e, Nusselt=%e" % (
solver.iteration, solver.sim_time, timestep, max_u, Nusselt))
except:
logger.error('Exception raised, triggering end of main loop.')
raise
finally:
solver.log_stats()
ax.plot(times_list, Re_list))
ax.set_title('\n'.join(wrap('Ra = 2e03, Boundary conditions: Stress-Free, Constant Temperature T(r=1) = 0')))
ax.set_ylabel('Re_max')
ax.set_xlabel('Time')
fig.tight_layout()
fig.savefig('stress-free_const-temp_ra-2000')

[Carla]

Hello again,

We have solved problem one: There was a difference (a factor of 1/2 in the temperature scale) in the nondimensionalisation of the equations as compared with that in Chandrasekhar. By multiplying both the conduction temperature profile (1-r**2) as well as the source temperature T_source with a factor of 1/2, this problem could be solved, and the system behaves as expected for the Rayleigh values stated in Chandrasekhar!

However, if someone could confirm to me if the implementation of the Nusselt number makes sense I'd still be grateful.

Cheers,

Carla

[Daniel Lecoanet]

Hi Carla,

The Nusselt number calculation looks good to me.

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/611bfa49-0720-42d9-8723-77b92481f74cn%40googlegroups.com.

[Carla]

Thanks!

Carla