00001 #include "creepbs.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 #include "therisomat.h"
00011
00012 creepbs::creepbs (void)
00013 {
00014 timeMax=Mp->timecon.endtime ();
00015 double m = log10(timeMax);
00016 nRetTime=7;
00017 allocv (nRetTime,retTime);
00018 allocv (nRetTime,ert);
00019 retTime[0]=1.0e-9;
00020 retTime[1]=1.0e-3;
00021 retTime[2]=1.0;
00022
00023 retTime[3]=pow(10.,0.25*m);
00024 retTime[4]=pow(10.,0.5*m);
00025 retTime[5]=pow(10.,0.75*m);
00026 retTime[6]=timeMax;
00027 type_h=1;
00028 type_temp=1;
00029 h_s = 0.4;
00030 desht=0.0;
00031 k_s = 1.0;
00032 k_d = 1.0;
00033 timemat=0.0;
00034 alfa=0.0;
00035 allocm(nRetTime, nRetTime, apom);
00036 ccTime = -1.0e30;
00037 }
00038
00039 creepbs::~creepbs (void)
00040 {
00041 }
00042
00043 void creepbs::creepinit (long ipp,double val,nonmechquant nmq)
00044 {
00045 long nc=Mm->ip[ipp].ncompstr;
00046
00047 if(nmq == rel_hum)
00048 Mm->ip[ipp].eqother[nc*(nRetTime+1)+1]=val;
00049
00050 if(nmq == temperature )
00051 Mm->ip[ipp].eqother[nc*(nRetTime+1)+2]=val;
00052 }
00053
00054 void creepbs::read (XFILE *in)
00055 {
00056 xfscanf (in,"%lf %lf %lf %lf %lf %lf %lf %lf %lf %ld %ld",
00057 &tb,&t_w,&fc,&wc,&sc,&gc,&c_s,&a1,&k_d,&type_h,&type_temp);
00058
00059 if (type_h==1){xfscanf (in,"%lf ",&h_s); }
00060 if (type_temp==1){xfscanf (in,"%lf ",&temp_s); }
00061
00062 }
00063
00064
00065
00066
00067
00068
00069
00070 double creepbs::approx (vector &areacoord,vector &nodval)
00071 {
00072 double f;
00073 scprd (areacoord,nodval,f);
00074 return f;
00075 }
00076
00077
00078
00079
00080
00081
00082
00083 void creepbs::inv_sym (matrix &a)
00084 {
00085 long i,j,k;
00086 double diag,an;
00087 for (k=0;k<a.n;k++){
00088 diag = a[k][k];
00089 if (diag<=0.0) { abort();}
00090 for (i=0;i<a.n;i++){
00091 an = a[i][k]/diag;
00092 if (i!=k) {
00093 for (j=i;j<a.n;j++){
00094 if(j!=k) {
00095 a[i][j] = a[i][j] - an*a[k][j];
00096 a[j][i] = a[i][j];
00097 }
00098 }
00099 }
00100 a[i][k] = an;
00101 a[k][i] = an;
00102 }
00103 a[k][k] = - 1.0/diag;
00104 }
00105
00106 for (i=0;i<a.n;i++){
00107 for (j=0;j<a.n;j++){
00108 a[i][j] = - a[i][j];
00109 }
00110 }
00111 }
00112
00113 void creepbs::updateval()
00114 {
00115 if(type_h!=1&&type_temp!=1) updatevalchange();
00116 else updatevalkons();
00117
00118 }
00119 void creepbs::updatevalchange()
00120 {
00121 long i,j,k;
00122 double pomt,pomt1,qq,delYi,delYj,m0,m;
00123
00124 if (ccTime != Mp->timecon.actualtime ())
00125 {
00126
00127 ddTime=Mp->timecon.actualforwtimeincr ();
00128
00129 ccTime=Mp->timecon.actualtime ();
00130 pomt=ccTime+ddTime/2.;
00131 qq=2./3.;
00132
00133 fillm(0.0, apom);
00134 if (timemat<=pomt) {
00135 pomt1=pomt+0.0001;
00136 m0 = log10(pomt1);
00137 m = log10(2.*timeMax)*0.999;
00138 k=0;
00139 while (pomt1<=2.*timeMax)
00140 {
00141 for (i=0;i<nRetTime;i++)
00142 {
00143 delYi=pow((pomt/retTime[i]),qq)-pow((pomt1/retTime[i]),qq);
00144 for (j=0;j<nRetTime;j++)
00145 {
00146 delYj= pow((pomt/retTime[j]),qq)-pow((pomt1/retTime[j]),qq);
00147 apom[i][j]=apom[i][j]+(1.-exp(delYi))*(1.-exp(delYj));
00148 }
00149 }
00150 k=k+1;
00151 pomt1=pow(10.,m0+(m-m0)/40*k);
00152 }
00153 inv_sym(apom);
00154 }
00155 fprintf(stdout, "\n ### inversion\n");
00156 }
00157 }
00158 void creepbs::updatevalkons()
00159 {
00160 long i,j,k;
00161 double pomt,pomt1,qq,jt,delYi,delYj,m0,m;
00162 vector bpom(nRetTime);
00163
00164 if (ccTime != Mp->timecon.actualtime ())
00165 {
00166
00167 ddTime=Mp->timecon.actualforwtimeincr ();
00168
00169 ccTime=Mp->timecon.actualtime ();
00170 pomt=ccTime+ddTime/2.;
00171 qq=2./3.;
00172 jt=0.0;
00173
00174 if(timemat<=pomt) {
00175 et=0.;
00176 for (i=0;i<nRetTime;i++){
00177 bpom[i]=0.;
00178 for (j=0;j<nRetTime;j++){
00179 apom[i][j]=0.;
00180 }
00181 }
00182
00183 pomt1=pomt+0.0001;
00184 m0 = log10(pomt1);
00185 m = log10(2.*timeMax)*0.999;
00186 k=0;
00187 while (pomt1<=2.*timeMax){
00188
00189 {b3_law(jt, pomt, pomt1);}
00190
00191
00192 for (i=0;i<nRetTime;i++){
00193 delYi=pow((pomt/retTime[i]),qq)-pow((pomt1/retTime[i]),qq);
00194 bpom[i]=bpom[i]+(1.-exp(delYi))*jt;
00195 for (j=0;j<nRetTime;j++){
00196 delYj= pow((pomt/retTime[j]),qq)-pow((pomt1/retTime[j]),qq);
00197 apom[i][j]=apom[i][j]+(1.-exp(delYi))*(1.-exp(delYj));
00198 }
00199 }
00200 k=k+1;
00201 pomt1=pow(10.,m0+(m-m0)/40*k);
00202
00203 }
00204
00205 inv_sym(apom);
00206 for (i=0;i<nRetTime;i++){
00207 ert[i]=0.;
00208 for (j=0;j<nRetTime;j++){
00209 ert[i]=ert[i]+apom[i][j]*bpom[j];
00210 }
00211 ert[i]=1./ert[i];
00212 if(ccTime<0.0001){
00213 delYi=pow((0.0001+ddTime)/retTime[i],qq)-pow(0.0001/retTime[i],qq);}
00214 else {
00215 delYi=pow((ccTime+ddTime)/retTime[i],qq)-pow(ccTime/retTime[i],qq);}
00216 et=et+(delYi-1.+exp(-delYi))/delYi/ert[i];
00217 }
00218 if(et>1.e-50) et=1./et;
00219 else{fprintf (Out,"\n\n Stiff is zero"); abort();}
00220
00221
00222 }
00223
00224 else {
00225
00226 et=1.;
00227 for (i=0;i<nRetTime;i++){ ert[i]=0.;}
00228 }
00229
00230 }
00231
00232 }
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243 void creepbs::nlstresses (long ipp)
00244 {
00245
00246 nc=Mm->ip[ipp].ncompstr;
00247
00248
00249
00250
00251
00252 ss=Mm->ip[ipp].ssst;
00253
00254
00255
00256 ddTime=Mp->timecon.actualforwtimeincr ();
00257
00258 ccTime=Mp->timecon.actualtime ();
00259
00260
00261 imat=Mm->ip[ipp].idm[Mm->ip[ipp].gemid()];
00262
00263
00264 if(type_h!=1){
00265 h_s=Mm->givenonmechq(rel_hum, ipp);
00266 }
00267
00268
00269 if(type_temp!=1){
00270 temp_s=Mm->givenonmechq(temperature, ipp);
00271 }
00272
00273 get_h (ipp);
00274 get_temp (ipp);
00275
00276 if (Mp->phase==1){
00277 phase1(ipp);
00278 }
00279
00280 if (Mp->phase==2){
00281 phase2(ipp);
00282 }
00283
00284 }
00285
00286
00287 void creepbs::phase1 (long ipp)
00288 {
00289
00290
00291 long i,j,ii;
00292 double qq,dj,endTime;
00293
00294 vector epscr(nc),depsshtem(nc),sig(nc);
00295 matrix d(nc,nc),screep(nc,nRetTime);
00296
00297
00298
00299
00300 endTime=ccTime+ddTime;
00301 get_desht (desht, ccTime, endTime);
00302 for (i=0; i<nc; i++){
00303 epscr[i]=0.0;
00304 }
00305 qq=2.0/3.0;
00306 for (i=0;i<nRetTime;i++){
00307 ii=nc*(i+1);
00308 dj=pow(endTime/retTime[i],qq)-pow(ccTime/retTime[i],qq);
00309 for (j=0; j<nc; j++){
00310 screep[j][i]=Mm->ip[ipp].eqother[ii+j];
00311 epscr[j]=epscr[j]+(1.-exp(-dj))*screep[j][i];
00312 }
00313 }
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333 epscr[0]= 0.0;
00334 epscr[1]= 0.0;
00335 if(ss==planestrain) epscr[2]=0.0;
00336 else if(ss==axisymm) epscr[2]=0.0;
00337 else if(ss==spacestress) epscr[2]=0.0;
00338
00339
00340
00341
00342
00343
00344
00345
00346 Mm->stiff_deps_creep (ipp, 0,epscr, nc,nc);
00347 }
00348
00349
00350
00351 void creepbs::phase2 (long ipp)
00352 {
00353
00354 long i,j,ii;
00355 double qq,dj, endTime;
00356 vector epscr(nc),deps(nc),depsshtem(nc),sig(nc);
00357 matrix d(nc,nc),screep(nc,nRetTime);
00358
00359 endTime=ccTime+ddTime;
00360 fprintf (Out,"\n\n phase 2 v case t prirustek dt");
00361 fprintf (Out," %e\t%e", ccTime, ddTime);
00362
00363 for (i=0; i<nc; i++){
00364 deps[i]=Mm->ip[ipp].strain[i]-Mm->ip[ipp].eqother[i];
00365 }
00366 qq=2.0/3.0;
00367 for (i=0;i<nRetTime;i++){
00368 ii=nc*(i+1);
00369 dj=pow(endTime/retTime[i],qq)-pow(ccTime/retTime[i],qq);
00370 for (j=0; j<nc; j++){
00371 screep[j][i]=Mm->ip[ipp].eqother[ii+j];
00372 deps[j]=deps[j]-(1.-exp(-dj))*screep[j][i];
00373 }
00374 }
00375 get_desht (desht, ccTime, endTime);
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388 deps[0]=0.0;
00389 deps[1]=0.0;
00390 if(ss==planestrain) deps[2]=0.0;
00391 else if(ss==axisymm) deps[2]=0.0;
00392 else if(ss==spacestress) deps[2]=0.0;
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404 Mm->matstiff(d,ipp);
00405 e0=Mm->eliso[imat].e;
00406 mi=Mm->eliso[imat].nu;
00407 mxv (d,deps,sig);
00408
00409 for (i=0; i<nc; i++){Mm->ip[ipp].eqother[i]=Mm->ip[ipp].strain[i];}
00410
00411 seps_time (screep,sig);
00412 for (i=0; i<nRetTime; i++){
00413 ii=nc*(i+1);
00414 for (j=0; j<nc; j++){
00415 Mm->ip[ipp].eqother[ii+j]=screep[j][i];
00416 }
00417 }
00418
00419
00420
00421
00422
00423
00424 Mm->ip[ipp].eqother[nc*(nRetTime+1)+0]+=desht;
00425 Mm->ip[ipp].eqother[nc*(nRetTime+1)+1]=h_s;
00426 Mm->ip[ipp].eqother[nc*(nRetTime+1)+2]=temp_s;
00427
00428
00429 }
00430
00431 void creepbs::get_h (long ipp)
00432 {
00433 h_slast=1.0;
00434 if(type_h !=1 ){h_slast=Mm->ip[ipp].eqother[nc*(nRetTime+1)+1]; }
00435 }
00436
00437
00438 void creepbs::get_temp (long ipp)
00439 {
00440 temp_slast=283.0;
00441 if(type_temp != 1 ){temp_slast=Mm->ip[ipp].eqother[nc*(nRetTime+1)+2]; }
00442 }
00443
00444
00445
00446
00447
00448
00449
00450
00451
00452
00453 void creepbs::matstiff (matrix &d, long ipp)
00454
00455
00456 {
00457 double qq=1.;
00458
00459 nc=Mm->ip[ipp].ncompstr;
00460 ss=Mm->ip[ipp].ssst;
00461
00462 imat=Mm->ip[ipp].idm[Mm->ip[ipp].gemid()];
00463
00464 Mm->elmatstiff (d,ipp);
00465
00466 e0=Mm->eliso[imat].e;
00467 mi=Mm->eliso[imat].nu;
00468
00469 if(type_h!=1 && type_temp!=1) matstiffchange(qq, ipp);
00470 else matstiffkons(qq);
00471
00472
00473 qq = 1.0;
00474
00475 cmulm( qq, d);
00476 }
00477 void creepbs::matstiffkons (double &qq)
00478
00479 {
00480
00481 qq=et/e0;
00482 }
00483
00484 void creepbs::matstiffchange (double &qq, long ipp)
00485
00486
00487 {
00488 long i,j,k;
00489 double jt,pomt,pomt1,q,delYi,m0,m;
00490 vector bpom(nRetTime);
00491
00492 if(Mm->ip[ipp].hmt & 1)
00493 {i=Mm->ip[ipp].idm[Mm->ip[ipp].nm-1]; alfa=Mm->tidilat[i].alpha;}
00494
00495 if(type_h!=1) {h_s=Mm->givenonmechq(rel_hum, ipp);}
00496 if(type_temp!=1) {temp_s=Mm->givenonmechq(temperature, ipp);}
00497 get_h (ipp);
00498 get_temp (ipp);
00499
00500 pomt=ccTime+ddTime/2.;
00501 q=2./3.;
00502 jt=0.0;
00503
00504 if(timemat<=pomt) {
00505 et=0.;
00506 for (i=0;i<nRetTime;i++){
00507 bpom[i]=0.;
00508
00509
00510
00511 }
00512
00513 pomt1=pomt+0.0001;
00514 m0 = log10(pomt1);
00515 m = log10(2.*timeMax)*0.999;
00516 k=0;
00517 while (pomt1<=2.*timeMax){
00518
00519 {b3_law(jt, pomt, pomt1);}
00520
00521
00522 for (i=0;i<nRetTime;i++){
00523 delYi=pow((pomt/retTime[i]),q)-pow((pomt1/retTime[i]),q);
00524 bpom[i]=bpom[i]+(1.-exp(delYi))*jt;
00525 for (j=0;j<nRetTime;j++){
00526
00527
00528 }
00529 }
00530 k=k+1;
00531 pomt1=pow(10.,m0+(m-m0)/40*k);
00532
00533 }
00534
00535
00536 for (i=0;i<nRetTime;i++){
00537 ert[i]=0.;
00538 for (j=0;j<nRetTime;j++){
00539 ert[i]=ert[i]+apom[i][j]*bpom[j];
00540 }
00541 ert[i]=1./ert[i];
00542 if(ccTime<0.0001){
00543 delYi=pow((0.0001+ddTime)/retTime[i],q)-pow(0.0001/retTime[i],q);}
00544 else {
00545 delYi=pow((ccTime+ddTime)/retTime[i],q)-pow(ccTime/retTime[i],q);}
00546 et=et+(delYi-1.+exp(-delYi))/delYi/ert[i];
00547 }
00548 if(et>1.e-50) et=1./et;
00549 else{fprintf (Out,"\n\n Stiff is zero"); abort();}
00550 }
00551
00552 else {
00553 et=1.;
00554 for (i=0;i<nRetTime;i++){ ert[i]=0.;}
00555 }
00556
00557
00558
00559
00560
00561 qq=et/e0;
00562
00563
00564 }
00565
00566
00567
00568
00569
00570
00571
00572
00573
00574 void creepbs::seps_time (matrix &screep,vector &sig)
00575 {
00576 int i;
00577 double q,dj,pomt,pp;
00578
00579 q=2./3.;
00580 if((ccTime+ddTime)==ddTime) {
00581 }
00582 else{
00583 pomt=ccTime+ddTime;
00584 if(ss==planestress){
00585 for (i=0;i<nRetTime;i++){
00586 dj=pow(pomt/retTime[i],q)-pow(ccTime/retTime[i],q);
00587 screep[0][i]=exp(-dj)*screep[0][i]+(sig[0]-mi*sig[1])*(1.-exp(-dj))/dj/ert[i];
00588 screep[1][i]=exp(-dj)*screep[1][i]+(sig[1]-mi*sig[0])*(1.-exp(-dj))/dj/ert[i];
00589 screep[2][i]=exp(-dj)*screep[2][i]+ sig[2]*(1.+mi) *(1.-exp(-dj))/dj/ert[i];
00590 screep[3][i]=0.0;
00591 }
00592 }
00593 else if(ss==planestrain){
00594 pp=1.+mi;
00595 for (i=0;i<nRetTime;i++){
00596 dj=pow(pomt/retTime[i],q)-pow(ccTime/retTime[i],q);
00597 screep[0][i]=exp(-dj)*screep[0][i]+ pp*(sig[0]*(1.-mi)-mi*sig[1])*(1.-exp(-dj))/dj/ert[i];
00598 screep[1][i]=exp(-dj)*screep[1][i]+ pp*(sig[1]*(1.-mi)-mi*sig[0])*(1.-exp(-dj))/dj/ert[i];
00599 screep[2][i]=exp(-dj)*screep[2][i]+ pp*sig[2]*2. *(1.-exp(-dj))/dj/ert[i];
00600 screep[3][i]=0.0;
00601 }
00602 }
00603 else if(ss==axisymm){
00604 for (i=0;i<nRetTime;i++){
00605 dj=pow(pomt/retTime[i],q)-pow(ccTime/retTime[i],q);
00606 screep[0][i]=exp(-dj)*screep[0][i]+(sig[0]-mi*sig[1]-mi*sig[2])*(1.-exp(-dj))/dj/ert[i];
00607 screep[1][i]=exp(-dj)*screep[1][i]+(sig[1]-mi*sig[0]-mi*sig[2])*(1.-exp(-dj))/dj/ert[i];
00608 screep[2][i]=exp(-dj)*screep[2][i]+(sig[2]-mi*sig[0]-mi*sig[1])*(1.-exp(-dj))/dj/ert[i];
00609 screep[3][i]=exp(-dj)*screep[3][i]+ sig[3]*(1.+mi) *(1.-exp(-dj))/dj/ert[i];
00610 }
00611 }
00612 else if(ss==spacestress){
00613 for (i=0;i<nRetTime;i++){
00614 dj=pow(pomt/retTime[i],q)-pow(ccTime/retTime[i],q);
00615 screep[0][i]=exp(-dj)*screep[0][i]+(sig[0]-mi*sig[1]-mi*sig[2])*(1.-exp(-dj))/dj/ert[i];
00616 screep[1][i]=exp(-dj)*screep[1][i]+(sig[1]-mi*sig[0]-mi*sig[2])*(1.-exp(-dj))/dj/ert[i];
00617 screep[2][i]=exp(-dj)*screep[2][i]+(sig[2]-mi*sig[0]-mi*sig[1])*(1.-exp(-dj))/dj/ert[i];
00618 screep[3][i]=exp(-dj)*screep[3][i]+ sig[3]*(1.+mi) *(1.-exp(-dj))/dj/ert[i];
00619 screep[4][i]=exp(-dj)*screep[4][i]+ sig[4]*(1.+mi) *(1.-exp(-dj))/dj/ert[i];
00620 screep[5][i]=exp(-dj)*screep[5][i]+ sig[5]*(1.+mi) *(1.-exp(-dj))/dj/ert[i];
00621 }
00622 }
00623 }
00624
00625 }
00626
00627
00628
00629
00630
00631
00632
00633
00634
00635
00636 void creepbs::get_desht (double &des_hn, double t0, double t)
00637 {
00638 double esn,kh,tau,ht,ht0,dh_s,dtemp;
00639
00640
00641
00642
00643
00644
00645
00646
00647
00648
00649
00650
00651
00652
00653 if (tb<=t0){
00654
00655
00656 if ( type_h ==1){
00657 if ( (t0>=t_w) ){
00658
00659 tau=k_s*k_s*k_d*k_d*85000./pow(t,0.08)/pow(fc,0.25);
00660 esn=-a1*1.2*(0.019*pow((wc*c_s),2.1)/pow(fc,0.28)+270.);
00661 ht0 =tanh(sqrt((t0-t_w)/tau));
00662 ht =tanh(sqrt((t -t_w)/tau));
00663 if(h_s<=0.98){kh=(1.0-pow(h_s,3.0));}
00664 else {kh=-10.0*(1.-h_s);}
00665 des_hn = esn*kh/1000000.0*(ht-ht0);
00666 }
00667 else {
00668 des_hn=0.0;
00669 }
00670 }
00671 else {
00672
00673 dh_s = h_s-h_slast;
00674 esn = a1*1.2*(0.019*pow((wc*c_s),2.1)/pow(fc,0.28)+270.);
00675 kh = 3.0*h_s*h_s;
00676 des_hn = esn*kh*dh_s*1.e-6;
00677 }
00678 if ( type_temp ==1 ){
00679
00680 }
00681 else {
00682
00683 dtemp = temp_s-temp_slast;
00684 esn=alfa*1.e6;
00685 des_hn=des_hn+esn*dtemp*1.e-6;
00686 }
00687 }
00688
00689 else {
00690 des_hn=0.0;
00691 }
00692
00693 des_hn = 0.0;
00694
00695 }
00696
00697
00698
00699
00700
00701
00702
00703
00704
00705 void creepbs::b3_law (double &jt, double t0, double t)
00706 {
00707 double n,c0,cd, ac,ag,m,x,y,q,q1,q2,q3,q4,q5,esn,kh,tau,ht,ht0,dh_s,dtemp;
00708
00709
00710
00711
00712
00713
00714
00715
00716
00717
00718
00719
00720
00721
00722
00723
00724
00725
00726
00727
00728 if (tb<t0){
00729 if (t_w<tb) t_w = tb;
00730 if (t_w<=0.0) t_w = 1.0;
00731 ac=sc+gc; ag=ac/gc; m=0.5;
00732
00733 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;
00734 if (x<=4){n=0.12;}
00735 else {n=0.12+(0.07*pow(x,6.0)/(5130.0+pow(x,6.0)));}
00736 q1=600000.0/4734.0/sqrt(fc);
00737 x=1.7*pow(t0,0.12)+8.0;
00738 y=(0.086*pow(t0,(2.0/9.0))+1.21*pow(t0,(4.0/9.0)));
00739 q=1.0/y/pow( ( 1.0+ pow( (pow(t0,m)/log(1.+pow((t-t0),n))/y),x ) ),(1./x));
00740 q2=185.4*sqrt(c_s)/pow(fc,0.9);
00741 q3=0.29*pow(wc,4.0)*q2;
00742 q4=20.3/pow(ac,0.7);
00743
00744 c0=q2*q+q3*log(1.0+pow((t-t0),n))+q4*log(t/t0);
00745
00746
00747 if ( type_h ==1){
00748 if ( (t0>=t_w) ){
00749
00750 tau=k_s*k_s*k_d*k_d*85000./pow(t,0.08)/pow(fc,0.25);
00751 esn=a1*1.2*(0.019*pow((wc*c_s),2.1)/pow(fc,0.28)+270.);
00752 ht0 =tanh(sqrt((t0-t_w)/tau));
00753 ht =tanh(sqrt((t -t_w)/tau));
00754 if(h_s<=0.98){kh=(1.0-pow(h_s,3.0));}
00755 else {kh=-10.0*(1.-h_s);}
00756
00757 q5 = 0.757/fc/pow(esn,0.6);
00758 cd = q5*sqrt( exp(-8+8.*(1.0-h_s)*ht)-exp(-8.+8.*(1.0-h_s)*ht0) );
00759 }
00760 else {
00761 cd = 0.0;
00762 }
00763 }
00764 else {
00765
00766 dh_s = h_s-h_slast;
00767 esn = a1*1.2*(0.019*pow((wc*c_s),2.1)/pow(fc,0.28)+270.);
00768 kh = 3.0*h_s*h_s;
00769 q5 = 0.35/fc;
00770 cd = fabs(q5*esn*kh*dh_s);
00771 }
00772 if ( type_temp ==1 ){
00773
00774 }
00775 else {
00776
00777 dtemp = temp_s-temp_slast;
00778 esn=alfa*1.e6;
00779 q5= 1.5/fc;
00780 cd = cd+fabs(q5*esn*dtemp);
00781 }
00782 }
00783
00784 else {
00785 q1 = 0.0;
00786 c0 = 0.0;
00787 cd = 0.0;
00788 }
00789
00790
00791
00792 jt=(q1+c0+cd)*1.e-6;
00793
00794
00795 }