Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Stratified plane Couette flow test #310

Closed
ali-ramadhan opened this issue Jul 10, 2019 · 33 comments · Fixed by #381
Closed

Stratified plane Couette flow test #310

ali-ramadhan opened this issue Jul 10, 2019 · 33 comments · Fixed by #381
Assignees
@ali-ramadhan
Labels
science 🌊 Sometimes addictive, sometimes allergenic

Comments

@ali-ramadhan
Copy link
Member

ali-ramadhan commented Jul 10, 2019 •
edited
Loading

We agreed that stratified plane Couette flow could be a test case against which we can verify the implementation of AMD in PR #309, although the published results [Vreugdenhil & Taylor (2018)] use a slightly modified version of AMD.

They report a number of LES runs, and since we don't have a vertically stretched grid, it might be easiest to try and reproduce the Pr = 0.7 case which they do with (Nx, Ny, Nz) = (64, 49, 64) grid points [for Ri=0 or 0.01, for Ri = 0.04 it's (64, 65, 64)]. I might suggest focusing on the Ri = 0 case as they also run a resolved DNS with (Nx, Ny, Nz) = (256, 129, 256) against which we might be able to compare (might not be possible, depending on vertical grid stretching).

Unfortunately, I don't think we can reproduce their results without a vertically stretched grid... If I understand the paper correctly, the grid is stretched according to

y_j = h*tanh(Sf*(2*(j-1)/(Ny-1))) / tanh(Sf)

so for h=1 you get grid spacings of Δy ~ 0.125 away from the wall and Δy ~ 0.0000035 adjacent to the wall. So a faithful reproduction would need 285,000+ vertical levels lol.

I don't fully understand how they got their values for the vertical grid cell size adjacent to the wall ∆y_w^+, but it seems like if h = 100,000 then our values agree (they get a spacing of Δy ~ 0.35 adjacent to the wall). Either way, the ratio between the thickest and thinnest spacings is ~35,000.

Not sure if there's still a way we could compare results in this case without a vertically stretched grid...

Absolute worst case scenario, I've been working on an implementation of a vertically stretched grid (see PRs #283 and #306) but it's a work in progress and might not be the best use of our efforts right now.

Reference: Catherine A. Vreugdenhil and John R. Taylor, Large-eddy simulations of stratified plane Couette flow using the anisotropic minimum-dissipation model, Physics of Fluids 30, 085104 (2018).

cc @glwagner @rafferrari
@rafferrari
Copy link

rafferrari commented Jul 10, 2019 via email

@ali-ramadhan
Copy link
Member Author

ali-ramadhan commented Jul 10, 2019 •
edited
Loading

Alternatively, we can compare with the original databases of stratified plane Couette flow DNS runs by Deusebio et al. (2015) who varied Re and Ri, and Zhou et al. (2017) who varied Pr and Ri.

That way we can design our own experiment, but maybe we still can't do anything rigorous without high-resolution near the walls.

Maybe we could keep increasing the vertical resolution and see if our solutions (the statistics) converge towards the resolved DNS solution? We could go from Nz = 50 to Nz = 2000 levels or something.

References:

@ali-ramadhan ali-ramadhan added the science 🌊 Sometimes addictive, sometimes allergenic label Jul 10, 2019
@ali-ramadhan ali-ramadhan self-assigned this Jul 10, 2019
@ali-ramadhan ali-ramadhan added this to the LES closures milestone Jul 10, 2019
@johncmarshall54
Copy link

johncmarshall54 commented Jul 10, 2019 via email

@glwagner
Copy link
Member

@ali-ramadhan thanks for looking into that!

I think a convergence test with coarse horizontal resolution and varying vertical resolution is probably a good idea, perhaps producing a plot like Fig 11 in Vreugdenhil and Taylor.

In my reading, ∆y_w^+ is the vertical spacing in 'wall units', which is the physical spacing divided by nu/u_\tau, where nu is viscosity and u_\tau is the friction velocity defined in eq 16.

In Table II, values for both ∆y_w^+ and ∆y_c^+ (grid spacing at wall and at domain center) are given. It looks like many of the runs may be feasible. Here's some code:

 y(j, Ny, Sf, h) = h / tanh(Sf) * tanh(Sf*(2*(j-1)/Ny - 1)) 
Δy(j, Ny, Sf, h) = y(j+1, Ny, Sf, h) - y(j, Ny, Sf, h)

N_unstretched(Ny_VT, Sf) = ceil(Int, 2 / Δy(2, Ny_VT, Sf, 1)) 

@show N_unstretched(49, 3)
@show N_unstretched(65, 2.5)
@show N_unstretched(97, 3)
@show N_unstretched(129, 2)

It looks like some of their runs are feasible, since

julia> N_unstretched(49, 3) # runs 2-4
1141

julia> N_unstretched(65, 2.5) # run 5
768

julia> N_unstretched(97, 3) # runs 7-11
2710

julia> N_unstretched(129, 2) # run 12
805

Runs B-G (which are included to demonstrate resolutions that are apparently inadequate) may also be feasible.

One concern about the suitability of stratified Couette flow is mentioned at the end of the paper...

Stratified plane Couette flow is a challenging test case
for the LES model because it has a linearly stable laminar
state which introduces requirement 1. The results of Abkar
and Moin13 suggest that AMD LES performs even better in
other stratified wall-bounded flows. Balancing all these concerns
is key to using the AMD model in LES of wall-bounded
stratified flow. Nevertheless the AMD model is able to capture
turbulent intermittency and mean and turbulent flowproperties
in stratified plane Couette flow.

@johncmarshall54 we are planning to implement stretched vertical grid functionality after the LES is implemented.

@ali-ramadhan
Copy link
Member Author

Ali, i thought the stretching was sorted now, after the hackathon. Without stretching the code will have limited use. John

Yes well, we know how to do a very fast Poisson solve in a vertically stretched grid but still working on integrating it (PR #306). The fast solver (which relies on cyclic reduction) may need some massaging to work with Neumann boundary conditions for pressure... The plan is to merge the vertically stretched grid with the finite volume operators.

@ali-ramadhan
Copy link
Member Author

@glwagner Never mind, I am an idiot. I divided by Ny-1 by mistake. The matching factor of 0.35 might have been pure coincidence.

We need closer to O(1000) vertical levels, which should be pretty doable.

@glwagner glwagner mentioned this issue Jul 25, 2019
Merged
2 tasks
@ali-ramadhan ali-ramadhan removed this from the LES closures milestone Aug 1, 2019
@ali-ramadhan
Copy link
Member Author

Temperature and velocity profiles look qualitatively correct which is promising but the sub-grid-scale viscosity and diffusivity don't match. At least the order of magnitude on each is roughly correct.

Most notably their version of AMD allows the SGS ν and κ to go below the background values, which seems essential for ν, although not sure what that physically means. Seems that the VerstappenAnisotropicMinimumDissipation closure in the new-closures branch allows for this so I'll try rerunning with it.

The SGS profiles seem kind of noisy but maybe I just need to bump up the resolution.

The SGS/turbulent Prandtl number is still too low but this should get fixed if the SGS ν and κ profiles are right.

The wall velocity Uw=1 so the x-axis is U/Uw.
plots_stratified_couette_flow_u_profiles

The wall temperature is Tw = 0.01 so this plot should match the profile from Vreugdenhil & Taylor.
plots_stratified_couette_flow_T_profiles

correct_profiles

plots_stratified_couette_flow_nu_profiles

plots_stratified_couette_flow_kappaT_profiles

correct_SGS_profiles

@glwagner
Copy link
Member

glwagner commented Aug 18, 2019 via email

@ali-ramadhan
Copy link
Member Author

ali-ramadhan commented Aug 18, 2019 •
edited
Loading

Thanks for the explanation @glwagner.

I didn't realize that the background diffusivity had to be subtracted from the nonlinear diffusivity. In that case the ν profile isn't super far off although both nonlinear diffusivity profiles (especially the κ profile) suggest that perhaps the boundary layer is not being resolved.

I'm rerunning the simulation but with the vertical resolution tripled. Hopefully this will improve the comparison.

@glwagner
Copy link
Member

@ali-ramadhan a few more comments:

  • I don't understand why the temperature profiles you've plotted have spikes in them, especially since they are horizontal averages. What's going on? The sub grid scale diffusivities also seem to have far too much variance.

  • Can you plot multiple equilibrated horizontally averaged solutions for differing model constant C? I think we may need to determine the correct value of this constant for our numerical method in order for our results to match. It looks like the SGS diffusivity values are too small, so this might help...

  • We can't be sure what resolution an Oceananigans.jl model needs to be to match Diablo results, since Diablo uses a spectral discretization in the horizontal. That's something to keep in mind. Note also that V&T take the spectral discretization into account when defining their 'filter widths' in their AMD implementation.

  • If our results matched Diablo without resolving the viscous boundary layer near the walls, I think this would not constitute a verification of Oceananigans (instead, it would indicate a problem with Oceananigans, or the presence of some kind of 'cancellation of errors'). We must resolve the viscous boundary layers near the walls in this experiment for the results to be valid.

For reference, the AMD viscosity is calculated here:

https://github.com/climate-machine/Oceananigans.jl/blob/352c018b6cb8801500bae61ac858ddc146bf3f48/src/turbulence_closures/anisotropic_minimum_dissipation.jl#L50

the diffusivities are similar.

@ali-ramadhan
Copy link
Member Author

Thanks for the comments! I've rerun the simulation with Ri = 0.04 so that the viscous boundary layer is resolved. Horizontal resolution may still be an issue but I can roughly double it. We may also have to figure out which C to use, but might be more useful to do when we know we have a run that looks reasonable?

Just posting vertical profiles from that run. They look better but non-linear diffusivity profiles are still off. The shape of the profiles looks different from the figure I posted above. The magnitudes also seem very low (they're barely above the background values). Magnitude could be a plotting error (if it is I can't see it...).

Discrepency could be because I'm doing the Pr = 0.7 case while the profiles posted in V&T figure 7 are for the Pr = 7 case. I assume things will be different, but I don't have the intuition to tell how big of a difference to expect.

Note: nu and kappaT profiles are nu_SGS/nu and kappa_SGS/kappa.

plots_stratified_couette_flow_u_profiles (1)
plots_stratified_couette_flow_T_profiles (1)
plots_stratified_couette_flow_nu_profiles (1)
plots_stratified_couette_flow_kappaT_profiles (1)
plots_stratified_couette_flow_Pr_profiles

@ali-ramadhan
Copy link
Member Author

Just reporting some dimensionless numbers for that run (it's number 5 from table 1 of V&T).

The friction Reynolds number is a little different between the two walls but is generally between 230-250. This does not match run 5's value of 183 but seems to match the Ri=0 values. If we're not resolving something, maybe we're effectively running at a lower Ri number? Not sure if it should be fluctuating this much?

image

The friction Nusselt number is quite different between the top and bottom walls, but generally around 5 and 7. Bottom wall matches table 1's value of 7.1 pretty closely. Again, not sure if it should fluctuate this much, and if so how to time average properly.

image

image

@ali-ramadhan ali-ramadhan mentioned this issue Aug 28, 2019
Merged
@ali-ramadhan
Copy link
Member Author

Just to document what's been done:

So I ran the Ri = 0.01, Pr = 0.7 case that Cat ran at three different horizontal resolutions and unfortunately the results are quite different from hers:

  1. U and T don't fully go to zero at the top and bottom. This could be because we only have 64 vertical levels, not enough to resolve the boundary layers, but otherwise the profiles look qualitatively correct.
  2. Our νSGS/ν and κSGS/κ are too low by a factor of ~2. They decrease as you increase the resolution which makes sense. κ profile shape looks similar to Cat's. ν profile shape should be flatter.
  3. Re is roughly a factor of 2 too large. Nu is too large by ~60%. (edited)
    I might try changing the Poincare "constant" to see if it makes a difference here. Currently C = 1/12, so I'll try 1/24 and 1/6. Might also try 128 vertical levels, but that's no longer a fair comparison.

Stratified_couette_flow_UT

Stratified_couette_flow_VK

Stratified_couette_flow_ReNu

@ali-ramadhan
Copy link
Member Author

To make sure it's not just low resolution being bad, I reran everything at 256x256x128 (to keep aspect ratio the same).
SGS diffusivity profiles look okay at C = 1/6 (although shapes still don't fully match Cat's) so it's also possible that I'm calculating Re and Nu incorrectly...

Stratified_couette_flow_UT_C
Stratified_couette_flow_VK_C
Stratified_couette_flow_ReNu_C

@ali-ramadhan
Copy link
Member Author

ali-ramadhan commented Sep 6, 2019 •
edited
Loading

So going up to 256x256x128 didn't help. Re and Nu don't actually change much.
Re and Nu are still too big by a factor of ~2...

Stratified_couette_flow_UT_C_N256
Stratified_couette_flow_VK_C_N256
Stratified_couette_flow_ReNu_C_N256

@ali-ramadhan
Copy link
Member Author

Apparently I'm just an idiot.

h = model.grid.Lz
Reτ = h * √uτ²⁺ / ν
Nu = (qw⁺ * h)/(κ * Θw)

The height of the domain is actually 2h based on the domain geometry described in Cat's Table 1. So I should have used h = model.grid.Lz / 2.

There's the factor of 2.

@glwagner
Copy link
Member

glwagner commented Sep 6, 2019

🤷‍♂️🧠🏃‍♂️🙉🦉

@glwagner
Copy link
Member

glwagner commented Sep 6, 2019

How does the flow field look?

@ali-ramadhan
Copy link
Member Author

Rerunning Ri = 0.01 and Ri = 0.04 to make sure both come out matching Cat's simulations.

Will also plot some flow fields from the 256x256x128 runs to see what they look like.

@rafferrari
Copy link

rafferrari commented Sep 6, 2019 via email

@glwagner
Copy link
Member

glwagner commented Sep 6, 2019

@rafferrari

√uτ²⁺ is a velocity scale (it's written oddly I think...), h is a length scale, and ν is ν.

@ali-ramadhan should correct me, but it looks like κ * Θw / h is a diffusive heat flux scaling across the domain, whereas qw⁺ is the wall temperature flux (not sure how its calculated).

@glwagner
Copy link
Member

glwagner commented Sep 6, 2019

The Nusselt number here is the ratio between the measured turbulent temperature flux and the temperature flux we expect in the laminar solution?

@rafferrari
Copy link

rafferrari commented Sep 6, 2019 via email

@ali-ramadhan
Copy link
Member Author

ali-ramadhan commented Sep 6, 2019 •
edited
Loading

@rafferrari we're using the frictional Re and Nu as defined by equation (20) of Vreugdenhil & Taylor (2018).

frictional Re = (friction velocity) * (domain height / 2) / ν

Re is still a velocity scale U times a length scale L divided by ν just with specific choices for U and L.

advective heat transfer = (wall heat flux qw) / (wall temperature Θw)
conductive heat transfer = κ / (domain height / 2)
frictional Nu = (advective heat transfer) / (conductive heat transfer)

So Nu is still the ratio of advective to conductive heat transfer, just with specific definitions again. Well, my "heat transfers" have units of m/s I guess so maybe they're not exactly heat transfer rates...

@glwagner Yeah it's written a little oddly as we can easily calculate the friction velocity squared (which I then take square roots of), but perhaps the script will be clearer if I just use uτ.

The superscript + and - indicate the top and bottom wall values, but as @glwagner has pointed out the superscript + is reserved for wall units so perhaps I need to improve my notation here as well.

@ali-ramadhan
Copy link
Member Author

I thought one of the advantages of Julia is that one can use easy to read notation. The notation for Re and Nu is anything but intuitive. Should we discuss how to improve on that?

Yes there's definitely room for improvement here.

@glwagner
Copy link
Member

glwagner commented Sep 6, 2019 •
edited
Loading

@ali-ramadhan how is qw calculated?

Are you sure your definition of Nu is correct? I would interpret it this way:

Laminar solution heat flux = Θw κ / h

(Because the wall-to-wall temperature difference is 2Θw, and the wall spacing is 2h, so the laminar temperature gradient would be Θw/h.)

Turbulent heat flux = qw.

Then

Nusselt number = Turbulent heat flux / laminar heat flux = qw * h / (Θw * κ).

I'll review the notation in the pull request.

@ali-ramadhan
Copy link
Member Author

@glwagner Quite possibly.

qw is calculated from equation (16) of Vreugdenhil & Taylor (2018)

qw = (friction velocity) * (friction temperature)

@glwagner
Copy link
Member

glwagner commented Sep 6, 2019 •
edited
Loading

Thanks @ali-ramadhan.

For the record, this means that the wall heat flux is define in a rational, sane manner, ie:

qw = κ | ∂y Θ |_wall |

aka, the (absolute value of the) heat flux out of the wall.

(In V & T, y is the wall-normal coordinate.)

@ali-ramadhan technically we should be using the total diffusivity (LES + molecular) to calculate qw (because these new coarse simulations do not resolve the wall boundary layer). Since the SGS diffusivity is quite small next to the wall (some fraction of the molecular diffusivity), this correction may not be huge.

@rafferrari
Copy link

rafferrari commented Sep 6, 2019 via email

@glwagner
Copy link
Member

glwagner commented Sep 6, 2019 •
edited
Loading

@rafferrari I agree its confusing for us --- it comes from engineering, where the vertical velocity is "u₂"...

These simulations run to steady-state, in which case the horizontally averaged total vertical temperature flux is constant at every level. Thus it is sufficient to calculate the heat flux at the wall.

Technically speaking, the Nusselt number is the ratio between turbulent heat flux and laminar heat flux. In situations with strong turbulence of course the heat flux is dominated by the advective component... ? In LES we need to be somewhat delicate, because the total heat flux includes "resolved-advective" and "unresolved-advective" (subgrid scale) components. Luckily in equilibrated problems these distinctions are moot.

@ali-ramadhan
Copy link
Member Author

Terrible notation in my humble opinion. w stands for vertical velocity. Why is the advective heat flux qw instead of w theta?

We should improve it before merging the test and polishing it off. @glwagner also suggested agreeing on some common notation for our papers.

Following the V&T paper it's q_w (subscript w for wall) but there is no subscript w in Unicode (or if there is, it's not available in Julia) so I ended up writing qw but maybe q_w, qʷ, or even q_wall would be better notation.

@glwagner
Copy link
Member

glwagner commented Sep 6, 2019

+1 for q_wall for the purposes of this verification example.

@ali-ramadhan
Copy link
Member Author

So I reran the simulations with Ri = 0.01 and Ri = 0.04 to compare with Cat's coarse uniform grid DIABLO simulations.

Some comments:

  1. Horizontally averaged velocity and temperature profiles look qualitatively correct. Oceananigans Ri = 0.01 velocity curve could be smoother but perhaps we just need higher resolution.
  2. ν and κ profiles are different. Oceananigans' magnitude is lower and shape is a bit different. Perhaps this is expected due to differences in the numerical method. Perhaps the Poincaré "constant" needs to be set to a suitable value for our finite volume numerics. Maybe we don't expect exact agreement here.
  3. Most importantly, the frictional Reynolds (Re) and Nusselt (Nu) numbers match more closely now. Re values are within ~10% of each other. Oceananigans Nu values are ~20% smaller. Probably higher resolution is required to match more closely with Diablo and DNS?
  4. Oceananigans uses z for the vertical direction while DIABLO uses y.

Stratified_Couette_flow_UT_Ri
Stratified_Couette_flow_VK_Ri
Stratified_Couette_flow_ReNu_Ri

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees

@ali-ramadhan ali-ramadhan

Labels
science 🌊 Sometimes addictive, sometimes allergenic
Projects
None yet
Milestone
No milestone
Development

Successfully merging a pull request may close this issue.

Stratified Couette flow verification CliMA/Oceananigans.jl
4 participants
@glwagner @ali-ramadhan @johncmarshall54 @rafferrari