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()
