简介

对于多相流模拟，Palabos中也是实现了很多，比如shanchen模型做的瑞丽-泰勒不平衡和两相混合器，还有helee模型做的双液滴碰撞，怎么说呢，我都跑过，但是，由于我的体系是气体和水，密度比极大，为了能顺利毕业，最终我选择了palabos里的自由表面流动模型来进行模拟，是通过原始案例中的溃坝，改造而成的，忽略了气相，只考虑液滴在复杂多孔介质中的渗透时，每一时间步的界面位置。选择这个模型的主要原因，因为改起来简单，模型封装的很彻底，接触角不需要我再查文献设置，很多的边界条件都被封装为一个单一的枚举值，通过该枚举值进行调用，很适合堆数据。

代码

#include "palabos3D.h"
#include "palabos3D.hh"using namespace plb;#define DESCRIPTOR descriptors::ForcedD3Q19Descriptortypedef double T;// Smagorinsky constant for LES model.
// const T cSmago = 0.14;T lx;
T ly;
T lz;const T rhoEmpty = T(1);
const T rho = T(1000);plint writeImagesIter   = 400;
plint getStatisticsIter = 20;plint maxIter;
plint N;
plint nx, ny, nz;
T delta_t, delta_x;
Array<T,3> externalForce;
T nuPhys, nuLB, tau, omega, surfaceTensionPhys, surfaceTensionLB, contactAngle;std::string fNameIn,fNameOut;void readGeometry(std::string fNameIn, std::string fNameOut, MultiScalarField3D<int>& geometry)
{const plint nx = geometry.getNx();const plint ny = geometry.getNy();const plint nz = geometry.getNz();Box3D sliceBox(0,0, 0,ny-1, 0,nz-1);std::auto_ptr<MultiScalarField3D<int> > slice = generateMultiScalarField<int>(geometry, sliceBox);plb_ifstream geometryFile(fNameIn.c_str());for (plint iX=0; iX<nx; ++iX) {if (!geometryFile.is_open()) {pcout << "Error: could not open geometry file " << fNameIn << std::endl;exit(EXIT_FAILURE);}geometryFile >> *slice;copy(*slice, slice->getBoundingBox(), geometry, Box3D(iX,iX, 0,ny-1, 0,nz-1));}/*{VtkImageOutput3D<T> vtkOut("porousMediumOrigin", 1.0);vtkOut.writeData<float>(*copyConvert<int,T>(geometry, geometry.getBoundingBox()), "tag", 1.0);}*//*{std::auto_ptr<MultiScalarField3D<T> > floatTags = copyConvert<int,T>(geometry, geometry.getBoundingBox());std::vector<T> isoLevels;isoLevels.push_back(0.5);typedef TriangleSet<T>::Triangle Triangle;std::vector<Triangle> triangles;Box3D domain = floatTags->getBoundingBox().enlarge(-1);domain.x0++;domain.x1--;isoSurfaceMarchingCube(triangles, *floatTags, isoLevels, domain);TriangleSet<T> set(triangles);std::string outDir = fNameOut + "/";set.writeBinarySTL(outDir + "porousMedium.stl");}*/
}void setupParameters() {delta_x = 0.001 / N;lx = nx*delta_x;ly = ny*delta_x;lz = nz*delta_x;// Gravity in lattice units.T gLB = 9.8 * delta_t * delta_t/delta_x;externalForce  = Array<T,3>(0., 0., -gLB);tau            = (nuPhys*DESCRIPTOR<T>::invCs2*delta_t)/(delta_x*delta_x) + 0.5;omega          = 1./tau;    nuLB           = (tau-0.5)*DESCRIPTOR<T>::cs2; // Viscosity in lattice units.surfaceTensionLB = (rhoEmpty / rho) * (delta_t*delta_t) / (delta_x*delta_x*delta_x) * surfaceTensionPhys;
}bool insideFluid(T x, T y, T z)
{Array<T,3> pos(x, y, z);Array<T,3> center(nx/2.0, ny/2.0, 0);T r = norm(pos-center);if (r <= nz/25) {return true;}return false;
}int initialFluidFlags(plint iX, plint iY, plint iZ) {if (insideFluid(iX, iY, iZ)) {return freeSurfaceFlag::fluid;}else {return freeSurfaceFlag::empty;}
}void porousMediaSetup(FreeSurfaceFields3D<T,DESCRIPTOR>& fields,MultiScalarField3D<int>& geometry)
{//integrateProcessingFunctional(new ShortenBounceBack3D<T,DESCRIPTOR>, fields.lattice.getBoundingBox(), fields.freeSurfaceArgs, 0);// Set all outer-wall cells to "wall" (here, bulk-cells are also set to "wall", but it// doesn't matter, because they are overwritten on the next line).setToConstant(fields.flag, fields.flag.getBoundingBox(), (int)freeSurfaceFlag::wall);setToFunction(fields.flag, fields.flag.getBoundingBox().enlarge(-1), initialFluidFlags);setToConstant(fields.flag, geometry, 1,fields.flag.getBoundingBox().enlarge(-1), (int)freeSurfaceFlag::wall);/*setToConstant(fields.flag, Box3D(6*nx/13, 7*nx/13, 6*ny/13, 7*ny/13, 0, 0), (int)freeSurfaceFlag::fluid);setToConstant(fields.flag, Box3D(5*nx/11, 6*nx/11, 5*ny/11, 6*ny/11, 1, nx/11), (int)freeSurfaceFlag::fluid);*/setToConstant(fields.flag, Box3D(9*nx/19, 10*nx/19, 9*ny/19, 10*ny/19, 0, nz/38), (int)freeSurfaceFlag::fluid);//pcout << "5" << std::endl;fields.periodicityToggle(2,false);fields.periodicityToggle(1,false);fields.periodicityToggle(0,false);fields.defaultInitialize();/*VtkImageOutput3D<T> vtkOut("porousMedium", 1.0);vtkOut.writeData<float>(*copyConvert<int,T>(fields.flag, fields.flag.getBoundingBox()), "tag", 1.0);*/
}void writeResults(MultiBlockLattice3D<T,DESCRIPTOR>& lattice, MultiScalarField3D<T>& volumeFraction, plint iT)
{static const plint nx = lattice.getNx();static const plint ny = lattice.getNy();static const plint nz = lattice.getNz();const plint imSize = 600;Box3D slice(0, nx-1, ny/2, ny/2, 0, nz-1);ImageWriter<T> imageWriter("leeloo");imageWriter.writeScaledGif(createFileName("u", iT, 6),*computeVelocityNorm(lattice, slice),imSize, imSize); imageWriter.writeScaledGif(createFileName("rho", iT, 6),*computeDensity(lattice, slice),imSize, imSize);imageWriter.writeScaledGif(createFileName("volumeFraction", iT, 6), *extractSubDomain(volumeFraction, slice),imSize, imSize);// Use a marching-cube algorithm to reconstruct the free surface and write an STL file./*std::vector<T> isoLevels;isoLevels.push_back((T) 0.5);typedef TriangleSet<T>::Triangle Triangle;std::vector<Triangle> triangles;isoSurfaceMarchingCube(triangles, volumeFraction, isoLevels, volumeFraction.getBoundingBox());TriangleSet<T>(triangles).writeBinarySTL(createFileName(fNameOut+"/interface", iT, 6)+".stl");*/VtkImageOutput3D<T> vtkOut(createFileName("volumeFraction", iT, 6), 1.);vtkOut.writeData<float>(volumeFraction, "vf", 1.);
}void writeStatistics(FreeSurfaceFields3D<T,DESCRIPTOR>& fields) {pcout << " -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- " << std::endl;T averageMass = freeSurfaceAverageMass<T,DESCRIPTOR>(fields.freeSurfaceArgs, fields.lattice.getBoundingBox());pcout << "Average Mass: " << averageMass  << std::endl;T averageDensity = freeSurfaceAverageDensity<T,DESCRIPTOR>(fields.freeSurfaceArgs, fields.lattice.getBoundingBox());pcout << "Average Density: " << std::setprecision(12) << averageDensity  << std::endl;T averageVolumeFraction = freeSurfaceAverageVolumeFraction<T,DESCRIPTOR>(fields.freeSurfaceArgs, fields.lattice.getBoundingBox());pcout << "Average Volume-Fraction: " << std::setprecision(12) << averageVolumeFraction  << std::endl;pcout << " -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- " << std::endl;
}int main(int argc, char **argv)
{plbInit(&argc, &argv);if (global::argc() != 12) {pcout << "Error missing some input parameter\n";}try {global::argv(1).read(fNameOut);global::directories().setOutputDir(fNameOut+"/");global::argv(2).read(nuPhys);global::argv(3).read(surfaceTensionPhys);global::argv(4).read(contactAngle);global::argv(5).read(N);global::argv(6).read(delta_t);global::argv(7).read(maxIter);global::argv(8).read(nx);global::argv(9).read(ny);global::argv(10).read(nz);global::argv(11).read(fNameIn);}catch(PlbIOException& except) {pcout << except.what() << std::endl;pcout << "The parameters for this program are :\n";pcout << "1. Output directory name.\n";pcout << "2. kinematic viscosity in physical Units (m^2/s) .\n";pcout << "3. Surface tension in physical units (N/m).\n";pcout << "4. Contact angle (in degrees).\n";pcout << "5. number of lattice nodes for lz .\n";pcout << "6. delta_t .\n"; pcout << "7. maxIter .\n";pcout << "Reasonable parameters on a parallel machine are: " << (std::string)global::argv(0) << "tmp 1.01e-6 0.0728 150.0 1000 1.e-8 80000 501 501 501 RawData.dat\n";exit (EXIT_FAILURE);}setupParameters();pcout << "Reading the geometry file." << std::endl;MultiScalarField3D<int> geometry(nx,ny,nz);readGeometry(fNameIn, fNameOut, geometry);MultiScalarField3D<int> inputgeometry(nx,ny,nz/2+1);copy(geometry, Box3D(0,nx-1,0,ny-1,0,nz/2), inputgeometry, inputgeometry.getBoundingBox());pcout << "delta_t= " << delta_t << std::endl;pcout << "delta_x= " << delta_x << std::endl;pcout << "delta_t*delta_t/delta_x= " << delta_t*delta_t/delta_x << std::endl;pcout << "externalForce= " << externalForce[2] << std::endl;pcout << "surfaceTensionLB= " << surfaceTensionLB << std::endl;pcout << "relaxation time= " << tau << std::endl;pcout << "omega= " << omega << std::endl;pcout << "kinematic viscosity physical units = " << nuPhys << std::endl;pcout << "kinematic viscosity lattice units= " << nuLB << std::endl;global::timer("initialization").start();//pcout << "1" << std::endl;SparseBlockStructure3D blockStructure(createRegularDistribution3D(nx, ny, nz/2+1));//pcout << "2" << std::endl;/*Dynamics<T,DESCRIPTOR>* dynamics= new SmagorinskyBGKdynamics<T,DESCRIPTOR>(omega, cSmago);*/Dynamics<T,DESCRIPTOR>* dynamics = new IncBGKdynamics<T,DESCRIPTOR>(omega);FreeSurfaceFields3D<T,DESCRIPTOR> fields( blockStructure, dynamics->clone(), rhoEmpty,surfaceTensionLB, contactAngle, externalForce );porousMediaSetup(fields, inputgeometry);std::vector<MultiBlock3D*> freeSurfaceArgs = aggregateFreeSurfaceParams(fields.lattice, fields.rhoBar, fields.j, fields.mass, fields.volumeFraction,fields.flag, fields.normal, fields.helperLists, fields.curvature, fields.outsideDensity);Array<T,3> injectionVelocity((T)0.,(T)0.,(T)0.1);pcout << "Time spent for setting up lattices: "<< global::timer("initialization").stop() << std::endl;T lastIterationTime = T();//pcout << "6" << std::endl;for (plint iT = 0; iT <= maxIter; ++iT) {global::timer("iteration").restart();T sum_of_mass_matrix = T();T lost_mass = T();if (iT % getStatisticsIter==0) {pcout << std::endl;pcout << "ITERATION = " << iT << std::endl;pcout << "Time of last iteration is " << lastIterationTime << " seconds" << std::endl;writeStatistics(fields);sum_of_mass_matrix = fields.lattice.getInternalStatistics().getSum(0);pcout << "Sum of mass matrix: " << sum_of_mass_matrix << std::endl;lost_mass = fields.lattice.getInternalStatistics().getSum(1);pcout << "Lost mass: " << lost_mass << std::endl;pcout << "Total mass: " << sum_of_mass_matrix + lost_mass << std::endl;pcout << "Interface cells: " << fields.lattice.getInternalStatistics().getIntSum(0) << std::endl;}if (iT % writeImagesIter == 0) {global::timer("images").start();writeResults(fields.lattice, fields.volumeFraction, iT);pcout << "Total time spent for writing images: "<< global::timer("images").stop() << std::endl;}                           // This includes the collision-streaming cycle, plus all free-surface operations.fields.lattice.executeInternalProcessors();/*applyProcessingFunctional (new PouringLiquid3D<T,DESCRIPTOR>(dynamics->clone(), injectionVelocity), Box3D(4*nx/9, 5*nx/9, 4*ny/9, 5*ny/9, 0, 0), freeSurfaceArgs );*/applyProcessingFunctional ( new PouringLiquid3D<T,DESCRIPTOR>(dynamics->clone(), injectionVelocity), Box3D(9*nx/19, 10*nx/19, 9*ny/19, 10*ny/19, 0, 0), freeSurfaceArgs );fields.lattice.evaluateStatistics();fields.lattice.incrementTime();lastIterationTime = global::timer("iteration").stop();}
}

说明

怎么说呢，这应该算是一只超级缝合怪，渗透案例的复杂结构读入，对格子描述符的尝试，与溃坝同案例文件夹下的vof两相流模型中的部分设置，看模型源代码找到的模型已封装好的速度边界条件的接口，嗯，这应该是我palabos中最费心的一个代码了。当然了，虽然我拿着个毕业，但实际上我能写的东西已经远远超过了当初写这个代码时的水平，分享了其实也不心疼。给小朋友们留一个可以参考的对象其实蛮好的，哈哈哈！
如果看不懂不要紧，以后我也会逐个解析案例。
无图无真相

完

毕业了——课题代码开源（三）使用Palabos的自由表面流模型仿复杂多孔介质中的液滴渗透相关推荐

毕业了——课题代码开源（二）使用Palabos的单相流多孔渗流代码
简介这是我依据Palabos案例解析(一)permeability.cpp案例的原始案例修改的适合我自己的代码,其实没咋动,主要就是改一下最后速度的统计的标准和数据输出的一些小东西,其实和源代码几乎 ...
＜2021SC@SDUSC＞开源游戏引擎Overload代码分析三(OvWindowing结束):OvWindowing——Dialogs
2021SC@SDUSC Overload代码分析三:OvWindowing--Dialogs 前言 Dialogs 一.FileDialog FileDialog.h FileDialog.cpp ...
10余万行C代码开源之后，我被震惊了。。。
10余万行C代码开源之后,我被震惊了... 7月12日,涛思团队对外宣布将研发了两年多的产品TDengine开源,10多万行C代码,包括最核心的存储引擎和计算引擎都上传到了GitHub上.上周末7月1 ...
多位博士毕业去了985/211/三四流高校，后来怎么样了？
除了985和部分211高校,我国有着数量更多的三四流高校,博士毕业或博士后出站想到高校求职,应当如何选择? 本文转载知乎讨论,希望给即将进入高校工作的青年科研人员提供一些参考意见. @Tensor 1 ...
【通知】深度学习之人脸图像算法核心代码开源和勘误汇总
许多同学早就在问人脸图像的书什么时候代码开源,8月底我们如约而至,今天开源核心的代码以及对第一次印刷的勘误汇总. 代码开源本次跟之前的两本书一样,代码开源在我们官方的GitHub项目中,按照有实战部 ...
语音合成论文与韩国小哥“撞车”后续：英伟达“赶紧”把代码开源了
乾明编辑整理量子位出品 | 公众号 QbitAI 前两天,量子位报道了韩国小哥语音合成论文与英伟达撞车一事. 在得知自己的论文与英伟达的论文"撞车"之后,韩国小哥赶紧在arX ...
iPhone在线音乐盒，代码开源
iPhone在线音乐盒,代码开源时间: 2010-02-03 09:50 来源: http://www.cocoachina.com/bbs/ 点击: 124次这是一个初中学生用一个月时间自学Ob ...
电赛校赛总结----一维板球系统【代码开源】
2022/4/21 搭建了整体的机械结构,最后因为经费问题,选择了用去年风力摆的架子去搭摄像头[openmv],看当年的国赛题,选择的是ov7670,但我们讨论后觉得还是openmv的识别比较好,, ...
基于STM32与OneNet平台的智能家居系统设计（代码开源含自制APP代码）
前言:本文为手把手教学的基础物联网开发设计,项目包含对下位机(MCU对外设数据读取与控制)和上位机(包含服务平台和APP端)的设计.下位机选取STM32作为MCU,外设有LED灯和DHT11温湿度传 ...

毕业了——课题代码开源（三）使用Palabos的自由表面流模型仿复杂多孔介质中的液滴渗透

简介

代码

说明

完

毕业了——课题代码开源（三）使用Palabos的自由表面流模型仿复杂多孔介质中的液滴渗透相关推荐

最新文章

热门文章