Commit 757294b9 authored by Christoph Rettinger's avatar Christoph Rettinger

Merge branch 'mr_new_coupling' into 'master'

New lbm-mesa_pd-coupling

See merge request !240
parents 48bd1980 595ebbb2
Pipeline #21442 passed with stages
in 520 minutes and 2 seconds
......@@ -4,6 +4,7 @@ add_subdirectory( DEM )
add_subdirectory( FieldCommunication )
add_subdirectory( MeshDistance )
add_subdirectory( CouetteFlow )
add_subdirectory( FluidParticleCoupling )
add_subdirectory( ForcesOnSphereNearPlaneInShearFlow )
add_subdirectory( GranularGas )
add_subdirectory( LennardJones )
......
waLBerla_link_files_to_builddir( "*.dat" )
waLBerla_python_file_generates(GeneratedLBM.py
GeneratedLBM.cpp GeneratedLBM.h
)
waLBerla_python_file_generates(GeneratedLBMWithForce.py
GeneratedLBMWithForce.cpp GeneratedLBMWithForce.h
)
waLBerla_add_executable ( NAME SphereWallCollision FILES SphereWallCollision.cpp GeneratedLBM.py
DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling mesa_pd postprocessing timeloop vtk )
waLBerla_add_executable ( NAME SettlingSphereTenCate FILES SettlingSphereInBox.cpp GeneratedLBM.py
DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling mesa_pd postprocessing timeloop vtk )
waLBerla_add_executable ( NAME SphereMovingWithPrescribedVelocity FILES SphereMovingWithPrescribedVelocity.cpp GeneratedLBM.py
DEPENDS blockforest boundary core domain_decomposition field lbm mesa_pd lbm_mesapd_coupling postprocessing timeloop vtk )
waLBerla_add_executable ( NAME LubricationForceEvaluation FILES LubricationForceEvaluation.cpp GeneratedLBM.py
DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling mesa_pd postprocessing timeloop vtk )
waLBerla_add_executable ( NAME DragForceSphere FILES DragForceSphere.cpp GeneratedLBMWithForce.py
DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling mesa_pd postprocessing timeloop vtk )
waLBerla_add_executable ( NAME ForcesOnSphereNearPlane FILES ForcesOnSphereNearPlane.cpp GeneratedLBM.py
DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling mesa_pd postprocessing timeloop vtk )
waLBerla_add_executable ( NAME ObliqueDryCollision FILES ObliqueDryCollision.cpp
DEPENDS blockforest core mesa_pd postprocessing )
waLBerla_add_executable ( NAME ObliqueWetCollision FILES ObliqueWetCollision.cpp GeneratedLBM.py
DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling mesa_pd postprocessing timeloop vtk )
\ No newline at end of file
This diff is collapsed.
import sympy as sp
import pystencils as ps
from lbmpy.creationfunctions import create_lb_method_from_existing, create_lb_ast, create_lb_method
from lbmpy_walberla import generate_lattice_model
from pystencils_walberla import CodeGeneration
from lbmpy.creationfunctions import create_lb_collision_rule
from lbmpy.moments import MOMENT_SYMBOLS, is_even, get_order
from lbmpy.stencils import get_stencil
from lbmpy.maxwellian_equilibrium import get_moments_of_discrete_maxwellian_equilibrium
from lbmpy.forcemodels import *
from collections import OrderedDict
from lbmpy.methods import create_with_discrete_maxwellian_eq_moments
with CodeGeneration() as ctx:
# methods: TRTlike, KBC, SRT, cumulant
generatedMethod = "TRTlike"
print("Generating " + generatedMethod + " LBM method")
#
# x,y,z = MOMENT_SYMBOLS
#
# e_sq = x**2 + y**2 + z**2
#
# moments = [0] * 19
#
# moments[0] = sp.sympify(1)
# moments[1] = 19 * e_sq - 30
# moments[2] = (21 * e_sq**2 - 53 * e_sq + 24) / 2
#
# moments[3] = x
# moments[5] = y
# moments[7] = z
#
# moments[4] = (5*e_sq - 9) * x
# moments[6] = (5*e_sq - 9) * y
# moments[8] = (5*e_sq - 9) * z
#
# moments[9] = 3 * x**2 - e_sq
# moments[11] = y**2 - z**2
#
# moments[13] = x * y
# moments[14] = y * z
# moments[15] = x * z
#
# moments[10] = (3* e_sq - 5)*(3*x**2 - e_sq)
# moments[12] = (3* e_sq - 5)*(y**2 - z**2)
# moments[16] = (y**2 - z**2) * x
# moments[17] = (z**2 - x**2) * y
# moments[18] = (x**2 - y**2) * z
#
# m_eq = get_moments_of_discrete_maxwellian_equilibrium(stencil, moments, order=2, c_s_sq=sp.Rational(1,3))
#
# rr_dict = OrderedDict(zip(moments, omega))
#
# forcing = sp.symbols("forcing_:%d" % 3)
# forcing=(sp.symbols("fx"),0,0)
# forcemodel=Luo(forcing) #None
# method = create_with_discrete_maxwellian_eq_moments(stencil, rr_dict, compressible=False, force_model=forcemodel)
# at this stage, we have the MRT model from the LBM book!
#this is the schiller / walberla MRT
if generatedMethod == "TRTlike":
stencil = get_stencil("D3Q19", 'walberla')
omega = sp.symbols("omega_:%d" % len(stencil))
method = create_lb_method(stencil=stencil, method='mrt', maxwellian_moments=False)
def modification_func(moment, eq, rate):
omegaVisc = sp.Symbol("omega_visc")
omegaBulk = ps.fields("omega_bulk: [3D]", layout='fzyx')# sp.Symbol("omega_bulk")
omegaMagic = sp.Symbol("omega_magic")
if get_order(moment) <= 1:
return moment, eq, 0
elif rate == omega[1]:
return moment, eq, omegaBulk.center_vector
elif rate == omega[2]:
return moment, eq, omegaBulk.center_vector
elif is_even(moment):
return moment, eq, omegaVisc
else:
return moment, eq, omegaMagic
my_method = create_lb_method_from_existing(method, modification_func)
# optimizations
collision_rule = create_lb_collision_rule(lb_method=my_method,
optimization={"cse_global": True,
"cse_pdfs": False,
#"split": True
} )
generate_lattice_model(ctx, 'GeneratedLBM', collision_rule)
elif generatedMethod == "KBC":
methodName = 'trt-kbc-n4'
#omega_field = ps.fields("omegaField(1): [3D]", layout='fzyx') omega_output_field=omega_field.center,
collision_rule = create_lb_collision_rule(method=methodName,entropic=True, stencil=get_stencil("D3Q27", 'walberla'), compressible=True,
optimization={"cse_global": True,
"cse_pdfs": False,
"split": True})
generate_lattice_model(ctx, 'GeneratedLBM', collision_rule)
elif generatedMethod == "SRT":
methodName = 'srt'
collision_rule = create_lb_collision_rule(method=methodName,stencil=get_stencil("D3Q19", 'walberla'),
optimization={"cse_global": True,
"cse_pdfs": False,
"split": True})
generate_lattice_model(ctx, 'GeneratedLBM', collision_rule)
elif generatedMethod == "cumulant":
x,y,z = MOMENT_SYMBOLS
cumulants = [0] * 27
cumulants[0] = sp.sympify(1) #000
cumulants[1] = x #100
cumulants[2] = y #010
cumulants[3] = z #001
cumulants[4] = x*y #110
cumulants[5] = x*z #101
cumulants[6] = y*z #011
cumulants[7] = x**2 - y**2 #200 - 020
cumulants[8] = x**2 - z**2 #200 - 002
cumulants[9] = x**2 + y**2 + z**2 #200 + 020 + 002
cumulants[10] = x*y**2 + x*z**2 #120 + 102
cumulants[11] = x**2*y + y*z**2 #210 + 012
cumulants[12] = x**2*z + y**2*z #201 + 021
cumulants[13] = x*y**2 - x*z**2 #120 - 102
cumulants[14] = x**2*y - y*z**2 #210 - 012
cumulants[15] = x**2*z - y**2*z #201 - 021
cumulants[16] = x*y*z #111
cumulants[17] = x**2*y**2 - 2*x**2*z**2 + y**2*z**2 # 220- 2*202 +022
cumulants[18] = x**2*y**2 + x**2*z**2 - 2*y**2*z**2 # 220 + 202 - 2*002
cumulants[19] = x**2*y**2 + x**2*z**2 + y**2*z**2 # 220 + 202 + 022
cumulants[20] = x**2*y*z # 211
cumulants[21] = x*y**2*z # 121
cumulants[22] = x*y*z**2 # 112
cumulants[23] = x**2*y**2*z # 221
cumulants[24] = x**2*y*z**2 # 212
cumulants[25] = x*y**2*z**2 # 122
cumulants[26] = x**2*y**2*z**2 # 222
from lbmpy.moments import get_order
def get_relaxation_rate(cumulant, omega):
if get_order(cumulant) <= 1:
return 0
elif get_order(cumulant) == 2 and cumulant != x**2 + y**2 + z**2:
return omega
else:
return 1
stencil = get_stencil("D3Q27", 'walberla')
omega = sp.Symbol("omega")
rr_dict = OrderedDict((c, get_relaxation_rate(c, omega))
for c in cumulants)
from lbmpy.methods import create_with_continuous_maxwellian_eq_moments
my_method = create_with_continuous_maxwellian_eq_moments(stencil, rr_dict, cumulant=True, compressible=True)
collision_rule = create_lb_collision_rule(lb_method=my_method,
optimization={"cse_global": True,
"cse_pdfs": False})
generate_lattice_model(ctx, 'GeneratedLBM', collision_rule)
else:
print("Invalid generated method string! " + generatedMethod)
import sympy as sp
import pystencils as ps
from lbmpy.creationfunctions import create_lb_method_from_existing, create_lb_ast, create_lb_method
from lbmpy_walberla import generate_lattice_model
from pystencils_walberla import CodeGeneration
from lbmpy.creationfunctions import create_lb_collision_rule
from lbmpy.moments import MOMENT_SYMBOLS, is_even, get_order
from lbmpy.stencils import get_stencil
from lbmpy.maxwellian_equilibrium import get_moments_of_discrete_maxwellian_equilibrium
from lbmpy.forcemodels import *
from collections import OrderedDict
from lbmpy.methods import create_with_discrete_maxwellian_eq_moments
with CodeGeneration() as ctx:
forcing=(sp.symbols("fx"),0,0)
forcemodel=Luo(forcing) #None
# methods: TRTlike, KBC, SRT
generatedMethod = "TRTlike"
if generatedMethod == "TRTlike":
stencil = get_stencil("D3Q19", 'walberla')
omega = sp.symbols("omega_:%d" % len(stencil))
methodWithForce = create_lb_method(stencil=stencil, method='mrt', maxwellian_moments=False,force_model=forcemodel)
def modification_func(moment, eq, rate):
omegaVisc = sp.Symbol("omega_visc")
omegaBulk = ps.fields("omega_bulk: [3D]", layout='fzyx')# = sp.Symbol("omega_bulk")
#omegaBulk = sp.Symbol("omega_bulk")
omegaMagic = sp.Symbol("omega_magic")
if get_order(moment) <= 1:
return moment, eq, 0
elif rate == omega[1]:
return moment, eq, omegaBulk.center_vector
elif rate == omega[2]:
return moment, eq, omegaBulk.center_vector
elif is_even(moment):
return moment, eq, omegaVisc
else:
return moment, eq, omegaMagic
my_methodWithForce = create_lb_method_from_existing(methodWithForce, modification_func)
collision_rule = create_lb_collision_rule(lb_method=my_methodWithForce)
# ,
# optimization={"cse_global": True,
# "cse_pdfs": False,
# "split": True,
# "vectorization":{'instruction_set': 'avx',
# 'nontemporal': True,
# 'assume_inner_stride_one': True}})
generate_lattice_model(ctx, 'GeneratedLBMWithForce', collision_rule)
elif generatedMethod == "SRT":
collision_rule = create_lb_collision_rule(method='srt',stencil=get_stencil("D3Q19", 'walberla'), force_model=forcemodel,
optimization={"cse_global": True,
"cse_pdfs": False,
#"split": True
})
generate_lattice_model(ctx, 'GeneratedLBMWithForce', collision_rule)
elif generatedMethod == "TRT":
collision_rule = create_lb_collision_rule(method='trt',stencil=get_stencil("D3Q19", 'walberla'), force_model=forcemodel, maxwellian_moments=False)
generate_lattice_model(ctx, 'GeneratedLBMWithForce', collision_rule)
elif generatedMethod == "KBC":
collision_rule = create_lb_collision_rule(method='trt-kbc-n4',entropic=True, stencil=get_stencil("D3Q27", 'walberla'), compressible=True, force_model=forcemodel,
optimization={"cse_global": True,
"cse_pdfs": False,
#"split": True
})
generate_lattice_model(ctx, 'GeneratedLBMWithForce', collision_rule)
else:
print("Invalid generated method string! " + generatedMethod)
//======================================================================================================================
//
// This file is part of waLBerla. waLBerla 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 3 of
// the License, or (at your option) any later version.
//
// waLBerla 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 waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
//
//! \file ObliqueDryCollision.cpp
//! \author Christoph Rettinger <christoph.rettinger@fau.de>
//
//======================================================================================================================
#include "mesa_pd/collision_detection/AnalyticContactDetection.h"
#include "mesa_pd/common/ParticleFunctions.h"
#include "mesa_pd/data/ParticleAccessorWithShape.h"
#include "mesa_pd/data/ParticleStorage.h"
#include "mesa_pd/data/ShapeStorage.h"
#include "mesa_pd/kernel/DoubleCast.h"
#include "mesa_pd/kernel/ExplicitEulerWithShape.h"
#include "mesa_pd/kernel/VelocityVerletWithShape.h"
#include "mesa_pd/kernel/LinearSpringDashpot.h"
#include "mesa_pd/mpi/ReduceContactHistory.h"
#include "core/Environment.h"
#include "core/logging/Logging.h"
#include <iostream>
namespace walberla {
namespace mesa_pd {
/*
* Simulates oblique sphere-wall collision and checks rebound angle, i.e. the tangential part of the collision model.
*
*/
int main( int argc, char ** argv )
{
walberla::mpi::Environment env(argc, argv);
WALBERLA_UNUSED(env);
walberla::mpi::MPIManager::instance()->useWorldComm();
real_t uNin_SI = real_t(1.7); // m/s
real_t diameter_SI = real_t(0.00318); // m
//real_t density_SI = real_t(2500); // kg/m**3, not used
real_t uNin = real_t(0.02);
real_t diameter = real_t(20);
real_t radius = real_t(0.5) * diameter;
real_t density = real_c(2.5);
// these values have actually no influence here and are just computed for completeness
real_t dx_SI = diameter_SI / diameter;
real_t dt_SI = uNin / uNin_SI * dx_SI;
real_t impactAngle = real_t(0);
real_t dt = real_t(0.1); // = (1 / #sub steps)
real_t frictionCoeff_s = real_t(0.8); // no influence
real_t frictionCoeff_d = real_t(0.12); // paper: 0.125+-0.007
std::string filename = "TangentialCollision.txt";
real_t collisionTime = real_t(80);
real_t nu = real_t(0.22); //Poissons ratio
bool useVelocityVerlet = false;
real_t restitutionCoeff = real_t(0.83);
for( int i = 1; i < argc; ++i )
{
if( std::strcmp( argv[i], "--impactAngle" ) == 0 ) { impactAngle = real_c( std::atof( argv[++i] ) ); continue; }
if( std::strcmp( argv[i], "--dt" ) == 0 ) { dt = real_c( std::atof( argv[++i] ) ); continue; }
if( std::strcmp( argv[i], "--staticFriction" ) == 0 ) { frictionCoeff_s = real_c( std::atof( argv[++i] ) ); continue; }
if( std::strcmp( argv[i], "--dynamicFriction" ) == 0 ) { frictionCoeff_d = real_c( std::atof( argv[++i] ) ); continue; }
if( std::strcmp( argv[i], "--collisionTime" ) == 0 ) { collisionTime = real_c( std::atof( argv[++i] ) ); continue; }
if( std::strcmp( argv[i], "--nu" ) == 0 ) { nu = real_c( std::atof( argv[++i] ) ); continue; }
if( std::strcmp( argv[i], "--filename" ) == 0 ) { filename = argv[++i]; continue; }
if( std::strcmp( argv[i], "--useVV" ) == 0 ) { useVelocityVerlet = true; continue; }
if( std::strcmp( argv[i], "--coefficientOfRestitution" ) == 0 ) { restitutionCoeff = real_c( std::atof( argv[++i] ) ); continue; }
}
WALBERLA_LOG_INFO_ON_ROOT("******************************************************");
WALBERLA_LOG_INFO_ON_ROOT("** NEW SIMULATION **");
WALBERLA_LOG_INFO_ON_ROOT("******************************************************");
WALBERLA_LOG_INFO_ON_ROOT("dt_SI = " << dt_SI);
WALBERLA_LOG_INFO_ON_ROOT("dx_SI = " << dx_SI);
WALBERLA_LOG_INFO_ON_ROOT("impactAngle = " << impactAngle);
WALBERLA_LOG_INFO_ON_ROOT("dt = " << dt);
WALBERLA_LOG_INFO_ON_ROOT("frictionCoeff_s = " << frictionCoeff_s);
WALBERLA_LOG_INFO_ON_ROOT("frictionCoeff_d = " << frictionCoeff_d);
WALBERLA_LOG_INFO_ON_ROOT("collisionTime = " << collisionTime);
WALBERLA_LOG_INFO_ON_ROOT("nu = " << nu);
WALBERLA_LOG_INFO_ON_ROOT("integrator = " << (useVelocityVerlet ? "Velocity Verlet" : "Explicit Euler"));
WALBERLA_LOG_INFO_ON_ROOT("restitutionCoeff = " << restitutionCoeff);
//init data structures
auto ps = walberla::make_shared<data::ParticleStorage>(2);
auto ss = walberla::make_shared<data::ShapeStorage>();
using ParticleAccessor_T = mesa_pd::data::ParticleAccessorWithShape;
auto accessor = walberla::make_shared<ParticleAccessor_T >(ps, ss);
auto sphereShape = ss->create<data::Sphere>( radius );
ss->shapes[sphereShape]->updateMassAndInertia(density);
const real_t particleMass = real_t(1) / ss->shapes[sphereShape]->getInvMass();
const real_t Mij = particleMass; // * particleMass / ( real_t(2) * particleMass ); // Mij = M for sphere-wall collision
const real_t lnDryResCoeff = std::log(restitutionCoeff);
// normal material parameters
const real_t stiffnessN = Mij * ( math::pi * math::pi + lnDryResCoeff * lnDryResCoeff ) / (collisionTime * collisionTime); // Costa et al., Eq. 18
const real_t dampingN = - real_t(2) * Mij * lnDryResCoeff / collisionTime; // Costa et al., Eq. 18
WALBERLA_LOG_INFO_ON_ROOT("normal: stiffness = " << stiffnessN << ", damping = " << dampingN);
// from Thornton et al
const real_t kappa = real_t(2) * ( real_t(1) - nu ) / ( real_t(2) - nu ) ;
const real_t stiffnessT = kappa * stiffnessN;
const real_t dampingT = std::sqrt(kappa) * dampingN;
WALBERLA_LOG_INFO_ON_ROOT("tangential: kappa = " << kappa << ", stiffness T = " << stiffnessT << ", damping T = " << dampingT);
real_t uTin = uNin * impactAngle;
// create sphere
data::Particle&& p = *ps->create();
p.setPosition(Vec3(0,0,2*radius));
p.setLinearVelocity(Vec3(uTin, 0., -uNin));
p.setType(0);
// create plane
data::Particle&& p0 = *ps->create(true);
p0.setPosition(Vec3(0,0,0));
p0.setShapeID(ss->create<data::HalfSpace>(Vector3<real_t>(0,0,1)));
p0.setType(0);
data::particle_flags::set(p0.getFlagsRef(), data::particle_flags::INFINITE);
data::particle_flags::set(p0.getFlagsRef(), data::particle_flags::FIXED);
// velocity verlet
kernel::VelocityVerletWithShapePreForceUpdate vvPreForce( dt );
kernel::VelocityVerletWithShapePostForceUpdate vvPostForce( dt );
// explicit euler
kernel::ExplicitEulerWithShape explEuler( dt );
// collision response
collision_detection::AnalyticContactDetection acd;
kernel::DoubleCast double_cast;
kernel::LinearSpringDashpot dem(1);
mpi::ReduceContactHistory rch;
dem.setStiffnessN(0,0,stiffnessN);
dem.setStiffnessT(0,0,stiffnessT);
dem.setDampingN(0,0,dampingN);
dem.setDampingT(0,0,dampingT);
dem.setFrictionCoefficientStatic(0,0,frictionCoeff_s);
dem.setFrictionCoefficientDynamic(0,0,frictionCoeff_d);
WALBERLA_LOG_DEVEL("begin: vel = " << p.getLinearVelocity() << ", contact vel: " << getVelocityAtWFPoint(0,*accessor,p.getPosition() + Vec3(0,0,-radius)) );
uint_t steps = 0;
real_t maxPenetration = real_t(0);
do
{
if(useVelocityVerlet) vvPreForce(0,*accessor);
real_t penetration;
if (double_cast(0, 1, *accessor, acd, *accessor ))
{
penetration = acd.getPenetrationDepth();
maxPenetration = std::max( maxPenetration, std::abs(penetration));
dem(acd.getIdx1(), acd.getIdx2(), *accessor, acd.getContactPoint(), acd.getContactNormal(), acd.getPenetrationDepth(), dt);
auto force = accessor->getForce(0);
auto torque = accessor->getTorque(0);
WALBERLA_LOG_INFO(steps << ": penetration = " << penetration << " || vel = " << accessor->getLinearVelocity(0) << " || force = " << force << ", torque = " << torque);
}
rch(*ps);
if(useVelocityVerlet) vvPostForce(0,*accessor);
else explEuler(0, *accessor);
++steps;
} while (double_cast(0, 1, *accessor, acd, *accessor ) || p.getLinearVelocity()[2] < 0);
real_t uTout = p.getLinearVelocity()[0] - radius * p.getAngularVelocity()[1];
WALBERLA_LOG_DEVEL("end: linear vel = " << p.getLinearVelocity() << ", angular vel = " << p.getAngularVelocity());
real_t reboundAngle = uTout / uNin;
WALBERLA_LOG_INFO_ON_ROOT("gamma_in = " << impactAngle);
WALBERLA_LOG_INFO_ON_ROOT("gamma_out = " << reboundAngle);
WALBERLA_LOG_INFO_ON_ROOT("Thornton: sliding should occur if " << real_t(2) * impactAngle / ( frictionCoeff_d * ( real_t(1) + restitutionCoeff)) << " >= " << real_t(7) - real_t(1) / kappa );
WALBERLA_LOG_INFO_ON_ROOT("Max penetration = " << maxPenetration << " -> " << maxPenetration / radius * 100. << "% of radius");
std::ofstream file;
file.open( filename.c_str(), std::ios::out | std::ios::app );
file << impactAngle << " " << reboundAngle << "\n";
file.close();
return EXIT_SUCCESS;
}
} //namespace mesa_pd
} //namespace walberla
int main( int argc, char ** argv )
{
return walberla::mesa_pd::main(argc, argv);
}
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
Setup
{
case 3;
impactRatio 0;
variant Acceleration; // Gravity or Acceleration
// control
applyOutflowBCAtTop true;
// LBM
diameter 30;
magicNumber 0.1875;
bulkViscRateFactor 100.;
useOmegaBulkAdaption true;
// RPD
numRPDSubCycles 10;
useLubricationCorrection true;
lubricationCutOffDistanceNormal 0.666667;
lubricationCutOffDistanceTangentialTranslational 0.5;
lubricationCutOffDistanceTangentialRotational 0.5;
lubricationMinimalGapSizeNonDim 0.002;
useVelocityVerlet true;
collisionTime 120;
// coupling
averageForceTorqueOverTwoTimeSteps true;
conserveMomentum false;
reconstructorType Grad; // Eq EAN Grad Ext
// IO properties
baseFolder vtk_out_ObliqueWetCollision;
fileIO true;
vtkIOFreq 0;
}
Variant_Gravity
{
gravitationalAcceleration_SI 0.25;
tau 0.502;
useFullGravityInNormalDirection false;
domainSize <20,10,10>; // multiples of diameter
numberOfBlocksPerDirection <8,5,5>;
initialSphereHeight 7; // multiple of diameter
}
Variant_Acceleration
{
uIn 0.02;
StTarget 100;
useStTargetInNormalDirection true;
domainSize <12,12,24>; // multiples of diameter
numberOfBlocksPerDirection <5,5,8>;
initialSphereHeight 22.5; // multiple of diameter
applyArtificialGravityAfterAccelerating true;
applyUInNormalDirection true;
accelerationFactor 25;
}
Case1
{
material Steel;
fluid Glycerol78;
}
Case2
{
material Glass;
fluid Glycerol45;
}
Case3
{
material Glass;
fluid Glycerol33;
}
Case4
{
material Steel;
fluid Glycerol37;
}
Mat_Steel
{
densitySphere_SI 7780; // kg/m**3
diameter_SI 12.7e-3; // m
restitutionCoeff 0.97;
frictionCoeff 0.02;
poissonsRatio 0.27; //Joseph 2001
}
Mat_Glass
{
densitySphere_SI 2540; // kg/m**3
diameter_SI 12.7e-3; // m
restitutionCoeff 0.97;
frictionCoeff 0.15;
poissonsRatio 0.23; //Joseph 2001
}
// from https://www.internetchemie.info/chemie-lexikon/stoffe/g/glycerol.php
Fluid_Water
{
densityFluid_SI 998; // kg/m**3
dynamicViscosityFluid_SI 1e-3; // Pa s
}
Fluid_Glycerol33
{
densityFluid_SI 1081; // kg/m**3
dynamicViscosityFluid_SI 2.5e-3; // Pa s
}
Fluid_Glycerol37
{
densityFluid_SI 1091; // kg/m**3