三维多孔介质数值模拟代码.zip_fluent udf_smell62s_三维udf_三维介质_多孔三维

// 3D_PorousMeida.cpp : Defines the entry point for the console application. //direction of growth D3Q26 //2014/1/12 copyright by luoluobian #include "stdafx.h" #include <iostream> #include <cmath> #include <cstdlib> #include <iomanip> #include <fstream> #include <sstream> #include <string> #include <ctime> using namespace std; //D3Q19 const int Q = 19; int e[Q][3] = { { 0, 0, 0 }, { 1, 0, 0 }, { -1, 0, 0 }, { 0, 1, 0 }, { 0, -1, 0 }, { 0, 0, 1 }, { 0, 0, -1 }, { 1, 1, 0 }, { -1, -1, 0 }, { 1, -1, 0 }, { -1, 1, 0 }, { 1, 0, 1 }, { -1, 0, -1 }, { 1, 0, -1 }, { -1, 0, 1 }, { 0, 1, 1 }, { 0, -1, -1 }, { 0, 1, -1 }, { 0, -1, 1 } }; double w[Q] = { 1. / 3, 1. / 18, 1. / 18, 1. / 18, 1. / 18, 1. / 18, 1. / 18, 1. / 36, 1. / 36, 1. / 36, 1. / 36, 1. / 36, 1. / 36, 1. / 36, 1. / 36, 1. / 36, 1. / 36, 1. / 36, 1. / 36 }; //定义物理尺寸 const int NX = 60; const int NY = 60; const int NZ = 60; const int M = NX*NY*NZ; //lattice boltzmann method double rho[NX + 1][NY + 1][NZ + 1], u[NX + 1][NY + 1][NZ + 1][3], u0[NX + 1][NX + 1][NZ + 1][3]; double velocity; double f[NX + 1][NY + 1][NZ + 1][Q], F[NX + 1][NY + 1][NZ + 1][Q], feqq[NX + 1][NY + 1][NZ + 1][Q], ftemp[NX + 1][NY + 1][NZ + 1][Q]; int i, j, k, q, ip, jp, kp, n; double c, Re, dx, dy, dz, LX, LY, LZ, dt, rho0, tau_f, niu, omega, error, temp; //QuartetStructureGenerationSet int arrgrid[NX + 1][NY + 1][NZ + 1], solidBoundaryNode[NX + 1][NY + 1][NZ + 1]; double rho_in = 1.0; double rho_out = 0.5; void initilization(); void QuartetStructureGenerationSet(); double feq(int q, double rho, double u[3]); void Evolution(); void outputdata(int m); void Error(); void boundary(); void Macroscopic(); int main() { initilization(); QuartetStructureGenerationSet(); for (n = 0;; n++) { Evolution(); Macroscopic(); boundary(); if (n % 100 == 0) { Error(); cout << "The " << n << " th computation result: The max relative error of uvw is: " << setiosflags(ios::scientific) << error << endl; if (n >= 500) { if (n % 500 == 0) outputdata(n); if (error<1.0e-8) break; } } } return 0; } void initilization() { dx = 1.0; dy = 1.0; dz = 1.0; dt = dx; c = dx / dt; LX = dx*NX; LY = dy*NY; LZ = dz*NZ; tau_f = 1.35; omega = 1.0 / tau_f; rho0 = 0.5; for (k = 0; k <= NZ; k++) { for (j = 0; j <= NY; j++) { for (i = 0; i <= NX; i++) { u[i][j][k][0] = 0.0; u[i][j][k][1] = 0.0; u[i][j][k][2] = 0.0; rho[i][j][k] = rho0; for (q = 0; q < Q; q++) { f[i][j][k][q] = feq(q, rho[i][j][k], u[i][j][k]); } } } } } void QuartetStructureGenerationSet() { double d1_6 = 0.02; double d7_18 = d1_6 / 2.; double d19_26 = d1_6 / 8.; double n = 0.8; double cdd = 0.001; double numtotal_need = (1 - n)*NX*NY*NZ; int numsoild = 0; int Tnumsoild; int soild[M + 1][3]; srand((unsigned)time(NULL)); //第一步，遍历所有网格，生成固相内核 for (i = 0; i <= NX; i++) { for (j = 0; j <= NY; j++) { for (k = 0; k <= NZ; k++) { arrgrid[i][j][k] = 0; soild[numsoild][0] = 0; soild[numsoild][1] = 0; soild[numsoild][2] = 0; if (((rand() % 1000) / 1000.0) < cdd) { arrgrid[i][j][k] = 1; soild[numsoild][0] = i; soild[numsoild][1] = j; soild[numsoild][2] = k; numsoild = numsoild + 1; } } } } Tnumsoild = numsoild; //第二步，从固相核向周围26个方向生长 while (Tnumsoild < numtotal_need) { for (int index_soild = 0; index_soild < Tnumsoild; index_soild++) { int index_i = soild[index_soild][0]; int index_j = soild[index_soild][1]; int index_k = soild[index_soild][2]; //1向右方向生长 if (index_j < NY - 1) { i = index_i; j = index_j + 1; k = index_k; if (arrgrid[i][j][k] == 0 && ((rand() % 1000) / 1000.0) < d1_6) { numsoild = numsoild + 1; arrgrid[i][j][k] = 1; soild[numsoild][0] = i; soild[numsoild][1] = j; soild[numsoild][2] = k; } } //2向后方向生长 if (index_i > 0) { i = index_i - 1; j = index_j; k = index_k; if (arrgrid[i][j][k] == 0 && ((rand() % 1000) / 1000.0) < d1_6) { numsoild = numsoild + 1; arrgrid[i][j][k] = 1; soild[numsoild][0] = i; soild[numsoild][1] = j; soild[numsoild][2] = k; } } //3向左方向生长 if (index_j > 0) { i = index_i; j = index_j - 1; k = index_k; if (arrgrid[i][j][k] == 0 && ((rand() % 1000) / 1000.0) < d1_6) { numsoild = numsoild + 1; arrgrid[i][j][k] = 1; soild[numsoild][0] = i; soild[numsoild][1] = j; soild[numsoild][2] = k; } } //4向前方向生长 if (index_i < NX - 1) { i = index_i + 1; j = index_j; k = index_k; if (arrgrid[i][j][k] == 0 && ((rand() % 1000) / 1000.0) < d1_6) { numsoild = numsoild + 1; arrgrid[i][j][k] = 1; soild[numsoild][0] = i; soild[numsoild][1] = j; soild[numsoild][2] = k; } } //5向上方向生长 if (index_k < NZ - 1) { i = index_i; j = index_j; k = index_k + 1; if (arrgrid[i][j][k] == 0 && ((rand() % 1000) / 1000.0) < d1_6) { numsoild = numsoild + 1; arrgrid[i][j][k] = 1; soild[numsoild][0] = i; soild[numsoild][1] = j; soild[numsoild][2] = k; } } //6向下方向生长 if (index_k > 0) { i = index_i; j = index_j; k = index_k - 1; if (arrgrid[i][j][k] == 0 && ((rand() % 1000) / 1000.0) < d1_6) { numsoild = numsoild + 1; arrgrid[i][j][k] = 1; soild[numsoild][0] = i; soild[numsoild][1] = j; soild[numsoild][2] = k; } } //7向水平方向右前生长 if (index_i < NX - 1 && index_j < NY - 1) { i = index_i + 1; j = index_j + 1; k = index_k; if (arrgrid[i][j][k] == 0 && ((rand() % 1000) / 1000.0) < d7_18) { numsoild = numsoild + 1; arrgrid[i][j][k] = 1; soild[numsoild][0] = i; soild[numsoild][1] = j; soild[numsoild][2] = k; } } //8向水平方向左前生长 if (index_i < NX - 1 && index_j > 0) { i = index_i + 1; j = index_j - 1; k = index_k; if (arrgrid[i][j][k] == 0 && ((rand() % 1000) / 1000.0) < d7_18) { numsoild = numsoild + 1; arrgrid[i][j][k] = 1; soild[numsoild][0] = i; soild[numsoild][1] = j; soild[numsoild][2] = k; } } //9向水平方向右后生长 if (index_i > 0 && index_j < NY - 1) { i = index_i - 1; j = index_j + 1; k = index_k; if (arrgrid[i][j][k] == 0 && ((rand() % 1000) / 1000.0) < d7_18) { numsoild = numsoild + 1; arrgrid[i][j][k] = 1; soild[numsoild][0] = i; soild[numsoild][1] = j; soild[numsoild][2] = k; } } //10向水平方向左后生长 if (index_i > 0 && index_j > 0) { i = index_i - 1; j = index_j - 1; k = index_k; if (arrgrid[i][j][k] == 0 && ((rand() % 1000) / 1000.0) < d7_18) { numsoild = numsoild + 1; arrgrid[i][j][k] = 1; soild[numsoild][0] = i; soild[numsoild][1] = j; soild[numsoild][2] = k; } } //11向右上方向生长 if (index_j < NY - 1 && index_k < NZ - 1) { i = index_i; j = index_j + 1; k = index_k + 1; if (arrgrid[i][j][k] == 0 && ((rand() % 1000) / 1000.0) < d7_18) { numsoild = numsoild + 1; arrgrid[i][j][k] = 1; soild[numsoild][0] = i; soild[numsoild][1] = j; soild[numsoild][2] = k; } } //12向右下方向生长 if (index_j <NY - 1 && index_k >0) { i = index_i; j = index_j + 1; k = index_k - 1; if (arrgrid[i][j][k] == 0 && ((rand() % 1000) / 1000.0) < d7_18) { numsoild = numsoild + 1; arrgrid[i][j][k] = 1; soild[numsoild][0] = i; soild[numsoild][1] = j; soild[numsoild][2] = k; } } //13向左上方向生长 if (index_j > 0 && index_k < NZ - 1) { i = index_i; j = index_j - 1; k = index_k + 1; if (arrgrid[i][j][k] == 0 && ((rand() % 1000) / 1000.0) < d7_18) { numsoild = numsoild + 1; arrgrid[i][j][k] = 1; soild[numsoild][0] = i; soild[numsoild][1] = j; soild[numsoild][2] = k; } } //14向左下方向生长 if (index_j > 0 && index_k > 0) { i = index_i; j = index_j - 1; k = index_k - 1; if (arrgrid[i][j][k] == 0 && ((rand() % 1000) / 1000.0) < d7_18) { nums