从零开始的CAD|CAE开发: LBM源码实现分享

起因:

上期我们写了流体仿真的经典案例: 通过LBM,模拟计算涡流的形成,

当时承诺: 只要验证通过,就把代码开源出来;

ok.验证通过了,那么我也就将代码全都贴出来

代码开源并贴出:

 public class LidDrivenCavityFlow : IDisposable { public LidDrivenCavityFlow(int width = 200, int height = 100) {  Width = width;  Height = height;  Initialize(); } public int Width { get; set; } = 256; public int Height { get; set; } = 256; public double Error { get; set; } = 0; public int TimeStep { get; set; } = 0; public bool IsRunning { get; set; } = false; private double dx = 1.0;//X方向的步长 private double dy = 1.0;//Y方向的步长 private double Lx = 1.0;//X方向的长度; private double Ly = 1.0;//Y方向的长度; private double dt = 1.0;//时间步长; private double c = 1.0;//格子的速度 // 模拟参数 public double rho0 = 1.0; //密度 public double U = 0.1;//顶盖的速度; public double Re = 1000; //雷诺数 private double niu = 0; //粘力-运动粘度系数v private double tau_f = 1.0;//松弛时间-无量纲 public double Omega { get; set; } = 1.0;//松弛时间 private int[,] e; // 离散速度; // LBM 计算变量 private double[,,] f; // 演化前的分布函数 private double[,,] F; // 演化后的分布函数; private double[,] rho; // 密度场 private double[,] ux0; // n层时的速度X private double[,] uy0; // n层时的速度Y private double[,] ux; //n+1层 x方向速度 private double[,] uy; //n+1层 y方向速度 // D2Q9 模型常量 private const int Q = 9; //权重系数; private readonly double[] w = { 4.0/9, 1.0/9, 1.0/9, 1.0/9, 1.0/9, 1.0/36, 1.0/36, 1.0/36, 1.0/36 }; private readonly int[] cx = { 0, 1, 0, -1, 0, 1, -1, -1, 1 }; private readonly int[] cy = { 0, 0, 1, 0, -1, 1, 1, -1, -1 }; // 可视化相关 public VisualizationMode VisualizationMode { get; set; } = VisualizationMode.Streamline; private List[] streamlines; private int streamlineCounter = 0; public event Action VisualizationUpdated; public event Action StepCompleted; public void Initialize() {  dx = 1.0;  dy = 1.0;  Lx = dx*Width;  Ly = dy*Height;//长度;  dt = dx;  c = dx / dt; //格子的速度;  // 计算粘度  niu = U * Width / Re;  tau_f = 3.0 * niu + 0.5;// 松弛步长  Omega = 1.0 / tau_f;// 松弛时间-无量纲  e = new int[Q, 2]; //离散速度;  for (int i = 0; i < Q; i++)  {  e[i, 0] = cx[i];  e[i, 1] = cy[i];  }  // 初始化数组  f = new double[Width + 1, Height + 1, Q];  F = new double[Width + 1, Height + 1, Q];  rho = new double[Width + 1, Height + 1];  ux = new double[Width + 1, Height + 1];  uy = new double[Width + 1, Height + 1];  ux0 = new double[Width + 1, Height + 1];  uy0 = new double[Width + 1, Height + 1];  TimeStep = 0;  Error = 0;  // 初始化为静止流体  for (int i = 0; i <= Width; i++)  {  for (int j = 0; j <= Height; j++)  {  rho[i, j] = rho0;  ux[i, j] = 0;  uy[i, j] = 0;  ux0[i, j] = 0;  uy0[i, j] = 0;  ux[i, Height] = U; //顶盖速度为U;  for (int k = 0; k < Q; k++)  { f[i, j, k] = feq(k, rho[i, j], ux[i, j], uy[i, j]);  }  }  }  InitializeStreamlines();  UpdateVisualization(); } public void Step() {  if (!IsRunning) return;  // 碰撞迁移步骤,去除边界.  for (int i = 1; i < Width; i++)  {  for (int j = 1; j < Height; j++)  {  for (int k = 0; k < Q; k++)  { int ip = i - e[k, 0]; int jp = j - e[k, 1]; F[i, j, k] = f[ip, jp, k] + (feq(k, rho[ip, jp], ux[ip, jp], uy[ip, jp]) - f[ip, jp, k]) / tau_f;  }  }  }  // 计算宏观量  for (int i = 1; i < Width; i++)  {  for (int j = 1; j < Height; j++)  {  ux0[i, j] = ux[i, j]; //记录上一次结果  uy0[i, j] = ux[i, j];  rho[i, j] = 0;  ux[i, j] = 0.0;  uy[i, j] = 0.0;  for (int k = 0; k < Q; k++)  { f[i, j, k] = F[i, j, k]; rho[i, j] += f[i, j, k]; ux[i, j] += e[k, 0]* f[i, j, k]; uy[i, j] += e[k, 1]* f[i, j, k];  }  ux[i, j] /= rho[i, j];  uy[i, j] /= rho[i, j];  }  }  // ➤ 左右边界处理 （非平衡外推法）  for (int j = 1; j < Height; j++)  {  for (int k = 0; k < Q; k++)  {  var NX = Width;  // 右边界（i = Width - 1）  rho[NX, j] = rho[NX-1, j];  f[NX, j, k] = feq(k, rho[NX, j], ux[NX, j], uy[NX, j]) +  (f[NX-1, j, k] - feq(k, rho[NX-1, j], ux[NX-1, j], uy[NX-1, j]));  // 左边界（i = 0）  rho[0, j] = rho[1, j];  f[0, j, k] = feq(k, rho[0, j], ux[0, j], uy[0, j]) + (f[1, j, k] - feq(k, rho[1, j], ux[1, j], uy[1, j]));  }  }  // ➤ 上下边界处理  for (int i = 0; i <= Width; i++) //注意这里  {  for (int k = 0; k < Q; k++)  {  var NY = Height;  // 下边界（j = 0）  rho[i, 0] = rho[i, 1];  f[i, 0, k] = feq(k, rho[i, 0], ux[i, 0], uy[i, 0]) + (f[i, 1, k] - feq(k, rho[i, 1], ux[i, 1], uy[i, 1]));  // 上边界（j = Height - 1）  rho[i, NY] = rho[i, NY-1];  ux[i, NY] = U; // 设置移动壁速度  uy[i, NY] = 0.0; // 假设只在 x 方向移动  f[i, NY, k] = feq(k, rho[i, NY], ux[i, NY], uy[i, NY]) +  (f[i, NY-1, k] - feq(k, rho[i, NY-1], ux[i, NY-1], uy[i, NY-1]));  }  }  // 计算误差  if (TimeStep % 100 == 0)  {  CalculateError();  }  TimeStep++;  UpdateVisualization();  StepCompleted?.Invoke(TimeStep); } public void OutputData(int step, string filePath = null) {  // 如果没有指定路径，使用默认文件名  string fileName = filePath ?? $\"cavity_{step}.dat\";  using (StreamWriter outFile = new StreamWriter(fileName))  {  // 写入文件头  outFile.WriteLine(\"Title=\\\"LBM顶盖驱动流模拟\\\"\");  outFile.WriteLine(\"VARIABLES=\\\"X\\\",\\\"Y\\\",\\\"U\\\",\\\"V\\\"\");  outFile.WriteLine($\"ZONE T=\\\"Box\\\", I={Width + 1}, J={Height + 1}, F=POINT\");  // 写入数据（注意行列顺序与C++保持一致）  for (int j = 0; j <= Height; j++)  {  for (int i = 0; i <= Width; i++)  { // 使用InvariantCulture确保小数点格式一致 string line = string.Format(System.Globalization.CultureInfo.InvariantCulture, \"{0:F8} {1:F8} {2:F8} {3:F8}\", (double)i / Width, // X坐标归一化 (double)j / Height, // Y坐标归一化 ux[i, j], // X方向速度 uy[i, j]); // Y方向速度 outFile.WriteLine(line);  }  }  }  Console.WriteLine($\"结果已写入文件: {fileName}\"); } private void CalculateError() {  double temp1 = 0.0, temp2 = 0.0;  for (int i = 1; i < Width; i++)  {  for (int j = 1; j  RenderDensity(),  VisualizationMode.Velocity => RenderVelocity(),  VisualizationMode.Vorticity => RenderVorticity(),  VisualizationMode.Pressure => RenderPressure(), // 新增压力场渲染  VisualizationMode.Combined => RenderCombined(),  VisualizationMode.Streamline => RenderStreamline(),  _ => RenderCombined(), // 默认渲染密度  };  VisualizationUpdated?.Invoke(fluidImage); } // 辅助方法：获取双线性插值的速度值（提高流线平滑度） private (double ux, double uy) GetInterpolatedVelocity(double x, double y) {  // 确保坐标在有效范围内  x = Math.Clamp(x, 0, Width);  y = Math.Clamp(y, 0, Height);  // 找到当前点所在的网格单元（i,j）  int i = (int)Math.Floor(x);  int j = (int)Math.Floor(y);  // 确保索引在有效范围内  i = Math.Clamp(i, 0, Width - 2); // 确保i+1不会越界  j = Math.Clamp(j, 0, Height - 2); // 确保j+1不会越界  // 四个角点的坐标  double x0 = i;  double x1 = i + 1;  double y0 = j;  double y1 = j + 1;  // 双线性插值权重  double dx = x - x0;  double dy = y - y0;  // 确保权重在[0,1]范围内  dx = Math.Clamp(dx, 0, 1);  dy = Math.Clamp(dy, 0, 1);  // 获取四个角点的速度值  double ux00 = ux[i, j];  double ux10 = ux[i + 1, j];  double ux01 = ux[i, j + 1];  double ux11 = ux[i + 1, j + 1];  double uy00 = uy[i, j];  double uy10 = uy[i + 1, j];  double uy01 = uy[i, j + 1];  double uy11 = uy[i + 1, j + 1];  // 双线性插值计算ux  double ux_val =  (1 - dx) * (1 - dy) * ux00 +  dx * (1 - dy) * ux10 +  (1 - dx) * dy * ux01 +  dx * dy * ux11;  // 双线性插值计算uy  double uy_val =  (1 - dx) * (1 - dy) * uy00 +  dx * (1 - dy) * uy10 +  (1 - dx) * dy * uy01 +  dx * dy * uy11;  return (ux_val, uy_val); } private Bitmap RenderPressure() {  Bitmap bmp = new Bitmap(Width + 1, Height + 1);  BitmapData bmpData = bmp.LockBits(  new Rectangle(0, 0, bmp.Width, bmp.Height),  ImageLockMode.WriteOnly,  PixelFormat.Format32bppArgb);  // 计算压力范围  double minPressure = double.MaxValue;  double maxPressure = double.MinValue;  for (int i = 0; i < Width; i++)  {  for (int j = 0; j < Height; j++)  {  double pressure = rho[i, j] / 3.0; // 压力计算  if (pressure  maxPressure) maxPressure = pressure;  }  }  double range = Math.Max(maxPressure - minPressure, 0.001);  unsafe  {  byte* ptr = (byte*)bmpData.Scan0;  int bytesPerPixel = 4;  int stride = bmpData.Stride;  Parallel.For(0, Height, j =>  {  byte* row = ptr + (j * stride);  for (int i = 0; i < Width; i++)  { double pressure = rho[i, j] / 3.0; double normalized = (pressure - minPressure) / range; // 使用蓝-白-红色谱表示压力 Color color; if (normalized < 0.5) { // 蓝色到白色渐变 int blue = 255; int green = (int)(normalized * 2 * 255); int red = (int)(normalized * 2 * 255); color = Color.FromArgb(red, green, blue); } else { // 白色到红色渐变 int red = 255; int green = (int)((1.0 - (normalized - 0.5) * 2) * 255); int blue = (int)((1.0 - (normalized - 0.5) * 2) * 255); color = Color.FromArgb(red, green, blue); } row[i * bytesPerPixel] = color.B; // Blue row[i * bytesPerPixel + 1] = color.G; // Green row[i * bytesPerPixel + 2] = color.R; // Red row[i * bytesPerPixel + 3] = 255; // Alpha  }  });  }  bmp.UnlockBits(bmpData);  return bmp; } #region 绘制密度. /// 
 /// 绘制密度 /// 
 ///  private Bitmap RenderDensity() {  Bitmap fluidImage = new Bitmap(Width, Height);  // 使用LockBits提高性能  BitmapData bmpData = fluidImage.LockBits(  new Rectangle(0, 0, Width, Height),  ImageLockMode.WriteOnly,  PixelFormat.Format32bppArgb);  // 计算密度范围（动态适应）  double minDensity = double.MaxValue;  double maxDensity = double.MinValue;  for (int x = 0; x < Width; x++)  {  for (int y = 0; y < Height; y++)  {  if (rho[x, y]  maxDensity) maxDensity = rho[x, y];  }  }  double densityRange = maxDensity - minDensity;  if (densityRange < 0.0001) densityRange = 0.0001; // 避免除零  unsafe  {  byte* ptr = (byte*)bmpData.Scan0;  int bytesPerPixel = 4;  int stride = bmpData.Stride;  for (int y = 0; y < Height; y++)  {  byte* row = ptr + (y * stride);  for (int x = 0; x < Width; x++)  { // 动态映射密度到灰度值 double normalized = (rho[x, y] - minDensity) / densityRange; byte value = (byte)(normalized * 255); // BGRA格式 row[x * bytesPerPixel] = value; // Blue row[x * bytesPerPixel + 1] = value; // Green row[x * bytesPerPixel + 2] = value; // Red row[x * bytesPerPixel + 3] = 255; // Alpha  }  }  }  fluidImage.UnlockBits(bmpData);  return fluidImage; } /// 
 /// 绘制密度（增强对比度黑白灰版） /// 
 ///  private Bitmap RenderDensityGama() {  if (rho == null || Width <= 0 || Height <= 0)  return new Bitmap(1, 1); // 返回空图像以防无效状态  Bitmap fluidImage = new Bitmap(Width, Height);  // 使用LockBits提高性能  BitmapData bmpData = fluidImage.LockBits(  new Rectangle(0, 0, Width, Height),  ImageLockMode.WriteOnly,  PixelFormat.Format32bppArgb);  // 初始化修正 - 使用安全的初始值  double minDensity = double.MaxValue;  double maxDensity = double.MinValue;  // 首次遍历找出大致范围  for (int x = 0; x < Width; x++)  {  for (int y = 0; y < Height; y++)  {  double density = rho[x, y];  if (!double.IsFinite(density)) continue; // 跳过非数值  if (density  maxDensity) maxDensity = density;  }  }  // 范围有效性检查  const double defaultRange = 1.0;  if (minDensity > maxDensity)  {  minDensity = 0;  maxDensity = defaultRange;  }  // 计算裁剪边界  double range = Math.Max(maxDensity - minDensity, defaultRange * 0.01); // 最小范围保护  double lowerBound = minDensity + range * 0.05;  double upperBound = maxDensity - range * 0.05;  // 确保边界有效  if (lowerBound >= upperBound)  {  lowerBound = minDensity;  upperBound = maxDensity;  }  // 二次遍历计算实际有效范围  double effectiveMin = double.MaxValue;  double effectiveMax = double.MinValue;  int validPixelCount = 0; // 有效像素计数器  for (int x = 0; x < Width; x++)  {  for (int y = 0; y < Height; y++)  {  double density = rho[x, y];  // 跳过无效值和极端值  if (!double.IsFinite(density)) continue;  if (density  upperBound) continue;  validPixelCount++;  if (density  effectiveMax) effectiveMax = density;  }  }  // 当有效像素不足时的处理  if (validPixelCount < 10) // 至少需要10个有效点  {  effectiveMin = minDensity;  effectiveMax = maxDensity;  }  // 确保有效范围不为零  double densityRange = effectiveMax - effectiveMin;  if (densityRange < double.Epsilon)  densityRange = defaultRange;  unsafe  {  byte* ptr = (byte*)bmpData.Scan0;  int bytesPerPixel = 4;  int stride = bmpData.Stride;  // 对比度增强参数（带保护）  const double gamma = 1.8;  double contrastBoost = Math.Clamp(1.2, 0.1, 5.0);  for (int y = 0; y < Height; y++)  {  byte* row = ptr + (y * stride);  for (int x = 0; x < Width; x++)  { double rawValue = rho[x, y]; // 处理非数值和极端值 if (!double.IsFinite(rawValue)) { // 异常值显示为品红色 row[x * bytesPerPixel] = 255; // Blue row[x * bytesPerPixel + 1] = 0; // Green row[x * bytesPerPixel + 2] = 255; // Red row[x * bytesPerPixel + 3] = 255; // Alpha continue; } // 安全映射密度到[0,1]范围 double clamped = Math.Clamp(rawValue, effectiveMin, effectiveMax); double normalized = (clamped - effectiveMin) / densityRange; normalized = Math.Clamp(normalized, 0.0, 1.0); // 二次保护 // 伽马校正增强对比度（带安全保护） double gammaCorrected = Math.Pow(normalized, gamma); gammaCorrected = Math.Clamp(gammaCorrected, 0.0, 1.0); // S曲线增强（带安全计算） double contrastEnhanced = gammaCorrected; if (gammaCorrected < 0.5) { double factor = Math.Clamp(gammaCorrected * 2, 0.0, 1.0); contrastEnhanced = Math.Pow(factor, contrastBoost) * 0.5; } else { double factor = Math.Clamp((1 - gammaCorrected) * 2, 0.0, 1.0); contrastEnhanced = 1 - Math.Pow(factor, contrastBoost) * 0.5; } byte value = (byte)(Math.Clamp(contrastEnhanced, 0.0, 1.0) * 255); // BGRA格式 row[x * bytesPerPixel] = value; // Blue row[x * bytesPerPixel + 1] = value; // Green row[x * bytesPerPixel + 2] = value; // Red row[x * bytesPerPixel + 3] = 255; // Alpha  }  }  }  fluidImage.UnlockBits(bmpData);  return fluidImage; }  #endregion private Bitmap RenderVelocity() {  Bitmap fluidImage = new Bitmap(Width, Height);  BitmapData bmpData = fluidImage.LockBits(  new Rectangle(0, 0, Width, Height),  ImageLockMode.WriteOnly,  PixelFormat.Format32bppArgb);  // 计算当前帧最大速度  double maxSpeed = 0.001;  for (int x = 0; x < Width; x++)  {  for (int y = 0; y  maxSpeed) maxSpeed = speed;  }  }  unsafe  {  byte* ptr = (byte*)bmpData.Scan0;  int bytesPerPixel = 4;  int stride = bmpData.Stride;  for (int y = 0; y < Height; y++)  {  byte* row = ptr + (y * stride);  for (int x = 0; x < Width; x++)  { double speed = Math.Sqrt(ux[x, y] * ux[x, y] + uy[x, y] * uy[x, y]); double normalizedSpeed = speed / maxSpeed; // 使用HSV色相表示速度方向 double angle = Math.Atan2(uy[x, y], ux[x, y]); double hue = (angle + Math.PI) / (2 * Math.PI); // 饱和度表示速度大小（非线性增强） double saturation = Math.Pow(normalizedSpeed, 0.5); // 平方根增强低速度可见性 Color color = HsvToRgb(hue, saturation, 1.0); row[x * bytesPerPixel] = color.B; // Blue row[x * bytesPerPixel + 1] = color.G; // Green row[x * bytesPerPixel + 2] = color.R; // Red row[x * bytesPerPixel + 3] = 255; // Alpha  }  }  }  fluidImage.UnlockBits(bmpData);  return fluidImage; }  private Bitmap RenderCombined() {  int scale = 4;  Bitmap fluidImage = new Bitmap(Width * scale, Height * scale);  using (Graphics g = Graphics.FromImage(fluidImage))  {  g.SmoothingMode = SmoothingMode.AntiAlias;  g.Clear(Color.Black);  // 绘制密度背景（使用纹理刷提高性能） using (TextureBrush brush = new TextureBrush(densityBmp))  {  brush.ScaleTransform(scale, scale);  g.FillRectangle(brush, 0, 0, fluidImage.Width, fluidImage.Height);  }  // 自适应箭头参数  int arrowCount = Math.Max(10, Width / 15); // 箭头数量自适应  int step = Math.Max(3, Width / arrowCount);  double minArrowSpeed = 0.01 * U; // 基于顶盖速度  for (int x = step / 2; x < Width; x += step)  {  for (int y = step / 2; y  minArrowSpeed) { Point start = new Point(x * scale + scale / 2, y * scale + scale / 2); // 箭头长度自适应（速度越大箭头越长） double arrowLength = scale * 2 + speed * scale * 5; Point end = new Point( (int)(start.X + ux[x, y] * arrowLength), (int)(start.Y + uy[x, y] * arrowLength)); // 箭头颜色根据速度大小变化 // 计算速度比例值，限制在0-255范围内 int sv = Math.Clamp((int)(255 * speed / U), 0, 255); Color arrowColor = Color.FromArgb( 255, // Alpha (完全不透明) sv, // Red 2, // Green 2 // Blue ); DrawArrow(g, start, end, arrowColor, Math.Max(1, scale / 4)); }  }  }  }  return fluidImage; } /// 
 /// 漩涡 /// 
 ///  private Bitmap RenderVorticity()// {  Bitmap fluidImage = new Bitmap(Width, Height);  BitmapData bmpData = fluidImage.LockBits(  new Rectangle(0, 0, Width, Height),  ImageLockMode.WriteOnly,  PixelFormat.Format32bppArgb);  double[,] vorticity = new double[Width, Height];  double maxVorticity = 0.001;  // 计算涡量场（内部点）  for (int x = 1; x < Width - 1; x++)  {  for (int y = 1; y  maxVorticity) maxVorticity = Math.Abs(vorticity[x, y]);  }  }  // 边界处理（置零）  for (int x = 0; x < Width; x++)  {  vorticity[x, 0] = 0;  vorticity[x, Height - 1] = 0;  }  for (int y = 0; y < Height; y++)  {  vorticity[0, y] = 0;  vorticity[Width - 1, y] = 0;  }  unsafe  {  byte* ptr = (byte*)bmpData.Scan0;  int bytesPerPixel = 4;  int stride = bmpData.Stride;  for (int y = 0; y < Height; y++)  {  byte* row = ptr + (y * stride);  for (int x = 0; x  0) ? (int)(v * 255) : 0; int b = (v < 0) ? (int)(-v * 255) : 0; int g = (int)((1 - Math.Abs(v)) * 200); // 动态绿色分量 r = Math.Clamp(r, 0, 255); g = Math.Clamp(g, 0, 255); b = Math.Clamp(b, 0, 255); row[x * bytesPerPixel] = (byte)b; // Blue row[x * bytesPerPixel + 1] = (byte)g; // Green row[x * bytesPerPixel + 2] = (byte)r; // Red row[x * bytesPerPixel + 3] = 255; // Alpha  }  }  }  fluidImage.UnlockBits(bmpData);  return fluidImage; } #region 流线 private Bitmap RenderStreamline2() {  // 每4步重置一次流线  if (streamlineCounter % 4 == 0 || streamlines == null)  {  InitializeStreamlines();  }  streamlineCounter++;  Bitmap fluidImage = new Bitmap(Width, Height);  using (Graphics g = Graphics.FromImage(fluidImage))  {  // 使用纯黑色背景  g.Clear(Color.Black);  using (Bitmap densityBmp = RenderDensityColor())  {  ColorMatrix matrix = new ColorMatrix();  matrix.Matrix33 = .5f; // 极低透明度  ImageAttributes attributes = new ImageAttributes();  attributes.SetColorMatrix(matrix, ColorMatrixFlag.Default, ColorAdjustType.Bitmap);  g.DrawImage( densityBmp, new Rectangle(0, 0, Width, Height), 0, 0, Width, Height, GraphicsUnit.Pixel, attributes);  }  // 使用抗锯齿保证线条平滑  g.SmoothingMode = SmoothingMode.AntiAlias;  // 绘制所有流线  for (int i = 0; i  1)  { // 使用朴素的细线（宽度0.5-1.0） float lineWidth = 0.3f + (i % 3) * 0.25f; // 使用固定颜色（白色）或根据速度动态着色 //Color lineColor = Color.White; // 或者根据速度动态着色（可选） Color lineColor = GetDynamicLineColor(line); using (Pen pen = new Pen(lineColor, lineWidth)) { // 绘制平滑曲线 g.DrawLines(pen, line.ToArray()); }  }  }  }  return fluidImage; } private StreamlineRenderer mSteamHelper { get; set; } private Bitmap RenderStreamline() {  // 根据速度场数据完成流线渲染  if (mSteamHelper == null)  mSteamHelper = new StreamlineRenderer();  // ⚠️ 这里建议确认 Width 和 Height 是否已经 +1，否则可能越界或空白边缘  mSteamHelper.SetDataParams(Width + 1, Height + 1, ux, uy);  // 创建最终的渲染图像  Bitmap fluidImage = new Bitmap(Width, Height);  using (Graphics g = Graphics.FromImage(fluidImage))  {  // 使用黑色背景清除画布  g.Clear(Color.Black);  // 渲染密度图（带低透明度叠加）  using (Bitmap densityBmp = RenderDensityColor())  {  ColorMatrix matrix = new ColorMatrix  { Matrix33 = 0.5f // 设置图像整体 alpha 为 50%  };  using (ImageAttributes attributes = new ImageAttributes())  { attributes.SetColorMatrix(matrix, ColorMatrixFlag.Default, ColorAdjustType.Bitmap); g.DrawImage( densityBmp, new Rectangle(0, 0, Width, Height), 0, 0, Width, Height, GraphicsUnit.Pixel, attributes);  }  }  // 生成流线图像并绘制在画布上（叠加方式）  mSteamHelper.RenderToBitmap(g);  }  return fluidImage; } // 绘制发光点（带光晕效果） private void DrawGlowingPoint(Graphics g, PointF point, int size, Color color) {  RectangleF rect = new RectangleF(point.X - size/2, point.Y - size/2, size, size);  // 绘制光晕  using (SolidBrush glowBrush = new SolidBrush(Color.FromArgb(100, Color.White)))  {  g.FillEllipse(glowBrush, point.X - size, point.Y - size, size*2, size*2);  }  // 绘制主点  using (SolidBrush brush = new SolidBrush(color))  {  g.FillEllipse(brush, rect);  } } // 动态计算流线颜色（基于平均速度） private Color GetDynamicLineColor(List line) {  if (line.Count = 0 && x = 0 && y  0.8) return Color.FromArgb(255, 255, 50, 50); // 高速: 亮红色  if (normalizedSpeed > 0.5) return Color.FromArgb(255, 255, 150, 50); // 中高速: 橙色  if (normalizedSpeed > 0.3) return Color.FromArgb(255, 255, 255, 100); // 中速: 亮黄色  return Color.FromArgb(255, 100, 200, 255); // 低速: 亮蓝色 } private void InitializeStreamlines() {  // 减少流线数量（4-8条）  int lineCount = Math.Max(10, Math.Min(8, Width / 40));  streamlines = new List[lineCount];  streamlineCounter = 0; // 重置流线计数器  Random rand = new Random();  for (int i = 0; i < lineCount; i++)  {  streamlines[i] = new List();  float startX, startY;  // 在主要涡旋区域放置起点（避开固体边界）  do  {  if (i % 2 == 0)  { // 主涡旋区域（靠近顶盖中心） startX = (float)(Width * (0.3 + 0.4 * rand.NextDouble())); // X: 30%-70% startY = (float)(Height * (0.6 + 0.3 * rand.NextDouble())); // Y: 60%-90%（靠近顶盖）  }  else  { // 次涡旋区域（靠近底部中心） startX = (float)(Width * (0.3 + 0.4 * rand.NextDouble())); // X: 30%-70% startY = (float)(Height * (0.1 + 0.2 * rand.NextDouble())); // Y: 10%-30%（靠近底部）  }  } while (IsSolidBoundary((int)startX, (int)startY)); // 确保起点不在固体边界  // 添加起点到流线  streamlines[i].Add(new PointF(startX, startY));  // 沿速度场追踪后续点（最多500步，防止无限循环）  int maxSteps = 600;  int currentStep = 0;  double currentX = startX;  double currentY = startY;  while (currentStep < maxSteps)  {  // 获取当前网格点的速度（双线性插值提高精度）  (double uxVal, double uyVal) = GetInterpolatedVelocity(currentX, currentY);  // 计算移动步长（速度越大，步长越小，避免跳过细节）  float stepSize = (float)Math.Max(0.1f, 1.0f / (Math.Sqrt(uxVal * uxVal + uyVal * uyVal) + 1e-3f));  // 计算下一步坐标  double nextX = currentX + uxVal * stepSize;  double nextY = currentY + uyVal * stepSize;  // 检查是否超出计算域边界  if (nextX  Width || nextY  Height) break;  // 添加新点到流线  streamlines[i].Add(new PointF((float)nextX, (float)nextY));  // 更新当前坐标  currentX = nextX;  currentY = nextY;  currentStep++;  }  } } #endregion // 辅助方法：判断是否为固体边界（顶盖、左右、底边界） private bool IsSolidBoundary(int x, int y) {  // 顶盖驱动流中，边界是固体壁面  return x == 0 || x == Width || y == 0 || y == Height; } private bool IsBoundary(int x, int y) {  return x == 0 || x == Width || y == 0 || y == Height; } private Color HsvToRgb(double h, double s, double v) {  h = Math.Max(0, Math.Min(1, h - Math.Floor(h)));  s = Math.Max(0, Math.Min(1, s));  v = Math.Max(0, Math.Min(1, v));  int hi = (int)(Math.Floor(h * 6)) % 6;  double f = h * 6 - Math.Floor(h * 6);  double p = v * (1 - s);  double q = v * (1 - f * s);  double t = v * (1 - (1 - f) * s);  double r, g, b;  switch (hi)  {  case 0: r = v; g = t; b = p; break;  case 1: r = q; g = v; b = p; break;  case 2: r = p; g = v; b = t; break;  case 3: r = p; g = q; b = v; break;  case 4: r = t; g = p; b = v; break;  default: r = v; g = p; b = q; break;  }  return Color.FromArgb(  (int)(r * 255),  (int)(g * 255),  (int)(b * 255)); } // 更健壮的HSV转RGB函数 private static Color ColorFromHSV(float hue, float saturation, float value) {  // 输入验证  if (float.IsNaN(hue) || float.IsInfinity(hue)) hue = 0;  if (float.IsNaN(saturation) || float.IsInfinity(saturation)) saturation = 0;  if (float.IsNaN(value) || float.IsInfinity(value)) value = 0;  // 参数边界限制  hue = Math.Clamp(hue % 360, 0, 360);  if (hue < 0) hue += 360; // 处理负值  saturation = Math.Clamp(saturation, 0, 1);  value = Math.Clamp(value, 0, 1);  // 核心计算（保证中间值范围安全）  float c = Math.Clamp(value * saturation, 0, 1);  float x = Math.Clamp(c * (1 - Math.Abs((hue / 60) % 2 - 1)), 0, 1);  float m = Math.Clamp(value - c, 0, 1);  float r = 0, g = 0, b = 0;  // 角度范围处理  float sector = hue / 60;  if (sector < 1)  {  r = c; g = x; b = 0;  }  else if (sector < 2)  {  r = x; g = c; b = 0;  }  else if (sector < 3)  {  r = 0; g = c; b = x;  }  else if (sector < 4)  {  r = 0; g = x; b = c;  }  else if (sector < 5)  {  r = x; g = 0; b = c;  }  else // 5 <= sector < 6  {  r = c; g = 0; b = x;  }  // 安全转换到RGB字节值  int red = Math.Clamp((int)((r + m) * 255), 0, 255);  int green = Math.Clamp((int)((g + m) * 255), 0, 255);  int blue = Math.Clamp((int)((b + m) * 255), 0, 255);  // 最后防线：确保颜色有效  if (red < 0 || green < 0 || blue  255 || green > 255 || blue > 255)  {  // 记录错误但不中断流程  //Debug.WriteLine($\"Invalid HSV conversion: h={hue}, s={saturation}, v={value} => r={r}, g={g}, b={b}\");  return Color.Magenta; // 视觉标记错误  }  return Color.FromArgb(red, green, blue); } #endregion public void Dispose() {  f = null;  F = null;  rho = null;  ux = null;  uy = null;  ux0 = null;  uy0 = null;  GC.Collect(); } public void Start() => IsRunning = true; public void Pause() => IsRunning = false; public void Reset() { Pause(); Initialize(); } public void SetGridSize(int width, int height) {  Width = Math.Clamp(width, 16, 512);  Height = Math.Clamp(height, 16, 512);  Initialize(); } public void SetVisualizationMode(VisualizationMode mode) {  VisualizationMode = mode;  UpdateVisualization(); } }

结尾:以上代码仅供学习参考;