LBM溃坝算例代码

/* This file is part of the Palabos library.
 *
 * Copyright （C） 2011-2017 FlowKit Sarl
 * Route d‘Oron 2
 * 1010 Lausanne Switzerland
 * E-mail contact: contact@flowkit.com
 *
 * The most recent release of Palabos can be downloaded at 
 * 
 *
 * The library Palabos is free software: you can redistribute it and/or
 * modify it under the terms of the GNU Affero General Public License as
 * published by the Free Software Foundation either version 3 of the
 * License or （at your option） any later version.
 *
 * The library 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 Affero General Public License for more details.
 *
 * You should have received a copy of the GNU Affero General Public License
 * along with this program.  If not see .
*/

/* The breaking dam free surface problem. This code demonstrates the basic usage of the
 * free surface module in Palabos. Surface tension and contact angles are optional. 
 */

#include “palabos3D.h“
#include “palabos3D.hh“

using namespace plb;

#define DEscriptOR descriptors::ForcedD3Q19Descriptor

typedef double T;


// Smagorinsky constant for LES model.
const T cSmago = 0.14;

// Physical dimensions of the system （in meters）.
const T lx = 3.22;
const T ly = 1.0;
const T lz = 1.0;

const T rhoEmpty = T（1）;
    
plint writeImagesIter   = 100;
plint getStatisticsIter = 20;

plint maxIter;
plint N;
plint nx ny nz;
T delta_t delta_x;
Array externalForce;
T nuPhys nuLB tau omega Bo surfaceTensionLB contactAngle;

std::string outDir;
plint obstacleCenterXYplane obstacleLength obstacleWidth obstacleHeight beginWaterReservoir waterReservoirHeight;
plint waterLevelOne waterLevelTwo waterLevelThree waterLevelFour;

void setupParameters（） {
    delta_x = lz / N;
    nx = util::roundToInt（lx / delta_x）;
    ny = util::roundToInt（ly / delta_x）;
    nz = util::roundToInt（lz / delta_x）;

    // Gravity in lattice units.
    T gLB = 9.8 * delta_t * delta_t/delta_x;
    externalForce = Array（0. 0. -gLB）;
    tau            = （nuPhys*DEscriptOR::invCs2*delta_t）/（delta_x*delta_x） + 0.5;
    omega          = 1./tau;    
    nuLB           = （tau-0.5）*DEscriptOR::cs2; // Viscosity in lattice units.
    
    surfaceTensionLB = rhoEmpty * gLB * N * N / Bo;

    obstacleCenterXYplane = util::roundToInt（0.744*N）;
    obstacleLength        = util::roundToInt（0.403*N）;
    obstacleWidth         = util::roundToInt（0.161*N）;
    obstacleHeight        = util::roundToInt（0.161*N）;
    beginWaterReservoir   = util::roundToInt（（0.744+1.248）*N）;
    waterReservoirHeight  = util::roundToInt（0.55*N）;
    
    waterLevelOne   = util::roundToInt（0.496*N）;
    waterLevelTwo   = util::roundToInt（2.*0.496*N）;
    waterLevelThree = util::roundToInt（3.*0.496*N）;
    waterLevelFour  = util::roundToInt（（3.*0.496 + 1.150）*N