起因:
上期我們寫了流體仿真的經典案例: 通過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();}}
