从零开始的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; //粘力-运动粘度系数vprivate 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层时的速度Xprivate double[,] uy0; // n层时的速度Yprivate 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<PointF>[] streamlines;private int streamlineCounter = 0;public event Action<Bitmap> VisualizationUpdated;public event Action<int> 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 < Height; j++){temp1 += (ux[i, j] - ux0[i, j]) * (ux[i, j] - ux0[i, j]) +(uy[i, j] - uy0[i, j]) * (uy[i, j] - uy0[i, j]);temp2 += ux[i, j] * ux[i, j] + uy[i, j] * uy[i, j];}}Error = Math.Sqrt(temp1) / (Math.Sqrt(temp2) + 1e-30);}private double feq(int k, double rho_val, double ux_val, double uy_val){double eu = cx[k] * ux_val + cy[k] * uy_val;double uv = ux_val * ux_val + uy_val * uy_val;return w[k] * rho_val * (1.0 + 3.0 * eu + 4.5 * eu * eu - 1.5 * uv);}#region 可视化方法public void UpdateVisualization(){if (streamlineCounter % 8 == 0){InitializeStreamlines();}streamlineCounter++;Bitmap fluidImage = VisualizationMode switch{VisualizationMode.Density => 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];// 双线性插值计算uxdouble ux_val =(1 - dx) * (1 - dy) * ux00 +dx * (1 - dy) * ux10 +(1 - dx) * dy * ux01 +dx * dy * ux11;// 双线性插值计算uydouble 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 < minPressure) minPressure = pressure;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; // Bluerow[i * bytesPerPixel + 1] = color.G; // Greenrow[i * bytesPerPixel + 2] = color.R; // Redrow[i * bytesPerPixel + 3] = 255; // Alpha}});}bmp.UnlockBits(bmpData);return bmp;}#region 绘制密度./// <summary>/// 绘制密度/// </summary>/// <returns></returns>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] < minDensity) minDensity = rho[x, 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; // Bluerow[x * bytesPerPixel + 1] = value; // Greenrow[x * bytesPerPixel + 2] = value; // Redrow[x * bytesPerPixel + 3] = 255; // Alpha}}}fluidImage.UnlockBits(bmpData);return fluidImage;}/// <summary>/// 绘制密度（增强对比度黑白灰版）/// </summary>/// <returns></returns>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 < minDensity) minDensity = density;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 < lowerBound || density > upperBound) continue;validPixelCount++;if (density < effectiveMin) effectiveMin = density;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; // Bluerow[x * bytesPerPixel + 1] = 0; // Greenrow[x * bytesPerPixel + 2] = 255; // Redrow[x * bytesPerPixel + 3] = 255; // Alphacontinue;}// 安全映射密度到[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; // Bluerow[x * bytesPerPixel + 1] = value; // Greenrow[x * bytesPerPixel + 2] = value; // Redrow[x * bytesPerPixel + 3] = 255; // Alpha}}}fluidImage.UnlockBits(bmpData);return fluidImage;}#endregionprivate 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 < Height; y++){double speed = Math.Sqrt(ux[x, y] * ux[x, y] + uy[x, y] * uy[x, y]);if (speed > 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; // Bluerow[x * bytesPerPixel + 1] = color.G; // Greenrow[x * bytesPerPixel + 2] = color.R; // Redrow[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 < Height; y += step){double speed = Math.Sqrt(ux[x, y] * ux[x, y] + uy[x, y] * uy[x, y]);//只绘制有显著速度的区域if (speed > 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, // Red2, // Green2 // Blue);DrawArrow(g, start, end, arrowColor, Math.Max(1, scale / 4));}}}}return fluidImage;}/// <summary>/// 漩涡/// </summary>/// <returns></returns>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 < Height - 1; y++){// 中心差分计算涡量double dudy = (uy[x, y + 1] - uy[x, y - 1]) / 2.0;double dvdx = (ux[x + 1, y] - ux[x - 1, y]) / 2.0;vorticity[x, y] = dvdx - dudy;if (Math.Abs(vorticity[x, 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 < Width; x++){double v = vorticity[x, y] / maxVorticity;int r = (v > 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; // Bluerow[x * bytesPerPixel + 1] = (byte)g; // Greenrow[x * bytesPerPixel + 2] = (byte)r; // Redrow[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 < streamlines.Length; i++){var line = streamlines[i];if (line.Count > 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<PointF> line){if (line.Count < 5) return Color.Red;double totalSpeed = 0;int count = 0;// 计算流线上点的平均速度foreach (PointF point in line){int x = (int)Math.Floor(point.X);int y = (int)Math.Floor(point.Y);if (x >= 0 && x < Width && y >= 0 && y < Height){double speed = Math.Sqrt(ux[x, y] * ux[x, y] + uy[x, y] * uy[x, y]);totalSpeed += speed;count++;}}if (count == 0) return Color.Red;double avgSpeed = totalSpeed / count;double normalizedSpeed = Math.Min(1.0, avgSpeed / U); // U是顶盖速度// 根据速度大小返回不同颜色if (normalizedSpeed > 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<PointF>[lineCount];streamlineCounter = 0; // 重置流线计数器Random rand = new Random();for (int i = 0; i < lineCount; i++){streamlines[i] = new List<PointF>();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 < 0 || nextX > Width || nextY < 0 || 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 < 0 || red > 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);}#endregionpublic 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();}}