Flow Over Cylinder¶

The simulation of a high Reynolds flow around a circular cylinder is used as validation for the HRRBGK collision operator implementation, the case is the same used in Jacob et al., 2018. The dynamic implementation of the tuning parameter is expected to provide more accurate results of velocity profiles for the HRRBGK in comparison to RRBGK.

[1]:

from nassu.cfg.model import ConfigScheme

filename = "validation/external_aero/02_flow_over_cylinder/02_flow_over_cylinder.nassu.yaml"

sim_cfgs = ConfigScheme.sim_cfgs_from_file_dct(filename)

The computational domain is periodic at \(y\)-direction, uses free slip BCs in \(z\)-direction, uniform velocity at inlet and Neumann with constant pressure at outlet.

Load experimental values for comparison:

[2]:

import pathlib

import pandas as pd

import nassu.viz as common

common.use_style()

comparison_folder = "validation/external_aero/02_flow_over_cylinder/reference"
files = [
    "ux_avg_upstream",
    "ux_rms_upstream",
    "ux_avg_102D",
    "ux_rms_102D",
    "uz_avg_102D",
    "uz_rms_102D",
    "ux_avg_154D",
    "ux_rms_154D",
    "uz_avg_154D",
    "uz_rms_154D",
    "ux_avg_202D",
    "ux_rms_202D",
    "uz_avg_202D",
    "uz_rms_202D",
]
get_filename_csv = lambda f: pathlib.Path(comparison_folder) / (f + ".csv")


df_vel = {f: pd.read_csv(get_filename_csv(f), delimiter=",") for f in files}

Results¶

[3]:

import numpy as np
import pandas as pd

from nassu.cfg.schemes.simul import SimulationConfigs

u_ref = 0.0585


def add_pos_to_points(df_points: pd.DataFrame, start_pos: float, end_pos: float):
    df_points["pos"] = np.linspace(start_pos, end_pos, num=len(df_points), endpoint=True)


def add_stats_to_points(df_points: pd.DataFrame, df_hs: dict[str, pd.DataFrame]):
    df_data = {k: df.drop(columns=["time_step"]) for k, df in df_hs.items()}
    for u_name in ["ux", "uy", "uz"]:
        df_points[f"{u_name}_avg"] = (df_data[u_name].mean() / u_ref).to_numpy(dtype=np.float32)
        df_points[f"{u_name}_rms"] = ((df_data[u_name].std() ** 2) / (u_ref**2)).to_numpy(
            np.float32
        )


def read_hs(sim_cfg: SimulationConfigs, line_name: str):
    line_ref = sim_cfg.output.exports["default_series"].series.lines[line_name]
    df_points = pd.read_csv(line_ref.points_filename)
    df_hs = {u: line_ref.read_full_data(u) for u in ["ux", "uy", "uz"]}
    return (line_ref, df_points, df_hs)

Read experimental lines and add data

[4]:

lines = {
    "upstream": {"start_pos": 0.5, "end_pos": 6.5},
    "line_102": {"start_pos": -1.5, "end_pos": 1.5},
    "line_154": {"start_pos": -1.5, "end_pos": 1.5},
    "line_202": {"start_pos": -1.5, "end_pos": 1.5},
}


def get_lines_from_sim(sim_cfg: SimulationConfigs):
    dct_sim_cfg = {}
    for line_name, dct_line in lines.items():
        hs, df_points, df_hs = read_hs(sim_cfg, line_name)
        # Filter time steps
        df_hs = {u: df[df["time_step"] > 6000] for u, df in df_hs.items()}
        # df_hs = df_hs[df_hs["time_step"] < 10000]

        add_pos_to_points(df_points, dct_line["start_pos"], dct_line["end_pos"])
        add_stats_to_points(df_points, df_hs)
        dct_sim_cfg[line_name] = {"hs": hs, "df_points": df_points}
    return dct_sim_cfg


sim_cfgs_lines = [get_lines_from_sim(sim_cfg) for sim_cfg in sim_cfgs.values()]

Set styles and legends for plots

[5]:

lg = [f"{s.models.LBM.vel_set} {s.models.LBM.coll_oper}" for s in sim_cfgs.values()]

sim_colors = [
    common.colors.sim,
    common.colors.green,
    common.colors.blue,
    common.colors.pink,
    common.colors.sky,
]

[6]:

def plot_line_values(ax, line_name: str, vel_name: str, exp_filename: str):
    for i, sim_lines in enumerate(sim_cfgs_lines):
        df_points = sim_lines[line_name]["df_points"]
        ax.plot(df_points["pos"], df_points[vel_name], color=sim_colors[i], linestyle="-")
    df_exp = df_vel[exp_filename]
    ax.plot(df_exp.iloc[:, 0], df_exp.iloc[:, 1], **common.markers.exp(shape="o"))
    ax.legend(lg + ["Exp."])

[7]:

import matplotlib.pyplot as plt

fig, ax = common.fig_double()

plot_line_values(ax[0], "upstream", "ux_avg", "ux_avg_upstream")
ax[0].set_ylabel("$u_{\mathrm{avg}}/U_{0}$")
ax[0].set_xlabel("$x/D$")


plot_line_values(ax[1], "upstream", "ux_rms", "ux_rms_upstream")
ax[1].set_ylabel("$u_{\mathrm{rms}}/U_{0}$")
ax[1].set_xlabel("$x/D$")

plt.tight_layout()
plt.show(fig)

<>:6: SyntaxWarning: invalid escape sequence '\m'
<>:11: SyntaxWarning: invalid escape sequence '\m'
<>:6: SyntaxWarning: invalid escape sequence '\m'
<>:11: SyntaxWarning: invalid escape sequence '\m'
/tmp/ipykernel_777446/3952554935.py:6: SyntaxWarning: invalid escape sequence '\m'
  ax[0].set_ylabel("$u_{\mathrm{avg}}/U_{0}$")
/tmp/ipykernel_777446/3952554935.py:11: SyntaxWarning: invalid escape sequence '\m'
  ax[1].set_ylabel("$u_{\mathrm{rms}}/U_{0}$")

../../../_images/validation_external_aero_02_flow_over_cylinder_02_flow_over_cylinder_12_1.png

[8]:

fig, ax = plt.subplots(3, 2)
fig.set_size_inches(12, 18)

for i, dist in enumerate([102, 154, 202]):
    plot_line_values(ax[i][0], f"line_{dist}", "ux_avg", f"ux_avg_{dist}D")
    ax[i][0].set_ylabel("$u_{\mathrm{avg}}/U_{0}$")
    ax[i][0].set_xlabel("$z/D$")
    ax[i][0].set_title(f"{dist / 100:.2f}D $u$")

    plot_line_values(ax[i][1], f"line_{dist}", "ux_rms", f"ux_rms_{dist}D")
    ax[i][1].set_ylabel("$u_{\mathrm{rms}}/U_{0}$")
    ax[i][1].set_xlabel("$z/D$")
    ax[i][1].set_title(f"{dist / 100:.2f}D $u$")

plt.tight_layout()
plt.show(fig)

<>:6: SyntaxWarning: invalid escape sequence '\m'
<>:11: SyntaxWarning: invalid escape sequence '\m'
<>:6: SyntaxWarning: invalid escape sequence '\m'
<>:11: SyntaxWarning: invalid escape sequence '\m'
/tmp/ipykernel_777446/1054197010.py:6: SyntaxWarning: invalid escape sequence '\m'
  ax[i][0].set_ylabel("$u_{\mathrm{avg}}/U_{0}$")
/tmp/ipykernel_777446/1054197010.py:11: SyntaxWarning: invalid escape sequence '\m'
  ax[i][1].set_ylabel("$u_{\mathrm{rms}}/U_{0}$")

../../../_images/validation_external_aero_02_flow_over_cylinder_02_flow_over_cylinder_13_1.png

[9]:

fig, ax = plt.subplots(3, 2)
fig.set_size_inches(12, 18)

for i, dist in enumerate([102, 154, 202]):
    plot_line_values(ax[i][0], f"line_{dist}", "uz_avg", f"uz_avg_{dist}D")
    ax[i][0].set_ylabel("$w_{\mathrm{avg}}/U_{0}$")
    ax[i][0].set_xlabel("$z/D$")
    ax[i][0].set_title(f"{dist / 100:.2f}D $w$")

    plot_line_values(ax[i][1], f"line_{dist}", "uz_rms", f"uz_rms_{dist}D")
    ax[i][1].set_ylabel("$w_{\mathrm{rms}}/U_{0}$")
    ax[i][1].set_xlabel("$z/D$")
    ax[i][1].set_title(f"{dist / 100:.2f}D $w$")

plt.tight_layout()
plt.show(fig)

<>:6: SyntaxWarning: invalid escape sequence '\m'
<>:11: SyntaxWarning: invalid escape sequence '\m'
<>:6: SyntaxWarning: invalid escape sequence '\m'
<>:11: SyntaxWarning: invalid escape sequence '\m'
/tmp/ipykernel_777446/2044914664.py:6: SyntaxWarning: invalid escape sequence '\m'
  ax[i][0].set_ylabel("$w_{\mathrm{avg}}/U_{0}$")
/tmp/ipykernel_777446/2044914664.py:11: SyntaxWarning: invalid escape sequence '\m'
  ax[i][1].set_ylabel("$w_{\mathrm{rms}}/U_{0}$")

../../../_images/validation_external_aero_02_flow_over_cylinder_02_flow_over_cylinder_14_1.png

Flow field¶

Instantaneous velocity magnitude on the spanwise mid-plane (plane_series.mid_span): vortex shedding behind the cylinder, one column per velocity set / collision operator variant.

[10]:

from nassu import viz

viz.enable_offscreen()

cfg0 = next(iter(sim_cfgs.values()))
body, geom = viz.read_body(cfg0, "cylinder")
view = viz.frame_body(geom, "y", slice_coord=4.0, downstream=2.0, half_extent=4.0)

panels = [
    viz.Panel(
        f"{cfg.models.LBM.vel_set} {cfg.models.LBM.coll_oper}",
        viz.PlaneSource.from_cfg(cfg, series="plane_series_roi", plane="mid_span_wake"),
        view,
    )
    for cfg in sim_cfgs.values()
]
step = panels[0].source.steps[-1]

# One streamwise (x-z) snapshot per case, coloured by velocity magnitude. A
# single side-by-side grid is too wide/short to read, so each case is rendered
# as its own (stacked) image.
for panel in panels:
    plotter = viz.render_grid(
        [panel],
        steps=[step],
        scalar="u_mag",
        cmap="viridis",
        clim=(0.0, 0.09),
        bar_title="|u|",
        bodies=[body],
    )
    plotter.show()

/tmp/ipykernel_777446/595302993.py:32: UserWarning: Using static image for notebook display.
Install trame for interactive backends: pip install "pyvista[jupyter]"
  plotter.show()

../../../_images/validation_external_aero_02_flow_over_cylinder_02_flow_over_cylinder_16_1.png

../../../_images/validation_external_aero_02_flow_over_cylinder_02_flow_over_cylinder_16_2.png

Density field¶

The same streamwise plane coloured by density. The flow is weakly compressible, so rho stays within about 1 % of the reference; a tight, data-driven colour range around the mean reveals the pressure/acoustic footprint of the shedding wake. One snapshot per case, as above.

[11]:

rho_all = np.concatenate(
    [np.asarray(panel.source.read(step).point_data["rho"]).ravel() for panel in panels]
)
rho_lo, rho_hi = np.percentile(rho_all, [1.0, 99.0])
for panel in panels:
    plotter = viz.render_grid(
        [panel],
        steps=[step],
        scalar="rho",
        cmap="coolwarm",
        clim=(float(rho_lo), float(rho_hi)),
        bar_title="rho",
        bodies=[body],
    )
    plotter.show()

/tmp/ipykernel_777446/1268032618.py:17: UserWarning: Using static image for notebook display.
Install trame for interactive backends: pip install "pyvista[jupyter]"
  plotter.show()

../../../_images/validation_external_aero_02_flow_over_cylinder_02_flow_over_cylinder_18_1.png

../../../_images/validation_external_aero_02_flow_over_cylinder_02_flow_over_cylinder_18_2.png

Version¶

[12]:

sim_cfg = next(iter(sim_cfgs.values()))
sim_info = sim_cfg.output.read_info()

nassu_commit = sim_info["commit"]
nassu_version = sim_info["version"]
print("Version:", nassu_version)
print("Commit hash:", nassu_commit)

Version: 2.0.0a0
Commit hash: aaa6995a09d88707a467e1a8e1de7fc7f229de5b

Configuration¶

[13]:

from IPython.display import Code

Code(filename=filename)

[13]:

variables:
  # Lattice freestream velocity. Used by the inlet BC, the init equation and
  # the Mach-derived NonEqTBL pressure-gradient multiplier below.
  u_inf_lattice: 0.0585
  # Lattice Mach number Ma_LBM = U_LBM / cs, cs = 1/sqrt(3). Driving the NonEqTBL
  # pressure-gradient multiplier from this value rescales the raw `dp/ds`
  # magnitude in lattice units back to the incompressible-equivalent expectation
  # of the TBL ODE - see `theory/wall_model/neq_pres_filter`.
  Ma_LBM: !math ${u_inf_lattice} * (3 ** 0.5)

simulations:
  - name: flowOverCylinder
    save_path: ./validation/external_aero/02_flow_over_cylinder/results

    n_steps: 12000

    report: {frequency: 100}

    domain:
      domain_size:
        x: 160
        y: 8
        z: 80
      block_size: 8
      bodies:
        cylinder:
          IBM:
            run: true
            cfg_use: cylinder_cfg
          geometry_path: fixture/stl/basic/cylinder_refined.stl
          small_triangles: "add"
          area: {min: 0.25, max: 1.0}
          transformation:
            scale: [1.0, 1.0, 1.0]
            translation: [23.75, -24, 38.75]
      refinement:
        static:
          default:
            volumes_refine:
              - start: [20, 0, 16]
                end: [160, 8, 64]
                lvl: 1
                is_abs: true
              - start: [20, 0, 24]
                end: [80, 8, 56]
                lvl: 2
                is_abs: true
              - start: [20, 0, 32]
                end: [48, 8, 48]
                lvl: 3
                is_abs: true
              - start: [20, 0, 36]
                end: [38, 8, 44]
                lvl: 4
                is_abs: true
              - start: [23, 0, 38]
                end: [30, 8, 42]
                lvl: 5
                is_abs: true

    data:
      export_IBM_nodes:
        ibm_exp:
          body_name: cylinder
          frequency: 5000
      body_nodes:
        # Per-source-triangle skin-friction primitives (issue #879): friction
        # velocity (`friction_u_tau`), `friction_y_plus` and the viscous-traction
        # vector, aggregated onto the cylinder source triangles and written as a
        # flat CSV/HDF time series. `nassu.viz.friction` derives tau_w / Cf and the
        # form/friction drag split from these (issue #876).
        friction:
          body_name: cylinder
          mode: triangle
          frequency: 1000
          quantities:
            - friction_u_tau
            - friction_y_plus
            - traction_x
            - traction_y
            - traction_z
            - area
            - rho_interp
      exports:
        default:
          macrs: [rho, u, omega_LES, sigma]
          interval:
            frequency: 3000
            lvl: 0
          target:
            volume: {}
          outputs:
            instantaneous: true
        default_series:
          macrs: ["rho", "u"]
          interval: {frequency: 10, lvl: 5}
          target:
            lines:
              upstream:
                dist: 0.25
                start_pos: [26.25, 4, 40]
                end_pos: [41.25, 4, 40]
              line_102:
                dist: 0.25
                start_pos: [27.65, 4, 36.25]
                end_pos: [27.65, 4, 43.75]
              line_154:
                dist: 0.25
                start_pos: [28.85, 4, 36.25]
                end_pos: [28.85, 4, 43.75]
              line_202:
                dist: 0.25
                start_pos: [30.05, 4, 36.25]
                end_pos: [30.05, 4, 43.75]
          outputs:
            instantaneous: true
        plane_series:
          macrs: ["rho", "u"]
          interval: {frequency: 1000, lvl: 0}
          target:
            planes:
              # Coarse spanwise mid-plane overview (one node apart, full
              # domain): debug-grade view of the whole flow.
              mid_span:
                axis: y
                axis_pos: 4
                dist: 1
          outputs:
            instantaneous: true
        plane_series_roi:
          macrs: ["rho", "u"]
          interval: {frequency: 100, lvl: 0}
          target:
            planes:
              # Fine plane restricted to the cylinder / near-wake region of
              # interest: vortex-shedding sequence, cheaper and
              # higher-frequency than the volume dump.
              mid_span_wake:
                axis: y
                axis_pos: 4
                min: [15, 28]
                max: [64, 52]
                dist: 0.1
          outputs:
            instantaneous: true
    models:
      precision:
        default: single

      LBM:
        tau: 0.5001125
        vel_set: !unroll [D3Q27, D3Q27]
        coll_oper: !unroll [RRBGK, HRRBGK]

      initialization:
        equations:
          rho: "1.0"
          ux: !sub "${u_inf_lattice}"
          uy: "0"
          uz: "0"
      engine:
        name: CUDA

      BC:
        periodic_dims: [false, true, false]
        BC_map:
          - pos: F
            BC: RegularizedNeumannOutlet
            rho: 1.0
            wall_normal: F
            order: 1

          - pos: B
            BC: RegularizedNeumannOutlet
            rho: 1.0
            wall_normal: B
            order: 1

          - pos: E
            BC: RegularizedNeumannOutlet
            rho: 1.0
            wall_normal: E
            order: 2

          - pos: W
            BC: UniformFlow
            wall_normal: W
            ux: !math ${u_inf_lattice}
            uy: 0
            uz: 0
            rho: 1
            order: 2

      IBM:
        dirac_delta: 3_points
        forces_accomodate_time: 0
        reset_forces: true
        body_cfgs:
          cylinder_cfg:
            n_iterations: 3
            forces_factor: 0.25
            wall_model:
              name: NonEqTBL
              dist_ref: 3.125
              dist_shell: 0.125
              start_step: 2000
              params:
                z0: 0.00155
                TDMA_max_error: 5e-06
                TDMA_max_iters: 30
                TDMA_min_div: 51
                TDMA_max_div: 51
                NeqWM_u_friction_floor: 1.0e-4
                # Magnitude-rescaling of dp/ds from the LBM weakly-compressible
                # signal to the incompressible-equivalent expectation of the
                # TBL ODE - see the theory page for the derivation.
                NeqWM_pres_grad_mult: !math ${Ma_LBM}

      multiblock:
        overlap_F2C: 2

      LES:
        model: Smagorinsky
        sgs_cte: 0.17