normmat.cpp

Finite element program for mechanical problem. It can solve various problem in solid problem

#include <stdlib.h>  
#include <stdio.h>


#include <math.h>


#include "normmat.h"


#include "mechmat.h"


#include "intpoints.h"


#include "global.h"


#include "matrix.h"








normmat::normmat (void)
{
  tcn=NULL;
  
  zero=1.0e-10;

  iter=1;

  tcn = new double [4];
  ki = new double [4];
  a1sp = new double [4];
  b1sp = new double [4];
  allocm (2,2,kn);
  allocm (4,4,kp);
  allocv (2,fkn);
  allocv (4,fkp);
  m = new double [4];
  mp = new double [4];
  mpp = new double [4];
  tcv = new double [4];
  tcn = new double [4];
  epscsp = new double [4];
  epss1sp = new double [4];			
  a1sp = new double [4];
  a2spp = new double [4];
  a2sp = new double [4];
  b1sp = new double [4];
  b2spp = new double [4];
  b2sp = new double [4];
  mr = new double [4];
  epscr = new double [4];
  epss1r = new double [4];
  epss1spp = new double [4];
  xr = new double [4];
  xsp = new double [4];
  ir1 = new double [4];
  ims1 = new double [4];
  ims2 = new double [4];
  ir2 = new double [4];
  ir3 = new double [4];
  ims3 = new double [4];
  ic1 = new double [4];
  ic2 = new double [4];
  icp = new double [4];
  c1 = new double [4];
  c2 = new double [4];
  c1p = new double [4];
  c2p = new double [4];
  kih = new double [4];
  kid = new double [4];
  ki = new double [4];

}





normmat::~normmat (void)


{


  delete [] tcn;


}





void normmat::read (XFILE *in)
{
  xfscanf (in,"%le %le %le %le %le %le  %le  %le %le %le %le %le %le %le %le",
	  &fck,&fy,&ft,&epsc1,&epscm,&epss1lim,&ash,&asd,&ch,&cd,&e0c,&h,&b,&fp,&fn);
  
  
  compute_fc ();
  compute_ecnom ();
  compute_epsuk ();
  compute_e0cn ();
  compute_k ();
  compute_a ();
  compute_d ();
  
}







//pevnost betonu v tlaku uvazovana pro vypocet


void normmat::compute_fc ()
{
  fc=-8000.0 + fck;
}





 //modul pruznosti betonu v tlaku prislusny 0,6 fc v kPa


void normmat::compute_ecnom ()
{
  ecnom = (1.6321156877891085*fc)/epsc1;
}





  //pretvoreni oceli pro stanoveni modulu pruznosti e0cn


void normmat::compute_epsuk ()
{
  epsuk = 0.05;
}


 


//modul pruznosti odpovidajici zplastizovani oceli v kPa
void normmat::compute_e0cn ()
{
  e0cn = (ft - fy)/(-(fy/e0c) + epsuk);
}


 


//soucinitel pro definovani pracovniho diagramu betonu
void normmat::compute_k ()
{
   k = (11*ecnom*epsc1)/(10.*fc);
}


 


//plocha prurezu


void normmat::compute_a ()
{
  a = b*h;
}


 


//ucinna vyska prurezu


void normmat::compute_d ()
{
  d = h - cd;
}





  
//svisla koncova sila ve stycniku i


void normmat::compute_z1 (long ipp,double lp)
{
  /*
	double wi,wfi,wj,wfj;

	wi = Mm->ip[ipp].strain[1]; 
	wfi = Mm->ip[ipp].strain[2];
	wj = Mm->ip[ipp].strain[4]; 
	wfj = Mm->ip[ipp].strain[5];

	z1 = ecnom*b*(h^3)/12.*((12.*wi)/pow(lp,3)- 
      (12.*wj)/pow(lp,3) - (6.*wfi)/pow(lp,2) - 
      (6.*wfj)/pow(lp,2))-(fp*lp)/2.;
      */
	
  z1 = Mm->ip[ipp].stress[1];
	//z1 = 0.;
} 


  
//momentova koncova sila ve stycniku i


void normmat::compute_z2 (long ipp,double lp)
{
/*
	double wi,wfi,wj,wfj;

	wi = Mm->ip[ipp].strain[1]; 
	wfi = Mm->ip[ipp].strain[2];
	wj = Mm->ip[ipp].strain[4]; 
	wfj = Mm->ip[ipp].strain[5];

	z2 = ecnom*b*(h^3)/12.*((-6.*wi)/pow(lp,2)+ 
      (6.*wj)/pow(lp,2) + (4.*wfi)/lp + 
      (2.*wfj)/lp)+(fp*pow(lp,2))/12.;
*/
	
  z2 = Mm->ip[ipp].stress[2];
	//z2=1000.;
} 





  
//vodorovna koncova sila ve stycniku i


void normmat::compute_z3p (long ipp,double lp)
{
	

	z3p = Mm->ip[ipp].stress[0];

	
} 






//vodorovna koncova sila ve stycniku i


void normmat::compute_z3 ()
{
	
/*
	double ui,uj;

	ui = Mm->ip[ipp].strain[0]; 
	uj = Mm->ip[ipp].strain[3]; 
	
	z3 = ecnom*a*(ui/lp - 
		uj/lp)+(fn*lp)/2.;
*/

	

  if (z3p <= 0.0)
    
    { z3 = 0.1; }
  
  else
    
    { z3 = z3p; }
  
  //z3 = Mm->ip[ipp].stress[0];
  //z3=2000;
} 





  
//vyztuz pro tahovou oblast od ucinku ohyboveho momentu as1


double normmat::compute_as1 (double lp)
{
  if (-(fp*pow(lp,2))/8. - (lp*z1)/2. - z2 <= 0){
    as1=asd;
  }
  else {
    as1=ash;
  }
  return as1;
}
 


//vyztuz pro tlakovou oblast od ucinku ohyboveho momentu as2


double normmat::compute_as2 (double lp)
{
  if (-(fp*pow(lp,2))/8. - (lp*z1)/2. - z2 <= 0){
   as2=ash;
 }
  else {
    as2=asd;
  }
  
  return as2;
}





 


//ohybovy moment na dannem elementu v polovine intervalu


void normmat::compute_m (double lp)
{

  m[0]=-(fp*pow(lp, 2))/128. - (lp*z1)/8. - z2;


  m[1]=(-9*fp*pow(lp, 2))/128. - (3*lp*z1)/8. - z2;


  m[2]=(-25*fp*pow(lp, 2))/128. - (5*lp*z1)/8. - z2;


  m[3]=(-49*fp*pow(lp, 2))/128. - (7*lp*z1)/8. - z2;


   


 }


 


//absolutni hodnota ohyboveho momentu na dannem elementu v polovine intervalu


void normmat::compute_mpp (double lp)
{
	if(m[0]==0.0){mpp[0]= 0.01;}
	else {mpp[0]= fabs(-(fp*pow(lp, 2))/128. - (lp*z1)/8. - z2);}

	if(m[1]==0.0){mpp[1]= 0.01;}
	else {mpp[1]= fabs((-9*fp*pow(lp, 2))/128. - (3*lp*z1)/8. - z2);}

	if(m[2]==0.0){mpp[2]= 0.01;}
	else {mpp[2]= fabs((-25*fp*pow(lp, 2))/128. - (5*lp*z1)/8. - z2);}

	if(m[3]==0.0){mpp[3]= 0.01;}
	else {mpp[3]= fabs((-49*fp*pow(lp, 2))/128. - (7*lp*z1)/8. - z2);}

  
}


 


//normalova sila na dannem elementu v polovine intervalu


double normmat::compute_n0 (double lp)


{


  double n0;


  


  


  if (-z3 + fn*lp/2 <= -0.1 ){


    n0=-z3 + fn*lp/2;}


  


  else{


    n0=-0.1;


  }


  


  return n0;


}






 
//pocet iteraci


long normmat::compute_iter ()


{


  long iter;	


  


  //iter = je rovna cislu iterace, ktera prave probiha pro danny zatezovaci stav, nulta iterace bude prvni vypocet, ktery bude stanoven linearne;


  iter = 1;


  return iter;


}






 
//funkce ktera rusi vliv betonove casti prurezu


void normmat::compute_tcv ()
{	
  long i;

  for (i=0;i<4;i++){
    
    
    //tcn[i]=Mm->ip.other[];
    
    
  }
  
  
  
  
  
  //compute_tcni (ipp,lp,tcn);


  


  //tcn = je funkce, ktera bude prenesena z predchoziho iterace a je rovna fci tcn, pro nultou iteraci je rovna (0,0,0,0);


  


  /*


  tcv[0]=tcn[0];


  tcv[1]=tcn[1];


  tcv[2]=tcn[2];


  tcv[3]=tcn[3];


  */





  tcv[0]=0.0;


  tcv[1]=0.0;


  tcv[2]=0.0;


  tcv[3]=0.0;


  


}


 


//VLIV NORMALOVE SILY


 


//pretvoreni od ucinku normalove sily


double normmat::compute_epsn ()


{


	double epsn;





	epsn =(-2*a*fc*epsc1 + 3*a*fc*k*epsc1 + as1*e0c*k*pow(epsc1, 2) + 


      as2*e0c*k*pow(epsc1, 2) - 


      sqrt(pow(epsc1, 


              2)*(-8*a*fc*(-1 + k)*k*n0 + 


                pow(a*fc*(-2 + 3*k) + (as1 + as2)*e0c*k*epsc1, 2))))/(4.*a*


      fc*(-1 + k));





	return epsn;


}


 


//maximalni tlakova normalova sila, kterou lze zatizit prurez


double normmat::compute_nmax ()


{	


	double nmax;





	nmax =as1*e0c*epsc1 + 


  as2*e0c*epsc1 + (a*fc*(-2*(-1 + k)*epsc1 + (-2 + 3*k)*epsc1))/(k*epsc1);





  return nmax;


}


 


//prvni clen pro linearizaci


double normmat::compute_an1 ()


{	


	double an1;





	an1 = sqrt(pow(epsc1, 


            2)*(-8*a*fc*(-1 + k)*k*n0 + 


              pow(a*fc*(-2 + 3*k) + (as1 + as2)*e0c*k*epsc1, 2)))/(k*


          pow(epsc1, 2));





  return an1;


}


 


//druhy clen pro linearizaci


double normmat::compute_an2 ()


{	


	double an2;





	an2 = n0 + (sqrt(pow(epsc1, 


                2)*(-8*a*fc*(-1 + k)*k*n0 + 


                  pow(a*fc*(-2 + 3*k) + (as1 + as2)*e0c*k*epsc1, 


                      2)))*(epsc1*(a*fc*(2 - 3*k) - (as1 + as2)*e0c*k*epsc1) +


               sqrt(pow(epsc1, 


                      2)*(-8*a*fc*(-1 + k)*k*n0 + 


                        pow(a*fc*(-2 + 3*k) + (as1 + as2)*e0c*k*epsc1, 


                            2)))))/(4.*a*fc*(-1 + k)*k*pow(epsc1, 2));





  return an2;


}


 


//prace na betonove casti prurezu pro prubeh normalove sily (-z3 + fn x) pres element delky lp 


double normmat::compute_in1 (double lp)


{


	double in1;	


	


	


  if (fn == 0. && z3>= 0.0){





    in1=1/(32.*pow(a, 2)*pow(fc, 2)*pow(-1 + k, 2)*k*pow(epsc1, 2))*(lp*


      pow(a*fc*(-2 + 3*k)*epsc1 + (as1 + as2)*e0c*k*pow(epsc1, 2) - 


          sqrt(pow(epsc1, 


                  2)*(8*a*fc*(-1 + k)*k*z3 + 


                    pow(a*fc*(-2 + 3*k) + (as1 + as2)*e0c*k*epsc1, 2))), 


        2)*(a*fc*(-2 + 3*k)*epsc1 - (as1 + as2)*e0c*k*pow(epsc1, 2) + 


          sqrt(pow(epsc1, 


                  2)*(8*a*fc*(-1 + k)*k*z3 + 


                    pow(a*fc*(-2 + 3*k) + (as1 + as2)*e0c*k*epsc1, 2)))));





  }


  else if (fn != 0. && z3>= 0.0){





    in1=-(1/(960.*pow(a, 3)*pow(fc, 3)*fn*pow(-1 + k, 3)*


            pow(k, 2)))*(-pow(a, 4)*pow(fc, 4)*


        pow(2 - 3*k, 


          4)*(sqrt(pow(epsc1, 


                    2)*(-8*a*fc*(-1 + k)*k*(fn*lp - z3) + 


                      pow(a*fc*(-2 + 3*k) + (as1 + as2)*e0c*k*epsc1, 2))) - 


            sqrt(pow(epsc1, 


                    2)*(8*a*fc*(-1 + k)*k*z3 + 


                      pow(a*fc*(-2 + 3*k) + (as1 + as2)*e0c*k*epsc1, 


                          2)))) + 


      9*pow(as1 + as2, 4)*pow(e0c, 4)*pow(k, 4)*