<jupyter_start><jupyter_text>RANS from DNS Our goal here, in general, is to find new valid equations that describe fluid flow. We will try to find new closures for the steady RANS equations based on direct numericalsimulation (DNS) of a boundary layer. Steady RANS equations with generic Reynolds stress effects$$(\vec{U} \cdot \nabla) \vec{U}+ \frac{1}{\rho} \nabla P - \nu \nabla^2 \vec{U}= \mathbf{R},$$where in this case $\mathbf{R}$ is simply the effects of the Reynolds stresses (i.e., the opposite of the gradient), not the Reynolds stresses themselves.Some ideas for what $\mathbf{R}$ could be:$$\mathbf{R} = A\nabla P^2 + B\nabla K + C \nabla \times \vec{U}     + D\nabla(\nabla \times \vec{U})^2    + E \vec{U}$$ Algorithm1. Pick terms (in addition to non-Reynolds stress Navier--Stokes terms).2. Create a random list of points in space that is at least as large as the number   of terms.3. At each point, acquire all data for all terms for all times.4. Average data at each point for all times.5. Solve for coefficients using a linear model.  Terms$$U \frac{\partial U}{\partial x} + V \frac{\partial U}{\partial y} + W \frac{\partial U}{\partial z}+ \frac{1}{\rho}\frac{\partial P}{\partial x} - \nu \left(     \frac{\partial^2 U}{\partial x^2}    + \frac{\partial^2 U}{\partial y^2}     + \frac{\partial^2 U}{\partial z^2} \right)$$$$= A \left( \frac{\partial U}{\partial x} \right)^2 + B \left( \frac{\partial U}{\partial y} \right)^2+ C \left( \frac{\partial U}{\partial z} \right)^2+ D \left( \frac{\partial P}{\partial x} \right)^2+ E \frac{\partial^2 P}{\partial x^2}+ F U \frac{\partial P}{\partial x}$$$$U \frac{\partial V}{\partial x} + V \frac{\partial V}{\partial y} + W \frac{\partial V}{\partial z}+ \frac{1}{\rho}\frac{\partial P}{\partial y} - \nu \left(     \frac{\partial^2 V}{\partial x^2}    + \frac{\partial^2 V}{\partial y^2}     + \frac{\partial^2 V}{\partial z^2} \right)$$$$= A \left( \frac{\partial V}{\partial x} \right)^2 + B \left( \frac{\partial V}{\partial y} \right)^2+ C \left( \frac{\partial V}{\partial z} \right)^2+ D \left( \frac{\partial P}{\partial y} \right)^2+ E \frac{\partial^2 P}{\partial y^2}+ F V \frac{\partial P}{\partial y}$$  Terms in index notationTo be general and consistent, since we don't have any x- or z-variation$$\frac{\partial U_i}{\partial t} + U_j \frac{\partial U_i}{\partial x_j} + \frac{1}{\rho}\frac{\partial P}{\partial x_i}- \nu \frac{\partial ^2 U_i}{\partial x_j x_j}=A \frac{(\partial U_i)^2}{\partial x_j \partial x_j}+ B \frac{\partial U_j U_j}{\partial x_i}+ C \frac{\partial P^2}{\partial x_i}+ D \left( \frac{\partial P}{\partial x_i} \right)^2+ E U_j \frac{\partial P}{\partial x_j}$$ What these coefficients describe:* $A$: The square of the velocity gradient* $B$: The gradient of kinetic energy* $C$: The gradient of squared pressure* $D$: The square of the pressure gradientOther possible quantities* Absolute distance from a solid boundary<jupyter_code>%load_ext autoreload
%autoreload 2

import matplotlib
matplotlib.use("nbAgg")

import h5py
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import os
import seaborn
import jhtdb
seaborn.set()

%matplotlib inline<jupyter_output><empty_output><jupyter_text>Using data from pyJHTDB1. Pick a bunch of points randomly throughout the domain, at least more than the number of terms we want to test.2. Add points in each direction for computing spatial derivatives.3. Get $\vec{u}$, $p$, and their gradients for all time at all points in the list.4. Calculate terms based on mean values.5. Use a regression model to determine coefficients on each term.6. Repeat this process to ensure the coefficients don't change?7. Run a RANS simulation with this new model and check the results against the mean profiles.<jupyter_code>df = jhtdb.read_stats()
# Let's check continuity
div = df["dudx"] + df["dvdy"] + df["dwdz"]
div.describe()<jupyter_output><empty_output><jupyter_text>The check above gives us an idea on how accurate these gradient calculationsare.<jupyter_code>df.columns
ax = (
    df.loc[df.index.get_level_values("x")[-1000]]
    .reset_index()
    .plot(x="u", y="y", legend=False, ylabel="$U$", xlabel="$y$")
)


# Check the momentum equation
from jhtdb import nu, rho

momx_no_res = (
    df.u * df.dudx
    + df.v * df.dudy
    + df.w * df.dudz
    + (1 / rho) * df.dpdx
    - nu * (df.d2udx2 + df.d2udy2 + df.d2udz2)
)

momx = momx_no_res + (df.duudx_fd + df.duvdy_fd)

momx.dropna().describe()<jupyter_output><empty_output><jupyter_text>Check the $y$-component of the momentum equation$$\frac{\partial V}{\partial t}+ U \frac{\partial V}{\partial x}+ V \frac{\partial V}{\partial y}+ W \frac{\partial V}{\partial z}= - \frac{1}{\rho} \frac{\partial P}{\partial y}+ \nu \left(     \frac{\partial^2 V}{\partial x^2}    + \frac{\partial^2 V}{\partial y^2}    + \frac{\partial^2 V}{\partial z^2}    \right)$$$$- \left(    \frac{\partial \overline{u'v'}}{\partial x}    + \frac{\partial \overline{v'v'}}{\partial y}    + \frac{\partial \overline{v'w'}}{\partial x}    \right)$$<jupyter_code># Check the momentum equation
from jhtdb import nu, rho

momy_no_res = (
    df.u * df.dvdx
    + df.v * df.dvdy
    + df.w * df.dvdz
    + (1 / rho) * df.dpdy
    - nu * (df.d2vdx2 + df.d2vdy2 + df.d2wdz2)
)

momy = momy_no_res + (df.duvdx + df.duvdy_fd)

momy_no_res.dropna().describe()

momy_no_res.plot.hist()
# Let's see how large the Reynold's stress residual is at different locations
df1 = df.dropna().reset_index()
df1.plot.scatter(
    x="x", y="y", color=momy_no_res.dropna().values, cmap="viridis"
)


# TODO: Compute a bunch of quantities and add to the data
# Mean kinetic energy
# Squared gradients
# Gradients multiplied by each other
# Gradients multiplied by mean values
# Compute the Reynolds stress residual as a target for an ML model
# Solve a linear regression for the coefficients of all derived terms
# Throw out terms with coefficients below a threshold
# TODO: Write as a RANS model for OpenFOAM and solve this same problem there
# First run a baseline case with a high Re kOmegaSST model
# How to handle wall functions?
# TODO: Check mean flow from OpenFOAM simulation matches DNS<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15069377218457566, 'fraction_numerical': 0.012262020006453695, 'mean_word_length': 3.979116465863454, 'pattern_counts': {'":': 0, '<': 17, '<?xml version=': 0, '>': 17, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '495143', 'license': 'MIT', 'n_tokens_mistral': 2114, 'n_tokens_neox': 2047, 'n_words': 734, 'path': 'notebook.ipynb', 'repo_name': 'petebachant/boundary-layer-turbulence-modeling', 'size': '6198'}