基于Lattice Boltzmann方法的C++ CFD仿真效果异常求助
2026-8-7
C++实现LBM仿真效果不佳求助
我尝试用C++实现自主CFD程序，参考YouTube上的Lattice Boltzmann方法视频，但仿真效果达不到视频中Python实现的LBM效果。使用SDL2在CPU端做可视化，不追求速度，只想要美观的仿真动画效果。
代码实现
Cell类代码
//cell class class cell { public: double Fi[nL] = {0,0,0,0,0,0,0,0,0}; double density = 0; double momentumX = 0; double momentumY = 0; double velocityX = 0; double velocityY = 0; double Fieq[nL] = {0,0,0,0,0,0,0,0,0}; //obstacle bool obstacle = false; void densityOperator() { for (int i = 0; i < nL; i++) { density += Fi[i]; } } void momentumOperator() { for (int i = 0; i < nL; i++) { momentumX += Fi[i] * cX[i]; momentumY += Fi[i] * cY[i]; } } void velocityOperator() { for (int i = 0; i < nL; i++) { if (density == 0) { density += 0.001; } velocityX += momentumX / density; // prolly very slow velocityY += momentumY / density; //velocityX += cX[i]; //velocityY += cY[i]; } } void FieqOperator() { for (int i = 0; i < nL; i++) { Fieq[i] = weights[i] * density * ( 1 + (cX[i] * velocityX + cY[i] * velocityY) / Cs + pow((cX[i] * velocityX + cY[i] * velocityY), 2) / (2 * pow(Cs, 4)) - (velocityX * velocityX + velocityY * velocityY) / (2 * pow(Cs, 2)) ); } } void FiOperator() { for (int i = 0; i < nL; i++) { Fi[i] = Fi[i] - (timestep / tau) * (Fi[i] - Fieq[i]); } } void addRightVelocity() { Fi[0] = 1.f; Fi[1] = 1.f; Fi[2] = 1.f; Fi[3] = 6.f; Fi[4] = 1.f; Fi[5] = 1.f; Fi[6] = 1.f; Fi[7] = 1.f; Fi[8] = 1.f; } };
索引转换函数
使用vector存储cell而非二维数组，通过以下函数将x,y坐标转换为一维索引：
int index(int x, int y) { return x * nY + y; }
全局变量定义
//box const int nX = 400; const int nY = 100; //viscosity float tau = 0.5; // 0.53 //time delta time per iteration float timestep = 1; //distance between cells float dist = 1000; //Speed of sound float Cs = 1 / sqrt(3) * (dist / timestep); //viscociti float v = pow(Cs, 2) * (tau - timestep / 2); // tau will need to be much smaller //time steps int nT = 3000; //lattice speeds and weights const int nL = 9; //Ci vector direction, discrete velocity int cX[9] = { 0, 0, 1, 1, 1, 0, -1, -1, -1 }; int cY[9] = { 0, 1, 1, 0, -1, -1, -1, 0 , 1 }; //weights, based on navier stokes float weights[9] = { 4 / 9, 1 / 9, 1 / 36, 1 / 9, 1 / 36, 1 / 9, 1 / 36, 1 / 4, 1 / 36 }; //opposite populations int cO[9] = { 0, 5, 6, 7, 8, 1, 2, 3, 4 };
主函数代码
int main() { //init vector cells for (int x = 0; x < nX; x++) { for (int y = 0; y < nY; y++) { cell cellUnit; cells.push_back(cellUnit); TempCells.push_back(cellUnit); } } //SDL //SDL //------------------------------------------------------------- SDL_Window* window = nullptr; SDL_Renderer* renderer = nullptr; SDL_Init(SDL_INIT_VIDEO); SDL_CreateWindowAndRenderer(nX* 3, nY * 3, 0, &window, &renderer); SDL_RenderSetScale(renderer, 3, 3); SDL_SetRenderDrawColor(renderer, 0, 0, 0, 255); SDL_RenderClear(renderer); //-------------------------------------------------------------// //Circle Object Gen for (int x = 0; x < nX; x++) { for (int y = 0; y < nY; y++) { //cicle position int circleX = 5; int circleY = 50; //circle radius float radius = 10; //distance bewtween cell and circle pos float distance = sqrt(pow(circleX - x, 2) + pow(circleY - y, 2)); if (distance < radius) { cells[index(x,y)].obstacle = true; } else { cells[index(x, y)].obstacle = false; } } } //add velocity for (int x = 0; x < nX; x++) { for (int y = 0; y < nY; y++) { cells[index(x, y)].addRightVelocity(); //random velocity for (int i = 0; i < nL; i++) { cells[index(x,y)].Fi[i] += (rand() % 200) / 100; } } } for (int t = 0; t < nT; t++) { //SDL //-------------------------------------------------------------- //clear renderer if (t % 20 == 0) { SDL_SetRenderDrawColor(renderer, 255, 255, 255, 255); SDL_RenderClear(renderer); } //-------------------------------------------------------------- //streaming: //because we will loop over the same populations we do not want to switch the same population twice for (int x = 0; x < nX; x++) { for (int y = 0; y < nY; y++) { if (x == 0) { cells[index(x, y)].Fi[3] += 0.4; } //for populations for (int i = 0; i < nL; i++) { //boundary //checs if cell is object or air if (cells[index(x, y)].obstacle == false) { //air //targetet cell int cellX = x + cX[i]; int cellY = y + cY[i]; //out of bounds check + rearange to other side if (cellX < 0) { //left to right cellX = nX; } if (cellX >= nX) { //right to left cellX = 0; continue; } if (cellY < 0) { //top to buttom cellY = nY; } if (cellY >= nY) { //bottom to top cellY = 0; } //if neighborinig cell is object --> collision with object if (cells[index(cellX, cellY)].obstacle == true) { //Boundary handling https://youtu.be/jfk4feD7rFQ?t=2821 TempCells[index(x,y)].Fi[cO[i]] = cells[index(x, y)].Fi[i]; } //if not then stream to neighbor air cell with oposite population TempCells[index(cellX, cellY)].Fi[cO[i]] = cells[index(x, y)].Fi[i]; } else { //wall //SDL GRAPICHS if (t % 20 == 0) { SDL_SetRenderDrawColor(renderer, 0, 0, 0, 255); SDL_RenderDrawPoint(renderer, x, y); } } } } } for (int x = 0; x < nX; x++) { for (int y = 0; y < nY; y++) { for (int i = 0; i < nL; i++) { cells[index(x, y)].Fi[i] = TempCells[index(x, y)].Fi[cO[i]]; } } } //collision: for (int x = 0; x < nX; x++) { for (int y = 0; y < nY; y++) { //density: cells[index(x, y)].densityOperator(); //momentum: cells[index(x, y)].momentumOperator(); //velocity: cells[index(x, y)].velocityOperator(); //Fieq + new new Fi: for (int i = 0; i < nL; i++) { cells[index(x, y)].FieqOperator(); } //SDL Graphics if (t % 20 == 0) { if (cells[index(x, y)].obstacle == false) { SDL_SetRenderDrawColor(renderer, cells[index(x, y)].density, cells[index(x, y)].density , 255 , 255); SDL_RenderDrawPoint(renderer, x, y); } } } } for (int x = 0; x < nX; x++) { for (int y = 0; y < nY; y++) { cells[index(x, y)].FiOperator(); } } //SDL Graphics if (t % 20 == 0 ) { SDL_RenderPresent(renderer); } } return 0; }
当前仿真效果
- 密度可视化（粉色表示高密度，超出255范围循环显示）
- 密度可视化
- 密度可视化
求助需求
代码较为杂乱且未优化，希望有C++实现LBM经验的开发者提供建议。仿真能运行，但无法达到目标视频中的美观动画效果。
编辑说明：已重置density、velocityX/Y、momentumX/Y等参数。
内容的提问来源于stack exchange，提问作者Christian
相关产品推荐
相关产品推荐
