AeroSim Docs
Flow Over Mounted Cube¶
The simulation of a wall-mounted cube is used for the validation of the combined features usage. The multiblock refined for a IBM body is used, the synthetic eddy method is applied at inlet to reproduce a category II velocity profile from EU standard. The turbulent flow around a surface-mounted cube is well estabilished in literature with results of experimental data available.
We seek to reproduce correct values of the average and standard deviation of the pressure in the cube surface.
[1]:
from nassu.cfg.model import ConfigScheme
filename = "./tests/validation/cases/08_flow_over_mounted_cube.nassu.yaml"
sim_cfgs = ConfigScheme.sim_cfgs_from_file_dct(filename)
The lateral boundaries and top of domain use free slip boundary conditions, for the outlet, a Neumann BC with fixed pressure is adopted.
Results¶
The results of the pressure coefficient are compared to the average of different laboratory results. An excellent agreement was achieved. The root mean squared pressure coefficent also presented great accordance with experimental data.
[2]:
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pathlib
from tests.validation.notebooks import common
common.use_style()
def get_experimental_profile_Cp(reynolds: float) -> dict[str, dict[str, pd.DataFrame]]:
files_tau: dict[float, dict[str, dict[str, str]]] = {
20000: {
"frontal": {
"Lim": "Wall_mounted_cube/Re_20000/pressure_coefficient_frontal_lim2009.csv",
"Holscher": "Wall_mounted_cube/Re_20000/pressure_coefficient_frontal_holscher1998.csv",
},
"side": {
"Lim": "Wall_mounted_cube/Re_20000/pressure_coefficient_side_lim2009.csv",
},
"rms": {
"Holscher": "Wall_mounted_cube/Re_20000/pressure_coefficient_rms_top_lim2009.csv",
},
"skew": {
"Holscher": "Wall_mounted_cube/Re_20000/pressure_coefficient_skew_top_lim2009.csv",
},
}
}
files_get = files_tau[reynolds]
vals_exp: dict[str, dict[str, pd.DataFrame]] = {}
for type_file, dct in files_get.items():
vals_exp[type_file] = {}
for name, f in dct.items():
filename = pathlib.Path("tests/validation/comparison") / f
# ([pos], [Cp])
df = pd.read_csv(filename, delimiter=",")
vals_exp[type_file][name] = df
return vals_exp
exp = get_experimental_profile_Cp(20000)
holscher_frontal = exp["frontal"]["Holscher"]
holscher_top_rms = exp["rms"]["Holscher"]
holscher_top_skew = exp["skew"]["Holscher"]
lim_side = exp["side"]["Lim"]
Calculation of reference velocity:
[3]:
sim_cfgs.keys()
[3]:
dict_keys([('cubeFlowMultilevelLES', 0), ('cubeFlowMultiLevelLESnoCube', 0)])
[4]:
sim_cfg = sim_cfgs["cubeFlowMultiLevelLESnoCube", 0]
line_ref = sim_cfg.output.series["velocity"].points["pitot"]
df_ref = line_ref.read_full_data("ux")
df_ref = df_ref[df_ref["time_step"] > 10000]
ux_avg = df_ref.mean()
u_ref = ux_avg[1]
u_ref
/tmp/ipykernel_1558770/689290571.py:9: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
u_ref = ux_avg[1]
[4]:
np.float32(0.050575547)
[5]:
import scipy
from scipy.ndimage import gaussian_filter
from scipy.stats import skew
[6]:
sim_cfg = sim_cfgs["cubeFlowMultilevelLES", 0]
def get_cp_line(line_names: list[str]) -> pd.DataFrame:
line_series = sim_cfg.output.series["pressure"].lines
pressure_series = sim_cfg.output.series["pressure"].points["pressure_point"]
df_pref = pressure_series.read_full_data("rho")
df_pref = df_pref[df_pref["time_step"] > 10000]
rho_ref = df_pref["0"].to_numpy().T
df_ret = pd.DataFrame({})
for i, line in enumerate(line_names):
df_hs = line_series[line].read_full_data("rho")
df_hs = df_hs[df_hs["time_step"] > 10000]
df_rho = df_hs.drop(columns=["time_step"])
df_rho["rho_ref"] = rho_ref
df_cp = (df_rho.subtract(df_rho["rho_ref"], axis=0)) / (1.5 * (u_ref**2))
df_cp = df_cp.drop(columns=["rho_ref"])
cp_avg = df_cp.mean().to_numpy().T
cp_rms = df_cp.std().to_numpy().T
cp_skew = skew(df_cp).T
df_cp = pd.DataFrame({"cp_avg": cp_avg, "cp_rms": cp_rms, "cp_skew": cp_skew})
df_cp["pos"] = np.linspace(i, i + 1, num=len(cp_avg), endpoint=True)
df_ret = pd.concat([df_ret, df_cp])
return df_ret
lines_cp_frontal = [f"long_line{i}" for i in range(1, 4)]
df_cp_frontal = get_cp_line(lines_cp_frontal)
lines_cp_side = [f"tr_line{i}" for i in range(1, 4)]
df_cp_side = get_cp_line(lines_cp_side)
[7]:
fig, ax = common.fig_triple()
ax[0].plot(
holscher_frontal["x/h"],
holscher_frontal["Cp"],
**common.markers.exp(shape="o"),
label="Holscher (1998)",
)
ax[0].plot(
df_cp_frontal["pos"],
df_cp_frontal["cp_avg"],
**common.markers.sim_line(linestyle="-"),
label="AeroSim",
)
ax[0].set_ylim(-1.5, 1.5)
ax[0].set_ylabel(r"$C_{p,\mathrm{avg}}$")
ax[0].set_xlabel(r"$x/L$")
ax[0].legend()
n_points = len(df_cp_frontal["pos"])
df_top = df_cp_frontal[(df_cp_frontal["pos"] >= 1.0) & (df_cp_frontal["pos"] <= 2.0)]
ax[1].plot(
holscher_top_rms["x/h"],
holscher_top_rms["Cp"],
**common.markers.exp(shape="o"),
label="Lim (2009)",
)
ax[1].plot(
df_top["pos"], df_top["cp_rms"], **common.markers.sim_line(linestyle="-"), label="AeroSim"
)
ax[1].set_ylabel(r"$C_{p,\mathrm{rms}}$")
ax[1].set_xlabel(r"$x/L$")
ax[1].legend()
ax[1].set_ylim(0, 0.6)
ax[1].set_xlim(1, 2)
ax[2].plot(
holscher_top_rms["x/h"],
holscher_top_skew["Cp"],
**common.markers.exp(shape="o"),
label="Lim (2009)",
)
ax[2].plot(
df_top["pos"], df_top["cp_skew"], **common.markers.sim_line(linestyle="-"), label="AeroSim"
)
ax[2].set_ylabel(r"$C_{p,\mathrm{skew}}$")
ax[2].set_xlabel(r"$x/L$")
ax[2].legend()
ax[2].set_ylim(-1.5, 1)
ax[2].set_xlim(1, 2)
plt.tight_layout()
plt.show(fig)
Good accuracy of results were obtained against experimental data of the average and root mean squared pressure coefficient towards longitudinal direction. The results from Holscher, 1998 consist of an average of multiple wind tunnel experiments of a turbulent flow over a surface-mounted cube.
[8]:
fig, ax = common.fig_single()
ax.plot(lim_side["x/h"], lim_side["Cp"], **common.markers.exp(shape="o"), label="Lim (2009)")
ax.plot(
df_cp_side["pos"],
df_cp_side["cp_avg"],
**common.markers.sim_line(linestyle="-"),
label="AeroSim",
)
ax.set_ylabel(r"$C_{p,\mathrm{avg}}$")
ax.set_xlabel(r"$x/L$")
ax.set_ylim(-1.5, 1)
ax.legend()
plt.tight_layout()
plt.show(fig)
Good agreement was also observed along the transversal perimeter of the surface-mounted cube.
Power Spectral Density¶
The power spectral density is also computed in order to check if there are anomalies in the velocity fluctuations.
[9]:
import pandas as pd
[10]:
sim_cfg = sim_cfgs["cubeFlowMultilevelLES", 0]
def read_df_spectrum(point_name: str):
spectrum_series = sim_cfg.output.spectrum[point_name]
df = pd.read_hdf(spectrum_series.data_filename)
df = df[df["time_step"] >= 10000]
return df
df_sa_pressure = read_df_spectrum("point_pressure_top")
df_sa_velocity_top = df_sa_pressure
df_sa_velocity_wake = read_df_spectrum("point_velocity_wake")
[11]:
Lu = 2.8
f = 1
rho_top = df_sa_pressure["rho"].tolist()
rho_top = np.array(rho_top, dtype=np.float32)
xf_1, yf_img = scipy.signal.periodogram(rho_top, 32, scaling="density")
avg = np.average(rho_top)
st = np.std(rho_top)
yf_img = yf_img * (u_ref / Lu) / (st * st)
xf_1 = xf_1 * (Lu / u_ref)
yf_1 = np.real(yf_img)
ux_top = df_sa_velocity_top["ux"].tolist()
ux_top = np.array(ux_top, dtype=np.float32)
xf_2, yf_img = scipy.signal.periodogram(ux_top, f, scaling="density")
avg = np.average(ux_top)
st = np.std(ux_top)
yf_img = yf_img * (u_ref / Lu) / (st * st)
xf_2 = xf_2 * (Lu / u_ref)
yf_2 = np.real(yf_img)
ux_wake = df_sa_velocity_wake["ux"].tolist()
ux_wake = np.array(ux_wake, dtype=np.float32)
xf_3, yf_img = scipy.signal.periodogram(ux_wake, f, scaling="density")
avg = np.average(ux_wake)
st = np.std(ux_wake)
yf_img = yf_img * (u_ref / Lu) / (st * st)
xf_3 = xf_3 * (Lu / u_ref)
yf_3 = np.real(yf_img)
[12]:
def filter_avg_data(data):
data = gaussian_filter(data, sigma=7)
return data
yf_avg_1 = filter_avg_data(yf_1)
yf_avg_2 = filter_avg_data(yf_2)
yf_avg_3 = filter_avg_data(yf_3)
[13]:
fig, ax = common.fig_triple()
ax[0].set_title("Density top")
ax[0].plot(xf_1, yf_1, color=common.colors.sim, alpha=0.7, label="AeroSim")
ax[0].set_ylabel(r"$f\,S_{uu}/ \sigma_{u}^{2}$")
ax[0].set_xlabel(r"$f\,L_{u}/u$")
ax[0].set_xscale("log")
ax[0].set_yscale("log")
ax[0].legend()
ax[1].set_title("Velocity top")
ax[1].plot(xf_2, yf_2, color=common.colors.sim, alpha=0.7, label="AeroSim")
ax[1].set_xlabel(r"$f\,L_{u}/u$")
ax[1].set_xscale("log")
ax[1].set_yscale("log")
ax[1].legend()
ax[2].set_title("Velocity wake")
ax[2].plot(xf_3, yf_3, color=common.colors.sim, alpha=0.7, label="AeroSim")
ax[2].set_xlabel(r"$f\,L_{u}/u$")
ax[2].set_xscale("log")
ax[2].set_yscale("log")
ax[2].legend()
plt.tight_layout()
plt.show(fig)
Version¶
[14]:
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: 1.6.45
Commit hash: a2287d74877f61304fc618ad5e2c77a25b3f7f39
Configuration¶
[15]:
from IPython.display import Code
Code(filename=filename)
[15]:
variables:
simul:
dev_time: 10000
stats_time: 30000
scale: !math 1/16
plane_height: 5.01
sigma_sem: 20
body_pos: 150
domain:
length: 600
width: 160
height: 96
var:
cube_length: 2.8
body_lvl: 5
aux:
offset: !math 2*(1/(2**${var.body_lvl}))
ref_delta: 0.02
simulations:
- name: cubeFlowMultilevelLES
save_path: ./tests/validation/results/08_flow_over_mounted_cube/
n_steps: !math "${simul.dev_time}+${simul.stats_time}"
report:
frequency: 1000
data:
divergence: { frequency: 1 }
monitors:
fields:
macrs_stats:
macrs: [rho, u]
stats: [min, max, mean]
interval: { start_step: 0, frequency: 10 }
instantaneous:
full_domain: { interval: { frequency: 0 }, macrs: [rho, u] }
export_IBM_nodes:
ibm_exp:
body_name: cube
frequency: 5000
probes:
historic_series:
pressure:
macrs: ["rho"]
interval: { frequency: 1, lvl: 0 }
points:
pressure_point:
pos:
!math [
"${simul.body_pos}",
"0.5*${domain.width}",
"0.5*${domain.height}-32",
]
lines:
long_line1:
dist: 0.0625
start_pos:
!math [
"${simul.body_pos} - 0.5*${var.cube_length} - ${aux.offset}",
"0.5*${domain.width}",
"${simul.plane_height}",
]
end_pos:
!math [
"${simul.body_pos} - 0.5*${var.cube_length} - ${aux.offset}",
"0.5*${domain.width}",
"${simul.plane_height} + ${var.cube_length}",
]
long_line2:
dist: 0.0625
start_pos:
!math [
"${simul.body_pos} - 0.5*${var.cube_length}",
"0.5*${domain.width}",
"${simul.plane_height} + ${var.cube_length} + ${aux.offset}",
]
end_pos:
!math [
"${simul.body_pos} + 0.5*${var.cube_length}",
"0.5*${domain.width}",
"${simul.plane_height} + ${var.cube_length} + ${aux.offset}",
]
long_line3:
dist: 0.0625
start_pos:
!math [
"${simul.body_pos} + 0.5*${var.cube_length} + ${aux.offset}",
"0.5*${domain.width}",
"${simul.plane_height} + ${var.cube_length}",
]
end_pos:
!math [
"${simul.body_pos} + 0.5*${var.cube_length} + ${aux.offset}",
"0.5*${domain.width}",
"${simul.plane_height}",
]
tr_line1:
dist: 0.0625
start_pos:
!math [
"${simul.body_pos}",
"0.5*${domain.width} - 0.5*${var.cube_length} - ${aux.offset}",
"${simul.plane_height}",
]
end_pos:
!math [
"${simul.body_pos}",
"0.5*${domain.width} - 0.5*${var.cube_length} - ${aux.offset}",
"${simul.plane_height} + ${var.cube_length}",
]
tr_line2:
dist: 0.0625
start_pos:
!math [
"${simul.body_pos}",
"0.5*${domain.width} - 0.5*${var.cube_length}",
"${simul.plane_height} + ${var.cube_length} + ${aux.offset}",
]
end_pos:
!math [
"${simul.body_pos}",
"0.5*${domain.width} + 0.5*${var.cube_length}",
"${simul.plane_height} + ${var.cube_length} + ${aux.offset}",
]
tr_line3:
dist: 0.0625
start_pos:
!math [
"${simul.body_pos}",
"0.5*${domain.width} + 0.5*${var.cube_length} + ${aux.offset}",
"${simul.plane_height} + ${var.cube_length}",
]
end_pos:
!math [
"${simul.body_pos}",
"0.5*${domain.width} + 0.5*${var.cube_length} + ${aux.offset}",
"${simul.plane_height}",
]
spectrum_analysis:
macrs: ["rho", "u"]
points:
point_pressure_top:
pos:
!math [
"${simul.body_pos}",
"0.5*${domain.width}",
"${simul.plane_height} + ${var.cube_length} + ${aux.offset}",
]
point_velocity_top:
pos:
!math [
"${simul.body_pos}",
"0.5*${domain.width}",
"${simul.plane_height} + 1.14*${var.cube_length}",
]
point_velocity_wake:
pos:
!math [
"${simul.body_pos} + 1.5*${var.cube_length}",
"0.5*${domain.width} + 0.5*${var.cube_length}",
"${simul.plane_height} + ${var.cube_length}",
]
domain:
domain_size:
x: !math "${domain.length}"
y: !math "${domain.width}"
z: !math "${domain.height}"
block_size: 8
bodies:
cube:
IBM:
cfg_use: body_wm
order: 2
lnas_path: fixture/lnas/basic/cube_no_floor.lnas
small_triangles: "add"
transformation:
scale:
!math [
"0.1*${var.cube_length}",
"0.1*${var.cube_length}",
"0.1*${var.cube_length}",
]
translation:
!math [
"${simul.body_pos} - 0.5*${var.cube_length}",
"0.5*${domain.width} - 0.5*${var.cube_length}",
"${simul.plane_height}",
]
full_plane:
IBM:
cfg_use: terrain_wm
lnas_path: fixture/lnas/wind_tunnel/full_plane.lnas
small_triangles: "add"
transformation:
translation: !math [0, 0, "${simul.plane_height}"]
obstacles_category_II:
IBM:
order: 1
lnas_path: fixture/lnas/wind_tunnel/category_II/plates_Nx160Ny70_6x2_spacing16x32_offset19y.lnas
volumes_limits:
body_transformed:
- start: !math [4, 20, 0]
end:
!math [
"${simul.body_pos} - 0.75*${var.cube_length}",
"${domain.width} - 20",
"${domain.height}",
]
small_triangles: "add"
transformation:
scale:
!math [
"64.0*${simul.scale}",
"64.0*${simul.scale}",
"32.0*${simul.scale}",
]
translation: !math [0, 0, "${simul.plane_height}"]
refinement:
static:
lvl1:
volumes_refine:
- start: [0.0, 55.0, 0.0]
end: [272.0, 105.0, 32.0]
lvl: 1
is_abs: true
lvl2:
volumes_refine:
- start: [0.0, 67.5, 0.0]
end: [208.0, 92.5, 20.0]
lvl: 2
is_abs: true
lvl3:
volumes_refine:
- start: [0.0, 73.75, 2.0]
end: [176.0, 86.25, 14.0]
lvl: 3
is_abs: true
lvl4:
volumes_refine:
- start: [140.0, 76.875, 4.0]
end: [160.0, 83.125, 11.0]
lvl: 4
is_abs: true
body_refinement:
bodies:
- body_name: cube
lvl: 5
normal_offsets: !range [-0.25, 0.751, 0.125]
transformation:
translation: !math [0, 0, "-${aux.ref_delta}"]
- body_name: cube
lvl: 5
normal_offsets: !range [-0.25, 0.751, 0.25]
models:
precision:
default: single
LBM:
tau: 0.50002
vel_set: D3Q27
coll_oper: RRBGK
initialization:
sem_field: true
engine:
name: CUDA
BC:
periodic_dims: [false, false, false]
BC_map:
- pos: E
BC: RegularizedNeumannOutlet
rho: 1.0
wall_normal: E
order: 2
- pos: F
BC: Neumann
wall_normal: F
order: 1
- pos: B
BC: RegularizedHWBB
wall_normal: B
order: 1
- pos: N
BC: Neumann
wall_normal: N
order: 0
- pos: S
BC: Neumann
wall_normal: S
order: 0
SEM:
eddies:
lengthscale:
{
x: !sub "${simul.sigma_sem}",
y: !sub "${simul.sigma_sem}",
z: !sub "${simul.sigma_sem}",
}
eddies_vol_density: 300
domain_limits_yz:
start: !math ["${simul.sigma_sem}", "-${simul.sigma_sem}"]
end:
!math [
"${domain.width}-${simul.sigma_sem}",
"${domain.height}-${simul.sigma_sem}",
]
profile:
csv_profile_data: "fixture/SEM/category_vprofile/profile_log_cat2_H150_Uh0.06.csv"
z_offset: !math "${simul.plane_height}"
length_mul: !math "${simul.scale}"
LES:
model: Smagorinsky
sgs_cte: 0.17
IBM:
forces_accomodate_time: 0
dirac_delta: "3_points"
body_cfgs:
default:
n_iterations: 3
forces_factor: 1.0
terrain_wm:
n_iterations: 3
forces_factor: 0.5
kinetic_energy_correction: false
wall_model:
name: EqLog
dist_ref: 3.125
dist_shell: 0.125
start_step: 5000
params:
z0: !math "0.05*${simul.scale}"
TDMA_max_error: 1e-04
TDMA_max_iters: 10
TDMA_max_div: 51
body_wm:
n_iterations: 3
forces_factor: 0.5
wall_model:
name: EqTBL
dist_ref: 3.125
dist_shell: 0.125
start_step: 10000
params:
z0: !math "0.05*${simul.scale}"
TDMA_max_error: 1e-04
TDMA_max_iters: 10
TDMA_min_div: 51
TDMA_max_div: 51
- name: cubeFlowMultiLevelLESnoCube
parent: cubeFlowMultilevelLES
run_simul: true
n_steps: !math "${simul.dev_time}+${simul.stats_time}"
data:
divergence: { frequency: 1 }
export_IBM_nodes: !not-inherit {}
probes: !not-inherit
historic_series:
velocity:
macrs: ["u"]
interval: { frequency: 1, lvl: 0 }
points:
pitot:
pos:
!math [
"${simul.body_pos}",
"0.5*${domain.width}",
"${simul.plane_height} + ${var.cube_length}",
]
lines:
velocity_profile:
dist: 0.0625
start_pos:
!math [
"${simul.body_pos}",
"0.5*${domain.width}",
"${simul.plane_height}",
]
end_pos:
!math [
"${simul.body_pos}",
"0.5*${domain.width}",
"${simul.plane_height} + ${var.cube_length}",
]
domain:
bodies: !not-inherit
full_plane:
IBM:
cfg_use: terrain_wm
lnas_path: fixture/lnas/wind_tunnel/full_plane.lnas
small_triangles: "add"
transformation:
translation: !math [0, 0, "${simul.plane_height}"]
obstacles_category_II:
IBM:
order: 1
lnas_path: fixture/lnas/wind_tunnel/category_II/plates_Nx160Ny70_6x2_spacing16x32_offset19y.lnas
volumes_limits:
body_transformed:
- start: !math [4, 20, 0]
end:
!math [
"${simul.body_pos} - 0.75*${var.cube_length}",
"${domain.width} - 20",
"${domain.height}",
]
small_triangles: "add"
transformation:
scale:
!math [
"64.0*${simul.scale}",
"64.0*${simul.scale}",
"32.0*${simul.scale}",
]
translation: !math [0, 0, "${simul.plane_height}"]
refinement: !not-inherit
static:
lvl1:
volumes_refine:
- start: [0.0, 55.0, 0.0]
end: [272.0, 105.0, 32.0]
lvl: 1
is_abs: true
lvl2:
volumes_refine:
- start: [0.0, 67.5, 0.0]
end: [208.0, 92.5, 20.0]
lvl: 2
is_abs: true
lvl3:
volumes_refine:
- start: [0.0, 73.75, 2.0]
end: [176.0, 86.25, 12.0]
lvl: 3
is_abs: true
lvl4:
volumes_refine:
- start: [140.0, 76.875, 4.0]
end: [160.0, 83.125, 8.0]
lvl: 4
is_abs: true