📄 tempf2d.h

📁 是一个结构试验设计的加工曲线的跟踪实现没有密码。第一次
💻 H
📖 第 1 页 / 共 5 页
字号:
上一页 1 2 3 45
	double a=0.0,b=0.0;
	for(loop=0; loop<4; loop++)        //初始化
		Pre[loop]=0.0;
	for(loop=0; loop<4; loop++)   
		for(loop1=0; loop1<4; loop1++)
			Kre[loop][loop1]=0.0;

	for(loop=0; loop<TotalLen; loop++)  
		KrG[loop]=0.0;
	for(loop=0; loop<TNode_num; loop++) 
		PrG[loop]=0.0;

	for(loop=0; loop<4; loop++) //把摄氏度化成绝对温度 炉内温度
		AireT[loop]=AireTemp[loop]*Defer+273.15;   

     //叠加截面底边单元边界条件	
	for(loop=0; loop<Tdiv1+divProt*2; loop++)   
	{
		//截面边界上温度  化成绝对温度
		BJTemp=(NTemp[*(TCElemNode_Info+4*loop+0)]+NTemp[*(TCElemNode_Info+4*loop+1)])/2.0+273.15;   
		//线性化后的近似辐射换热系数
		Ratio=Hr[0]*Ste_Bol*(AireT[0]*AireT[0]+BJTemp*BJTemp)*(AireT[0]+BJTemp);
	
		a=fabs(*(Node_Coor+*(TCElemNode_Info+4*loop+1)*2)-*(Node_Coor+*(TCElemNode_Info+4*loop+0)*2))/2.0;

		Kre[0][0]=Kre[1][1]=2.0*a*Ratio/3.0;   
		Kre[0][1]=Kre[1][0]=1.0*a*Ratio/3.0;
		Kre[0][2]=Kre[0][3]=Kre[1][2]=Kre[1][3]=0.0;
		Kre[2][0]=Kre[2][1]=Kre[2][2]=Kre[2][3]=0.0;   
		Kre[3][0]=Kre[3][1]=Kre[3][2]=Kre[3][3]=0.0;

		Pre[0]=Pre[1]=AireTemp[0]*Defer*a*Ratio; 
		Pre[2]=Pre[3]=0.0;

		for(loop1=0; loop1<4; loop1++)	
		{
			iRow=*(TCElemNode_Info+4*loop+loop1);
			for(loop2=0; loop2<=loop1; loop2++)
			{
				iCol=*(TCElemNode_Info+4*loop+loop2);
				if(iRow>=iCol)
					  ibuf=DiagAddress[iRow]-(iRow-iCol);
				else 
					ibuf=DiagAddress[iCol]-(iCol-iRow);
				KrG[ibuf] +=Kre[loop1][loop2];					
			}
			PrG[iRow] +=Pre[loop1];
		}
	}

	//叠加截面右边单元边界条件
	for(loop=Tdiv1+divProt*2-1; loop<(Tdiv1+divProt*2)*(Tdiv2+divProt*2); )
	{
		BJTemp=(NTemp[*(TCElemNode_Info+4*loop+1)]+NTemp[*(TCElemNode_Info+4*loop+3)])/2.0+273.15;
		//线性化后的近似辐射换热系数
		Ratio=Hr[1]*Ste_Bol*(AireT[1]*AireT[1]+BJTemp*BJTemp)*(AireT[1]+BJTemp);

		a=fabs(*(Node_Coor+*(TCElemNode_Info+4*loop+1)*2)-*(Node_Coor+*(TCElemNode_Info+4*loop+0)*2))/2.0;
		b=fabs(*(Node_Coor+*(TCElemNode_Info+4*loop+2)*2+1)-*(Node_Coor+*(TCElemNode_Info+4*loop+0)*2+1))/2.0;

		Kre[1][1]=Kre[3][3]=2.0*b*Ratio/3.0;   
		Kre[1][3]=Kre[3][1]=1.0*b*Ratio/3.0;  
		Kre[0][0]=Kre[0][1]=Kre[0][2]=Kre[0][3]=0.0;		  
		Kre[2][0]=Kre[2][1]=Kre[2][2]=Kre[2][3]=0.0;
		Kre[1][0]=Kre[1][2]=Kre[3][0]=Kre[3][2]=0.0; 

		Pre[1]=Pre[3]=AireTemp[1]*Defer*b*Ratio; 
		Pre[0]=Pre[2]=0.0;
		
		for(loop1=0; loop1<4; loop1++)
		{
			iRow=*(TCElemNode_Info+4*loop+loop1);
			for(loop2=0; loop2<=loop1; loop2++)
			{
				iCol=*(TCElemNode_Info+4*loop+loop2);
				if(iRow>=iCol) 
				   ibuf=DiagAddress[iRow]-(iRow-iCol);
				else 
					ibuf=DiagAddress[iCol]-(iCol-iRow);
				KrG[ibuf] +=Kre[loop1][loop2];
			}
			PrG[iRow] +=Pre[loop1];
		}
		loop+=Tdiv1+divProt*2;
	}

	//叠加截面上边单元边界条件
	for(loop=(Tdiv2+divProt*2-1)*(Tdiv1+divProt*2); loop<(Tdiv1+divProt*2)*(Tdiv2+divProt*2); loop++)
	{
		BJTemp=(NTemp[*(TCElemNode_Info+4*loop+2)]+NTemp[*(TCElemNode_Info+4*loop+3)])/2.0+273.15;
		//线性化后的近似辐射换热系数
		Ratio=Hr[2]*Ste_Bol*(AireT[2]*AireT[2]+BJTemp*BJTemp)*(AireT[2]+BJTemp);

		a=fabs(*(Node_Coor+*(TCElemNode_Info+4*loop+1)*2)-*(Node_Coor+*(TCElemNode_Info+4*loop+0)*2))/2.0;
		b=fabs(*(Node_Coor+*(TCElemNode_Info+4*loop+2)*2+1)-*(Node_Coor+*(TCElemNode_Info+4*loop+0)*2+1))/2.0;

		Kre[2][2]=Kre[3][3]=2.0*a*Ratio/3.0;         //
		Kre[2][3]=Kre[3][2]=1.0*a*Ratio/3.0;
		Kre[0][0]=Kre[0][1]=Kre[0][2]=Kre[0][3]=0.0;
		Kre[1][0]=Kre[1][1]=Kre[1][2]=Kre[1][3]=0.0; 		  
		Kre[2][0]=Kre[2][1]=Kre[3][0]=Kre[3][1]=0.0;
		
		Pre[2]=Pre[3]=AireTemp[2]*Defer*a*Ratio;    //
		Pre[0]=Pre[1]=0.0;

		for(loop1=0; loop1<4; loop1++)
		{
			iRow=*(TCElemNode_Info+4*loop+loop1);
			for(loop2=0; loop2<=loop1; loop2++)
			{
				iCol=*(TCElemNode_Info+4*loop+loop2);
				if(iRow>=iCol) 
				   ibuf=DiagAddress[iRow]-(iRow-iCol);
				else 
					ibuf=DiagAddress[iCol]-(iCol-iRow);
				KrG[ibuf] +=Kre[loop1][loop2];					
			}
			PrG[iRow] +=Pre[loop1];
		}
	}

	//叠加截面左边单元边界条件
	for(loop=0; loop<=(Tdiv2+divProt*2-1)*(Tdiv1+divProt*2); )
	{
		BJTemp=(NTemp[*(TCElemNode_Info+4*loop)]+NTemp[*(TCElemNode_Info+4*loop+2)])/2.0+273.15;
		//线性化后的近似辐射换热系数
		Ratio=Hr[3]*Ste_Bol*(AireT[3]*AireT[3]+BJTemp*BJTemp)*(AireT[3]+BJTemp);

		a=fabs(*(Node_Coor+*(TCElemNode_Info+4*loop+1)*2)-*(Node_Coor+*(TCElemNode_Info+4*loop+0)*2))/2.0;
		b=fabs(*(Node_Coor+*(TCElemNode_Info+4*loop+2)*2+1)-*(Node_Coor+*(TCElemNode_Info+4*loop+0)*2+1))/2.0;

		Kre[0][0]=Kre[2][2]=2.0*b*Ratio/3.0;   
		Kre[0][2]=Kre[2][0]=1.0*b*Ratio/3.0;
		Kre[0][1]=Kre[0][3]=Kre[2][1]=Kre[2][3]=0.0;
		Kre[1][0]=Kre[1][1]=Kre[1][2]=Kre[1][3]=0.0; 		  
		Kre[3][0]=Kre[3][1]=Kre[3][2]=Kre[3][3]=0.0;
		
		Pre[0]=Pre[2]=AireTemp[3]*Defer*b*Ratio; 
		Pre[1]=Pre[3]=0.0;
			
		for(loop1=0; loop1<4; loop1++)
		{
			iRow=*(TCElemNode_Info+4*loop+loop1);
			for(loop2=0; loop2<=loop1; loop2++)
			{
				iCol=*(TCElemNode_Info+4*loop+loop2);
				if(iRow>=iCol)
				   ibuf=DiagAddress[iRow]-(iRow-iCol);
				else 
					ibuf=DiagAddress[iCol]-(iCol-iRow);
				KrG[ibuf] +=Kre[loop1][loop2];
			}
			PrG[iRow] +=Pre[loop1];
		}
		loop+=Tdiv1+divProt*2;
	}
}
//********************************一维矩阵乘法***********************//
void Multiply_1D(double *B,double *C,double *Temp0)
{
	int loop,loop1;
	for(loop=0; loop<TNode_num; loop++)
		B[loop]=0.0;
	for(loop=0; loop<TNode_num; loop++)    //C矩阵乘以Temp0矩阵
	{
		for(loop1=0; loop1<TNode_num; loop1++)
		{
			if(loop>loop1){
				if((loop-loop1)<=(DiagAddress[loop]-DiagAddress[loop-1]-1))
					B[loop] +=C[DiagAddress[loop]-(loop-loop1)]*Temp0[loop1];
			}
			else if(loop<loop1){
				if((loop1-loop)<=(DiagAddress[loop1]-DiagAddress[loop1-1]-1))
					B[loop] +=C[DiagAddress[loop1]-(loop1-loop)]*Temp0[loop1];
			}
			else
				B[loop] +=C[DiagAddress[loop1]]*Temp0[loop1];
		}
	}
}
//**************************高斯消元解方程组(Ax=pz)******************//
void Gauss_1D(double *A,double *pz,double *x)
{
	int loop,loop1,loop2;
	double dbuf=0.0;

	for(loop=0; loop<TNode_num; loop++)
		x[loop]=0.0;
	for(loop=0; loop<TNode_num-1; loop++)
	{
		for(loop1=loop+1; loop1<TNode_num; loop1++)
		{
			if(DiagAddress[loop1]-(loop1-loop)>DiagAddress[loop1-1])
			{
				if(A[DiagAddress[loop]]==0.0 && A[DiagAddress[loop1]-(loop1-loop)]==0.0){
					cout<<"前推出现0/0型"<<endl;
					exit(0);
				}
				dbuf=A[DiagAddress[loop1]-(loop1-loop)]/A[DiagAddress[loop]];
				for(loop2=loop1; loop2<TNode_num; loop2++)
				{
					if(DiagAddress[loop2]-(loop2-loop)>DiagAddress[loop2-1])
						A[DiagAddress[loop2]-(loop2-loop1)] -=A[DiagAddress[loop2]-(loop2-loop)]*dbuf;
				}
			pz[loop1] -=dbuf*pz[loop];		
			}
		}
	}
	pz[TNode_num-1] /=A[DiagAddress[TNode_num-1]];
	if(A[DiagAddress[TNode_num-1]]==0.0 && pz[TNode_num-1]==0.0){
		cout<<"前推出现0/0型"<<endl;
		exit(0);
	}
	
	for(loop=TNode_num-2; loop>=0; loop--)
	{
		for(loop1=loop+1; loop1<TNode_num; loop1++)
		{
			if(DiagAddress[loop1]-(loop1-loop)>DiagAddress[loop1-1])
				pz[loop] -=A[DiagAddress[loop1]-(loop1-loop)]*pz[loop1];
		}
		if(A[DiagAddress[loop]]==0.0 && pz[loop]==0.0){
			cout<<"回代出现0/0型"<<endl;
			exit(0);
		}

		pz[loop] /=A[DiagAddress[loop]];
	}
	for(loop=0;loop<TNode_num; loop++)
		x[loop]=pz[loop];

}
//*******************************提取指定点的温度*************************//
void Cal_Temp_TiQu(int Num_TiQu,double *CorX_TiQu,double *CorY_TiQu,double *Temp,double *Temp_TiQu)
{
	int loop,loop1;
	int PointElemNum=0;            
	double xyytemp1=0.0,xyytemp2=0.0,xyytemp3=0.0,xyytemp4=0.0;

	for(loop=0; loop<Num_TiQu; loop++)
		Temp_TiQu[loop]=0.0;

	for(loop=0; loop<Num_TiQu; loop++)
	{   
		//点的坐标所对应的温度场分析时的单元号
		for(loop1=0;loop1<TConElem_num;loop1++)
		{
			if(CorX_TiQu[loop]>=*(Node_Coor+*(TCElemNode_Info+4*loop1)*2)-0.0000000000000001  &&  CorX_TiQu[loop]<=*(Node_Coor+*(TCElemNode_Info+4*loop1+1)*2)+0.0000000000000001)
			{
				if(CorY_TiQu[loop]>=*(Node_Coor+*(TCElemNode_Info+4*loop1)*2+1)-0.0000000000000001  &&  CorY_TiQu[loop]<=*(Node_Coor+*(TCElemNode_Info+4*loop1+2)*2+1)+0.0000000000000001)
				{
					PointElemNum=loop1;                 
				}			  
			}
		}
		//温度场分析单元的a和b 过书式(6-20) P89
	    xyytemp3=fabs(*(Node_Coor+*(TCElemNode_Info+4*PointElemNum+1)*2)-*(Node_Coor+*(TCElemNode_Info+4*PointElemNum)*2));
		xyytemp4=fabs(*(Node_Coor+*(TCElemNode_Info+4*PointElemNum+2)*2+1)-*(Node_Coor+*(TCElemNode_Info+4*PointElemNum)*2+1));
		//节点相对单元中心的相对坐标x和y 
		xyytemp1=CorX_TiQu[loop]-(*(Node_Coor+*(TCElemNode_Info+4*PointElemNum+1)*2)+*(Node_Coor+*(TCElemNode_Info+4*PointElemNum)*2))/2.0;
		xyytemp2=CorY_TiQu[loop]-(*(Node_Coor+*(TCElemNode_Info+4*PointElemNum+2)*2+1)+*(Node_Coor+*(TCElemNode_Info+4*PointElemNum)*2+1))/2.0;
	    
		Temp_TiQu[loop]=*(Temp+*(TCElemNode_Info+PointElemNum*4))*1.0/4.0*(1-xyytemp1/(xyytemp3/2.0))*(1-xyytemp2/(xyytemp4/2.0));
        Temp_TiQu[loop] +=*(Temp+*(TCElemNode_Info+PointElemNum*4+1))*1.0/4.0*(1+xyytemp1/(xyytemp3/2.0))*(1-xyytemp2/(xyytemp4/2.0));
        Temp_TiQu[loop] +=*(Temp+*(TCElemNode_Info+PointElemNum*4+2))*1.0/4.0*(1-xyytemp1/(xyytemp3/2.0))*(1+xyytemp2/(xyytemp4/2.0));
        Temp_TiQu[loop] +=*(Temp+*(TCElemNode_Info+PointElemNum*4+3))*1.0/4.0*(1+xyytemp1/(xyytemp3/2.0))*(1+xyytemp2/(xyytemp4/2.0));
	}
}
//****************判断点是否处于明火中,已经爆裂****************//
void Cal_Spall_Ind(double CorX,double CorY,int &Spall_Ind)
{	
	Spall_Ind=0;//Spall_Ind=1:说明该点已经处于爆裂区中
	if(Fire0==0 && Fire1==1 && Fire2==0  && Fire3==0)//右
	{
		if(Sectwid-CorX<Spall_thick) 
			Spall_Ind=1;
	}
	else if(Fire0==0 && Fire1==1 && Fire2==1  && Fire3==0)//右 上
	{
		if(Sectwid-CorX<Spall_thick || Sectlen-CorY<Spall_thick) 
			Spall_Ind=1;
	}
	else if(Fire0==1 && Fire1==1 && Fire2==1  && Fire3==0)//右 上 下
	{
		if(CorY<Spall_thick || Sectlen-CorY<Spall_thick || Sectwid-CorX<Spall_thick) 
			Spall_Ind=1;
	}
	else if(Fire0==1 && Fire1==1 && Fire2==1  && Fire3==1)
	{
		if(CorX<Spall_thick || CorY<Spall_thick || Sectwid-CorX<Spall_thick || Sectlen-CorY<Spall_thick)
			Spall_Ind=1;
	}
	else if(Fire0==0 && Fire1==1 && Fire2==0  && Fire3==1)//左 右
	{  
		if(CorX<Spall_thick || Sectwid-CorX<Spall_thick) 
			Spall_Ind=1;
	}
	else{
		cout<<"受火形式(Face)不对,请检查";
		exit(0);
	}
}
//****************混凝土的比热容******************************//
double CalcuCSpecialHeat(double Temp,int Index)
{
	double SpecialHeat=0.0;
	double SpecialHeatPeakT=0.0;               //CWaterPeakT时最大热容

	if(CWaterContent>=0.001 && Water_Method==1)   //Lie-韩林海考虑含水量方法 
	{	
		if(Index==1)      //过镇海建议的公式
		{
			if(Temp<20.0)
				SpecialHeat=913.2*(1.0-CWaterContent)*PConcrete+Cw*PWater*CWaterContent;
			else if(Temp>=20.0 && Temp<100.0)
				SpecialHeat=(900.0+80.0*Temp/120.0-4.0*Temp*Temp/120.0/120.0)*PConcrete*(1.0-CWaterContent)+Cw*PWater*CWaterContent;
			else if(Temp>=100.0 && Temp<=1200.0)
				SpecialHeat=(900.0+80.0*Temp/120.0-4.0*Temp*Temp/120.0/120.0)*PConcrete;
			else
				SpecialHeat=1300.0*PConcrete;
		}
		else if(Index==2) //T.T Lie建议的公式(1991年)
		{
			if(Temp<0.0)
				SpecialHeat=1.7*1E6*(1.0-CWaterContent)+Cw*PWater*CWaterContent;
			else if(Temp>=0.0 && Temp<100.0)
				SpecialHeat=(0.005*Temp+1.7)*1E6*(1.0-CWaterContent)+Cw*PWater*CWaterContent;
			else if(Temp>=100.0 && Temp<=200.0)
				SpecialHeat=(0.005*Temp+1.7)*1E6;
			else if(Temp>200.0 && Temp<=400.0)
				SpecialHeat=2.7*1E6;
			else if(Temp>400.0 && Temp<=500.0)
				SpecialHeat=(0.013*Temp-2.5)*1E6;
			else if(Temp>500.0 && Temp<=600.0)
				SpecialHeat=(-0.013*Temp+10.5)*1E6;
			else
				SpecialHeat=2.7*1E6;
		}
		else if(Index==3)  //T.T Lie建议的公式(1994年)
		{
			if(Temp<100.0)
				SpecialHeat=2.566*1.0E6*(1.0-CWaterContent)+Cw*PWater*CWaterContent; 
			else if(Temp>=100.0 && Temp<=400.0)
				SpecialHeat=2.566*1.0E6;
			else if(Temp>400.0 && Temp<=410.0)
				SpecialHeat=(0.1765*Temp-68.034)*1.0E6;
			else if(Temp>410.0 && Temp<=445.0)
				SpecialHeat=(-0.05043*Temp+25.00671)*1.0E6;
			else if(Temp>445.0 && Temp<=500.0)
				SpecialHeat=2.566*1.0E6;
			else if(Temp>500.0 && Temp<=635.0)
				SpecialHeat=(0.01603*Temp-5.44881)*1.0E6;
			else if(Temp>635.0 && Temp<=715.0)
				SpecialHeat=(0.16635*Temp-100.90225)*1.0E6;
			else if(Temp>715.0 && Temp<=785.0)
				SpecialHeat=(-0.22103*Temp+176.07343)*1.0E6;
			else if(Temp>785.0)
				SpecialHeat=2.56
上一页 1 2 3 45
⌨️ 快捷键说明

复制代码 Ctrl + C
搜索代码 Ctrl + F
全屏模式 F11
切换主题 Ctrl + Shift + D
显示快捷键 ?
增大字号 Ctrl + =
减小字号 Ctrl + -