#include "stdafx.h"
#include "HDuctHydraulic.h"
#include "math.h"

#define ZERO_AIR 1.0e-6

HDuctHydraulic::HDuctHydraulic()
{
	m_dFlux = 0.0;      // 风量
	m_dRough = 0.15;	// 内壁粗糙度
	m_dTem = 20.0;		// 温度

	m_nPipeShape = 1;	// 类型
	m_dWidth = 200.0;   // 管宽
	m_dHeight = 200.0;  // 管高
	m_dDia = 200.0;     // 直径
	
	m_dLength = 0.0;
	m_dVel = 0.0;
	m_dRe = 0.0;
	m_dFricCoef = 0.0;
	m_dR = 0.0;
	m_dFricLoss = 0.0;
	m_dLocCoef = 0.0;
	m_dOtherLocLoss = 0.0;
	m_dLocLoss = 0.0;
	m_dTotalLoss = 0.0;
	
	m_dDensity = 1.205;
	m_dMoveVis = 15.7e-6;
	m_csPipeName = _T("");
	
	m_dHeightStart = 0.0;
	m_dHeightEnd = 0.0;
	m_dVelDia = 0.0;
	m_dFluxDia = 0.0;
	m_dOutAirDensity = 0.0;
}

HDuctHydraulic::~HDuctHydraulic()
{
	
}

double HDuctHydraulic::GetVel() const
{
	return m_dVel;
}
double HDuctHydraulic::GetR() const
{
	return m_dR;
}
double HDuctHydraulic::GetFricLoss() const
{
	return m_dFricLoss;
}
void HDuctHydraulic::SetFlux(double dFlux)
{
	if (dFlux < 0.0)
		return;
	m_dFlux = dFlux;
}
void HDuctHydraulic::SetShape(BOOL bIsCircle)
{
	if (bIsCircle)
	{
		m_nPipeShape = 0;
	}
	else
	{
		m_nPipeShape = 1;
	}
}
void HDuctHydraulic::SetDia(double dDia)
{
	if (dDia < 0.0)
		return;
	m_dDia = dDia;
}
void HDuctHydraulic::SetWidth(double dWidth)
{
	if (dWidth < 0.0)
		return;
	m_dWidth = dWidth;
}
void HDuctHydraulic::SetHeight(double dHeight)
{
	if (dHeight < 0.0)
		return;
	m_dHeight = dHeight;
}
void HDuctHydraulic::SetLength(double dLength)
{
	if (dLength < 0.0)
		return;
	m_dLength = dLength;
}

void HDuctHydraulic::SetRough(double dRough)
{
	if (dRough < 0.0)
		return;
	m_dRough = dRough;
}
void HDuctHydraulic::SetTemperature(double dTem)
{
	m_dTem = dTem;
	m_dDensity = CalDesity(dTem);
	m_dMoveVis = CalMoveVis(dTem);
}
void HDuctHydraulic::SetDensity(double dDensity)
{
	m_dDensity = dDensity;
}
void HDuctHydraulic::SetMoveVis(double dMoveVis)
{
	m_dMoveVis = dMoveVis;
}

// 计算沿程阻力(计算前设置流量，形状，管径，长度，密度，黏度，粗糙度)
void HDuctHydraulic::CalFricLoss()
{
	if (m_nPipeShape == 0)
	{
		m_dVel = CalVelocity(m_dFlux, m_dDia);				// 计算流速
		m_dRe = CalRe(m_dVel, m_dDia, m_dMoveVis);			// 计算Re
		m_dFricCoef = CalFricCoef(m_dRough, m_dDia, m_dRe);	// 计算摩擦系数
		m_dR = CalR1(m_dVel, m_dFricCoef, m_dDensity, m_dDia);	// 计算比摩阻	
		m_dFricLoss = CalFricLoss(m_dFricCoef, m_dDia, m_dLength, m_dDensity, m_dVel);
	}
	else
	{
		m_dVel = CalVelocity(m_dFlux, m_dWidth, m_dHeight);
		m_dVelDia = CalVelDia(m_dWidth, m_dHeight);
		m_dFluxDia = CalFluxDia(m_dWidth, m_dHeight);
		m_dRe = CalRe(m_dVel, m_dVelDia, m_dMoveVis);			// 计算Re
		m_dFricCoef = CalFricCoef(m_dRough, m_dVelDia, m_dRe);	// 计算摩擦系数
		m_dR = CalR1(m_dVel, m_dFricCoef, m_dDensity, m_dVelDia);	// 计算比摩阻
		m_dFricLoss = CalFricLoss(m_dFricCoef, m_dVelDia, m_dLength, m_dDensity, m_dVel);
	}
}

// 根据流量(m^3/h)和管径(mm)计算圆管流速(m/s)
double HDuctHydraulic::CalVelocity(double dFlux, double dDia) const
{
	if(dDia < ZERO_AIR)
		return 0.0;
	return dFlux / (9.0e-4 * PI * pow(dDia, 2));
}

// 根据流量(m^3/h)、管宽(mm)和管高(mm)计算矩形管流速(m/s)
double HDuctHydraulic::CalVelocity(double dFlux, double dWidth, double dHeight) const
{
	if((dWidth < ZERO_AIR) || (dHeight < ZERO_AIR))
		return 0.0;
	return dFlux / (3.6e-3 * dWidth * dHeight);
}

// 根据流速(m/s)、摩擦系数、密度(kg/m^3)和管径(mm)计算比摩阻
double HDuctHydraulic::CalR1(double dVelocity, double dFricCoef, double dDensity, double dDia) const
{
	if ((dVelocity < 0.0) || (dFricCoef < ZERO_AIR) || (dDensity < ZERO_AIR) 
		|| (dDia < ZERO_AIR))
		return 0.0;
	return dFricCoef * dDensity * pow(dVelocity, 2) / (0.002 * dDia);
}

// 流速(m/s)、流速当量直径(mm)和黏度(m^2/s)计算雷诺数
double HDuctHydraulic::CalRe(double dVelocity, double dVelDia, double dMoveVis) const
{
	if ((dVelocity < 0.0) || (dVelDia < ZERO_AIR) || (dMoveVis < 1.0e-10))
		return ZERO_AIR;
	return dVelocity * 0.001 * dVelDia / dMoveVis;
}

// 根据粗糙度(mm)、流速当量直径(mm)和雷诺数计算摩擦系数, 用柯列勃洛克-怀特公式计算
double HDuctHydraulic::CalFricCoef(double dRough, double dCalDia, double dRe) const
{
	if ((dRough < 0.0) || (dCalDia < ZERO_AIR) || (dRe < ZERO_AIR))
		return ZERO_AIR;
	
	if (dRe < 2000)
		return 64 / dRe;
	
	double dA, dB;	// 系数a, b
	dA = dRough / (3.7 * dCalDia);
	dB = 2.51 / dRe;
	double dNewX, dOldX;	// 用于迭代计算的新旧数据 =1/(FricCoef)^0.5
	dNewX = 5.0;	// 初始值为5
	
	double e = 1.0;		// 两次迭代之差
	double t = 0.0;
	int nFlag = 0;		// 控制最大迭代次数
	while ((e > 0.0001) && (nFlag < 100)) 
	{
		dOldX = dNewX;
		t = dA + dB * dOldX;
		if (t < ZERO_AIR)
			return 0.0;	// 认为无解
		dNewX = dOldX - (dOldX + 2 * log10(t)) * (t) * log(10.0) / ((t) * log(10.0) + 2 * dB);
		e = fabs(dNewX - dOldX);		
		nFlag++;
	}
	if (nFlag >= 100)
		return 0.0;		// 认为无解
	return 1 / pow(dNewX, 2);	
}

// 根据摩擦系数、管径(mm)、管长(m)、密度(kg/m^3)和流速(m/s)计算沿程阻力(Pa)
double HDuctHydraulic::CalFricLoss(double dFricCoef, double dDia, double dLength,
				   double dDensity, double dVelocity) const
{
	if ((dFricCoef < ZERO_AIR) || (dDia < ZERO_AIR) || (dLength < 0)
		|| (dDensity < ZERO_AIR) || (dVelocity < 0))
		return 0.0;
	return dFricCoef * dLength * dDensity * pow(dVelocity, 2) / (0.002 * dDia);
}

// 根据局部阻力系数、密度(kg/m^3)、流速(m/s)和其它局部阻力(Pa)计算局部阻力(Pa)
double HDuctHydraulic::CalLocLoss(double dLocCoef, double dDensity, double dVelocity) const
{
	if ((dLocCoef < 0) || (dDensity < ZERO_AIR) || (dVelocity < 0))
		return 0.0;
	return dLocCoef * dDensity * pow(dVelocity, 2) / 2;
}

// 根据沿程阻力(Pa)、局部阻力(Pa)、起点标高(m)、终点标高(m)、管内密度(kg/m^3)和管外密度(kg/m^3)计算总阻力(Pa)
double HDuctHydraulic::CalTotalLoss(double dFricLoss, double dLocLoss, double dStHeight,
					double dEndHeight, double dDensity, double dOutAirDensity) const
{
	if ((dFricLoss < 0.0) || (dFricLoss < 0.0) 
		|| (dDensity < ZERO_AIR) || (dOutAirDensity < ZERO_AIR))
		return 0.0;	
	return dFricLoss + dLocLoss + (dEndHeight - dStHeight)
		* (dDensity - dOutAirDensity) * 9.8;
}

// 根据管宽(mm)和管高(mm)计算流速当量直径(mm)
double HDuctHydraulic::CalVelDia(double dWidth, double dHeight) const
{
	if ((dWidth < ZERO_AIR) || (dHeight < ZERO_AIR))
		return 0.0;	
	return 2 * dWidth * dHeight / (dWidth + dHeight);
}

// 根据管宽(mm)和管高(mm)计算流量当量直径(mm)
double HDuctHydraulic::CalFluxDia(double dWidth, double dHeight) const
{
	if ((dWidth < ZERO_AIR) || (dHeight < ZERO_AIR))
		return 0.0;	
	return 1.3 * pow(dWidth * dHeight, 0.625) / pow(dWidth + dHeight, 0.25);	
}

// 根据管径(mm)计算圆形断面积(m2)
double HDuctHydraulic::CalSecArea(double dDia) const
{
	if (dDia < ZERO_AIR)
		return 0.0;	
	return PI * pow(0.001 * dDia, 2) / 4;
}

// 根据管宽(mm)和管高(mm)计算矩形断面积(m2)
double HDuctHydraulic::CalSecArea(double dWidth, double dHeight) const
{
	if ((dWidth < ZERO_AIR) || (dHeight < ZERO_AIR))
		return 0.0;	
	return 1.0e-6 * dWidth * dHeight;
}

/*----------------------------------------------------------------------------*
* 功能描述: 用插值法从一维表格中取值
* 备    注:	dAbscissa为表头；共m列数据；dX表示列值;表头升序
-----------------------------------------------------------------------------*/
double CalTable(double dAbscissa[], double *pTable, double dX, int m)
{
	if (dX < dAbscissa[0])
		dX = dAbscissa[0];
	else if (dX > dAbscissa[m - 1])
		dX = dAbscissa[m - 1];
	
	for(int i = 0; i < m; i++)
	{
		if (dAbscissa[i] == dAbscissa[i - 1])	// 与前一个表头相等
			continue;
		
		if (dX == dAbscissa[i])					// 等于某个表头
			return pTable[i];
		else if (dX < dAbscissa[i])				// 小于某个表头，用查值
			return (pTable[i] - pTable[i - 1])*(dX - dAbscissa[i - 1]) 
			/ (dAbscissa[i] - dAbscissa[i - 1]) + pTable[i - 1];
	}
	return 0.0;
}
double HDuctHydraulic::CalDesity(double dTemperature) const
{
	double dDensity[15] = {1.293, 1.270, 1.248, 1.226, 1.205,
		1.185, 1.165, 1.146, 1.128, 1.093,
		1.060, 1.029, 1.000, 0.973, 0.947};
	double dTem[15] = {0, 5, 10, 15, 20,
		25, 30, 35, 40, 50,
		60, 70, 80, 90, 100};
	
	return CalTable(dTem, (double*)dDensity, dTemperature, 15);
}
double HDuctHydraulic::CalMoveVis(double dTemperature) const
{
	double dMovVis[18] = {13.7, 14.7, 15.7, 16.6, 17.6,
		18.6, 19.6, 20.5, 21.7, 22.9,
		23.6, 26.2, 28.5, 30.6, 33.2, 
		35.8, 42.8, 49.9};
	double dTem[18] = {0, 10, 20, 30, 40,
		50, 60, 70, 80, 90,
		100, 120, 140, 160, 180,
		200, 250, 300};
	
	return CalTable(dTem, (double*)dMovVis, dTemperature, 18) * 1.0e-6;
}