Skip to content

17_Kelvin_Helmholtz: add VOF Kelvin-Helmholtz instability tutorial #4

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
florian-simvia wants to merge 1 commit into
base: main
Choose a base branch
from
Open

17_Kelvin_Helmholtz: add VOF Kelvin-Helmholtz instability tutorial #4

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 5 additions & 0 deletions 17_Kelvin_Helmholtz/CASE/DATA/run.cfg
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,5 @@
[job_defaults]

n_procs: 1
n_threads: 1

205 changes: 205 additions & 0 deletions 17_Kelvin_Helmholtz/CASE/DATA/setup.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,205 @@
<?xml version="1.0" encoding="utf-8"?><Code_Saturne_GUI case="CASE" solver_version="9.1" study="code_saturne" version="2.0">
<additional_scalars/>
<analysis_control>
<output>
<mesh id="-1" label="Fluid domain" type="cells">
<all_variables status="on"/>
<location>all[]</location>
<writer id="-1"/>
</mesh>
<mesh id="-2" label="Boundary" type="boundary_faces">
<all_variables status="on"/>
<location>all[]</location>
<writer id="-1"/>
</mesh>
<probe_format choice="CSV"/>
<probe_recording_frequency>1</probe_recording_frequency>
<probes_interpolation choice=""/>
<probes_snap choice=""/>
<writer id="-1" label="results">
<directory name="postprocessing"/>
<format name="ensight" options="separate_meshes"/>
<frequency period="time_value">0.1</frequency>
<output_at_end status="on"/>
<time_dependency choice="fixed_mesh"/>
</writer>
</output>
<profiles/>
<time_averages/>
<time_parameters>
<maximum_time>10</maximum_time>
<property name="courant_number" label="CourantNb"/>
<property name="fourier_number" label="FourierNb"/>
<time_passing>0</time_passing>
<time_step_ref>0.004</time_step_ref>
</time_parameters>
</analysis_control>
<boundary_conditions>
<boundary label="X0" name="1" nature="wall">X0</boundary>
<boundary label="X1" name="2" nature="wall">X1</boundary>
<boundary label="Y0" name="3" nature="symmetry">Y0</boundary>
<boundary label="Y1" name="4" nature="symmetry">Y1</boundary>
<boundary label="Z0" name="5" nature="symmetry">Z0</boundary>
<boundary label="Z1" name="6" nature="symmetry">Z1</boundary>
<symmetry label="Y0" field_id="none"/>
<symmetry label="Y1" field_id="none"/>
<symmetry label="Z0" field_id="none"/>
<symmetry label="Z1" field_id="none"/>
<wall label="X1" field_id="none">
<velocity_pressure choice="off"/>
</wall>
<wall label="X0" field_id="none">
<velocity_pressure choice="off"/>
</wall>
</boundary_conditions>
<calculation_management>
<start_restart>
<frozen_field status="off"/>
</start_restart>
</calculation_management>
<lagrangian model="off"/>
<numerical_parameters>
<velocity_pressure_algo choice="simplec"/>
</numerical_parameters>
<physical_properties>
<fluid_properties>
<material choice="user_material"/>
<method choice="user_properties"/>
<property name="density" choice="predefined_law" label="Density">
<initial_value>1.17862</initial_value>
<value_0>1</value_0>
<value_1>2</value_1>
</property>
<property name="dynamic_diffusion" choice="constant" label="DiffDyn">
<initial_value>0.01</initial_value>
<listing_printing status="off"/>
<postprocessing_recording status="off"/>
</property>
<property name="molecular_viscosity" choice="predefined_law" label="LamVisc">
<initial_value>1.83e-05</initial_value>
<value_0>0.0001</value_0>
<value_1>0.0001</value_1>
</property>
<property name="surface_tension" choice="constant" label="SurfTen">
<initial_value>0</initial_value>
<listing_printing status="off"/>
<postprocessing_recording status="off"/>
</property>
<reference_pressure>101325</reference_pressure>
<reference_temperature>293.15</reference_temperature>
</fluid_properties>
<gravity>
<gravity_x>0</gravity_x>
<gravity_y>0</gravity_y>
<gravity_z>0</gravity_z>
</gravity>
<notebook/>
<omega>
<omega_x>0</omega_x>
<omega_y>0</omega_y>
<omega_z>0</omega_z>
</omega>
<time_tables/>
</physical_properties>
<solution_domain>
<extrusion/>
<faces_cutting status="off"/>
<joining/>
<mesh_cartesian>
<x_direction ncells="128" min="0.0" max="1.0" prog="1.0" law="constant"/>
<y_direction ncells="64" min="0.25" max="0.75" prog="1.0" law="constant"/>
<z_direction ncells="1" min="0.0" max="1.0" prog="1.0" law="constant"/>
</mesh_cartesian>
<mesh_origin choice="mesh_cartesian"/>
<mesh_smoothing status="off"/>
<meshes_list/>
<periodicity>
<face_periodicity mode="translation" name="1">
<fraction>0.1</fraction>
<plane>25</plane>
<selector>X0 or X1</selector>
<translation>
<translation_x>1</translation_x>
<translation_y>0</translation_y>
<translation_z>0</translation_z>
</translation>
<verbosity>1</verbosity>
<visualization>1</visualization>
</face_periodicity>
</periodicity>
<thin_walls/>
<volumic_conditions>
<zone label="all_cells" id="1" initialization="on" physical_properties="on">all[]</zone>
</volumic_conditions>
</solution_domain>
<thermophysical_models>
<ale_method/>
<atmospheric_flows model="off">
<large_scale_meteo status="off"/>
<read_meteo_data status="off"/>
</atmospheric_flows>
<compressible_model model="off"/>
<conjugate_heat_transfer>
<external_coupling>
<syrthes_instances/>
</external_coupling>
</conjugate_heat_transfer>
<gas_combustion model="off">
<thermodynamical_pressure status="off"/>
</gas_combustion>
<hgn_model model="no_mass_transfer">
<initialization>
<formula zone_id="1"/>
</initialization>
<variable name="void_fraction" type="model" label="void_fraction">
<rhs_reconstruction>1</rhs_reconstruction>
</variable>
</hgn_model>
<immersed_boundaries/>
<internal_coupling>
<coupled_scalars/>
<solid_zones/>
</internal_coupling>
<porosities/>
<radiative_transfer model="off"/>
<reference_values>
<length/>
</reference_values>
<source_terms/>
<thermal_scalar model="off"/>
<turbomachinery model="off">
<joining/>
</turbomachinery>
<turbulence model="off">
<gravity_terms status="on"/>
<initialization zone_id="1" choice="reference_value"/>
<reference_velocity>1</reference_velocity>
<wall_function>0</wall_function>
</turbulence>
<velocity_pressure>
<initialization>
<formula zone_id="1">velocity[0] = 0.;
velocity[1] = 0.;
velocity[2] = 0.;</formula>
</initialization>
<property name="stress" label="Stress" support="boundary"/>
<property name="stress_normal" label="Stress, normal" support="boundary">
<postprocessing_recording status="off"/>
</property>
<property name="stress_tangential" label="Stress, tangential" support="boundary">
<postprocessing_recording status="off"/>
</property>
<property name="total_pressure" label="total_pressure"/>
<property name="yplus" label="Yplus" support="boundary"/>
<variable name="pressure" label="Pressure">
<rhs_reconstruction>2</rhs_reconstruction>
</variable>
<variable name="velocity" dimension="3" label="Velocity">
<rhs_reconstruction>1</rhs_reconstruction>
</variable>
</velocity_pressure>
</thermophysical_models>
<user_functions>
<calculator/>
</user_functions>
</Code_Saturne_GUI>
158 changes: 158 additions & 0 deletions 17_Kelvin_Helmholtz/CASE/SRC/cs_user_initialization.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,158 @@
/*============================================================================
* User initialization for the Kelvin-Helmholtz instability (VOF).
*
* A shear layer between two fluid phases with a sinusoidal perturbation
* of the interface triggers the Kelvin-Helmholtz rolling vortex instability.
*============================================================================*/

/* code_saturne version 9.1 */

/*
This file is part of code_saturne, a general-purpose CFD tool.

Copyright (C) 1998-2025 EDF S.A.

This program is free software; you can redistribute it and/or modify it under
the terms of the GNU General Public License as published by the Free Software
Foundation; either version 2 of the License, or (at your option) any later
version.

This program is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
details.

You should have received a copy of the GNU General Public License along with
this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/

/*----------------------------------------------------------------------------*/

#include "cs_defs.h"

/*----------------------------------------------------------------------------
* Standard C library headers
*----------------------------------------------------------------------------*/

#include <assert.h>
#include <math.h>

#if defined(HAVE_MPI)
#include <mpi.h>
#endif

/*----------------------------------------------------------------------------
* PLE library headers
*----------------------------------------------------------------------------*/

#include <ple_coupling.h>

/*----------------------------------------------------------------------------
* Local headers
*----------------------------------------------------------------------------*/

#include "cs_headers.h"

/*----------------------------------------------------------------------------*/

BEGIN_C_DECLS

/*----------------------------------------------------------------------------*/
/*!
* \file cs_user_initialization.cpp
*
* \brief Kelvin-Helmholtz instability initialization (VOF).
*/
/*----------------------------------------------------------------------------*/

/*============================================================================
* User function definitions
*============================================================================*/

/*----------------------------------------------------------------------------*/
/*!
* \brief Define initial conditions for the Kelvin-Helmholtz problem.
*
* Sets a hyperbolic tangent shear velocity profile and a smoothed
* sinusoidal VOF interface between two fluid phases.
*
* \param[in, out] domain pointer to a cs_domain_t structure
*/
/*----------------------------------------------------------------------------*/

void
cs_user_initialization(cs_domain_t *domain)
{
const cs_lnum_t n_cells = domain->mesh->n_cells;

/* Do not reinitialize on restart */
if (domain->time_step->nt_prev > 0)
return;

/* ---- Parameters (adapt if needed) ---- */

const cs_real_t Lx = 1.0; /* domain length in x */
const cs_real_t y0 = 0.5; /* mean interface position */
const cs_real_t eps = 0.03; /* interface perturbation amplitude */
const cs_real_t U0 = 0.5; /* shear velocity magnitude */
const cs_real_t delta = 0.01; /* shear layer thickness */
const cs_real_t thickness = 0.01; /* interface smoothing thickness */

/* ---- Get mesh coordinates (v9.x: array [n_cells][3]) ---- */

const cs_real_t (*cell_cen)[3] =
domain->mesh_quantities->cell_cen;

/* ---- Get fields ---- */

/* Volume fraction field (VOF) */
cs_field_t *vf = cs_field_by_name_try("void_fraction");
if (vf == nullptr)
vf = cs_field_by_name_try("volume_fraction");
if (vf == nullptr)
vf = cs_field_by_name_try("alpha");

if (vf == nullptr)
bft_error(__FILE__, __LINE__, 0,
"VOF fraction field not found (tried: void_fraction, volume_fraction, alpha).");

/* Velocity field */
cs_real_t *vel = CS_F_(vel)->val;

/* ---- Initialization loop ---- */

for (cs_lnum_t c_id = 0; c_id < n_cells; c_id++) {

const cs_real_t x = cell_cen[c_id][0];
const cs_real_t y = cell_cen[c_id][1];

/* 1) Wavy interface */
const cs_real_t y_int =
y0 + eps * sin(4.0 * M_PI * x / Lx);

/* 2) Smoothed VOF using tanh */
const cs_real_t s = (y - y_int) / thickness;

cs_real_t alpha =
0.5 * (1.0 + tanh(s));

/* Clamp safety */
if (alpha < 0.0) alpha = 0.0;
if (alpha > 1.0) alpha = 1.0;

vf->val[c_id] = alpha;

/* 3) Smooth shear velocity profile */
const cs_real_t r = (y - y0) / delta;

const cs_real_t u =
U0 * tanh(r);

vel[3*c_id + 0] = u; /* u */
vel[3*c_id + 1] = 0.0; /* v */
vel[3*c_id + 2] = 0.0; /* w */
}
}

END_C_DECLS
Loading