659 lines
26 KiB
C#
659 lines
26 KiB
C#
using MathNet.Numerics.Interpolation;
|
||
using System;
|
||
using System.Collections.Generic;
|
||
using System.Linq;
|
||
|
||
namespace ConductivityApp.GBStandard
|
||
{
|
||
/// <summary>
|
||
/// GB/T 32064-2015 标准计算器（严格遵循国标）
|
||
/// 瞬态平面热源测试法计算导热系数和热扩散系数
|
||
/// </summary>
|
||
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 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; }
|
||
}
|
||
#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;
|
||
#endregion
|
||
|
||
#region 构造函数
|
||
public GB32064Calculator(BridgeConfig bridge, ProbeConfig probe, TestConfig test)
|
||
{
|
||
_bridge = bridge;
|
||
_probe = probe;
|
||
_test = test;
|
||
|
||
ValidateConfigurations();
|
||
InitializeTauDTable();
|
||
}
|
||
#endregion
|
||
|
||
#region 配置验证
|
||
private void ValidateConfigurations()
|
||
{
|
||
var errors = new List<string>();
|
||
|
||
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)
|
||
{
|
||
throw new ArgumentException($"配置验证失败:\n{string.Join("\n", errors)}");
|
||
}
|
||
}
|
||
#endregion
|
||
|
||
#region 核心计算方法
|
||
|
||
public CalculationResult Calculate(double[] timeArray, double[] deltaUArray, double currentMA = 120.0)
|
||
{
|
||
try
|
||
{
|
||
// 1. 数据验证
|
||
if (!ValidateInputData(timeArray, deltaUArray))
|
||
{
|
||
return InvalidResult("输入数据验证失败");
|
||
}
|
||
|
||
double I0 = currentMA / 1000.0; // mA转A
|
||
|
||
// 2. 计算ΔT(t) - 国标公式(3)
|
||
double[] deltaT = CalculateTemperatureIncrements(deltaUArray, I0);
|
||
|
||
// 3. 迭代求解热扩散系数α和校正时间tc
|
||
var (alpha, tc) = IterateThermalDiffusivity(timeArray, deltaT);
|
||
|
||
// 4. 验证τ_max²是否符合国标要求
|
||
double r = _probe.RadiusMM / 1000.0;
|
||
double tmax = timeArray.Max();
|
||
double tauMaxSquared = CalculateTauMaxSquared(tmax, tc, alpha, r);
|
||
|
||
if (tauMaxSquared < TAU_MAX_SQUARED_LOWER || tauMaxSquared > TAU_MAX_SQUARED_UPPER)
|
||
{
|
||
return InvalidResult($"τ_max²({tauMaxSquared:F3})不在国标要求范围[{TAU_MAX_SQUARED_LOWER}, {TAU_MAX_SQUARED_UPPER}]内，测试无效");
|
||
}
|
||
|
||
// 5. 计算导热系数λ - 国标公式(4)
|
||
double lambda = CalculateThermalConductivity(timeArray, deltaT, alpha, tc);
|
||
|
||
// 6. 计算比热容Cp - Cp = λ / (ρ × α)
|
||
double cp = CalculateSpecificHeatCapacity(lambda, alpha);
|
||
|
||
// 7. 验证测试结果有效性
|
||
var validation = ValidateTestResults(timeArray, alpha, tc, tauMaxSquared);
|
||
|
||
// 8. 计算R²
|
||
double rSquared = CalculateRSquared(timeArray, deltaT, alpha, tc, lambda);
|
||
|
||
return new CalculationResult
|
||
{
|
||
IsValid = validation.IsValid,
|
||
ThermalConductivity = lambda,
|
||
ThermalDiffusivity = alpha,
|
||
SpecificHeatCapacity = cp,
|
||
CorrectionTime = tc,
|
||
RSquared = rSquared,
|
||
ValidationMessage = string.Join("\n", validation.Messages)
|
||
};
|
||
}
|
||
catch (Exception ex)
|
||
{
|
||
return InvalidResult($"计算失败: {ex.Message}");
|
||
}
|
||
}
|
||
|
||
/// <summary>
|
||
/// 计算温度增量ΔT(t) - 国标公式(3)
|
||
/// </summary>
|
||
private double[] CalculateTemperatureIncrements(double[] deltaUArray, double I0)
|
||
{
|
||
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)
|
||
{
|
||
throw new InvalidOperationException($"分母接近0，无法计算ΔT(t)");
|
||
}
|
||
|
||
deltaT[i] = numerator / denominator;
|
||
}
|
||
|
||
return deltaT;
|
||
}
|
||
|
||
/// <summary>
|
||
/// 迭代求解热扩散系数α和校正时间tc
|
||
/// </summary>
|
||
private (double alpha, double tc) IterateThermalDiffusivity(double[] timeArray, double[] deltaT)
|
||
{
|
||
// 初始猜测值
|
||
double bestAlpha = 1e-6; // 常见建筑材料的热扩散系数
|
||
//double bestAlpha = 3.58e-7; // 常见建筑材料的热扩散系数
|
||
double bestTc = timeArray.Max() * 0.01; // 测试时间的1%作为初始校正时间
|
||
double bestError = double.MaxValue;
|
||
|
||
// 使用最小二乘法迭代优化
|
||
int maxIterations = 200;
|
||
double convergenceThreshold = 1e-8;
|
||
|
||
for (int 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;
|
||
if (improvement > 0)
|
||
{
|
||
bestAlpha = bestCandidate.alpha;
|
||
bestTc = bestCandidate.tc;
|
||
bestError = bestCandidate.error;
|
||
}
|
||
|
||
// 收敛条件
|
||
if (Math.Abs(improvement) < convergenceThreshold && iteration > 20)
|
||
{
|
||
break;
|
||
}
|
||
}
|
||
|
||
return (bestAlpha, bestTc);
|
||
}
|
||
|
||
/// <summary>
|
||
/// 计算拟合误差
|
||
/// </summary>
|
||
private double CalculateFitError(double[] timeArray, double[] deltaT, double alpha, double tc)
|
||
{
|
||
double totalError = 0;
|
||
int count = 0;
|
||
double r = _probe.RadiusMM / 1000.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);
|
||
// 理论ΔT = [P0 / (π^(3/2) × r × λ)] × D(τ)
|
||
// 这里假设λ=1，因为我们只关心相对误差
|
||
double theoreticalDeltaT = _test.P0 * dTau / (PI_POW_1_5 * r);
|
||
double error = Math.Pow(deltaT[i] - theoreticalDeltaT, 2);
|
||
totalError += error;
|
||
count++;
|
||
}
|
||
catch
|
||
{
|
||
continue;
|
||
}
|
||
}
|
||
}
|
||
}
|
||
|
||
return count > 0 ? Math.Sqrt(totalError / count) : double.MaxValue;
|
||
}
|
||
|
||
/// <summary>
|
||
/// 计算导热系数λ - 国标公式(4)，使用国标范围内的数据
|
||
/// </summary>
|
||
private double CalculateThermalConductivity(double[] timeArray, double[] deltaT, double alpha, double tc)
|
||
{
|
||
double r = _probe.RadiusMM / 1000.0;
|
||
var validPoints = new List<(double dTau, double deltaT)>();
|
||
|
||
// 收集有效数据点（严格按国标τ范围筛选）
|
||
int validCount = 0;
|
||
int totalCount = 0;
|
||
|
||
for (int i = 0; i < timeArray.Length; i++)
|
||
{
|
||
totalCount++;
|
||
|
||
// 只使用校正后的数据（t > t_c）
|
||
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);
|
||
validPoints.Add((dTau, deltaT[i]));
|
||
validCount++;
|
||
}
|
||
catch
|
||
{
|
||
continue;
|
||
}
|
||
}
|
||
}
|
||
}
|
||
|
||
if (validPoints.Count < 10)
|
||
{
|
||
throw new InvalidOperationException(
|
||
$"有效数据点不足({validPoints.Count})，无法计算导热系数\n" +
|
||
$"国标要求：在τ范围[{TAU_MAX_SQUARED_LOWER}, {TAU_MAX_SQUARED_UPPER}]内应有足够数据点");
|
||
}
|
||
|
||
// 线性回归：ΔT(τ) = [P₀ / (π^(3/2) × r × λ)] × D(τ)
|
||
double[] dTauArray = validPoints.Select(p => p.dTau).ToArray();
|
||
double[] deltaTArray = validPoints.Select(p => p.deltaT).ToArray();
|
||
|
||
var (slope, intercept) = LinearRegression(dTauArray, deltaTArray);
|
||
|
||
// 检查斜率是否合理
|
||
if (Math.Abs(slope) < 1e-12)
|
||
{
|
||
throw new InvalidOperationException("回归斜率接近0，计算结果不可靠");
|
||
}
|
||
|
||
// 计算导热系数 λ = P₀ / (π^(3/2) × r × slope) - 国标公式(4)
|
||
double lambda = _test.P0 / (PI_POW_1_5 * r * slope);
|
||
|
||
return lambda;
|
||
}
|
||
|
||
/// <summary>
|
||
/// 计算无量纲时间τ - 国标公式(5)
|
||
/// τ = √[(t - t_c) / (r²/a)]
|
||
/// 等价于：τ = √[(t - t_c) * a / r²]
|
||
/// </summary>
|
||
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;
|
||
}
|
||
|
||
/// <summary>
|
||
/// 计算比热容Cp - Cp = λ / (ρ × α)
|
||
/// </summary>
|
||
private double CalculateSpecificHeatCapacity(double lambda, double alpha)
|
||
{
|
||
if (alpha <= 0 || _test.SampleDensity <= 0)
|
||
{
|
||
throw new InvalidOperationException("热扩散系数或密度无效，无法计算比热容");
|
||
}
|
||
|
||
double cp = lambda / (_test.SampleDensity * alpha);
|
||
return cp;
|
||
}
|
||
|
||
/// <summary>
|
||
/// 验证测试结果 - 严格遵循国标5.3.4节要求
|
||
/// </summary>
|
||
private ValidationResult ValidateTestResults(double[] timeArray, double alpha, double tc, double tauMaxSquared)
|
||
{
|
||
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;
|
||
result.Messages.Add($"τ_max²({tauMaxSquared:F3}) < {TAU_MAX_SQUARED_LOWER}，测试时间不足，结果无效");
|
||
}
|
||
else if (tauMaxSquared > TAU_MAX_SQUARED_UPPER)
|
||
{
|
||
result.IsValid = false;
|
||
result.Messages.Add($"τ_max²({tauMaxSquared:F3}) > {TAU_MAX_SQUARED_UPPER}，测试时间过长，结果无效");
|
||
}
|
||
else
|
||
{
|
||
result.Messages.Add($"τ_max²({tauMaxSquared:F3})满足国标要求[{TAU_MAX_SQUARED_LOWER}, {TAU_MAX_SQUARED_UPPER}]");
|
||
}
|
||
|
||
// 2. 计算探测深度 ΔP_probe = 2√(a·t_max) （国标5.3.4）
|
||
double probingDepth = 2 * Math.Sqrt(alpha * tmax);
|
||
double probingDepthRatio = probingDepth / r;
|
||
|
||
// 3. 验证探测深度是否满足国标要求
|
||
if (probingDepthRatio < PROBING_DEPTH_RATIO_LOWER)
|
||
{
|
||
result.IsValid = false;
|
||
result.Messages.Add($"探测深度/探头半径({probingDepthRatio:F2}) < {PROBING_DEPTH_RATIO_LOWER}，样品厚度可能不足");
|
||
}
|
||
else if (probingDepthRatio > PROBING_DEPTH_RATIO_UPPER)
|
||
{
|
||
result.IsValid = false;
|
||
result.Messages.Add($"探测深度/探头半径({probingDepthRatio:F2}) > {PROBING_DEPTH_RATIO_UPPER}，可能超出测试范围");
|
||
}
|
||
else
|
||
{
|
||
result.Messages.Add($"探测深度/探头半径({probingDepthRatio:F2})满足国标要求[{PROBING_DEPTH_RATIO_LOWER}, {PROBING_DEPTH_RATIO_UPPER}]");
|
||
}
|
||
|
||
// 4. 验证数据点数量（国标5.3.3.3要求采集次数大于100次）
|
||
if (timeArray.Length < MIN_DATA_COUNT)
|
||
{
|
||
result.Messages.Add($"警告：采集次数({timeArray.Length})少于国标要求的{MIN_DATA_COUNT}次");
|
||
}
|
||
else
|
||
{
|
||
result.Messages.Add($"采集次数({timeArray.Length})满足国标要求(≥{MIN_DATA_COUNT}次)");
|
||
}
|
||
|
||
// 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)
|
||
{
|
||
result.Messages.Add($"警告：数据采集间隔({minInterval:F3}s)小于国标要求的{MIN_TIME_INTERVAL}s");
|
||
}
|
||
else
|
||
{
|
||
result.Messages.Add($"数据采集间隔({minInterval:F3}s)满足国标要求(≥{MIN_TIME_INTERVAL}s)");
|
||
}
|
||
}
|
||
|
||
// 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节 测试结果有效性验证");
|
||
|
||
return result;
|
||
}
|
||
|
||
/// <summary>
|
||
/// 计算τ_max²
|
||
/// </summary>
|
||
private double CalculateTauMaxSquared(double tmax, double tc, double alpha, double r)
|
||
{
|
||
// 国标公式：τ_max² = (t_max - t_c) × a / r²
|
||
if (r <= 0) throw new ArgumentException("探头半径必须大于0");
|
||
if (tmax <= tc) return 0;
|
||
|
||
double tauMaxSquared = (tmax - tc) * alpha / (r * r);
|
||
return tauMaxSquared;
|
||
}
|
||
|
||
/// <summary>
|
||
/// 计算R²拟合优度
|
||
/// </summary>
|
||
private double CalculateRSquared(double[] timeArray, double[] deltaT, double alpha, double tc, double lambda)
|
||
{
|
||
var observed = new List<double>();
|
||
var predicted = new List<double>();
|
||
double r = _probe.RadiusMM / 1000.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);
|
||
}
|
||
catch
|
||
{
|
||
continue;
|
||
}
|
||
}
|
||
}
|
||
}
|
||
|
||
if (observed.Count < 2) 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 (Math.Abs(ssTotal) < 1e-12) return 0;
|
||
|
||
double rSquared = 1 - (ssResidual / ssTotal);
|
||
return rSquared;
|
||
}
|
||
#endregion
|
||
|
||
#region 辅助方法
|
||
private CalculationResult InvalidResult(string message)
|
||
{
|
||
return new CalculationResult
|
||
{
|
||
IsValid = false,
|
||
ValidationMessage = message,
|
||
ThermalConductivity = 0,
|
||
ThermalDiffusivity = 0,
|
||
SpecificHeatCapacity = 0,
|
||
CorrectionTime = 0,
|
||
RSquared = 0
|
||
};
|
||
}
|
||
|
||
private bool ValidateInputData(double[] timeArray, double[] deltaUArray)
|
||
{
|
||
if (timeArray == null || deltaUArray == null)
|
||
{
|
||
throw new ArgumentNullException("时间和电压数据不能为空");
|
||
}
|
||
|
||
if (timeArray.Length != deltaUArray.Length)
|
||
{
|
||
throw new ArgumentException("时间和电压数据长度不一致");
|
||
}
|
||
|
||
if (timeArray.Length < 20)
|
||
{
|
||
throw new ArgumentException($"数据点数量({timeArray.Length})不足，至少需要20个点");
|
||
}
|
||
|
||
// 检查时间单调递增
|
||
for (int i = 1; i < timeArray.Length; i++)
|
||
{
|
||
if (timeArray[i] <= timeArray[i - 1])
|
||
{
|
||
throw new ArgumentException($"时间数据必须严格单调递增，第{i}点不符合要求");
|
||
}
|
||
}
|
||
|
||
return true;
|
||
}
|
||
|
||
private (double slope, double intercept) LinearRegression(double[] x, double[] y)
|
||
{
|
||
if (x.Length != y.Length || x.Length < 2)
|
||
throw new ArgumentException("线性回归需要至少2个数据点且x、y长度相等");
|
||
|
||
double xAvg = x.Average();
|
||
double yAvg = y.Average();
|
||
|
||
double numerator = 0;
|
||
double denominator = 0;
|
||
|
||
for (int i = 0; i < x.Length; i++)
|
||
{
|
||
numerator += (x[i] - xAvg) * (y[i] - yAvg);
|
||
denominator += (x[i] - xAvg) * (x[i] - xAvg);
|
||
}
|
||
|
||
if (Math.Abs(denominator) < 1e-12)
|
||
throw new InvalidOperationException("数据方差太小，无法进行线性回归");
|
||
|
||
double slope = numerator / denominator;
|
||
double intercept = yAvg - slope * xAvg;
|
||
|
||
return (slope, intercept);
|
||
}
|
||
#endregion
|
||
|
||
#region τ-D(τ)插值表初始化
|
||
private void InitializeTauDTable()
|
||
{
|
||
// 使用国标理论公式计算 D(τ) = (√(τ²+1) - 1) / τ
|
||
var tauList = new List<double>();
|
||
var dList = new List<double>();
|
||
|
||
// 根据国标要求，τ范围通常为0.01-3.0
|
||
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);
|
||
}
|
||
#endregion
|
||
|
||
#region 内部类
|
||
private class ValidationResult
|
||
{
|
||
public bool IsValid { get; set; }
|
||
public List<string> Messages { get; } = new List<string>();
|
||
}
|
||
#endregion
|
||
|
||
#region 土壤测试辅助方法（但不影响核心计算）
|
||
/// <summary>
|
||
/// 用于土壤测试的参数设置建议（仅建议，不影响计算）
|
||
/// </summary>
|
||
public static class SoilTestRecommendations
|
||
{
|
||
/// <summary>
|
||
/// 根据土壤类型获取建议的测试参数
|
||
/// </summary>
|
||
/// <param name="soilType">土壤类型（干土、湿土、冻土等）</param>
|
||
/// <returns>建议的测试时间范围</returns>
|
||
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
|
||
}
|
||
} |