t255.cpp

此程序为研究异形柱最不利荷载角的C++程序

sita=sita1+dsita;
cossita=cos(sita);sinsita=sin(sita);
a=0;w=1;aa=0;
dr=0.1;
do{
eicm=0;eism=0;

r=r1+aa*w*dr;

np=0;nc=0;
ns1=0;nsc=0;

for(i=0;i<z0;i++)
{ri=r-danc[i][0]*cossita-danc[i][1]*sinsita;
 eic[i]=ri*fai;
 if(eic[i]<eicm)
	 eicm=eic[i];
 { 		  
       nc=nc+strc1(eic[i],fc)*ds;
    //  mcx=mcx+strc1(eic[i],fc)*ds*(danc[i][1]-y0);
    //  mcy=mcy+strc1(eic[i],fc)*ds*(danc[i][0]-x0);
 }
}
  for(i=z0;i<ik;i++)
         {ri=r-danc[i][0]*cossita-danc[i][1]*sinsita;
          eic[i]=ri*fai;
         
          nc=nc+strc(eic[i],z,kc,fc)*ds;
   //       mcx=mcx+strc(eic[i],z,kc,fc)*ds*(danc[i][1]-y0);
   //       mcy=mcy+strc(eic[i],z,kc,fc)*ds*(danc[i][0]-x0);
}


for(i=0;i<10;i++)
{
     ri=r-dans[i][0]*cossita-dans[i][1]*sinsita;
     eis[i]=fai*ri;
     if(eis[i]>eism)
		 eism=eis[i];
         ns1=ns1+strs(eis[i],fy)*as1;
         nsc=nsc+strc(eis[i],z,kc,fc)*as1;
  //       msx1=msx1+strs(eis[i],fy)*as1*(dans[i][1]-y0);
  //       msy1=msy1+strs(eis[i],fy)*as1*(dans[i][0]-x0);
 //        mscx=mscx+strc(eis[i],z,fc)*as1*(dans[i][1]-y0);
 //        mscy=mscy+strc(eis[i],z,fc)*as1*(dans[i][0]-x0);
}
  
for(i=10;i<numdans;i++)
{
     ri=r-dans[i][0]*cossita-dans[i][1]*sinsita;
     eis[i]=fai*ri;
     if(eis[i]>eism)
		 eism=eis[i];
         ns1=ns1+strs(eis[i],fy)*as3;
         nsc=nsc+strc(eis[i],z,kc,fc)*as3;
  //       msx1=msx1+strs(eis[i],fy)*as3*(dans[i][1]-y0);
 //        msy1=msy1+strs(eis[i],fy)*as3*(dans[i][0]-x0);
 //        mscx=mscx+strc(eis[i],z,fc)*as3*(dans[i][1]-y0);
//         mscy=mscy+strc(eis[i],z,fc)*as3*(dans[i][0]-x0);
}
np=nc+ns1-nsc;
//msx=msx1-mscx;
//msy=msy1-mscy;
a++;
aa=(a+1);     //原为(a+1)/2
if((np+n)>0)
{ w=-1;
}else w=1;    //w=-w

if(fabs((np+n)/n)>0.4)
{
 dr=1;
} else
dr=0.1;
 
nnnn++;nnn++;
}while((fabs((np+n)/n))>=0.01);

r1=r;

//k++;
//u=(k+1)/2;
//q=-q;
 
mcx=0;mcy=0;
msx=0;msy=0;mscx=0;mscy=0;
msx1=0;msy1=0;
mx=0;my=0;
for(i=0;i<z0;i++)
{
// ri=r-danc[i][0]*cossita-danc[i][1]*sinsita;
// eic=ri*fai;
  
       mcx=mcx+strc1(eic[i],fc)*ds*(danc[i][1]-y0);
       mcy=mcy+strc1(eic[i],fc)*ds*(danc[i][0]-x0);

}
  for(i=z0;i<ik;i++)
        {
//	     ri=r-danc[i][0]*cossita-danc[i][1]*sinsita;
  //      eic=ri*fai;
         
        
          mcx=mcx+strc(eic[i],z,kc,fc)*ds*(danc[i][1]-y0);
          mcy=mcy+strc(eic[i],z,kc,fc)*ds*(danc[i][0]-x0);
}


for(i=0;i<10;i++)
{
//     ri=r-dans[i][0]*cossita-dans[i][1]*sinsita;
//     eis[i]=fai*ri;

         msx1=msx1+strs(eis[i],fy)*as1*(dans[i][1]-y0);
         msy1=msy1+strs(eis[i],fy)*as1*(dans[i][0]-x0);
         mscx=mscx+strc(eis[i],z,kc,fc)*as1*(dans[i][1]-y0);
         mscy=mscy+strc(eis[i],z,kc,fc)*as1*(dans[i][0]-x0);
}
  
for(i=10;i<numdans;i++)
{
//     ri=r-dans[i][0]*cossita-dans[i][1]*sinsita;
//    eis[i]=fai*ri;


         msx1=msx1+strs(eis[i],fy)*as3*(dans[i][1]-y0);
         msy1=msy1+strs(eis[i],fy)*as3*(dans[i][0]-x0);
         mscx=mscx+strc(eis[i],z,kc,fc)*as3*(dans[i][1]-y0);
         mscy=mscy+strc(eis[i],z,kc,fc)*as3*(dans[i][0]-x0);
}     
 
   
 mx=mcx+msx1-mscx;
     my=mcy+msy1-mscy;
  shi=atan2((-mx),(-my));
  wucha=shi-alf;
	  dsita=-0.5*wucha;
 cout<<dsita<<"  "<<sita<<"    "<<shi<<"   "<<nnnn<<"   "; 
 sita1=sita;
 }while(fabs(wucha)>=0.01);                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                 



m=sqrt(mx*mx+my*my)/1000000.0;
fai1=1000*fai;
//**************************
//求混凝土达到0.0033时的曲率
for(i=0;i<ik;i++)
{
//ri=r-danc[i][0]*cossita-danc[i][1]*sinsita;
//eic=ri*fai;
if(eic[i]<-0.0033) 
{
 fac[kkk]=fai;
 ang2[kkk]=sita1;
 dist2[kkk]=r;
// cout<<eic<<"   "<<danc[i][0]<<"   "<<danc[i][1]<<"  ";
 kkk++;
 break;
}
}



//***********************
//求钢筋达到屈服时的曲率
for(i=0;i<numdans;i++)
{
ri=r-dans[i][0]*cossita-dans[i][1]*sinsita;
eis[i]=ri*fai;
if(eis[i]>fy/ES)
{fag[kk]=fai;
ang1[kk]=sita;
dist1[kk]=r;
//cout<<eis[i]<<"  "<<dans[i][0]<<"  "<<dans[i][1]<<"  ";

kk++;
break;
}
}
//**********************
//计算纵筋压曲时的end值
for(i=0;i<numdans;i++)
{
ri=r-dans[i][0]*cossita-dans[i][1]*sinsita;
eis[i]=ri*fai;
if(eis[i]<=-esb) 
{end=1;
cout<<end<<"   "<<eis[i]<<"  "<<dans[i][0]<<"   "<<dans[i][1]<<"  ";
break;}
if(eis[i]>=0.01) 
{end=2;
cout<<end<<"   "<<eis[i]<<"  "<<dans[i][0]<<"   "<<dans[i][1]<<"  ";
break;}
}
for(i=0;i<z0;i++)
{
if(eic[i]<=-0.0033) 
{end=3;
break;
}
}

    





    tt++;
    mfai[tt][0]=m;
    mfai[tt][1]=fai1;
    mfai[tt][2]=eicm;
    mfai[tt][3]=eism;
    mfai[tt][4]=r;
    mfai[tt][5]=sita1;	
   
/*	for(i=1;i<=tt;i++)
	{      
    if(mfai[i][0]>mmax[km])
    mmax[km]=mfai[i][0]; 
	}*/
   if(m>mmax[km])
   {mmax[km]=m; 
   mmax1[km]=fabs(mx/1000);
   shi1[km]=atan2(fabs(my),fabs(mx));
   }
    } while((mfai[tt][0]>=0.95*mmax[km])&&(end==0));
	//*****************************
	//跳出循环～～～～～～～～～～·

    angle=sita1;
	dist=r;
	anga=ang1[1];
	dista=dist1[1];
	angb=ang2[1];
	distb=dist2[1];
	//释放堆内存空间
	delete[]ang1;
	delete[]ang2;
	delete[]dist1;
	delete[]dist2;
	
	

     
	    num2=num3=num4=0;

	
	//*******************************
	//计算混凝土压碎的单元，1表示压碎
	for(i=0;i<ik;i++)
	{
//	ri=r-danc[i][0]*cossita-danc[i][1]*sinsita;
//	eic=ri*fai;
	if(eic[i]<-0.0033)
		{
		num2++;
	 conc[i]=1;
		}
	 else
		 conc[i]=0;
	}

	//*****************************************
	//计算达到极限曲率时钢筋屈服的单元，1表示屈服
	for(i=0;i<numdans;i++)
	{
		ri=r-dans[i][0]*cossita-dans[i][1]*sinsita;
		eis[i]=ri*fai;
		//求解极限状态时各钢筋单元的应变
		strs3[i]=eis[i];
		if(eis[i]>fy/ES)
		{
			num3++;
			stel[i]=1;
		}
			else 
				stel[i]=0;
	//**********************************************
	//计算达到极限曲率时纵筋压曲的钢筋单元，1表示压曲
			if(eis[i]<-esb)
		{
				num4++;
			stelp[i]=1;
		}
			else
				stelp[i]=0;
	}
    //***************************
	// 求解屈服曲率
   if(fac[1]>=fag[1])
  {facg=fag[1]*1000;
   quf='g';
  }
  else
  {facg=fac[1]*1000;
  quf='h';
  }
   miu=fai1/facg;
   //***********************************
   //求解钢筋开始屈服时各钢筋单元的应变
       for(i=0;i<numdans;i++)
		{ri=dista-dans[i][0]*cos(anga)-dans[i][1]*sin(anga);
          strs1[i]=fag[1]*ri;
        }
   //***********************************
   //求解混凝土压应变达到0.0033时各钢筋单元的应变
      for(i=0;i<numdans;i++)
		{ri=distb-dans[i][0]*cos(angb)-dans[i][1]*sin(angb);
          strs2[i]=fac[1]*ri;
		 }
    alf=alf+0.4;
     km++; 
 }while(alf<1.0*pi);


	  ofstream out,out1,out2,out3;
out.open("second-270.dat");
out1.open("mmax.dat");
out2.open("mmax1.dat");
out3.open("shi1.dat");
for(i=0;i<km;i++)
{out1<<mmax[i]<<endl;
out2<<mmax1[i]<<endl;
out3<<shi1[i]<<endl;
}
out<<setw(10)<<"轴比压n＝" <<zhou<<" 轴力 "<<n<<endl;
out<<"esb="<<esb<<endl;
out<<"配箍率p＝" <<p<<endl;
out<<"配筋率ρ="<<rou<<endl;
out<<"荷载角alf="<<alf<<endl;
out<<"end="<<end<<endl;
out<<num2<<" "<<num3<<" "<<num4<<endl;
out<<"屈服曲率："<<facg<<"  "<<"极限曲率："<<fai1<<"\n";
out<<"最大承载力："<<mmax[km-1]<<"截面延性："<<miu<<"\n";
out<<"最先屈服单元是："<<quf<<endl;
out<<setw(15)<<'M'<<setw(15)<<"fai"<<endl;
//****************************************
for(i=1;i<=tt;i++)
out<<setw(15)<<mfai[i][0]<<setw(15)<<mfai[i][1]<<endl;
out<<setw(15)<<"eicm"<<setw(15)<<"eism"<<endl;
for(i=1;i<=tt;i++) 
out<<setw(15)<<mfai[i][2]<<setw(15)<<mfai[i][3]<<endl;
out<<setw(15)<<'R'<<setw(15)<<"sita"<<endl;
for(i=1;i<=tt;i++)
out<<setw(15)<<mfai[i][4]<<setw(15)<<mfai[i][5]<<endl;
 for(i=0;i<ik;i++)
	{ out<<conc[i]<<" ";
    if(i%10==0)
		out<<endl;}
 out<<endl;
     for(i=0;i<numdans;i++)
	 out<<stel[i]<<"    "<<stelp[i]<<endl;
//******************************
out<<"混凝土应变达到0.0033时的曲率为:"<<endl;
out<<fac[1]*1000<<endl;
out<<"R="<<distb<<endl;
out<<"中和轴法线角度："<<angb<<endl;
//**********************************
out<<"钢筋屈服时截面的曲率为:"<<endl;
out<<fag[1]*1000<<endl;
out<<"R="<<dista<<endl;
out<<"中和轴法线角度："<<anga<<endl;

//******************************
out<<"钢筋屈服时各钢筋单元的应变为:" <<endl;
for(i=0;i<numdans;i++)
out<<strs1[i]<<endl;
out<<"第一批混凝土单元屈服时各钢筋单元的应变为:" <<endl;
for(i=0;i<numdans;i++)
out<<strs2[i]<<endl;
out<<"极限破坏时各钢筋单元的应变为:" <<endl;
for(i=0;i<numdans;i++)
out<<strs3[i]<<endl;
 
out<<"ts"<<ts<<"x0"<<x0<<"y0"<<y0<<"n"<<nnn<<endl;
out1.close();
out.close();
}