00001 #include "creepbbeam.h"
00002 #include "global.h"
00003 #include "math.h"
00004 #include "globmat.h"
00005 #include "genfile.h"
00006 #include "adaptivity.h"
00007 #include "intpoints.h"
00008 #include "stdlib.h"
00009 #include "elastisomat.h"
00010
00011 creepbbeam::creepbbeam (void)
00012 {
00013 nRetTime=7;
00014 allocv (nRetTime,retTime);
00015 allocv (nRetTime,ert);
00016 retTime[0]=1.0e-9;
00017 retTime[1]=1.0e-3;
00018 retTime[2]=1.0;
00019 retTime[3]=14.0;
00020 retTime[4]=100.0;
00021 retTime[5]=800.0;
00022 retTime[6]=3000.0;
00023 type_h=1;
00024 type_temp=1;
00025 timeMax=Mp->timecon.endtime ();
00026 h_s = 0.4;
00027 temp_s=20;
00028 esht=0.0;
00029 nc=1;
00030 k_s = 1.0;
00031 r_s = 0.8;
00032 k_d = 0.12;
00033 timemat=0.0;
00034 }
00035
00036 creepbbeam::~creepbbeam (void)
00037 {
00038 }
00039
00040 void creepbbeam::creepinit (long mie)
00041 {
00042 }
00043
00044 void creepbbeam::read (FILE *in)
00045 {
00046 fscanf (in,"%lf %lf %lf %lf %lf %lf %lf %lf %ld %ld",
00047 &tb,&t_w,&fc,&wc,&sc,&gc,&c_s,&a1,&type_h,&type_temp);
00048
00049 if (type_h==1){fscanf (in,"%lf ",&h_s); }
00050 if (type_temp==1){fscanf (in,"%lf ",&temp_s); }
00051
00052 }
00053
00054
00055
00056
00057
00058
00059
00060 double creepbbeam::approx (vector &areacoord,vector &nodval)
00061 {
00062 double f;
00063 scprd (areacoord,nodval,f);
00064 return f;
00065 }
00066
00067
00068
00069
00070
00071
00072
00073 void creepbbeam::inv_sym (matrix &a)
00074 {
00075 long i,j,k;
00076 double diag,an;
00077 for (k=0;k<a.n;k++){
00078 diag = a[k][k];
00079 if (diag<=0.0) { abort();}
00080 for (i=0;i<a.n;i++){
00081 an = a[i][k]/diag;
00082 if (i!=k) {
00083 for (j=i;j<a.n;j++){
00084 if(j!=k) {
00085 a[i][j] = a[i][j] - an*a[k][j];
00086 a[j][i] = a[i][j];
00087 }
00088 }
00089 }
00090 a[i][k] = an;
00091 a[k][i] = an;
00092 }
00093 a[k][k] = - 1.0/diag;
00094 }
00095
00096 for (i=0;i<a.n;i++){
00097 for (j=0;j<a.n;j++){
00098 a[i][j] = - a[i][j];
00099 }
00100 }
00101 }
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111 void creepbbeam::nlstresses (long ipp)
00112 {
00113
00114
00115
00116 ddTime=Mp->timecon.actualforwtimeincr ();
00117
00118 ccTime=Mp->timecon.actualtime ();
00119 imat=Mm->ip[ipp].idm[Mm->ip[ipp].gemid()];
00120
00121
00122 get_h (ipp);
00123 get_temp (ipp);
00124 if (Mp->phase==1){
00125 phase1(ipp);
00126 }
00127
00128 if (Mp->phase==2){
00129 phase2(ipp);
00130 }
00131
00132 }
00133
00134 void creepbbeam::phase1 (long ipp)
00135 {
00136
00137
00138 long i,j,ii;
00139 double qq,dj;
00140 vector epscr(nc),deps(nc),sig(nc);
00141 matrix d(nc,nc),screep(nc,nRetTime);
00142
00143 Mm->matstiff(d,ipp);
00144 e0=Mm->eliso[imat].e;
00145 mi=Mm->eliso[imat].nu;
00146 for (i=0; i<nc; i++){epscr[i]=0.0;}
00147 qq=2.0/3.0;
00148 for (i=0;i<nRetTime;i++){
00149 ii=nc*(i+1);
00150 dj=pow((ccTime+ddTime)/retTime[i],qq)-pow(ccTime/retTime[i],qq);
00151 for (j=0; j<nc; j++){
00152 screep[j][i]=Mm->ip[ipp].other[ii+j];
00153 epscr[j]=epscr[j]+(1.-exp(-dj))*screep[j][i];
00154 }
00155 }
00156
00157
00158
00159
00160
00161
00162
00163 epscr[0]=epscr[0]+deps0;
00164 epscr[1]=epscr[1]+deps0;
00165 if(epscr.n==6) epscr[2]=epscr[2]+deps0;
00166
00167
00168
00169
00170 Mm->stiff_deps_creep (ipp, epscr, nc,nc);
00171 }
00172
00173
00174
00175 void creepbbeam::phase2 (long ipp)
00176 {
00177
00178 long i,j,ii;
00179 double qq,dj;
00180 vector epscr(nc),deps(nc),sig(nc);
00181 matrix d(nc,nc),screep(nc,nRetTime);
00182
00183 fprintf (Out,"\n\n phase 2 v case t prirustek dt");
00184 fprintf (Out," %e\t%e", ccTime, ddTime);
00185
00186 for (i=0; i<nc; i++){deps[i]=Mm->ip[ipp].strain[i]-Mm->ip[ipp].other[i];}
00187 qq=2.0/3.0;
00188 for (i=0;i<nRetTime;i++){
00189 ii=nc*(i+1);
00190 dj=pow((ccTime+ddTime)/retTime[i],qq)-pow(ccTime/retTime[i],qq);
00191 for (j=0; j<nc; j++){
00192 screep[j][i]=Mm->ip[ipp].other[ii+j];
00193 deps[j]=deps[j]-(1.-exp(-dj))*screep[j][i];
00194 }
00195 }
00196
00197 deps[0]=deps[0]-deps0;
00198 deps[1]=deps[1]-deps0;
00199 if(deps.n==6) deps[2]=deps[2]-deps0;
00200
00201
00202
00203
00204
00205 fprintf (Out,"\n\n delta eps");
00206 fprintf (Out,"\n %e %e %e",deps[0],deps[1],deps[2]);
00207
00208 Mm->matstiff(d,ipp);
00209 e0=Mm->eliso[imat].e;
00210 mi=Mm->eliso[imat].nu;
00211 mxv (d,deps,sig);
00212
00213 for (i=0; i<nc; i++){Mm->ip[ipp].other[i]=Mm->ip[ipp].strain[i];}
00214
00215 seps_time (screep,sig);
00216 for (i=0; i<nRetTime; i++){
00217 ii=nc*(i+1);
00218 for (j=0; j<nc; j++){
00219 Mm->ip[ipp].other[ii+j]=screep[j][i];
00220 }
00221 }
00222
00223
00224
00225
00226 Mm->ip[ipp].other[nc*(nRetTime+1)+1]=esht;
00227 Mm->ip[ipp].other[nc*(nRetTime+1)+2]=h_s;
00228 Mm->ip[ipp].other[nc*(nRetTime+1)+3]=temp_s;
00229 }
00230
00231
00232 void creepbbeam::get_h (long ipp)
00233 {
00234 h_slast=0.0;
00235 if(type_h != 1){h_slast=Mm->ip[ipp].other[nc*(nRetTime+1)+2]; }
00236 }
00237
00238
00239 void creepbbeam::get_temp (long ipp)
00240 {
00241 templast=0.0;
00242 if(type_temp != 1){templast=Mm->ip[ipp].other[nc*(nRetTime+1)+2]; }
00243 }
00244
00245
00246
00247
00248
00249
00250
00251
00252 void creepbbeam::matstiff (matrix &d, long ipp)
00253
00254
00255 {
00256 long i,j,k;
00257 double jt,et,pomt,pomt1,qq,delYi,delYj;
00258 vector bpom(nRetTime);
00259 matrix apom(nRetTime,nRetTime);
00260
00261
00262
00263
00264 Mm->elmatstiff (d,ipp);
00265
00266
00267 e0=Mm->eliso[imat].e;
00268 mi=Mm->eliso[imat].nu;
00269
00270
00271 ddTime=Mp->timecon.actualforwtimeincr ();
00272
00273 ccTime=Mp->timecon.actualtime ();
00274 pomt=ccTime+ddTime/2.;
00275 qq=2./3.;
00276 jt=0.0;
00277
00278
00279
00280
00281 if(timemat<=pomt) {
00282
00283 et=0.;
00284 for (i=0;i<nRetTime;i++){
00285 bpom[i]=0.;
00286 for (j=0;j<nRetTime;j++){
00287 apom[i][j]=0.;
00288 }
00289 }
00290 k=0;
00291 pomt1=pomt+0.0001;
00292 while (pomt1<=timeMax*5.){
00293
00294 b3_law(jt,esht, pomt, pomt1);
00295
00296
00297 for (i=0;i<nRetTime;i++){
00298 delYi=pow((pomt/retTime[i]),qq)-pow((pomt1/retTime[i]),qq);
00299 bpom[i]=bpom[i]+(1.-exp(delYi))*jt;
00300 for (j=0;j<nRetTime;j++){
00301 delYj= pow((pomt/retTime[j]),qq)-pow((pomt1/retTime[j]),qq);
00302 apom[i][j]=apom[i][j]+(1.-exp(delYi))*(1.-exp(delYj));
00303 }
00304 }
00305
00306 k=k+1;
00307 pomt1=pomt1+20*(k);
00308 }
00309
00310 inv_sym(apom);
00311 for (i=0;i<nRetTime;i++){
00312
00313 ert[i]=0.;
00314 for (j=0;j<nRetTime;j++){
00315
00316 ert[i]=ert[i]+apom[i][j]*bpom[j];
00317 }
00318
00319 ert[i]=1./ert[i];
00320 if(ccTime<0.0001){
00321 delYi=pow((0.0001+ddTime)/retTime[i],qq)-pow(0.0001/retTime[i],qq);}
00322 else {
00323 delYi=pow((ccTime+ddTime)/retTime[i],qq)-pow(ccTime/retTime[i],qq);}
00324
00325
00326 et=et+(delYi-1.+exp(-delYi))/delYi/ert[i];
00327 }
00328
00329 et=1./et;
00330
00331
00332
00333
00334 pomt1=pomt+0.0001;
00335 b3_law(jt,esht, pomt, pomt1);
00336 }
00337
00338 else {
00339
00340 et=e0/100000.;
00341 for (i=0;i<nRetTime;i++){ ert[i]=0.;}
00342 }
00343
00344
00345
00346
00347
00348
00349 qq=et/e0;
00350
00351 cmulm( qq, d);
00352 }
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362 void creepbbeam::seps_time (matrix &screep,vector &sig)
00363 {
00364 int i;
00365 double qq,dj,pomt;
00366
00367 qq=2./3.;
00368 if((ccTime+ddTime)==ddTime) {
00369 }
00370 else{
00371 pomt=ccTime+ddTime;
00372 if(sig.n==4){
00373 for (i=0;i<nRetTime;i++){
00374 dj=pow(pomt/retTime[i],qq)-pow(ccTime/retTime[i],qq);
00375 screep[0][i]=exp(-dj)*screep[0][i]+(sig[0]-mi*sig[1])*(1.-exp(-dj))/dj/ert[i];
00376 screep[1][i]=exp(-dj)*screep[1][i]+(sig[1]-mi*sig[0])*(1.-exp(-dj))/dj/ert[i];
00377 screep[2][i]=exp(-dj)*screep[2][i]+ sig[2]*(1.+mi) *(1.-exp(-dj))/dj/ert[i];
00378 screep[3][i]=0.0;
00379 }
00380 }
00381 else if(sig.n==6){
00382 for (i=0;i<nRetTime;i++){
00383 dj=pow(pomt/retTime[i],qq)-pow(ccTime/retTime[i],qq);
00384 screep[0][i]=exp(-dj)*screep[0][i]+(sig[0]-mi*sig[1]-mi*sig[2])*(1.-exp(-dj))/dj/ert[i];
00385 screep[1][i]=exp(-dj)*screep[1][i]+(sig[1]-mi*sig[0]-mi*sig[2])*(1.-exp(-dj))/dj/ert[i];
00386 screep[2][i]=exp(-dj)*screep[2][i]+(sig[2]-mi*sig[0]-mi*sig[1])*(1.-exp(-dj))/dj/ert[i];
00387 screep[3][i]=exp(-dj)*screep[3][i]+ sig[3]*(1.+mi) *(1.-exp(-dj))/dj/ert[i];
00388 screep[4][i]=exp(-dj)*screep[4][i]+ sig[4]*(1.+mi) *(1.-exp(-dj))/dj/ert[i];
00389 screep[5][i]=exp(-dj)*screep[5][i]+ sig[5]*(1.+mi) *(1.-exp(-dj))/dj/ert[i];
00390 }
00391 }
00392 }
00393
00394 }
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404 void creepbbeam::b2_law (double &jt,double &des_hn, double t0, double t )
00405 {
00406 double n, c, c0,cd, ac,ag,m,x,fi,esn,esn1,tau,sd,st,st0, kt,kh,et1,et2;
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431 if (tb<t0){
00432 if (t_w<tb) t_w = tb;
00433 if (t_w<=0.0) t_w = 1.0;
00434 ac=sc+gc; ag=ac/gc; m=-0.28-47.541025/fc/fc;
00435
00436 x=(2.1*ac/pow(sc,1.4)+0.1*pow((fc/6895.),1.5)*pow(wc,(1./3.))*pow(ag,2.2))*a1-4.;
00437 if (x<=4){n=0.12;}
00438 else {n=0.12+(0.07*pow(x,6.0)/(5130.0+pow(x,6.0)));}
00439
00440 fi=0.5*pow(10.,(3.*n))/(pow(28.,m)+0.025/wc);
00441 c0=fi*(pow((t0-tb),m)+0.025/wc)*pow((t-t0),n);
00442
00443
00444 x=( (1.25*sqrt(ac)+0.5*pow((gc/sc),2))*pow(((1.+sc)/wc),(1/3.))*sqrt(fc/6895.) )-12;
00445 if (x>=0.0){esn=0.0121;}
00446 else {esn=0.0121-0.0088/(390./pow(x,4)+1.);}
00447
00448 c=(0.000000125*wc*c_s-0.000012);
00449 if (c<0.000007) c=0.000007;
00450 if (c>0.000021) c=0.000021;
00451 kt=(temp_s+273.)/296*exp(5000/296.-5000./(temp_s+273.));
00452 c=c*kt*(0.05+sqrt(6.3/t_w));
00453 tau=0.267*(k_s*k_d)*(k_s*k_d)/c;
00454 if (5<=2){
00455
00456
00457 et1=e0/(1.+pow(0.1,(3.*n))*fi*( pow(607.,m) +0.025/wc));
00458 et2=e0/(1.+pow(0.1,(3.*n))*fi*(pow((t_w+tau),m)+0.025/wc));
00459 kh=1-pow(h_s,3.)+pow((1-h_s),5.);
00460 st= 1./(1.+pow(pow(tau/(t-t_w),r_s),0.5*r_s) );
00461 st0= 1./(1.+pow(pow(tau/(t0-t_w),r_s),0.5*r_s) );
00462 esn1=esn*et1/et2*st*kh/1000000.0;
00463
00464 x=56000*pow((ac*fc/6895.),0.3)*pow(gc,1.3)*pow((wc/esn),1.5)-0.85;
00465 if (x>=0.0){fi=0.008+0.027*(1./(1.+0.7/pow(x,1.4))); }
00466 else fi=0.008;
00467
00468 fi = 1./sqrt(1.+(t-t0)/10./tau)*fi;
00469 sd = 1./pow((1.+10.*tau/(t-t_w)),(2.6*n - 6.*n*n));
00470
00471 cd = fi*pow(t_w,(m/2.))*(1.0 - pow(h_s,1.5))*esn1*sd;
00472 des_hn=-esn*et1/et2*(st0-st)*kh/1000000.0;
00473 }
00474 else {
00475 cd = 0.0;
00476 des_hn=0.0;
00477 }
00478 }
00479 else {
00480 cd = 0.0;
00481 des_hn=0.0;
00482 }
00483
00484
00485 jt=(1+c0+cd)/(1.5*e0);
00486
00487
00488 }
00489
00490
00491 void creepbbeam::b3_law (double &jt,double &des_hn, double t0, double t)
00492 {
00493 double n,c0,cd, ac,ag,m,x,y,q,q1,q2,q3,q4,q5,esn,kh,tau,ht,ht0,dh_s,dtemp;
00494 short experiment=0;
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504
00505
00506
00507
00508
00509
00510
00511
00512
00513
00514
00515
00516 if (tb<t0){
00517 if (t_w<tb) t_w = tb;
00518 if (t_w<=0.0) t_w = 1.0;
00519 ac=sc+gc; ag=ac/gc; m=0.5;
00520
00521 x=(2.1*ac/pow(sc,1.4)+0.1*pow((fc/6895),1.5)*pow(wc,(1.0/3.0))*pow(ag,2.2))*a1-4;
00522 if (x<=4){n=0.12;}
00523 else {n=0.12+(0.07*pow(x,6.0)/(5130.0+pow(x,6.0)));}
00524 q1=600000.0/4734.0/sqrt(fc);
00525 x=1.7*pow(t0,0.12)+8.0;
00526 y=(0.086*pow(t0,(2.0/9.0))+1.21*pow(t0,(4.0/9.0)));
00527 q=1.0/y/pow( ( 1.0+ pow( (pow(t0,m)/log(1.+pow((t-t0),n))/y),x ) ),(1./x));
00528 q2=185.4*sqrt(c_s)/pow(fc,0.9);
00529 q3=0.29*pow(wc,4.0)*q2;
00530 q4=20.3/pow(ac,0.7);
00531
00532 c0=q2*q+q3*log(1.0+pow((t-t0),n))+q4*log(t/t0);
00533
00534
00535 if ( type_h ==1){
00536 if ( (t0>=t_w) ){
00537
00538 tau=k_s*k_s*k_d*k_d*85000./pow(t,0.08)/pow(fc,0.25);
00539 esn=-a1*1.2*(0.019*pow((wc*c_s),2.1)/pow(fc,0.28)+270.);
00540 ht0 =tanh(sqrt((t0-t_w)/tau));
00541 ht =tanh(sqrt((t -t_w)/tau));
00542 if(h_s<=0.98){kh=(1.0-pow(h_s,3.0));}
00543 else {kh=-10*(1.-h_s);}
00544 des_hn = esn*kh/1000000.0*(ht-ht0);
00545
00546 q5 = 0.757/fc/pow(fabs(esn),0.6);
00547 cd = q5*sqrt( -exp(-8+8.*(1.0-h_s)*ht)+exp(-8.+8.*(1.0-h_s)*ht0) );
00548 }
00549 else {
00550 cd = 0.0;
00551 des_hn=0.0;
00552 }
00553 }
00554 else {
00555
00556 dh_s = h_s-h_slast;
00557 esn = -a1*1.2*(0.019*pow((wc*c_s),2.1)/pow(fc,0.28)+270.);
00558 if(h_s<=0.98){kh=pow(h_s,3.0)-pow(h_slast,3.0);}
00559 else {kh=-10.0*dh_s;}
00560
00561 des_hn = esn*kh/1000000.0;
00562
00563 q5 = 0.757/fc/pow(fabs(esn),0.6);
00564 cd = -q5*dh_s;
00565 }
00566 if ( type_temp ==1 ){
00567
00568 }
00569 else {
00570
00571 dtemp = temp_s-templast;
00572 esn=alfa;
00573
00574 des_hn=des_hn+esn*dtemp;
00575
00576 q5= 1.5/fc;
00577 cd = cd+q5*esn*dtemp;
00578 }
00579 }
00580
00581 else {
00582 q1 = 0.0;
00583 c0 = 0.0;
00584 cd = 0.0;
00585 des_hn=0.0;
00586 }
00587
00588
00589 jt=(q1+c0+cd)*1.e-6;
00590
00591
00592 }