From 94d3e79a8617f88dc0219cfdeedfa3147833719d Mon Sep 17 00:00:00 2001 From: Adrian Kummerlaender Date: Mon, 24 Jun 2019 14:43:36 +0200 Subject: Initialize at openlb-1-3 --- src/dynamics/powerLawBGKdynamics.hh | 157 ++++++++++++++++++++++++++++++++++++ 1 file changed, 157 insertions(+) create mode 100644 src/dynamics/powerLawBGKdynamics.hh (limited to 'src/dynamics/powerLawBGKdynamics.hh') diff --git a/src/dynamics/powerLawBGKdynamics.hh b/src/dynamics/powerLawBGKdynamics.hh new file mode 100644 index 0000000..9c2f0c7 --- /dev/null +++ b/src/dynamics/powerLawBGKdynamics.hh @@ -0,0 +1,157 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2012, 2015 Mathias J. Krause, Vojtech Cvrcekt, Davide Dapelo + * E-mail contact: info@openlb.net + * The most recent release of OpenLB can be downloaded at + * + * + * 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. +*/ + +/** \file + * BGK Dynamics with adjusted omega -- generic implementation. + * Strain rate similar to "J.Boyd, J. Buick and S.Green: A second-order accurate lattice Boltzmann non-Newtonian flow model" + * Power Law similar to "Huidan Yu, Sharath S. Girimaji, Li-Shi Luo - DNS and LES of decaying isotropic turbulence with and without frame rotation using lattice Boltzmann method" + */ +#ifndef POWER_LAW_BGK_DYNAMICS_HH +#define POWER_LAW_BGK_DYNAMICS_HH + +#include "dynOmegaLatticeDescriptors.h" +#include "powerLawBGKdynamics.h" +#include "core/cell.h" +#include "core/util.h" +#include "lbHelpers.h" +#include "latticeDescriptors.h" +#include "math.h" + +namespace olb { + +////////////////////// Class PowerLawDynamics ////////////////////////// + +template +PowerLawDynamics::PowerLawDynamics (T m, T n , T nuMin, T nuMax) + : _m(m), + _n(n) +{ + _omegaMin = 2./(nuMax*2.*descriptors::invCs2() + 1.); + _omegaMax = 2./(nuMin*2.*descriptors::invCs2() + 1.); +} + +template +T PowerLawDynamics::computeOmegaPL( Cell& cell, T omega0, + T rho, T pi[util::TensorVal::n] ) +{ + T pre2 = pow(descriptors::invCs2()/2.* omega0/rho,2.); // strain rate tensor prefactor + T gamma = sqrt(2.*pre2*PiNeqNormSqr(cell)); // shear rate + + T nuNew = _m*pow(gamma,_n-1.); //nu for non-Newtonian fluid + T newOmega = 2./(nuNew*2.*descriptors::invCs2() + 1.); + + if (newOmega>_omegaMax) { + newOmega = _omegaMax; + } + if (newOmega<_omegaMin) { + newOmega = _omegaMin; + } + + return newOmega; +} + + +////////////////////// Class PowerLawBGKdynamics ////////////////////////// + +/** \param vs2_ speed of sound + * \param momenta_ a Momenta object to know how to compute velocity momenta + * \param momenta_ a Momenta object to know how to compute velocity momenta + */ +template +PowerLawBGKdynamics::PowerLawBGKdynamics ( + T omega, Momenta& momenta, T m, T n , T nuMin, T nuMax) + : BGKdynamics(omega, momenta), + PowerLawDynamics(m, n, nuMin, nuMax) +{ } + +template +void PowerLawBGKdynamics::collide ( + Cell& cell, + LatticeStatistics& statistics ) +{ + T rho, u[DESCRIPTOR::d], pi[util::TensorVal::n]; + this->_momenta.computeAllMomenta(cell, rho, u, pi); + + // Computation of the power-law omega. + // An external is used in place of BGKdynamics::_omega to keep generality and flexibility. + const auto oldOmega = cell.template getField(); + const auto newOmega = this->computeOmegaPL(cell, oldOmega, rho, pi); + + const T uSqr = lbHelpers::bgkCollision(cell, rho, u, newOmega); + cell.template setField(newOmega); + statistics.incrementStats(rho, uSqr); +} + +template +T PowerLawBGKdynamics::PiNeqNormSqr(Cell& cell ) +{ + return lbHelpers::computePiNeqNormSqr(cell); +} + + +////////////////////// Class ForcedPowerLawBGKdynamics ////////////////////////// + +/** \param vs2_ speed of sound + * \param momenta_ a Momenta object to know how to compute velocity momenta + * \param momenta_ a Momenta object to know how to compute velocity momenta + */ +template +PowerLawForcedBGKdynamics::PowerLawForcedBGKdynamics ( + T omega, Momenta& momenta, T m, T n , T nuMin, T nuMax) + : ForcedBGKdynamics(omega, momenta), + PowerLawDynamics(m, n, nuMin, nuMax) +{ } + +template +void PowerLawForcedBGKdynamics::collide ( + Cell& cell, + LatticeStatistics& statistics ) +{ + T rho, u[DESCRIPTOR::d], pi[util::TensorVal::n]; + this->_momenta.computeAllMomenta(cell, rho, u, pi); + + // Computation of the power-law omega. + // An external is used in place of BGKdynamics::_omega to keep generality and flexibility. + T oldOmega = cell.template getFieldPointer()[0]; + T newOmega = this->computeOmegaPL(cell, oldOmega, rho, pi); + + T* force = cell.template getFieldPointer(); + for (int iVel=0; iVel::bgkCollision(cell, rho, u, this->getOmega()); + lbHelpers::addExternalForce(cell, u, newOmega, rho); + cell.template getFieldPointer()[0] = newOmega; // updating omega + statistics.incrementStats(rho, uSqr); +} + +template +T PowerLawForcedBGKdynamics::PiNeqNormSqr(Cell& cell ) +{ + return lbHelpers::computeForcedPiNeqNormSqr(cell); +} + +} + +#endif -- cgit v1.2.3