using iTextSharp.text.pdf.parser.clipper; using MathNet.Numerics.Interpolation; using System; using System.Collections.Generic; using System.Linq; using System.Text; namespace ConductivityApp.GBStandard { /// /// GB/T 32064-2015 标准计算器(严格遵循国标) /// 瞬态平面热源测试法计算导热系数和热扩散系数 /// public class GB32064Calculator { #region 结构定义 public struct BridgeConfig { public double Rs { get; set; } // 串联电阻 (Ω) public double RL { get; set; } // 引线电阻 (Ω) } public struct ProbeConfig { public double R0 { get; set; } // 初始电阻 (Ω) public double Alpha { get; set; } // 电阻温度系数 (1/K) public double RadiusMM { get; set; } // 探头半径 (mm) public int CoilCount { get; set; } // 探头环数 } public struct TestConfig { public double P0 { get; set; } // 输出功率 (W) public double SampleDensity { get; set; } // 样品密度 (kg/m³) public string cmbSampleType { get; set; } } public struct CalculationResult { public bool IsValid { get; set; } public double ThermalConductivity { get; set; } // λ - W/(m·K) public double ThermalDiffusivity { get; set; } // α - m²/s public double SpecificHeatCapacity { get; set; } // Cp - J/(kg·K) public double CorrectionTime { get; set; } // tc - s public double RSquared { get; set; } public string ValidationMessage { get; set; } public string CalculationLog { get; set; } // 新增:计算过程日志 } #endregion #region 常量定义(根据国标GB/T 32064-2015) private const double TAU_MAX_SQUARED_LOWER = 0.3; // τ_max²的最小值(国标5.3.4要求) private const double TAU_MAX_SQUARED_UPPER = 1.0; // τ_max²的最大值(国标5.3.4要求) private const double PROBING_DEPTH_RATIO_LOWER = 1.1; // 探测深度/探头半径最小值 private const double PROBING_DEPTH_RATIO_UPPER = 2.0; // 探测深度/探头半径最大值 private const double MIN_DATA_COUNT = 100; // 最小数据点数(国标5.3.3.3) private const double MIN_TIME_INTERVAL = 0.1; // 最小时间间隔(s)(国标5.3.3.3) private double PI_POW_1_5 = Math.Pow(Math.PI, 1.5); // π^(3/2) #endregion #region 私有字段 private readonly BridgeConfig _bridge; private readonly ProbeConfig _probe; private readonly TestConfig _test; private IInterpolation _tauDInterpolation; private StringBuilder _logBuilder; // 新增:日志记录器 #endregion #region 构造函数 public GB32064Calculator(BridgeConfig bridge, ProbeConfig probe, TestConfig test) { _bridge = bridge; _probe = probe; _test = test; _logBuilder = new StringBuilder(); LogInfo("GB32064Calculator 初始化开始"); LogInfo($"探头配置: R0={_probe.R0}Ω, α={_probe.Alpha}/K, 半径={_probe.RadiusMM}mm, 环数={_probe.CoilCount}"); LogInfo($"电桥配置: Rs={_bridge.Rs}Ω, RL={_bridge.RL}Ω"); LogInfo($"测试配置: P0={_test.P0}W, 密度={_test.SampleDensity}kg/m³"); ValidateConfigurations(); InitializeTauDTable(); LogInfo("GB32064Calculator 初始化完成"); } #endregion #region 可配置参数 private int _maxIterations = 300; private double _rSquaredTarget = 0.95; private double _convergenceThreshold = 1e-8; private int _maxStagnationIterations = 30; /// /// 设置迭代参数 /// public void SetIterationParameters(int maxIterations = 300, double rSquaredTarget = 0.95, double convergenceThreshold = 1e-8, int maxStagnation = 30) { _maxIterations = maxIterations; _rSquaredTarget = rSquaredTarget; _convergenceThreshold = convergenceThreshold; _maxStagnationIterations = maxStagnation; LogInfo($"设置迭代参数: 最大迭代={_maxIterations}, R²目标={_rSquaredTarget:F3}"); } #endregion #region 配置验证 private void ValidateConfigurations() { LogInfo("开始配置验证"); var errors = new List(); if (_probe.R0 <= 0) errors.Add("探头初始电阻R0必须大于0"); if (_probe.Alpha <= 0) errors.Add("探头温度系数α必须大于0"); if (_probe.RadiusMM <= 0) errors.Add("探头半径必须大于0"); if (_probe.CoilCount < 1) errors.Add("探头环数应大于0"); if (_bridge.Rs <= 0) errors.Add("串联电阻Rs必须大于0"); if (_bridge.RL < 0) errors.Add("引线电阻RL不能为负"); if (_test.P0 <= 0) errors.Add("输出功率P0必须大于0"); if (_test.SampleDensity <= 0) errors.Add("样品密度必须大于0"); if (errors.Count > 0) { LogError($"配置验证失败,发现{errors.Count}个错误"); throw new ArgumentException($"配置验证失败:\n{string.Join("\n", errors)}"); } else { LogInfo("配置验证通过"); } } #endregion #region 核心计算方法 private double GetDefaultAlphaBySoilType() { switch (_test.cmbSampleType) { case "干土": return 0.8e-6; // 干土的α更小,修正为0.3e-6 case "湿土": return 2.0e-6; // 湿土典型值 case "冻土": return 2.5e-6; // 冻土典型值 default: return 1.6e-6; // 通用默认值 } } public CalculationResult Calculate(double[] timeArray, double[] deltaUArray, double currentMA = 120.0) { LogInfo("=".PadRight(80, '=')); LogInfo("开始计算导热系数和热扩散系数"); LogInfo($"输入数据: 时间点数={timeArray?.Length}, 电压点数={deltaUArray?.Length}, 电流={currentMA}mA"); try { // 1. 数据验证 LogInfo("步骤1: 数据验证"); if (!ValidateInputData(timeArray, deltaUArray)) { LogError("输入数据验证失败"); return InvalidResult("输入数据验证失败"); } LogInfo($"数据验证通过: 时间范围[{timeArray.Min():F4}s, {timeArray.Max():F4}s], 共{timeArray.Length}个点"); double I0 = currentMA / 1000.0; // mA转A LogInfo($"电流转换: {currentMA}mA -> {I0:F4}A"); // 2. 计算ΔT(t) - 国标公式(3) LogInfo("步骤2: 计算温度增量ΔT(t) - 国标公式(3)"); double[] deltaT = CalculateTemperatureIncrements(deltaUArray, I0); LogInfo($"ΔT计算完成: 平均值={deltaT.Average():E6}K, 最大值={deltaT.Max():E6}K, 最小值={deltaT.Min():E6}K"); // 3. 迭代求解热扩散系数α和校正时间tc LogInfo("步骤3: 迭代求解热扩散系数α和校正时间tc"); var (alpha, tc, iterationLog) = IterateThermalDiffusivity(timeArray, deltaT); LogInfo(iterationLog); // 记录迭代过程 LogInfo($"迭代完成: α={alpha:E6} m²/s, tc={tc:F4}s"); // 4. 验证τ_max²是否符合国标要求 LogInfo("步骤4: 验证τ_max²是否符合国标要求"); double r = _probe.RadiusMM / 1000.0; double tmax = timeArray.Max(); double tauMaxSquared = CalculateTauMaxSquared(tmax, tc, alpha, r); LogInfo($"τ_max²计算: tmax={tmax:F2}s, tc={tc:F4}s, α={alpha:E6} m²/s, r={r:F6}m"); LogInfo($"τ_max²结果: {tauMaxSquared:F4} (要求范围: [{TAU_MAX_SQUARED_LOWER}, {TAU_MAX_SQUARED_UPPER}])"); if (tauMaxSquared < TAU_MAX_SQUARED_LOWER || tauMaxSquared > TAU_MAX_SQUARED_UPPER) { LogError($"τ_max²({tauMaxSquared:F3})不在国标要求范围内"); return InvalidResult($"τ_max²({tauMaxSquared:F3})不在国标要求范围[{TAU_MAX_SQUARED_LOWER}, {TAU_MAX_SQUARED_UPPER}]内,测试无效"); } LogInfo("τ_max²验证通过"); // 5. 计算导热系数λ - 国标公式(4) LogInfo("步骤5: 计算导热系数λ - 国标公式(4)"); double lambda = CalculateThermalConductivity(timeArray, deltaT, alpha, tc); LogInfo($"导热系数计算结果: λ={lambda:F6} W/(m·K)"); // 6. 计算比热容Cp - Cp = λ / (ρ × α) LogInfo("步骤6: 计算比热容Cp"); double cp = CalculateSpecificHeatCapacity(lambda, alpha); LogInfo($"比热容计算结果: Cp={cp:F2} J/(kg·K)"); // 7. 验证测试结果有效性 LogInfo("步骤7: 验证测试结果有效性"); var validation = ValidateTestResults(timeArray, alpha, tc, tauMaxSquared); LogInfo($"测试结果验证: {(validation.IsValid ? "通过" : "未通过")}"); // 8. 计算R² LogInfo("步骤8: 计算拟合优度R²"); double rSquared = CalculateRSquared(timeArray, deltaT, alpha, tc, lambda); LogInfo($"R²计算结果: {rSquared:F6}"); // 生成最终日志 LogInfo("=".PadRight(80, '=')); LogInfo("计算完成"); LogInfo($"最终结果: λ={lambda:F6} W/(m·K), α={alpha:E6} m²/s, Cp={cp:F2} J/(kg·K)"); LogInfo($"校正参数: tc={tc:F4}s, R²={rSquared:F6}"); return new CalculationResult { IsValid = validation.IsValid, ThermalConductivity = lambda, ThermalDiffusivity = alpha, SpecificHeatCapacity = cp, CorrectionTime = tc, RSquared = rSquared, ValidationMessage = string.Join("\n", validation.Messages), CalculationLog = _logBuilder.ToString() // 保存完整日志 }; } catch (Exception ex) { LogError($"计算失败: {ex.Message}"); LogError($"异常堆栈: {ex.StackTrace}"); return InvalidResult($"计算失败: {ex.Message}"); } } /// /// 计算温度增量ΔT(t) - 国标公式(3) /// private double[] CalculateTemperatureIncrements(double[] deltaUArray, double I0) { LogDebug($"开始计算ΔT(t),共{deltaUArray.Length}个点"); double[] deltaT = new double[deltaUArray.Length]; for (int i = 0; i < deltaUArray.Length; i++) { double deltaU = deltaUArray[i]; // 国标公式(3): ΔT(t) = (Rs + RL + R0) × ΔU(t) / [(I0 × Rs - ΔU(t)) × α × R0] double numerator = (_bridge.Rs + _bridge.RL + _probe.R0) * deltaU; double denominator = (I0 * _bridge.Rs - deltaU) * _probe.Alpha * _probe.R0; if (Math.Abs(denominator) < 1e-12) { LogError($"第{i}点计算ΔT(t)时分母接近0: denominator={denominator:E6}"); throw new InvalidOperationException($"分母接近0,无法计算ΔT(t)"); } deltaT[i] = numerator / denominator; // 每100个点记录一次进度 if (i % 100 == 0 && i > 0) { LogDebug($"已计算{i}个ΔT点,当前值={deltaT[i]:E6}K"); } } LogDebug($"ΔT(t)计算完成,最后10个值: {string.Join(", ", deltaT.Skip(Math.Max(0, deltaT.Length - 10)).Select(d => d.ToString("E3")))}"); return deltaT; } /// /// 迭代求解热扩散系数α和校正时间tc /// private (double alpha, double tc, string iterationLog) IterateThermalDiffusivity(double[] timeArray, double[] deltaT) { StringBuilder iterationLog = new StringBuilder(); iterationLog.AppendLine("迭代过程记录:"); iterationLog.AppendLine("迭代 | α | tc | 误差 | 改进 | 状态"); iterationLog.AppendLine("-".PadRight(80, '-')); // 初始猜测值 double bestAlpha = GetDefaultAlphaBySoilType(); // 常见建筑材料的热扩散系数 double bestTc = timeArray.Max() * 0.01; // 测试时间的1%作为初始校正时间 double bestError = double.MaxValue; LogInfo($"迭代初始值: α={bestAlpha:E6} m²/s, tc={bestTc:F4}s"); // 使用最小二乘法迭代优化 int maxIterations = 500; double convergenceThreshold = GetDefaultAlphaBySoilType(); convergenceThreshold = 1e-8; int iteration; bool converged = false; for (iteration = 0; iteration < maxIterations; iteration++) { double alphaStep = bestAlpha * 0.1 / (iteration + 1); double tcStep = 0.001 / (iteration + 1); var candidates = new List<(double alpha, double tc, double error)> { (bestAlpha, bestTc, CalculateFitError(timeArray, deltaT, bestAlpha, bestTc)), (bestAlpha + alphaStep, bestTc, CalculateFitError(timeArray, deltaT, bestAlpha + alphaStep, bestTc)), (bestAlpha - alphaStep, bestTc, CalculateFitError(timeArray, deltaT, bestAlpha - alphaStep, bestTc)), (bestAlpha, bestTc + tcStep, CalculateFitError(timeArray, deltaT, bestAlpha, bestTc + tcStep)), (bestAlpha, bestTc - tcStep, CalculateFitError(timeArray, deltaT, bestAlpha, bestTc - tcStep)) }; var bestCandidate = candidates.OrderBy(c => c.error).First(); double improvement = bestError - bestCandidate.error; string status = "搜索"; // 收敛条件 if (Math.Abs(improvement) < convergenceThreshold && iteration > 20) { status = "收敛"; converged = true; } iterationLog.AppendLine($"{iteration + 1,3} | {bestCandidate.alpha:E6} | {bestCandidate.tc:F6} | {bestCandidate.error:E6} | {improvement:E6} | {status}"); if (improvement > 0) { bestAlpha = bestCandidate.alpha; bestTc = bestCandidate.tc; bestError = bestCandidate.error; } // 每10次迭代记录一次详细状态 if (iteration % 10 == 0) { LogDebug($"迭代{iteration}: α={bestAlpha:E6}, tc={bestTc:F6}, 误差={bestError:E6}, 改进={improvement:E6}"); } // 收敛后停止迭代 if (converged && iteration > 20) { LogInfo($"第{iteration + 1}次迭代达到收敛条件,停止迭代"); break; } } iterationLog.AppendLine("-".PadRight(80, '-')); iterationLog.AppendLine($"总迭代次数: {iteration + 1}次"); iterationLog.AppendLine($"是否收敛: {(converged ? "是" : "否")}"); iterationLog.AppendLine($"最终误差: {bestError:E6}"); iterationLog.AppendLine($"步长调整: α步长={bestAlpha * 0.1 / (iteration + 1):E6}, tc步长={0.001 / (iteration + 1):E6}"); if (!converged && iteration >= maxIterations) { LogWarning($"达到最大迭代次数({maxIterations})仍未收敛"); iterationLog.AppendLine($"警告: 达到最大迭代次数仍未完全收敛"); } return (bestAlpha, bestTc, iterationLog.ToString()); } //private (double alpha, double tc, string iterationLog) IterateThermalDiffusivity(double[] timeArray, double[] deltaT) //{ // StringBuilder iterationLog = new StringBuilder(); // iterationLog.AppendLine("迭代过程记录:"); // iterationLog.AppendLine("迭代 | α | tc | 误差 | 改进 | 状态"); // iterationLog.AppendLine("-".PadRight(80, '-')); // // 强制从更大的α值开始搜索 // double r = _probe.RadiusMM / 1000.0; // double tmax = timeArray.Max(); // // 基于τ²=0.8估算初始α(之前可能估算偏小) // double initialAlpha = 0.8 * r * r / (tmax * 0.99); // τ²=0.8,假设tc≈0 // // 关键修改:强制初始α不小于1e-6 // initialAlpha = Math.Max(initialAlpha, 2.0e-6); // initialAlpha = Math.Min(initialAlpha, 3.0e-6); // double bestAlpha = initialAlpha; // double bestTc = tmax * 0.01; // double bestError = CalculateFitError(timeArray, deltaT, bestAlpha, bestTc); // LogInfo($"强制初始α: α={bestAlpha:E6} m²/s (基于τ²=0.8估算), tc={bestTc:F4}s"); // int maxIterations = 200; // double convergenceThreshold = 1e-8; // int iteration; // bool converged = false; // // 记录最佳解的改进历史 // List errorHistory = new List { bestError }; // for (iteration = 0; iteration < maxIterations; iteration++) // { // // 动态步长:初期探索大,后期精细调整 // double alphaStep = bestAlpha * (0.3 / (1 + iteration * 0.05)); // 衰减更慢 // double tcStep = 0.01 / (1 + iteration * 0.05); // // 确保最小步长 // alphaStep = Math.Max(alphaStep, bestAlpha * 0.05); // tcStep = Math.Max(tcStep, 0.001); // // 生成候选解 // var candidates = new List<(double alpha, double tc, double error)> //{ // (bestAlpha, bestTc, bestError), // (bestAlpha + alphaStep, bestTc, CalculateFitError(timeArray, deltaT, bestAlpha + alphaStep, bestTc)), // (bestAlpha - alphaStep, bestTc, CalculateFitError(timeArray, deltaT, bestAlpha - alphaStep, bestTc)), // (bestAlpha, bestTc + tcStep, CalculateFitError(timeArray, deltaT, bestAlpha, bestTc + tcStep)), // (bestAlpha, bestTc - tcStep, CalculateFitError(timeArray, deltaT, bestAlpha, bestTc - tcStep)), // // 强制探索更大的α值 // (bestAlpha * 1.5, bestTc, CalculateFitError(timeArray, deltaT, bestAlpha * 1.5, bestTc)), // (bestAlpha * 0.67, bestTc, CalculateFitError(timeArray, deltaT, bestAlpha * 0.67, bestTc)), //}; // var bestCandidate = candidates.OrderBy(c => c.error).First(); // double improvement = bestError - bestCandidate.error; // string status = "搜索"; // // 收敛条件:连续5次改进小于阈值 // errorHistory.Add(bestCandidate.error); // if (errorHistory.Count > 5) // { // double recentImprovement = errorHistory[errorHistory.Count - 6] - errorHistory[errorHistory.Count - 1]; // if (Math.Abs(recentImprovement) < convergenceThreshold * 5 && iteration > 20) // { // status = "收敛"; // converged = true; // } // } // iterationLog.AppendLine($"{iteration + 1,3} | {bestCandidate.alpha:E6} | {bestCandidate.tc:F6} | {bestCandidate.error:E6} | {improvement:E6} | {status}"); // if (improvement > 0 || (improvement == 0 && bestCandidate.error < bestError)) // { // bestAlpha = bestCandidate.alpha; // bestTc = bestCandidate.tc; // bestError = bestCandidate.error; // } // // 每10次迭代记录一次 // if (iteration % 10 == 0) // { // LogDebug($"迭代{iteration}: α={bestAlpha:E6}, tc={bestTc:F6}, 误差={bestError:E6}"); // // 如果连续多次没有改进,尝试跳跃 // if (iteration > 30 && errorHistory.Count > 10) // { // double lastImprovement = errorHistory[errorHistory.Count - 11] - errorHistory[errorHistory.Count - 1]; // if (Math.Abs(lastImprovement) < 1e-6) // { // // 随机扰动 // double randomFactor = 0.8 + (new Random().NextDouble() * 0.4); // bestAlpha *= randomFactor; // LogDebug($"迭代{iteration}: 停滞检测,随机扰动α乘以{randomFactor:F3}"); // } // } // } // if (converged && iteration > 20) // { // LogInfo($"第{iteration + 1}次迭代达到收敛条件"); // break; // } // } // iterationLog.AppendLine("-".PadRight(80, '-')); // iterationLog.AppendLine($"总迭代次数: {iteration + 1}次"); // iterationLog.AppendLine($"是否收敛: {(converged ? "是" : "否")}"); // iterationLog.AppendLine($"最终误差: {bestError:E6}"); // // 验证α是否合理 // if (bestAlpha < 5e-7) // { // LogWarning($"警告:热扩散系数α={bestAlpha:E6}可能过小!建议:"); // LogWarning($" 1. 检查样品是否真的导热很差"); // LogWarning($" 2. 检查探头参数R0、α是否正确"); // LogWarning($" 3. 尝试增加测试时间以获得更大的τ²"); // } // return (bestAlpha, bestTc, iterationLog.ToString()); //} /// /// 计算拟合误差 - 使用过原点回归的残差平方和 /// private double CalculateFitError(double[] timeArray, double[] deltaT, double alpha, double tc) { double r = _probe.RadiusMM / 1000.0; var dTauList = new List(); var deltaTList = new List(); // 收集有效数据点 for (int i = 0; i < timeArray.Length; i++) { if (timeArray[i] > tc) { double tau = CalculateTau(timeArray[i], tc, alpha, r); // 检查τ是否在合理范围内 - 修复:使用正确的τ范围 [√0.3, 1] = [0.5477, 1] if (tau >= Math.Sqrt(TAU_MAX_SQUARED_LOWER) && tau <= Math.Sqrt(TAU_MAX_SQUARED_UPPER)) { // 边界检查 if (tau < 0.01 || tau > 3.0) continue; try { double dTau = _tauDInterpolation.Interpolate(tau); dTauList.Add(dTau); deltaTList.Add(deltaT[i]); } catch (Exception ex) { LogDebug($"第{i}点插值失败: τ={tau:F4}, 错误={ex.Message}"); continue; } } } } if (dTauList.Count < 10) { LogDebug($"误差计算: 有效数据点不足({dTauList.Count}),返回最大误差"); return double.MaxValue; } // 过原点线性回归: ΔT = k * D(τ) double numerator = 0; double denominator = 0; for (int i = 0; i < dTauList.Count; i++) { numerator += dTauList[i] * deltaTList[i]; denominator += dTauList[i] * dTauList[i]; } if (Math.Abs(denominator) < 1e-12) { LogDebug("误差计算: 分母接近0,返回最大误差"); return double.MaxValue; } double slope = numerator / denominator; // k = Σ(D*ΔT) / Σ(D²) // 计算均方根误差 (RMSE) double totalError = 0; for (int i = 0; i < dTauList.Count; i++) { double predicted = slope * dTauList[i]; // 过原点预测值 double error = deltaTList[i] - predicted; totalError += error * error; } double rmse = Math.Sqrt(totalError / dTauList.Count); double relativeRmse = rmse / Math.Abs(deltaTList.Average()); // 相对误差 LogDebug($"误差计算: 使用{dTauList.Count}个点, 斜率={slope:E6}, RMSE={rmse:E6}, 相对误差={relativeRmse:P2}"); return relativeRmse; // 返回相对误差 } ///// ///// 计算拟合误差 - 增加对α变化的敏感性 ///// //private double CalculateFitError(double[] timeArray, double[] deltaT, double alpha, double tc) //{ // double r = _probe.RadiusMM / 1000.0; // var dTauList = new List(); // var deltaTList = new List(); // // 收集有效数据点 // for (int i = 0; i < timeArray.Length; i++) // { // if (timeArray[i] > tc) // { // double tau = CalculateTau(timeArray[i], tc, alpha, r); // // 检查τ是否在国标范围内 [√0.3, 1] = [0.5477, 1] // if (tau >= Math.Sqrt(TAU_MAX_SQUARED_LOWER) && tau <= Math.Sqrt(TAU_MAX_SQUARED_UPPER)) // { // if (tau < 0.01 || tau > 3.0) continue; // try // { // double dTau = _tauDInterpolation.Interpolate(tau); // dTauList.Add(dTau); // deltaTList.Add(deltaT[i]); // } // catch { continue; } // } // } // } // if (dTauList.Count < 10) // { // return double.MaxValue; // 返回最大误差 // } // // 过原点线性回归 // double numerator = 0, denominator = 0; // for (int i = 0; i < dTauList.Count; i++) // { // numerator += dTauList[i] * deltaTList[i]; // denominator += dTauList[i] * dTauList[i]; // } // if (Math.Abs(denominator) < 1e-12) return double.MaxValue; // double slope = numerator / denominator; // // 关键修改:增加惩罚项,鼓励α在合理范围内 // double rmse = 0; // for (int i = 0; i < dTauList.Count; i++) // { // double predicted = slope * dTauList[i]; // double error = deltaTList[i] - predicted; // rmse += error * error; // } // rmse = Math.Sqrt(rmse / dTauList.Count); // // 修改1:对α过小的情况增加惩罚 // double alphaPenalty = 0; // if (alpha < 5e-7) // 如果α小于0.5e-6,增加惩罚 // { // alphaPenalty = (5e-7 - alpha) * 100; // 惩罚系数 // } // // 修改2:对斜率合理性增加检查(根据物理意义) // double expectedSlopeMin = _test.P0 / (PI_POW_1_5 * r * 3.0); // λ=3时的最小斜率 // double expectedSlopeMax = _test.P0 / (PI_POW_1_5 * r * 0.1); // λ=0.1时的最大斜率 // double slopePenalty = 0; // if (slope < expectedSlopeMin || slope > expectedSlopeMax) // { // slopePenalty = 1.0; // 显著惩罚 // } // // 综合误差 = RMSE + 惩罚项 // double totalError = rmse + alphaPenalty + slopePenalty; // LogDebug($"误差计算: α={alpha:E6}, 使用{dTauList.Count}点, RMSE={rmse:E3}, α惩罚={alphaPenalty:E3}, 斜率惩罚={slopePenalty:E3}, 总误差={totalError:E3}"); // return totalError; //} /// /// 计算导热系数λ - 国标公式(4),使用国标范围内的数据 /// private double CalculateThermalConductivity(double[] timeArray, double[] deltaT, double alpha, double tc) { LogInfo("开始计算导热系数λ"); double r = _probe.RadiusMM / 1000.0; var validPoints = new List<(double dTau, double deltaT)>(); // 收集有效数据点(严格按国标τ范围筛选) int validCount = 0; int totalCount = 0; int tGreaterTcCount = 0; LogDebug("筛选有效数据点用于线性回归:"); for (int i = 0; i < timeArray.Length; i++) { totalCount++; // 只使用校正后的数据(t > t_c) if (timeArray[i] > tc) { tGreaterTcCount++; double tau = CalculateTau(timeArray[i], tc, alpha, r); // 只使用τ在国标要求范围内的数据点 if (tau >= Math.Sqrt(TAU_MAX_SQUARED_LOWER) && tau <= Math.Sqrt(TAU_MAX_SQUARED_UPPER)) { try { double dTau = _tauDInterpolation.Interpolate(tau); validPoints.Add((dTau, deltaT[i])); validCount++; // 每50个有效点记录一次 if (validCount % 50 == 0) { LogDebug($"已找到{validCount}个有效点,当前点: t={timeArray[i]:F4}s, τ={tau:F4}, D(τ)={dTau:F6}, ΔT={deltaT[i]:E6}K"); } } catch (Exception ex) { LogDebug($"第{i}点插值失败: τ={tau:F4}, 错误={ex.Message}"); continue; } } } } LogInfo($"数据点统计: 总数={totalCount}, t>tc点数={tGreaterTcCount}, τ有效点数={validCount}"); if (validPoints.Count < 10) { LogError($"有效数据点不足({validPoints.Count}),无法计算导热系数"); throw new InvalidOperationException( $"有效数据点不足({validPoints.Count}),无法计算导热系数\n" + $"国标要求:在τ范围[{TAU_MAX_SQUARED_LOWER}, {TAU_MAX_SQUARED_UPPER}]内应有足够数据点"); } LogInfo($"使用{validPoints.Count}个有效数据点进行线性回归"); // 线性回归:ΔT(τ) = [P₀ / (π^(3/2) × r × λ)] × D(τ) double[] dTauArray = validPoints.Select(p => p.dTau).ToArray(); double[] deltaTArray = validPoints.Select(p => p.deltaT).ToArray(); LogDebug($"D(τ)范围: [{dTauArray.Min():F6}, {dTauArray.Max():F6}]"); LogDebug($"ΔT范围: [{deltaTArray.Min():E6}, {deltaTArray.Max():E6}]"); var (slope, intercept) = LinearRegression(dTauArray, deltaTArray, true); LogInfo($"线性回归结果: 斜率={slope:E6}, 截距={intercept:E6}, R={CalculateCorrelation(dTauArray, deltaTArray):F6}"); // 检查斜率是否合理 if (Math.Abs(slope) < 1e-12) { LogError($"回归斜率接近0 ({slope:E6}),计算结果不可靠"); throw new InvalidOperationException("回归斜率接近0,计算结果不可靠"); } // 计算导热系数 λ = P₀ / (π^(3/2) × r × slope) - 国标公式(4) double lambda = _test.P0 / (PI_POW_1_5 * r * slope); LogInfo($"导热系数计算: P0={_test.P0}W, r={r}m, π^(3/2)={PI_POW_1_5:F6}, slope={slope:E6}"); LogInfo($"最终导热系数: λ={lambda:F6} W/(m·K)"); return lambda; } /// /// 计算无量纲时间τ - 国标公式(5) /// τ = √[(t - t_c) / (r²/a)] /// 等价于:τ = √[(t - t_c) * a / r²] /// private double CalculateTau(double time, double tc, double alpha, double r) { if (time <= tc) { return 0; } // 国标公式(5): τ = √[(t - t_c) / (r²/a)] double theta = (r * r) / alpha; // 特征时间 θ = r²/a double tau = Math.Sqrt((time - tc) / theta); return tau; } /// /// 计算比热容Cp - Cp = λ / (ρ × α) /// private double CalculateSpecificHeatCapacity(double lambda, double alpha) { LogInfo($"开始计算比热容: λ={lambda:F6}, α={alpha:E6}, ρ={_test.SampleDensity}"); if (alpha <= 0 || _test.SampleDensity <= 0) { LogError($"参数无效: α={alpha:E6}, ρ={_test.SampleDensity}"); throw new InvalidOperationException("热扩散系数或密度无效,无法计算比热容"); } double cp = lambda / (_test.SampleDensity * alpha); LogInfo($"比热容计算: Cp = {lambda:F6} / ({_test.SampleDensity} × {alpha:E6}) = {cp:F2} J/(kg·K)"); return cp; } /// /// 验证测试结果 - 严格遵循国标5.3.4节要求 /// private ValidationResult ValidateTestResults(double[] timeArray, double alpha, double tc, double tauMaxSquared) { LogInfo("开始测试结果验证"); var result = new ValidationResult { IsValid = true }; double r = _probe.RadiusMM / 1000.0; // 转换为米 double tmax = timeArray.Max(); // 1. 验证τ_max²是否满足国标要求 if (tauMaxSquared < TAU_MAX_SQUARED_LOWER) { result.IsValid = false; string msg = $"τ_max²({tauMaxSquared:F3}) < {TAU_MAX_SQUARED_LOWER},测试时间不足,结果无效"; LogError(msg); result.Messages.Add(msg); } else if (tauMaxSquared > TAU_MAX_SQUARED_UPPER) { result.IsValid = false; string msg = $"τ_max²({tauMaxSquared:F3}) > {TAU_MAX_SQUARED_UPPER},测试时间过长,结果无效"; LogError(msg); result.Messages.Add(msg); } else { string msg = $"τ_max²({tauMaxSquared:F3})满足国标要求[{TAU_MAX_SQUARED_LOWER}, {TAU_MAX_SQUARED_UPPER}]"; LogInfo(msg); result.Messages.Add(msg); } // 2. 计算探测深度 ΔP_probe = 2√(a·t_max) (国标5.3.4) double probingDepth = 2 * Math.Sqrt(alpha * tmax); double probingDepthRatio = probingDepth / r; LogInfo($"探测深度计算: 2×√({alpha:E6}×{tmax:F2}) = {probingDepth * 1000:F2}mm, 探头半径={r * 1000:F2}mm, 比值={probingDepthRatio:F2}"); // 3. 验证探测深度是否满足国标要求 if (probingDepthRatio < PROBING_DEPTH_RATIO_LOWER) { result.IsValid = false; string msg = $"探测深度/探头半径({probingDepthRatio:F2}) < {PROBING_DEPTH_RATIO_LOWER},样品厚度可能不足"; LogError(msg); result.Messages.Add(msg); } else if (probingDepthRatio > PROBING_DEPTH_RATIO_UPPER) { result.IsValid = false; string msg = $"探测深度/探头半径({probingDepthRatio:F2}) > {PROBING_DEPTH_RATIO_UPPER},可能超出测试范围"; LogError(msg); result.Messages.Add(msg); } else { string msg = $"探测深度/探头半径({probingDepthRatio:F2})满足国标要求[{PROBING_DEPTH_RATIO_LOWER}, {PROBING_DEPTH_RATIO_UPPER}]"; LogInfo(msg); result.Messages.Add(msg); } // 4. 验证数据点数量(国标5.3.3.3要求采集次数大于100次) if (timeArray.Length < MIN_DATA_COUNT) { string msg = $"警告:采集次数({timeArray.Length})少于国标要求的{MIN_DATA_COUNT}次"; LogWarning(msg); result.Messages.Add(msg); } else { string msg = $"采集次数({timeArray.Length})满足国标要求(≥{MIN_DATA_COUNT}次)"; LogInfo(msg); result.Messages.Add(msg); } // 5. 验证数据采集间隔(国标5.3.3.3要求不小于0.1s) if (timeArray.Length > 1) { double minInterval = timeArray[1] - timeArray[0]; for (int i = 2; i < timeArray.Length; i++) { double interval = timeArray[i] - timeArray[i - 1]; if (interval < minInterval) minInterval = interval; } if (minInterval < MIN_TIME_INTERVAL) { string msg = $"警告:数据采集间隔({minInterval:F3}s)小于国标要求的{MIN_TIME_INTERVAL}s"; LogWarning(msg); result.Messages.Add(msg); } else { string msg = $"数据采集间隔({minInterval:F3}s)满足国标要求(≥{MIN_TIME_INTERVAL}s)"; LogInfo(msg); result.Messages.Add(msg); } } // 6. 生成详细验证报告 result.Messages.Insert(0, $"【国标GB/T 32064-2015测试结果有效性验证】"); result.Messages.Insert(1, $"测试总时间: {tmax:F2}s,校正时间: {tc:F4}s"); result.Messages.Insert(2, $"热扩散系数α: {alpha:E6} m²/s"); result.Messages.Insert(3, $"探头半径r: {r * 1000:F2}mm,探测深度: {probingDepth * 1000:F2}mm"); result.Messages.Insert(4, $"验证依据:5.3.4节 测试结果有效性验证"); LogInfo($"测试结果验证完成: {(result.IsValid ? "有效" : "无效")}"); return result; } /// /// 计算τ_max² /// private double CalculateTauMaxSquared(double tmax, double tc, double alpha, double r) { LogDebug($"计算τ_max²: tmax={tmax:F4}, tc={tc:F4}, α={alpha:E6}, r={r:F6}"); // 国标公式:τ_max² = (t_max - t_c) × a / r² if (r <= 0) throw new ArgumentException("探头半径必须大于0"); if (tmax <= tc) { LogWarning($"tmax({tmax:F4}) <= tc({tc:F4}),返回0"); return 0; } //double tauMaxSquared = (tmax - tc) * alpha / (r * r); double tauMax = Math.Sqrt((tmax - tc) * alpha / (r * r)); double tauMaxSquared = tauMax * tauMax; // τ² = τ×τ LogDebug($"τ_max²计算结果: ({tmax}-{tc})×{alpha:E6}/{r * r:E6} = {tauMaxSquared:F4}"); return tauMaxSquared; } /// /// 计算R²拟合优度 /// private double CalculateRSquared(double[] timeArray, double[] deltaT, double alpha, double tc, double lambda) { LogInfo("开始计算R²拟合优度"); var observed = new List(); var predicted = new List(); double r = _probe.RadiusMM / 1000.0; int validCount = 0; for (int i = 0; i < timeArray.Length; i++) { if (timeArray[i] > tc) { double tau = CalculateTau(timeArray[i], tc, alpha, r); if (tau >= TAU_MAX_SQUARED_LOWER && tau <= TAU_MAX_SQUARED_UPPER) { try { double dTau = _tauDInterpolation.Interpolate(tau); double theoretical = _test.P0 * dTau / (PI_POW_1_5 * r * lambda); observed.Add(deltaT[i]); predicted.Add(theoretical); validCount++; } catch { continue; } } } } LogInfo($"R²计算使用{validCount}个有效数据点"); if (observed.Count < 10) { LogWarning($"R²计算数据点不足({observed.Count}),返回0"); return 0; } double meanObserved = observed.Average(); double ssTotal = observed.Sum(o => Math.Pow(o - meanObserved, 2)); double ssResidual = observed.Zip(predicted, (o, p) => Math.Pow(o - p, 2)).Sum(); if (ssResidual > ssTotal) { LogWarning($"SSR({ssResidual:E6}) > SST({ssTotal:E6}),强制设为0"); return 0; } LogDebug($"R²计算: 平均值={meanObserved:E6}, SST={ssTotal:E6}, SSR={ssResidual:E6}"); if (Math.Abs(ssTotal) < 1e-12) { LogWarning("SST接近0,返回0"); return 0; } double rSquared = 1 - (ssResidual / ssTotal); LogInfo($"R²计算结果: {rSquared:F6} (SST={ssTotal:E6}, SSR={ssResidual:E6})"); return rSquared; } #endregion #region 辅助方法 private CalculationResult InvalidResult(string message) { LogError($"返回无效结果: {message}"); return new CalculationResult { IsValid = false, ValidationMessage = message, CalculationLog = _logBuilder.ToString(), ThermalConductivity = 0, ThermalDiffusivity = 0, SpecificHeatCapacity = 0, CorrectionTime = 0, RSquared = 0 }; } private bool ValidateInputData(double[] timeArray, double[] deltaUArray) { LogDebug("开始验证输入数据"); if (timeArray == null || deltaUArray == null) { LogError("时间和电压数据不能为空"); throw new ArgumentNullException("时间和电压数据不能为空"); } if (timeArray.Length != deltaUArray.Length) { LogError($"数据长度不一致: 时间{timeArray.Length}点, 电压{deltaUArray.Length}点"); throw new ArgumentException("时间和电压数据长度不一致"); } if (timeArray.Length < 20) { LogError($"数据点数量({timeArray.Length})不足,至少需要20个点"); throw new ArgumentException($"数据点数量({timeArray.Length})不足,至少需要20个点"); } // 检查时间单调递增 for (int i = 1; i < timeArray.Length; i++) { if (timeArray[i] <= timeArray[i - 1]) { LogError($"时间数据非单调递增: 第{i - 1}点={timeArray[i - 1]}, 第{i}点={timeArray[i]}"); throw new ArgumentException($"时间数据必须严格单调递增,第{i}点不符合要求"); } } // 统计信息 double timeSpan = timeArray[timeArray.Length - 1] - timeArray[0]; double avgInterval = timeSpan / (timeArray.Length - 1); LogInfo($"数据验证通过: 时间范围[{timeArray[0]:F4}s, {timeArray[timeArray.Length - 1]:F4}s], 时长={timeSpan:F2}s"); LogInfo($"数据统计: 点数={timeArray.Length}, 平均间隔={avgInterval:F4}s"); LogInfo($"电压范围: [{deltaUArray.Min():E6}V, {deltaUArray.Max():E6}V], 平均值={deltaUArray.Average():E6}V"); return true; } private (double slope, double intercept) LinearRegression(double[] x, double[] y, bool forceOrigin = true) { LogDebug($"线性回归: 输入{x.Length}个点"); if (x.Length != y.Length || x.Length < 2) { LogError($"线性回归参数错误: x长度={x.Length}, y长度={y.Length}"); throw new ArgumentException("线性回归需要至少2个数据点且x、y长度相等"); } if (forceOrigin) { // 过原点回归:y = kx double numerator = 0; double denominator = 0; for (int i = 0; i < x.Length; i++) { numerator += x[i] * y[i]; denominator += x[i] * x[i]; } if (Math.Abs(denominator) < 1e-12) { LogError("数据方差太小,无法进行线性回归"); throw new InvalidOperationException("数据方差太小,无法进行线性回归"); } double slope = numerator / denominator; double intercept = 0; LogDebug($"过原点回归: 斜率={slope:E6}"); return (slope, intercept); } else { double xAvg = x.Average(); double yAvg = y.Average(); LogDebug($"平均值: x̄={xAvg:E6}, ȳ={yAvg:E6}"); double numerator = 0; double denominator = 0; for (int i = 0; i < x.Length; i++) { double xDiff = x[i] - xAvg; double yDiff = y[i] - yAvg; numerator += xDiff * yDiff; denominator += xDiff * xDiff; } LogDebug($"回归计算: 分子={numerator:E6}, 分母={denominator:E6}"); if (Math.Abs(denominator) < 1e-12) { LogError("数据方差太小,无法进行线性回归"); throw new InvalidOperationException("数据方差太小,无法进行线性回归"); } double slope = numerator / denominator; double intercept = yAvg - slope * xAvg; LogDebug($"回归结果: 斜率={slope:E6}, 截距={intercept:E6}"); return (slope, intercept); } } /// /// 计算相关系数 /// private double CalculateCorrelation(double[] x, double[] y) { if (x.Length != y.Length || x.Length < 2) return 0; double xAvg = x.Average(); double yAvg = y.Average(); double numerator = 0; double xSumSq = 0; double ySumSq = 0; for (int i = 0; i < x.Length; i++) { double xDiff = x[i] - xAvg; double yDiff = y[i] - yAvg; numerator += xDiff * yDiff; xSumSq += xDiff * xDiff; ySumSq += yDiff * yDiff; } if (xSumSq == 0 || ySumSq == 0) return 0; return numerator / Math.Sqrt(xSumSq * ySumSq); } #endregion #region 日志记录方法 private void LogInfo(string message) { string logEntry = $"[INFO] {DateTime.Now:HH:mm:ss.fff} - {message}"; _logBuilder.AppendLine(logEntry); // 这里也可以输出到控制台或文件 // Console.WriteLine(logEntry); } private void LogDebug(string message) { string logEntry = $"[DEBUG] {DateTime.Now:HH:mm:ss.fff} - {message}"; _logBuilder.AppendLine(logEntry); // 调试信息可根据需要启用或禁用 // Console.WriteLine(logEntry); } private void LogWarning(string message) { string logEntry = $"[WARNING] {DateTime.Now:HH:mm:ss.fff} - {message}"; _logBuilder.AppendLine(logEntry); // Console.WriteLine(logEntry); } private void LogError(string message) { string logEntry = $"[ERROR] {DateTime.Now:HH:mm:ss.fff} - {message}"; _logBuilder.AppendLine(logEntry); // Console.WriteLine(logEntry); } #endregion #region τ-D(τ)插值表初始化 private void InitializeTauDTable() { LogInfo("初始化τ-D(τ)插值表"); // 根据国标GB/T 32064-2015,瞬态平面热源法的D(τ)函数 // D(τ) = [√(τ²+1) - 1] / τ 或 D(τ) = [√(τ²+1) - τ] / [√(τ²+1) + τ] // 需要确认哪个是正确的公式! var tauList = new List(); var dList = new List(); LogDebug("生成τ-D(τ)数据点"); // 先计算并记录两种公式的差异 LogDebug("τ-D(τ)公式验证:"); for (double tau = 0.5; tau <= 2.5; tau += 0.5) { double d1 = (Math.Sqrt(tau * tau + 1.0) - 1.0) / tau; double d2 = (Math.Sqrt(tau * tau + 1.0) - tau) / (Math.Sqrt(tau * tau + 1.0) + tau); LogDebug($"τ={tau:F2}: 公式1(D1)={d1:F6}, 公式2(D2)={d2:F6}, 比值={d1 / d2:F3}"); } for (double tau = 0.01; tau <= 3.0; tau += 0.01) { tauList.Add(tau); double d = (Math.Sqrt(tau * tau + 1.0) - 1.0) / tau; dList.Add(d); } double[] tauValues = tauList.ToArray(); double[] dValues = dList.ToArray(); _tauDInterpolation = CubicSpline.InterpolateAkimaSorted(tauValues, dValues); LogInfo($"τ-D(τ)插值表初始化完成,共{tauValues.Length}个点,τ范围[{tauValues[0]:F2}, {tauValues[tauValues.Length - 1]:F2}]"); } #endregion #region 内部类 private class ValidationResult { public bool IsValid { get; set; } public List Messages { get; } = new List(); } #endregion #region 土壤测试辅助方法(但不影响核心计算) /// /// 用于土壤测试的参数设置建议(仅建议,不影响计算) /// public static class SoilTestRecommendations { /// /// 根据土壤类型获取建议的测试参数 /// /// 土壤类型(干土、湿土、冻土等) /// 建议的测试时间范围 public static (double minTime, double maxTime) GetRecommendedTestTime(string soilType) { // 土壤热扩散系数典型范围:1e-7 到 1e-6 m²/s // 使用典型探头半径6.4mm计算 double r = 0.0064; // 米 // 根据τ_max² = (t * a) / r²,反推时间 // t = τ_max² * r² / a // 对于导热系数较低的土壤(干土、冻土) if (soilType.Contains("冻土") || soilType.Contains("干土")) { double alpha = 0.5e-6; // 较低的热扩散系数 double t_min = TAU_MAX_SQUARED_LOWER * r * r / alpha; double t_max = TAU_MAX_SQUARED_UPPER * r * r / alpha; return (t_min, t_max); } // 对于导热系数较高的土壤(湿土) else if (soilType.Contains("湿土")) { double alpha = 1.2e-6; // 较高的热扩散系数 double t_min = TAU_MAX_SQUARED_LOWER * r * r / alpha; double t_max = TAU_MAX_SQUARED_UPPER * r * r / alpha; return (t_min, t_max); } // 默认值 else { double alpha = 1e-6; // 典型热扩散系数 double t_min = TAU_MAX_SQUARED_LOWER * r * r / alpha; double t_max = TAU_MAX_SQUARED_UPPER * r * r / alpha; return (t_min, t_max); } } } #endregion } }