normmat.cpp

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

    in2=0.;





  }


  return in2;


 }


  


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


double normmat::compute_kinh ()


{


	double kinh;





	kinh = in1 + in2;


	


  return kinh;


 }


 


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


double normmat::compute_kind (double lp)


{





  double kind;


  


  


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


    





    kind=1.0/(16.*pow(a, 2)*an1*pow(fc, 2)*pow(-1 + k, 2))*(as1*e0c*lp*(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))))*(an1*(as1 + as2)*e0c*k*pow(epsc1, 2) + 


																						  a*fc*(4*an2*(-1 + k) + an1*(-2 + 3*k)*epsc1) - an1*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){


    





    kind=1/(96.*pow(a, 3)*an1*pow(fc, 3)*fn*pow(-1 + k, 3)*k)*(as1*


      e0c*(an1*pow(as1 + as2, 3)*pow(e0c, 3)*pow(k, 3)*


            pow(epsc1, 


              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)))) + 


          pow(a, 3)*


            pow(fc, 


              3)*(-2 + 3*k)*(2*an2*(-1 + k) + 


                an1*(-2 + 3*k)*epsc1)*(12*fn*(-1 + k)*k*lp*


                  epsc1 + (-2 + 


                      3*k)*(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))))) + 


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


            pow(epsc1, 


              2)*(2*an2*(as1 + as2)*


                  e0c*(-1 + 


                      k)*(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)))) + 


                an1*(3*as1*e0c*


                        epsc1*(4*fn*(-1 + k)*k*lp*


                              epsc1 + (-2 + 


                                  3*k)*(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))))) + 


                      3*as2*e0c*


                        epsc1*(4*fn*(-1 + k)*k*lp*


                              epsc1 + (-2 + 


                                  3*k)*(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))))) - 


                      8*(-1 + k)*(fn*lp*


                            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)))


+z3*(-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))))))) + 


          pow(a, 2)*pow(fc, 2)*


            k*(an1*epsc1*(-24*pow(fn, 2)*pow(-1 + k, 2)*k*pow(lp, 2)*


                        epsc1 + (-2 + 3*k)*(8*(-1 + k)*z3 + 


                            3*(as1 + as2)*e0c*(-2 + 3*k)*


                              epsc1)*(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)))) + 


                      8*fn*(-1 + k)*


                        lp*(3*pow(k, 2)*


                              epsc1*(2*z3 + 3*(as1 + as2)*e0c*epsc1) + 


                            2*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))) - 


                            3*k*(2*z3*epsc1 + 2*as1*e0c*pow(epsc1, 2) + 


                                  2*as2*e0c*pow(epsc1, 2) + 


                                  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)))))) + 


                4*an2*(-1 + 


                      k)*(as1*e0c*


                        epsc1*(6*fn*(-1 + k)*k*lp*


                              epsc1 + (-2 + 


                                  3*k)*(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))))) + 


                      as2*e0c*


                        epsc1*(6*fn*(-1 + k)*k*lp*


                              epsc1 + (-2 + 


                                  3*k)*(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))))) - 


                      4*(-1 + k)*(fn*lp*


                            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)))


+z3*(-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)))))))));


					    





  }





  else{





    kind=1.0;





  }


  return kind;





}


 


//pomer prace vnitrnich sil a pomocne prace vnitrnich sil na dannem elementu od ucinku normalove sily


double normmat::compute_kin ()


{


	double kin;





	kin = kinh/kind;


	


  return kin;


 }


  


//TLACENA VYZTUZ + NORMALOVA SILA + x = h


 


//pretvoreni betonu, kdy tlacena plocha x=h


double normmat::compute_epschr1 ()


{	


  double epsc,disk,pom,epschr1;	


  


  disk =pow(6*b*fc*pow(h, 2)*epsc1 - 9*b*fc*pow(h, 2)*k*epsc1 + 


        6*as2*ch*e0c*k*pow(epsc1, 2) + 6*as1*d*e0c*k*pow(epsc1, 2) - 


        6*as1*e0c*h*k*pow(epsc1, 2) - 6*as2*e0c*h*k*pow(epsc1, 2) - 


        8*b*fc*pow(h, 2)*epsn + 8*b*fc*pow(h, 2)*k*epsn, 2) - 


  4*b*fc*pow(h, 


      2)*(-4 + 4*k)*(6*h*k*n0*pow(epsc1, 2) + 6*b*fc*pow(h, 2)*epsc1*epsn - 


        9*b*fc*pow(h, 2)*k*epsc1*epsn + 6*as2*ch*e0c*k*pow(epsc1, 2)*epsn + 


        6*as1*d*e0c*k*pow(epsc1, 2)*epsn - 6*as1*e0c*h*k*pow(epsc1, 2)*epsn - 


        6*as2*e0c*h*k*pow(epsc1, 2)*epsn - 4*b*fc*pow(h, 2)*pow(epsn, 2) + 


        4*b*fc*pow(h, 2)*k*pow(epsn, 2));





  if (disk >= 0.0 ){





    epsc=1/(2.*b*fc*pow(h, 2)*(-4 + 4*k))*(-6*b*fc*pow(h, 2)*epsc1 + 


      9*b*fc*pow(h, 2)*k*epsc1 - 6*as2*ch*e0c*k*pow(epsc1, 2) - 


      6*as1*d*e0c*k*pow(epsc1, 2) + 6*as1*e0c*h*k*pow(epsc1, 2) + 


      6*as2*e0c*h*k*pow(epsc1, 2) + 8*b*fc*pow(h, 2)*epsn - 


      8*b*fc*pow(h, 2)*k*epsn - 


      sqrt(pow(6*b*fc*pow(h, 2)*epsc1 - 9*b*fc*pow(h, 2)*k*epsc1 + 


                  6*as2*ch*e0c*k*pow(epsc1, 2) + 


                  6*as1*d*e0c*k*pow(epsc1, 2) - 6*as1*e0c*h*k*pow(epsc1, 2) - 


                  6*as2*e0c*h*k*pow(epsc1, 2) - 8*b*fc*pow(h, 2)*epsn + 


                  8*b*fc*pow(h, 2)*k*epsn, 2) - 


            4*b*fc*pow(h, 


                2)*(-4 + 4*k)*(6*h*k*n0*pow(epsc1, 2) + 


                  6*b*fc*pow(h, 2)*epsc1*epsn - 


                  9*b*fc*pow(h, 2)*k*epsc1*epsn + 


                  6*as2*ch*e0c*k*pow(epsc1, 2)*epsn + 


                  6*as1*d*e0c*k*pow(epsc1, 2)*epsn - 


                  6*as1*e0c*h*k*pow(epsc1, 2)*epsn - 


                  6*as2*e0c*h*k*pow(epsc1, 2)*epsn - 


                  4*b*fc*pow(h, 2)*pow(epsn, 2) + 


                  4*b*fc*pow(h, 2)*k*pow(epsn, 2))));





  }


  else{


    epsc=-1.0;


  }


   pom =(-(fy/e0c) - epsn)-((-ch + h)*epsc)/h;


   


   if (pom >= 0.0 ){


    epschr1=-1.0;


  }


  else{


    epschr1=epsc;


  }


  return epschr1;


}


 


//hodnota ohyboveho momentu, kdy tlacena plocha x=h a nedochazi k teceni oceli As1 a As2


double normmat::compute_mhr1 ()


{


  double mhr1;	


  


  if(epschr1 < 0 && epschr1 != -1)


    { 


    mhr1 = 0.08333333333333333*(6*as1*


        e0c*(-2*cd + h)*(epsn + 


            1/(6.*(as2*ch + as1*d)*e0c*k*pow(epsc1, 2))*(6*as2*ch*e0c*k*


                    epschr1*pow(epsc1, 2) + 


                  3*d*k*pow(epsc1, 


                      2)*(n0 + as1*e0c*(epschr1 - epsn) - 


                        as2*e0c*(epschr1 + epsn)) + 


                  sqrt(3.)*


                    sqrt(k*pow(epsc1, 


                          2)*(3*k*pow(epsc1, 2)*


                              pow(-2*as2*ch*e0c*epschr1 + 


                                  d*(-n0 + 


                                        e0c*(as1*(-epschr1 + epsn) + 


                                        as2*(epschr1 + epsn))), 2) - 


                            2*(as2*ch + as1*d)*


                              e0c*(b*pow(d, 2)*fc*


                                    pow(epschr1 + epsn, 


                                      2)*(4*(-1 + k)*epschr1 + (6 - 9*k)*


                                        epsc1 + 4*(-1 + k)*epsn) +