OpenFoam turbulence model crashing when setting fluid velocity using preCICE

Hi,

I am attempting to simmulate the evaporation of nitrogen in a tank. I use openFoam to simmulate the gas, and precice to model the heat flux into the gas from the side-walls, and to set the inlet velocity of gas (wich emulates the evaporation of liquified nitrogen).

For now i am setting the heat flux into the gas from the walls using HTC boundary conditions, and setting the inflow from the evaporating liquid. The inflow velocity is set as a constant + a contribution based on the heat flux at the fluid-gas boundary. To validate the precice results, we have made a simulation setting the same inflow and heat transfer boundary conditions using codedFixedValue and externalWallHeatFluxTemperature boundary conditions, respectively in openFoam. See the figure below for a schematic of the setup

Whithout turbulence, the preCICE and refrence simulation give very similar results, as seen from the plot of the total heat flux accross the fluid-gas boundary:

and i think the difference is due to some slightly different constants. However, when i turn on turbulence the refrence simmulation works fine, while the precice simulation crashes after 4 timesteps, with the follwoing error message

#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::sigFpe::sigHandler(int) at ??:?
#2  ? in /lib/x86_64-linux-gnu/libc.so.6
#3  Foam::divide(Foam::Field<double>&, Foam::UList<double> const&, Foam::UList<double> const&) at ??:?
#4  void Foam::divide<Foam::fvPatchField, Foam::volMesh>(Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&) in ~/openfoam/OpenFOAM-v2206/platforms/linux64GccDPInt32Opt/bin/buoyantPimpleFoam
#5  Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > Foam::operator/<Foam::fvPatchField, Foam::volMesh>(Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > const&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&) at ??:?
#6  Foam::RASModels::realizableKE<Foam::EddyDiffusivity<Foam::ThermalDiffusivity<Foam::CompressibleTurbulenceModel<Foam::fluidThermo> > > >::correct() at ??:?
#7  ? in ~/openfoam/OpenFOAM-v2206/platforms/linux64GccDPInt32Opt/bin/buoyantPimpleFoam
#8  ? in /lib/x86_64-linux-gnu/libc.so.6
#9  __libc_start_main in /lib/x86_64-linux-gnu/libc.so.6
#10  ? in ~/openfoam/OpenFOAM-v2206/platforms/linux64GccDPInt32Opt/bin/buoyantPimpleFoam
Floating point exception (core dumped)
/home/eirik/code_dir/tank-insulation-coupling-2D-modifiedInlet/lh2-tank-gas-phase-2D-precice/post_processing.py:73: RuntimeWarning: invalid value encountered in scalar divide
  return np.trapz(values[mask], times[mask])/(times[mask][-1] - times[mask][0])
/home/eirik/code_dir/tank-insulation-coupling-2D-modifiedInlet
removing precice-run

However, the inflow velocity in the first timesteps looks ok:

Do you have any idea what might be causing the turbulence model to crash? now i am using realizeableKE, but i have also tired kEpsilon, kOmega and kOmegaSST with the same result.

Best regards,
Eirik Høydalsvik

I attach the precice-config, preciceDict and the python model for insulation and heat flow using precice python API (insulation-htc.py):

precice-config:

<?xml version="1.0"?>
<precice-configuration>
  <log>
    <sink
      filter="%Severity% > debug and %Rank% = 0"
      format="---[precice] %ColorizedSeverity% %Message%"
      enabled="true" />
  </log>

  <solver-interface dimensions="3">
    <data:scalar name="Temperature" />
    <data:scalar name="Heat-Flux-wall" />
    <data:scalar name="Heat-Flux-bottom" />
    <data:vector name="Velocity" />

    <mesh name="Fluid-Mesh">
      <use-data name="Temperature" />
      <use-data name="Heat-Flux-wall" />
      <use-data name="Heat-Flux-bottom" />
      <use-data name="Velocity" />
    </mesh>

    <mesh name="Solid-Mesh">
      <use-data name="Temperature" />
      <use-data name="Heat-Flux-wall" />
      <use-data name="Heat-Flux-bottom" />
      <use-data name="Velocity" />
    </mesh>

    <participant name="Fluid">
      <use-mesh name="Fluid-Mesh" provide="yes" />
      <use-mesh name="Solid-Mesh" from="Solid" />
      <read-data name="Heat-Flux-wall" mesh="Fluid-Mesh" />
      <read-data name="Velocity" mesh="Fluid-Mesh" />
      <write-data name="Temperature" mesh="Fluid-Mesh" />
      <write-data name="Heat-Flux-bottom" mesh="Fluid-Mesh" />
      <mapping:nearest-neighbor
        direction="read"
        from="Solid-Mesh"
        to="Fluid-Mesh"
        constraint="consistent" />
    </participant>

    <participant name="Solid">
      <export:vtk directory="preCICE-output" every-n-time-windows="1000"/>
      <use-mesh name="Fluid-Mesh" from="Fluid" />
      <use-mesh name="Solid-Mesh" provide="yes" />
      <mapping:nearest-neighbor
        direction="read"
        from="Fluid-Mesh"
        to="Solid-Mesh"
        constraint="consistent" />
      <read-data name="Temperature" mesh="Solid-Mesh" />
      <read-data name="Heat-Flux-bottom" mesh="Solid-Mesh" />
      <write-data name="Heat-Flux-wall" mesh="Solid-Mesh" />
      <write-data name="Velocity" mesh="Solid-Mesh" />
    </participant>

    <m2n:sockets from="Fluid" to="Solid" exchange-directory="../" />

     <coupling-scheme:serial-explicit>
      <time-window-size value="-1" method="first-participant"/>
      <participants first="Fluid" second="Solid" />
      <exchange data="Temperature" mesh="Fluid-Mesh" from="Fluid" to="Solid" />
      <exchange data="Heat-Flux-wall" mesh="Solid-Mesh" from="Solid" to="Fluid" />
      <exchange data="Heat-Flux-bottom" mesh="Fluid-Mesh" from="Fluid" to="Solid" />
      <exchange data="Velocity" mesh="Solid-Mesh" from="Solid" to="Fluid" initialize="yes"/>
  </coupling-scheme:serial-explicit>
  </solver-interface>
</precice-configuration>

preCICEdict

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      preciceDict;
}

preciceConfig "../precice-config.xml";

participant Fluid;

modules (CHT FF);

interfaces
{
  Interface1
  {
    mesh              Fluid-Mesh;
    patches           (walls);
    
    readData
    (
      Heat-Flux-wall
    );
    
    writeData
    (
      Temperature
    );
  };
  bottomInterface
  {
    mesh Fluid-Mesh;
    patches (bottom);

  readData
  (
    Velocity
  );

  writeData
  (
    Heat-Flux-bottom
  );
  }
};

insulation-htc.py

import numpy as np
import precice
import matplotlib.pyplot as plt
import pandas as pd
from pytank.case_tanks import NASADemo1m3Perlite
# TODO
# 1. rydd opp i denne filen
# 2. Finn grunn til avvik ved laminær strømning
# 3. Les om openfoam turbulens crash

# ------------------------------------------------------------
# Setting up tank geometry
# ------------------------------------------------------------

R = 8.31446261815324  # ideal gas constant J/(mol*K)

TANK = NASADemo1m3Perlite()
T_L = 77.2  # LN2 temperature (K)
M_N2 = 0.0280134  # Molar mass of N2 (kg/mol)
P_GAS = 100000  # Pa
RHO_N2 = P_GAS * M_N2 / (R * T_L)  # how to calc???


R_INNER = 0.6203504908994001  # inner radius of tank
H_FILL = 0.329494  # distance from origin to gas-liquid interface
A = TANK.get_internal_area()
# TODO: regn ut Q_L fra pytank
FILL_RATIO = 0.86  # ration liquid to total volume
A_CARGO = TANK.get_wetted_area(FILL_RATIO)

# setting up vector to recive temperature data from Fluid
N_THETA = 30
T_AMB = 295  # Kelvin
T_AMB_ARR = np.ones(N_THETA)*T_AMB

# calculation of heat transfer coefficient
Q_WALL_SAT = abs(
    TANK.calc_heat_flow(T_L, T_AMB) / A
)
HTC = Q_WALL_SAT / (
    T_AMB - T_L
)

THETA_BOTTOM = np.arccos(H_FILL/R_INNER)
# scale factor for 2d-wedge
WEDGE_FRACTION = 1/360
A_WEDGE = np.pi * (R_INNER * np.sin(THETA_BOTTOM))**2 * WEDGE_FRACTION

Q_L = Q_WALL_SAT * A_CARGO * WEDGE_FRACTION
DH_VAP = 199323.0
# constant mass flux from liquid to gas, caused by heat flow from wall to liquid
mdot_L = Q_L / DH_VAP

# TODO trenger bare 1 phi nå som vi ser på wedge
phi = np.array([0])
theta = np.linspace(0, THETA_BOTTOM, N_THETA)

x = R_INNER * np.outer(np.sin(theta), np.cos(phi)).flatten()
y = R_INNER * np.outer(np.cos(theta), np.ones(1)).flatten()
z = np.zeros_like(x)

# ------------------------------------------------------------
# Setting up precice
# ------------------------------------------------------------


participant_name = "Solid"
# Definerer hvilken participant dette er, samt hvordan fil de skal "kommunisere" gjennom
config_file_name = "../precice-config.xml"
solver_process_index = 0
solver_process_size = 1

interface = precice.Interface(
    participant_name, config_file_name, solver_process_index, solver_process_size)

mesh_name = "Solid-Mesh"  # Definerer meshet.
mesh_id = interface.get_mesh_id(mesh_name)

# wall interface
temp_name = "Temperature"
# De parametrene som skal brukes, se precice-config.xml hvordan de kommuniserer.
heatFluxWall_name = "Heat-Flux-wall"
temp_id = interface.get_data_id(temp_name, mesh_id)
heatFluxWall_id = interface.get_data_id(heatFluxWall_name, mesh_id)

verticesWall = np.array([x, y, z]).T
vertexWall_ids = interface.set_mesh_vertices(mesh_id, verticesWall)

# bottom interface
N_BOTTOM = 1000
velocity_name = "Velocity"
heatFluxBottom_name = "Heat-Flux-bottom"
velocity_id = interface.get_data_id(velocity_name, mesh_id)
heatFluxBottom_id = interface.get_data_id(heatFluxBottom_name, mesh_id)


dx_bottom = R_INNER * np.sin(THETA_BOTTOM) / N_BOTTOM
x_bottom = np.arange(dx_bottom/2, R_INNER *
                     np.sin(THETA_BOTTOM) - dx_bottom / 2, dx_bottom)
# x_bottom = np.linspace(0, np.cos(THETA_BOTTOM), N_BOTTOM)
y_bottom = np.ones_like(x_bottom)*H_FILL
z_bottom = np.zeros_like(x_bottom)
vertexBottom = np.array([x_bottom, y_bottom, z_bottom]).T
vertexBottom_ids = interface.set_mesh_vertices(mesh_id, vertexBottom)
# sf = dx_bottom * x_bottom * WEDGE_FRACTION * 2 * np.pi


# ------------------------------------------------------------
# Running simulation
# ------------------------------------------------------------

dt = 1e-8
t = 0
precice_dt = interface.initialize()

# initialize bottom velocity with contribution from insulation to liquid
if interface.is_action_required(precice.action_write_initial_data()):
    print("initializing velocity on bottom interface")
    ubottom_y = 1 / (RHO_N2) * mdot_L/A_WEDGE * np.ones_like(x_bottom)
    print("ubottom_y", ubottom_y)
    ubottom = np.array(
        [np.zeros_like(ubottom_y), ubottom_y, np.zeros_like(ubottom_y)]).T
    interface.write_block_vector_data(velocity_id, vertexBottom_ids, ubottom)
    interface.mark_action_fulfilled(precice.action_write_initial_data())

interface.initialize_data()
while interface.is_coupling_ongoing():
    # print("start loop")
    if interface.is_action_required(precice.action_write_iteration_checkpoint()):
        interface.mark_action_fulfilled(
            precice.action_write_iteration_checkpoint())
        t_checkpoint = t

    dt = np.minimum(dt, precice_dt)
    print("dt", dt, "precice")
    t += dt
    T = interface.read_block_scalar_data(
        temp_id, vertexWall_ids)  # leser temperatur inn i T-array

    # read velocity and write mass flux from/to precice
    qbottom = interface.read_block_scalar_data(
        heatFluxBottom_id, vertexBottom_ids)  # read heat flux from bottom

    # Usikker på fortegn hær
    mflux_LG = qbottom/DH_VAP
    # consider calculating density based on precice
    # ubottom_y = np.ones_like(qbottom) * -1.6515750578e-07
    ubottom_y = 1 / (RHO_N2) * (mdot_L/A_WEDGE + mflux_LG)
    ubottom = np.array(
        [np.zeros_like(ubottom_y), ubottom_y, np.zeros_like(ubottom_y)]).T
    # // vectorField U = -(sf/Af)*mdot_L/(rhopatch*gSum(Af)) - (sf/Af)*mflux_LG/(rhopatch);

    # print("sf / Af")
    # print(sf/A_WEDGE)
    # print("mdot_L")
    # print(mdot_L)
    # print("rhopatch")
    # print(RHO_N2 * (sf/A_WEDGE))
    print("RHON2")
    print(RHO_N2)
    print("mflux LG")
    print(mflux_LG)
    print("qbottom")
    print(qbottom)
    print("ubottom")
    print(ubottom)

    interface.write_block_vector_data(
        velocity_id, vertexBottom_ids, ubottom)

    # calculate heat-flux on wall and write to precice
    qwall = -HTC * (T - T_AMB_ARR)

    interface.write_block_scalar_data(heatFluxWall_id, vertexWall_ids, qwall)
    if interface.is_action_required(precice.action_read_iteration_checkpoint()):
        t = t_checkpoint
        interface.mark_action_fulfilled(
            precice.action_read_iteration_checkpoint())

    precice_dt = interface.advance(dt)  # Sender/mottar, går til neste tid
    dt = 0.001

interface.finalize()

I think i figured out my own problem:

I noticed that the edge of the fluid patch got zero velocity:

(this plot corresponds to bottom right corner in the schematic of the last post)

This was caused by a slight inaccuracy in the definition of the fluid boundary mesh, which meant that the nearest-neighbor mapping interpolated values from the wall to the fluid-gas boundary. Since i did not set the velocity values on the wall, it defaulted to zero. This abrupt shift to zero velocity caused the turbulence model to crash.

Am i using preCICE in the correct way? It seems dangerous to me that values from the different OpenFOAM patches are interpolated to each other.

Is it possible to have several meshes in each solver? If so, would having a separate mesh for the wall and bottom have fixed my problem?

Eirik