00001 #include "scaldam.h"
00002 #include "matrix.h"
00003 #include "vector.h"
00004 #include "elastisomat.h"
00005 #include "global.h"
00006 #include "intpoints.h"
00007 #include "tensor.h"
00008 #include "vecttens.h"
00009 #include <math.h>
00010
00011 #define nijac 20
00012 #define limit 1.0e-8
00013
00014
00015
00016
00017
00018
00019
00020
00021 scaldam::scaldam (void)
00022 {
00023 ft=0.0; uf = 0.0; c = 0.0; k = 0.0; cde = corr_on;
00024 }
00025
00026
00027
00028
00029
00030
00031
00032
00033 scaldam::~scaldam (void)
00034 {
00035
00036 }
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050 void scaldam::read (XFILE *in)
00051 {
00052 xfscanf (in,"%k%lf %k%lf %k%m %k%m","ft",&ft,"uf",&uf,"norm_type",¶mf_type_kwdset,(int *)(&ftype),
00053 "cor_dis_energy", &corr_disip_en_kwdset, (int *)(&cde));
00054 if ((ftype == norenergy) || (ftype == norposenergy))
00055 xfscanf(in, "%k%lf", "c", &c);
00056 if (ftype == vonmises)
00057 xfscanf(in, "%k%lf", "k", &k);
00058
00059 if (cde != corr_off)
00060 sra.read (in);
00061
00062
00063
00064
00065
00066 }
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080 void scaldam::print (FILE *out)
00081 {
00082 fprintf (out,"%le %le %d %d ",ft,uf,ftype,cde);
00083 if ((ftype == norenergy) || (ftype == norposenergy))
00084 fprintf(out, " %e ",c);
00085 if (ftype == vonmises)
00086 fprintf(out, " %e ",k);
00087
00088 if (cde != corr_off)
00089 sra.print (out);
00090 }
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108 void scaldam::damfuncpar(long ipp, vector &eps, vector &kappa)
00109 {
00110 double epseq;
00111 vector poseps;
00112 vector princsig;
00113 vector princeps;
00114 vector posprincsig;
00115 vector tmp;
00116 vector sig;
00117 matrix pvect;
00118 matrix sigt;
00119 matrix epst;
00120 matrix dev;
00121 matrix d;
00122 double i1e, j2e, a, b, e, nu, temp;
00123 long ncomp;
00124
00125 switch (ftype)
00126 {
00127 case norstrain:
00128
00129 scprd(eps, eps, epseq);
00130 kappa[0] = sqrt(epseq);
00131 break;
00132 case norenergy:
00133
00134 ncomp=Mm->ip[ipp].ncompstr;
00135 allocm(ncomp, ncomp, d);
00136 elmatstiff (d,ipp);
00137 vxm(eps, d, tmp);
00138 scprd(tmp, eps, epseq);
00139 kappa[0] = sqrt(epseq)*c;
00140 break;
00141 case norposstrain:
00142
00143 allocv(eps.n,poseps);
00144 extractposv (eps, poseps);
00145 scprd(poseps, poseps, epseq);
00146 kappa[0] = sqrt(epseq);
00147 break;
00148 case norposenergy:
00149
00150 ncomp=Mm->ip[ipp].ncompstr;
00151 allocm(ncomp, ncomp, d);
00152 allocv(ncomp, tmp);
00153 elmatstiff (d,ipp);
00154 allocv(eps.n,poseps);
00155 extractposv (eps, poseps);
00156 vxm(poseps, d, tmp);
00157 scprd(tmp, poseps, epseq);
00158 kappa[0] = sqrt(epseq)*c;
00159 break;
00160 case norrankine:
00161
00162 {
00163 e = Mm->give_actual_ym(ipp);
00164 ncomp=Mm->ip[ipp].ncompstr;
00165 allocv(ncomp, sig);
00166 allocm(ncomp, ncomp, d);
00167 allocm(3, 3, sigt);
00168 allocm(3, 3, pvect);
00169 allocv(3, princsig);
00170 elmatstiff (d,ipp);
00171 mxv(d, eps, sig);
00172 vector_tensor(sig, sigt, Mm->ip[ipp].ssst, stress);
00173 princ_val (sigt,princsig,pvect,nijac,limit,Mp->zero,3,1);
00174
00175
00176 if (princsig[2] > 0.0)
00177 kappa[0] = princsig[2] / e;
00178 else
00179 kappa[0] = 0.0;
00180 break;
00181 }
00182 case norrankinesmooth:
00183
00184 {
00185 e = Mm->give_actual_ym(ipp);
00186 ncomp=Mm->ip[ipp].ncompstr;
00187 allocv(ncomp, sig);
00188 allocm(ncomp, ncomp, d);
00189 allocm(3, 3, sigt);
00190 allocm(3, 3, pvect);
00191 allocv(3, princsig);
00192 allocv(3, posprincsig);
00193 elmatstiff (d,ipp);
00194 mxv(d, eps, sig);
00195 vector_tensor(sig, sigt, Mm->ip[ipp].ssst, stress);
00196 princ_val (sigt,princsig,pvect,nijac,limit,Mp->zero,3,1);
00197 extractposv (princsig, posprincsig);
00198 scprd(posprincsig, posprincsig, epseq);
00199 kappa[0] = sqrt(epseq) / e;
00200 break;
00201 }
00202 case normazar:
00203
00204 {
00205 allocm(3, 3, epst);
00206 allocm(3, 3, pvect);
00207 allocv(3, princeps);
00208 allocv(3,poseps);
00209 fillv(0.0, princeps);
00210 vector_tensor(eps, epst, Mm->ip[ipp].ssst, strain);
00211 princ_val (epst,princeps,pvect,nijac,limit,Mp->zero,3,1);
00212 extractposv (princeps, poseps);
00213 scprd(poseps, poseps, epseq);
00214 kappa[0] = sqrt(epseq);
00215 break;
00216 }
00217 case vonmises:
00218
00219 ncomp=Mm->ip[ipp].ncompstr;
00220 allocm(3, 3, epst);
00221 allocm(3, 3, dev);
00222 vector_tensor(eps, epst, Mm->ip[ipp].ssst, strain);
00223 i1e = first_invar(epst);
00224 deviator (epst,dev);
00225 j2e = second_invar (dev);
00226 e = Mm->give_actual_ym(ipp);
00227 nu = Mm->give_actual_nu(ipp);
00228 a = (k-1.0)/(2.0*k)/(1.0-2.0*nu);
00229 b = 3.0/k/(1+nu)/(1.0+nu);
00230 temp = a*a*i1e*i1e + b*j2e;
00231 kappa[0] = a*i1e+sqrt(temp);
00232
00233
00234
00235
00236
00237 break;
00238 default:
00239 print_err("unknown damage parameter function type", __FILE__, __LINE__, __func__);
00240 }
00241 return;
00242 }
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258 double scaldam::damfunction(long ipp, vector &kappa, vector &omegao)
00259 {
00260 double err,omega = 0.0, dtmp, tmp, h, e, uf_a, aft;
00261 long i,ni, eid;
00262 double indam;
00263
00264 ni=sra.give_ni ();
00265 err=sra.give_err ();
00266
00267 e = Mm->give_actual_ym(ipp);
00268
00269 aft = Mm->give_actual_ft(ipp);
00270
00271 indam = ft/e;
00272
00273 uf_a = uf*(aft/ft);
00274
00275 if (kappa[0] < indam)
00276
00277 return 0.0;
00278 if ((Mm->ip[ipp].hmt & 2) > 0)
00279 {
00280 omega = 1.0-(aft/e/kappa[0])*exp(-(kappa[0]-aft/e)/(uf_a-aft/e));
00281 if ((omega > 1.0) || (omega < 0.0))
00282 {
00283 print_err("computed omega=%le is out of prescribed range <0,1>\n"
00284 " elem %ld, ip %ld:\n"
00285 " kappa=%le, ft=%le, e=%le, uf=%le",
00286 __FILE__, __LINE__, __func__, omega, Mm->elip[ipp]+1, ipp+1, kappa[0], ft, e, uf_a);
00287 abort();
00288 }
00289 return omega;
00290 }
00291
00292 eid = Mm->elip[ipp];
00293 switch (Mt->give_dimension(eid))
00294 {
00295 case 1:
00296 h = Mt->give_length(eid);
00297 break;
00298 case 2:
00299 h = sqrt(Mt->give_area(eid));
00300 break;
00301 case 3:
00302 h = pow(Mt->give_volume(eid), 1.0/3.0);
00303 break;
00304 default:
00305 print_err("unknown dimension of element", __FILE__, __LINE__, __func__);
00306 }
00307 if (cde == corr_off)
00308
00309 {
00310 omega = 1.0-(ft/e/kappa[0])*exp(-(kappa[0]-ft/e)/(uf_a-ft/e));
00311 if ((omega > 1.0) || (omega < 0.0))
00312 {
00313 print_err("computed omega=%le is out of prescribed range <0,1>\n"
00314 " elem %ld, ip %ld:\n"
00315 " kappa=%le, ft=%le, e=%le, uf=%le",
00316 __FILE__, __LINE__, __func__, omega, Mm->elip[ipp]+1, ipp+1, kappa[0], ft, e, uf_a);
00317 abort();
00318 }
00319 return omega;
00320 }
00321
00322 omega = 1.0;
00323 for (i = 0; i < ni; i++)
00324 {
00325 dtmp = -e*kappa[0]+ft/uf_a*h*kappa[0]*exp(-h*omega*kappa[0]/uf_a);
00326 tmp = (1.0-omega)*e*kappa[0]-ft*exp(-h*omega*kappa[0]/uf_a);
00327 if (fabs(tmp) < err*ft)
00328 {
00329 omega += - tmp/dtmp;
00330 break;
00331 }
00332 omega += - tmp/dtmp;
00333 }
00334 if ((omega > 1.0) || (omega < 0.0) || (fabs(tmp) > ft*err))
00335
00336
00337 {
00338 omega = omegao[0];
00339 for (i = 0; i < ni; i++)
00340 {
00341 dtmp = -e*kappa[0]+ft/uf_a*h*kappa[0]*exp(-h*omega*kappa[0]/uf_a);
00342 tmp = (1.0-omega)*e*kappa[0]-ft*exp(-h*omega*kappa[0]/uf_a);
00343 if (fabs(tmp) < err*ft)
00344 {
00345 omega += - tmp/dtmp;
00346 break;
00347 }
00348 omega += - tmp/dtmp;
00349 }
00350 }
00351
00352 if (fabs(tmp) > ft*err)
00353 {
00354 print_err("iteration of omega doesn't converge with required error\n"
00355 " elem %ld, ip %ld:\n"
00356 " kappa=%le, ft=%le, e=%le, uf=%le, h=%le",
00357 __FILE__, __LINE__, __func__, Mm->elip[ipp]+1, ipp+1, kappa[0], ft, e, uf_a, h);
00358 abort();
00359 }
00360
00361 if ((omega > 1.0) || (omega < 0.0))
00362 {
00363 print_err("computed omega=%le is out of prescribed range <0,1>\n"
00364 " elem %ld, ip %ld:\n"
00365 " kappa=%le, ft=%le, e=%le, uf=%le, h=%le",
00366 __FILE__, __LINE__, __func__, omega, Mm->elip[ipp]+1, ipp+1, kappa[0], ft, e, uf_a, h);
00367 abort();
00368 }
00369
00370 return omega;
00371 }
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385 void scaldam::elmatstiff (matrix &d,long ipp)
00386 {
00387 double e, e0;
00388 long idem = Mm->ip[ipp].gemid();
00389 long idoem=Mm->givencompeqother(ipp,0);
00390 long tmp=Mm->givencompeqother(ipp,idem);
00391 idoem -= tmp;
00392
00393 Mm->elmatstiff (d,ipp);
00394
00395 e = Mm->give_actual_ym(ipp);
00396 e0 = Mm->give_actual_ym(ipp,idem,idoem);
00397
00398 cmulm(e/e0, d);
00399 }
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414 void scaldam::matstiff (matrix &d,long ipp,long ido)
00415 {
00416 double dp;
00417
00418 switch (Mp->nlman->stmat)
00419 {
00420 case initial_stiff:
00421 elmatstiff (d,ipp);
00422 break;
00423 case tangent_stiff:
00424 case secant_stiff:
00425 case incr_tangent_stiff:
00426 case ijth_tangent_stiff:
00427 elmatstiff (d,ipp);
00428 dp=Mm->ip[ipp].eqother[ido+1];
00429 if (dp > 0.999999)
00430 dp = 0.999999;
00431 cmulm (1.0-dp,d);
00432 break;
00433 default:
00434 print_err("unknown type of stifness matrix is required", __FILE__, __LINE__, __func__);
00435 }
00436 }
00437
00438
00439
00440
00441
00442
00443
00444
00445
00446
00447
00448
00449
00450
00451
00452 void scaldam::nlstresses (long ipp, long im, long ido)
00453 {
00454 long i,ncomp=Mm->ip[ipp].ncompstr;
00455 long eid;
00456 double dp, e, nu, h, u_crack;
00457 vector epsn(ncomp),sigma(ncomp), kappa(1), omegao(1);
00458 matrix d(ncomp, ncomp);
00459
00460
00461 e = Mm->give_actual_ym(ipp);
00462 if (Mm->ip[ipp].ssst == planestress)
00463 {
00464 nu = Mm->give_actual_nu(ipp);
00465 Mm->ip[ipp].strain[3] = -nu / (1.0 - nu) * (Mm->ip[ipp].strain[0]+Mm->ip[ipp].strain[1]);
00466 }
00467
00468 for (i=0;i<ncomp;i++){
00469 epsn[i] = Mm->ip[ipp].strain[i];
00470 }
00471 kappa[0] = Mm->ip[ipp].eqother[ido+0];
00472 omegao[0] = Mm->ip[ipp].eqother[ido+1];
00473
00474 elmatstiff(d,ipp);
00475
00476 dp = Mm->scal_dam_sol (ipp, im, ido, epsn, kappa, omegao, sigma,d);
00477
00478
00479 for (i=0;i<ncomp;i++){
00480 Mm->ip[ipp].stress[i]=sigma[i];
00481 }
00482
00483 eid = Mm->elip[ipp];
00484 switch (Mt->give_dimension(eid))
00485 {
00486 case 1:
00487 h = Mt->give_length(eid);
00488 break;
00489 case 2:
00490 h = sqrt(Mt->give_area(eid));
00491 break;
00492 case 3:
00493 h = pow(Mt->give_volume(eid), 1.0/3.0);
00494 break;
00495 default:
00496 print_err("unknown dimension of element", __FILE__, __LINE__, __func__);
00497 }
00498 u_crack = dp*kappa[0]*h;
00499
00500 Mm->ip[ipp].other[ido+0]=kappa[0];
00501 Mm->ip[ipp].other[ido+1]=dp;
00502 Mm->ip[ipp].other[ido+2]=u_crack;
00503 }
00504
00505
00506
00507
00508
00509
00510
00511
00512
00513
00514
00515
00516
00517
00518
00519 void scaldam::updateval (long ipp,long im,long ido)
00520 {
00521 long i,n = Mm->givencompeqother(ipp,im);
00522
00523 for (i=0;i<n;i++){
00524 Mm->ip[ipp].eqother[ido+i]=Mm->ip[ipp].other[ido+i];
00525 }
00526
00527 }
00528
00529
00530
00531
00532
00533
00534
00535
00536
00537
00538
00539
00540
00541
00542 double scaldam::give_actual_ft (long ipp, long im, long ido)
00543 {
00544
00545
00546
00547
00548
00549
00550
00551 return ft;
00552 }
00553
00554
00555
00556
00557
00558
00559
00560
00561
00562
00563
00564
00565
00566 double scaldam::epsefunction (long ipp)
00567 {
00568 double e = Mm->give_actual_ym(ipp);
00569 double ft = Mm->give_actual_ft(ipp);
00570
00571 return (ft /e) ;
00572 }
00573
00574
00575
00576
00577
00578
00579
00580
00581
00582
00583
00584
00585
00586 double scaldam::givedamage (long ipp, long ido)
00587 {
00588 return Mm->ip[ipp].eqother[ido+1];
00589 }
00590
00591
00592
00593
00594
00595
00596
00597
00598
00599
00600
00601 double scaldam::give_proczonelength (long ipp, long ido)
00602 {
00603 long eid;
00604 double h;
00605
00606
00607 eid = Mm->elip[ipp];
00608 switch (Mt->give_dimension(eid))
00609 {
00610 case 1:
00611 h = Mt->give_length(eid);
00612 break;
00613 case 2:
00614 h = sqrt(Mt->give_area(eid));
00615 break;
00616 case 3:
00617 h = pow(Mt->give_volume(eid), 1.0/3.0);
00618 break;
00619 default:
00620 print_err("unknown dimension of element", __FILE__, __LINE__, __func__);
00621 }
00622
00623 return h;
00624 }
00625
00626
00627
00628
00629
00630
00631
00632
00633
00634
00635
00636
00637 double scaldam::give_crackwidth (long ipp, long ido)
00638 {
00639 return Mm->ip[ipp].eqother[ido+2];
00640 }
00641
00642
00643
00644
00645
00646
00647
00648
00649
00650
00651
00652 void scaldam::changeparam (atsel &atm,vector &val)
00653 {
00654 long i;
00655
00656 for (i=0;i<atm.num;i++){
00657 switch (atm.atrib[i]){
00658 case 0:{
00659 ft=val[i];
00660 break;
00661 }
00662 case 1:{
00663 uf=val[i];
00664 break;
00665 }
00666 default:{
00667 print_err("wrong index of changed atribute",__FILE__,__LINE__,__func__);
00668 }
00669 }
00670 }
00671 }
00672