Periodic Heat Exchanger

Periodic Heat Exchanger#

image info

import matplotlib.pyplot as plt
import numpy as np

import xcompact3d_toolbox as x3d

Parameters#

Numerical precision

Use np.float64 if Xcompact3d was compiled with the flag -DDOUBLE_PREC, use np.float32 otherwise.

x3d.param["mytype"] = np.float64

Xcompact3d’s parameters

For more information about them, checkout the API reference.

prm = x3d.Parameters(
    filename="input.i3d",
    # BasicParam
    itype=12,
    p_row=0,
    p_col=0,
    nx=128,
    ny=129,
    nz=8,
    xlx=6.0,
    yly=6.0,
    zlz=0.375,
    nclx1=0,
    nclxn=0,
    ncly1=2,
    nclyn=2,
    nclz1=0,
    nclzn=0,
    iin=1,
    re=300.0,
    init_noise=0.0125,
    inflow_noise=0.0,
    dt=0.0025,
    ifirst=1,
    ilast=90000,
    ilesmod=1,
    iibm=2,
    gravx=0.0,
    gravy=-1.0,
    gravz=0.0,
    # NumOptions
    nu0nu=4.0,
    cnu=0.44,
    # InOutParam
    irestart=0,
    icheckpoint=45000,
    ioutput=500,
    iprocessing=100,
    # LESModel
    jles=4,
    # ScalarParam
    numscalar=1,
    nclxS1=0,
    nclxSn=0,
    nclyS1=2,
    nclySn=2,
    nclzS1=0,
    nclzSn=0,
    sc=[1.0],
    ri=[-0.25],
    uset=[0.0],
    cp=[1.0],
)

Setup#

Initial Condition#

Velocity profile: Since the boundary conditions for velocity at the top and at the bottom are no-slip in this case (ncly2=nclyn=2), the inflow profile for streamwise velocity must be zero near walls:

# This function gives the shape
fun = -((ds.y - prm.yly / 2.0) ** 2.0)

# This attribute will be shown in the figure
fun.attrs["long_name"] = r"Inflow Profile - f($x_2$)"

# Now, let's adjust its magnitude
fun -= fun.isel(y=0)
fun /= fun.x3d.simpson("y") / prm.yly

fun.plot()

[<matplotlib.lines.Line2D at 0x7f10d5736420>]

../_images/c14a4e566b91a6db83c8b2804a1f895edd0385948523b929a1afb63688624296.png

Now, let’s make sure that the dimensionless averaged velocity is unitary:

fun.x3d.simpson("y") / prm.yly

<xarray.DataArray 'y' ()> Size: 8B
array(1.)

A random noise will be applied at the three velocity components, we can create a modulation function mod to control were it will be applied. In this case, we will concentrate the noise near the center region and make it zero were \(y=0\) and \(y=L_y\). The domain is periodic in \(z\) nclz1=nclzn=0, so there is no need to make mod functions of \(z\). The functions looks like:

\[ \text{mod} = \exp\left(-0.5 (y - 0.5 L_y)^2 \right) . \]

See the code:

# Random noise with fixed seed,
# important for reproducibility, development and debugging
if prm.iin == 2:
    np.random.seed(seed=67)

mod = np.exp(-0.5 * (ds.y - prm.yly * 0.5) ** 2.0)

# This attribute will be shown in the figure
mod.attrs["long_name"] = r"Noise Modulation - f($x_2$)"

mod.plot()

[<matplotlib.lines.Line2D at 0x7f10d584c290>]

../_images/79a2ac61394c70ed534949cafaae2cb28b17a816d23630c14995dfb752cb4c1f.png

Now we reset velocity fields ds[key] *= 0.0, just to guarantee consistency in the case of multiple executions of this cell.

We then add a random number array with the right shape, multiply by the noise amplitude at the initial condition init_noise and multiply again by our modulation function mod, defined previously. Finally, we add the streamwise profile fun to ux and make the plots for reference, I’m adding extra options just to exemplify how easily we can slice the spanwise coordinate and produce multiple plots:

for key in "ux uy uz".split():
    #
    print(ds[key].attrs["name"])
    #
    ds[key] *= 0.0
    ds[key] += prm.init_noise * (np.random.random(ds[key].shape) - 0.5)
    ds[key] *= mod
    #
    if key == "ux":
        ds[key] += fun
    #
    ds[key].sel(z=slice(None, None, ds.z.size // 3)).plot(x="x", y="y", col="z", col_wrap=2, aspect=1, size=3)
    plt.show()
    #

plt.close("all")

Initial Condition for Streamwise Velocity

../_images/a495680312ce015f755c3bc51d6a0f75c3cd049eee9f07e1b35ff661a225d158.png

Initial Condition for Vertical Velocity

../_images/026266d0e7b0ab4256a2f40903a4dea48cc99ecc8b3077f292e0565056a4f3d4.png

Initial Condition for Spanwise Velocity

../_images/82118dd72ebbc97093927b39b7459d3068bd22db7467ed72d80eb12db4b67df2.png

For temperature, let’s start with one everywhere:

ds["phi"] *= 0.0
ds["phi"] += 1.0

Flow rate control#

The Sandbox flow configuration is prepared with a forcing term when the flow is periodic in the streamwise direction \(x_1\), in order to compensate viscous dissipation and keep the flow rate constant.

It is done with the help of a personalized volumetric integral operator (vol_frc), then, the streamwise velocity will be corrected at each simulated time-step as:

I = sum(vol_frc * ux)

ux = ux / I

For the composed trapezoidal rule in a uniform grid, the integral operator in the vertical direction can be computed as:

\[ vol_{frc} = \Delta y ~ [1/2, 1, \dots, 1, \dots, 1, 1/2] \]

For a unitary averaged velocity, vol_frc must be divided by the domain’s height. Besides that, for streamwise and spanwise averaged values, vol_frc must be divided by nx and nz.

Finally, vol_frc can be coded as:

ds["vol_frc"] *= 0.0

ds["vol_frc"] += prm.dy / prm.yly / prm.nx / prm.nz

ds["vol_frc"][{"y": 0}] *= 0.5
ds["vol_frc"][{"y": -1}] *= 0.5

ds.vol_frc.isel(z=0).plot(x="x", y="y", aspect=1, size=5)

<matplotlib.collections.QuadMesh at 0x7f10d603bb90>

../_images/527711e2f7ec4d13ff7ffbc76dfda4def499206a9938c5e76cb65e98a9dae458.png

Now, let’s make sure that sum(vol_frc * ux) is near to one:

(ds.vol_frc * ds.ux).sum(["x", "y", "z"])

<xarray.DataArray ()> Size: 8B
array(0.9999334)

And the last step, we can remove the solid geometry from our integral operator using the code:

ds["vol_frc"] = ds.vol_frc.where(epsi["epsi"] == False, 0.0)

ds.vol_frc.isel(z=0).plot(x="x", y="y", aspect=1, size=5)

<matplotlib.collections.QuadMesh at 0x7f10d6093b90>

../_images/c1444acafbf677ebb6daca492031fa22c78621f5f13ad8ce332e46f433accdc6.png

Writing to the disc#

is as simple as:

prm.dataset.write(ds)

prm.write()

Running the Simulation#

It was just to show the capabilities of xcompact3d_toolbox.sandbox, keep in mind the aspects of numerical stability of our Navier-Stokes solver. It is up to the user to find the right set of numerical and physical parameters.

The Sandbox Flow Configuration is still in prerelease, it can be found at fschuch/Xcompact3d.

Make sure that the compiling flags and options at Makefile are what you expect. Then, compile the main code at the root folder with make.

And finally, we are good to go:

mpirun -n [number of cores] ./xcompact3d |tee log.out