00001 #include "ortodam.h"
00002 #include "bmatrix.h"
00003 #include "matrix.h"
00004 #include "vector.h"
00005 #include "elastisomat.h"
00006 #include "global.h"
00007 #include "intpoints.h"
00008 #include "tensor.h"
00009 #include "vecttens.h"
00010 #include <math.h>
00011 #include <stdlib.h>
00012
00013 #define nijac 200
00014 #define limit 1.0e-15 // zero limit for principal values of tensors (Jacobi method)
00015 #define qlimit 1.0e-8 // zero limit for switch from quadratuc Bezier to pure quadrat
00016 #ifndef M_PI
00017 #define M_PI 3.14159265358979323846
00018 #endif
00019 #define sqr(x) (x)*(x)
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029 ortodam::ortodam (void)
00030 {
00031 cde = corr_off;
00032 fat = fatigue_off;
00033 ac = 0.0; bc = 0.0; y0c = 0.0; at = 0.0; bt = 0.0; y0t = 0.0;
00034 gfc = 0.0; gft = 0.0;
00035 ft = uft = ult = 0.0;
00036 fc = ufc = ulc = 0.0;
00037 pqt = pqc = 0;
00038 betac = betat = 0.0;
00039 }
00040
00041
00042
00043
00044
00045
00046
00047
00048 ortodam::~ortodam (void)
00049 {
00050
00051 }
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066 void ortodam::read (XFILE *in)
00067 {
00068 xfscanf (in, "%k%m", "dam_evol_f", &dam_evolfunc_kwdset, (int *)&damevf);
00069 switch(damevf)
00070 {
00071 case brittle:
00072 xfscanf (in,"%k%lf %k%lf", "f_t", &ft, "u_ft", &uft);
00073 xfscanf (in,"%k%lf %k%lf", "f_c", &fc, "u_fc", &ufc);
00074 xfscanf (in,"%k%m", "corr_dissip", &corr_disip_en_kwdset, (int *)(&cde));
00075 if (cde != corr_off)
00076 sra.read (in);
00077 break;
00078 case quasi_brittle:
00079 xfscanf (in, "%k%m", "corr_dissip", &corr_disip_en_kwdset, (int *)(&cde));
00080 switch (cde)
00081 {
00082 case corr_off:
00083 xfscanf (in, "%k%lf %k%lf %k%lf", "A_c", &ac, "B_c", &bc, "Y_0c", &y0c);
00084 xfscanf (in, "%k%lf %k%lf %k%lf", "A_t", &at, "B_t", &bt, "Y_0t", &y0t);
00085 break;
00086 case corr_on:
00087 xfscanf (in, "%k%lf %k%lf %k%lf", "G_fc", &gfc, "B_c", &bc, "Y_0c", &y0c);
00088 xfscanf (in, "%k%lf %k%lf %k%lf", "G_ft", &gft, "B_t", &bt, "Y_0t", &y0t);
00089 break;
00090 default:
00091 print_err(" unknown type of correction of disipated energy is required", __FILE__, __LINE__,__func__);
00092 }
00093 break;
00094 case quadbezier:
00095 xfscanf (in,"%k%lf %k%lf %k%lf", "f_t", &ft, "u_ft", &uft, "u_limt", &ult);
00096 if (fabs(uft/ult-0.5) < qlimit)
00097 pqt = 1;
00098 xfscanf (in,"%k%lf %k%lf %k%lf", "f_c", &fc, "u_fc", &ufc, "u_limc", &ulc);
00099 if (fabs(ufc/ulc-0.5) < qlimit)
00100 pqc = 1;
00101 xfscanf (in,"%k%m", "corr_dissip", &corr_disip_en_kwdset, (int *)(&cde));
00102 if (cde != corr_off)
00103 sra.read (in);
00104 break;
00105 default:
00106 print_err(" unknown type of damage evolution function is required", __FILE__, __LINE__,__func__);
00107 }
00108 xfscanf (in,"%k%m", "dam_eq_strain", ¶mf_type_kwdset, (int *)(&dameqstr));
00109 xfscanf (in,"%k%m", "fatigue", &fatigue_flag_kwdset, (int *)(&fat));
00110 if (fat == fatigue_on)
00111 xfscanf (in,"%k%lf %k%lf", "beta_t", &betat, "beta_c", &betac);
00112 }
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128 double ortodam::dam_eq_strain(double nu, long i, vector &peps)
00129 {
00130 double kappa = 0.0;
00131 switch (dameqstr)
00132 {
00133 case norstrain:
00134
00135 if (peps(i) > 0.0)
00136 kappa = peps(i);
00137
00138 if (peps(i) < 0.0)
00139 kappa = fabs(peps(i));
00140 break;
00141 case norenergy:
00142
00143 if (peps(i) > 0.0)
00144 {
00145 kappa = (1-nu)*peps(i)*peps(i);
00146 if (peps((i+1)%3) > 0.0)
00147 kappa += nu*peps(i)*peps((i+1)%3);
00148 if (peps((i+2)%3) > 0.0)
00149 kappa += nu*peps(i)*peps((i+2)%3);
00150 kappa /= (1-nu);
00151
00152
00153 kappa = sqrt(kappa);
00154 }
00155
00156 if (peps(i) < 0.0)
00157 {
00158 kappa = (1-nu)*peps(i)*peps(i);
00159 if (peps((i+1)%3) < 0.0)
00160 kappa += nu*peps(i)*peps((i+1)%3);
00161 if (peps((i+2)%3) < 0.0)
00162 kappa += nu*peps(i)*peps((i+2)%3);
00163 kappa /= (1-nu);
00164
00165
00166 kappa = sqrt(kappa);
00167 }
00168 break;
00169 case norrankine:
00170
00171 kappa = (1-nu)*peps(i) + nu*(peps((i+1)%3) + peps((i+2)%3));
00172 kappa /= (1-2*nu)*(1+nu);
00173 kappa = kappa;
00174 break;
00175 default:
00176 print_err(" unknown type of equivalent strain norm is required", __FILE__, __LINE__,__func__);
00177 }
00178 return kappa;
00179 }
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203 void ortodam::loadfunc(long ipp, double nu, vector &peps, vector &damt, vector &damc,
00204 double aat, double aac, vector &lft, vector &lfc)
00205 {
00206 long i;
00207 double e, af, auf, aul, kappa, tmp1, tmp2, a, b, h;
00208 vector pyt(3), pyc(3);
00209
00210 nu = Mm->give_actual_nu(ipp);
00211 switch(damevf)
00212 {
00213 case brittle:
00214 for (i=0; i<3; i++)
00215 {
00216
00217 if (peps(i) >= 0.0)
00218 {
00219 kappa = dam_eq_strain(nu, i, peps);
00220 e = Mm->give_actual_ym(ipp);
00221 af = Mm->give_actual_ft(ipp);
00222 if (kappa <= af/e)
00223 lft(i) = -1.0;
00224 else
00225 {
00226 e = Mm->give_actual_ym(ipp);
00227
00228 af = Mm->give_actual_ft(ipp);
00229
00230 auf = uft*(af/ft);
00231 tmp1 = -(kappa-af/e)/(auf-af/e);
00232 lft(i) = 1.0 - af/(e*kappa)*exp(tmp1) - damt(i);
00233 }
00234 }
00235
00236 if (peps(i) < 0.0)
00237 {
00238 kappa = dam_eq_strain(nu, i, peps);
00239 e = Mm->give_actual_ym(ipp);
00240 af = Mm->give_actual_fc(ipp);
00241 if (kappa <= af/e)
00242 lfc(i) = -1.0;
00243 else
00244 {
00245 e = Mm->give_actual_ym(ipp);
00246
00247 af = Mm->give_actual_fc(ipp);
00248
00249 auf = ufc*(af/fc);
00250 tmp1 = -(kappa-af/e)/(auf-af/e);
00251 lfc(i) = 1.0 - af/(e*kappa)*exp(tmp1) - damc(i);
00252 }
00253 }
00254 }
00255 break;
00256 case quasi_brittle:
00257 for (i=0; i<3; i++)
00258 {
00259
00260 if (peps(i) >= 0.0)
00261 {
00262 kappa = dam_eq_strain(nu, i, peps);
00263 if (kappa <= y0t)
00264 lft(i) = -1.0;
00265 else
00266 {
00267
00268 tmp1 = kappa - y0t;
00269 tmp2 = log(tmp1);
00270 tmp1 = exp(bt*tmp2);
00271
00272 lft(i) = (1.0 - damt(i))*(1.0+aat*tmp1) - 1.0;
00273 }
00274 }
00275
00276 if (peps(i) < 0.0)
00277 {
00278 kappa = dam_eq_strain(nu, i, peps);
00279 if (kappa <= y0c)
00280 lfc(i) = -1.0;
00281 else
00282 {
00283
00284 tmp1 = kappa - y0c;
00285 tmp2 = log(tmp1);
00286 tmp1 = exp(bc*tmp2);
00287
00288 lfc(i) = (1.0 - damc(i))*(1.0+aac*tmp1) - 1.0;
00289 }
00290 }
00291 }
00292 break;
00293 case quadbezier:
00294 h = 1.0;
00295 for (i=0; i<3; i++)
00296 {
00297
00298 if (peps(i) >= 0.0)
00299 {
00300 kappa = dam_eq_strain(nu, i, peps);
00301 if (kappa <= y0t)
00302 lft(i) = -1.0;
00303 else
00304 {
00305 e = Mm->give_actual_ym(ipp);
00306
00307 af = Mm->give_actual_ft(ipp);
00308
00309 auf = uft*(af/ft);
00310
00311 aul = ult*(af/ft);
00312 a = 4.0*auf*auf;
00313 b = (aul-2.0*auf);
00314 if (pqt)
00315 lft(i) = e*kappa-damt(i)*e*kappa-af*(1.0-damt(i)*kappa*h/auf+sqr(damt(i)*kappa*h)*(aul-auf)/(auf*aul*aul));
00316 else
00317
00318 lft(i) = e*kappa-damt(i)*e*kappa-af*(1.0-auf*(2.0*aul-3.0*auf)/(b*b)-2.0*(aul-auf)/(b*b)*sqrt(a+b*damt(i)*h*kappa)+damt(i)*h*kappa/b);
00319 }
00320 }
00321
00322 if (peps(i) < 0.0)
00323 {
00324 kappa = dam_eq_strain(nu, i, peps);
00325 if (kappa <= y0c)
00326 lfc(i) = -1.0;
00327 else
00328 {
00329 e = Mm->give_actual_ym(ipp);
00330
00331 af = Mm->give_actual_fc(ipp);
00332
00333 auf = ufc*(af/fc);
00334
00335 aul = ult*(af/ft);
00336 a = 4.0*auf*auf;
00337 b = (aul-2.0*auf);
00338 if (pqc)
00339 lfc(i) = e*kappa-damc(i)*e*kappa-af*(1.0-damc(i)*kappa*h/auf+sqr(damc(i)*kappa*h)*(aul-auf)/(auf*aul*aul));
00340 else
00341
00342 lfc(i) = e*kappa-damc(i)*e*kappa-af*(1.0-auf*(2.0*aul-3.0*auf)/(b*b)-2.0*(aul-auf)/(b*b)*sqrt(a+b*damc(i)*h*kappa)+damc(i)*h*kappa/b);
00343 }
00344 }
00345 }
00346 break;
00347 default:
00348 print_err(" unknown type of damage evolution function is required", __FILE__, __LINE__,__func__);
00349 }
00350 return;
00351 }
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371 double ortodam::brittle_damage(long ipp, double y, double e, double f, double uf, double omegao)
00372 {
00373 double err,omega = 0.0, dtmp, tmp, h, indam, kappa;
00374 long i,ni, eid;
00375
00376 ni=sra.give_ni ();
00377 err=sra.give_err ();
00378
00379 indam = f/e;
00380
00381 kappa = y;
00382 if (kappa < indam)
00383
00384 return 0.0;
00385 if ((Mm->ip[ipp].hmt & 2) > 0)
00386 {
00387 omega = 1.0-(f/(e*kappa))*exp(-(kappa-f/e)/(uf-f/e));
00388 return omega;
00389 }
00390
00391 eid = Mm->elip[ipp];
00392 switch (Mt->give_dimension(eid))
00393 {
00394 case 1:
00395 h = Mt->give_length(eid);
00396 break;
00397 case 2:
00398 h = sqrt(Mt->give_area(eid));
00399 break;
00400 case 3:
00401 h = pow(Mt->give_volume(eid), 1.0/3.0);
00402 break;
00403 default:
00404 print_err("unknown dimension of element", __FILE__, __LINE__, __func__);
00405 }
00406 if (cde == corr_off)
00407
00408 {
00409 omega = 1.0-(f/(e*kappa))*exp(-(kappa-f/e)/(uf-f/e));
00410 return omega;
00411 }
00412
00413 omega = 1.0;
00414 for (i = 0; i < ni; i++)
00415 {
00416 dtmp = -e*kappa+f/uf*h*kappa*exp(-h*omega*kappa/uf);
00417 tmp = (1.0-omega)*e*kappa-f*exp(-h*omega*kappa/uf);
00418 if (fabs(tmp) < ft*err)
00419 {
00420 omega += - tmp/dtmp;
00421 break;
00422 }
00423 omega += - tmp/dtmp;
00424 }
00425 if ((omega > 1.0) || (omega < 0.0) || (fabs(tmp) > ft*err))
00426
00427
00428 {
00429 omega = omegao;
00430 for (i = 0; i < ni; i++)
00431 {
00432 dtmp = -e*kappa+ft/uf*h*kappa*exp(-h*omega*kappa/uf);
00433 tmp = (1.0-omega)*e*kappa-ft*exp(-h*omega*kappa/uf);
00434 if (fabs(tmp) < ft*err)
00435 {
00436 omega += - tmp/dtmp;
00437 break;
00438 }
00439 omega += - tmp/dtmp;
00440 }
00441 }
00442
00443 if (fabs(tmp) > ft*err)
00444 {
00445 print_err("iteration of omega doesn't converge with required error\n"
00446 " elem %ld, ip %ld:\n"
00447 " kappa=%le, ft=%le, e=%le, uf=%le, h=%le",
00448 __FILE__, __LINE__, __func__, Mm->elip[ipp]+1, ipp+1, kappa, ft, e, uf, h);
00449 abort();
00450 }
00451
00452 if ((omega > 1.0) || (omega < 0.0))
00453 {
00454 print_err("computed omega=%le is out of prescribed range <0,1>\n"
00455 " elem %ld, ip %ld:\n"
00456 " kappa=%le, ft=%le, e=%le, uf=%le, h=%le",
00457 __FILE__, __LINE__, __func__, Mm->elip[ipp]+1, ipp+1, kappa, ft, e, uf, h);
00458 abort();
00459 }
00460
00461 return omega;
00462 }
00463
00464
00465
00466
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484 double ortodam::qbezier_damage(long ipp, double y, double e, double f, double uf, double ul,double omegao,long pq)
00485 {
00486 double err,omega = 0.0, dtmp, tmp, h, indam, kappa, a, b;
00487 long i,ni, eid;
00488
00489 ni=sra.give_ni ();
00490 err=sra.give_err ();
00491
00492 indam = f/e;
00493
00494 kappa = y;
00495 if (kappa-indam <= Mp->zero)
00496
00497 return 0.0;
00498 a = uf*uf;
00499 b = (ul-2.0*uf);
00500 if ((Mm->ip[ipp].hmt & 2) > 0)
00501 {
00502 if (kappa >= ul)
00503 return 1.0;
00504 omega = 1.0-2.0*uf*(ul-uf)/(b*b)-2.0*(ul-uf)*sqrt(a+b*(kappa-f/e))/(b*b)+(kappa-f/e)/b;
00505 omega = 1.0-(f/(e*(kappa-f/e)))*omega;
00506 return omega;
00507 }
00508
00509 eid = Mm->elip[ipp];
00510 switch (Mt->give_dimension(eid))
00511 {
00512 case 1:
00513 h = Mt->give_length(eid);
00514 break;
00515 case 2:
00516 h = sqrt(Mt->give_area(eid));
00517 break;
00518 case 3:
00519 h = pow(Mt->give_volume(eid), 1.0/3.0);
00520 break;
00521 default:
00522 print_err("unknown dimension of element", __FILE__, __LINE__, __func__);
00523 }
00524 if (cde == corr_off)
00525
00526 {
00527 if (kappa >= ul)
00528 return 1.0;
00529 if (pq)
00530 omega = (1.0-omega)*e*kappa-f*(1.0-(kappa-f/e)/uf+(ul-uf)/(uf*ul*ul)*(kappa-f/e)*(kappa-f/e));
00531 else
00532 {
00533 omega = 1.0-2.0*uf*(ul-uf)/(b*b)-2.0*(ul-uf)*sqrt(a+b*(kappa-f/e))/(b*b)+(kappa-f/e)/b;
00534 omega = 1.0-(f/(e*(kappa-f/e)))*omega;
00535 }
00536 return omega;
00537 }
00538 if (omega < 0.5)
00539 omega = 0.5;
00540 else
00541 omega = omegao;
00542 if (omega < 0.0)
00543 print_err("previous iteration of omega returned value < 0", __FILE__, __LINE__, __func__);
00544
00545 if (kappa*omegao*h >= ul)
00546 return 1.0;
00547 for (i = 0; i < ni; i++)
00548 {
00549 if (pq)
00550 {
00551 dtmp = 2.0*f*(uf-ul)*sqr(kappa*h)*omega/(uf*ul*ul)+kappa*h*f/uf-e*kappa;
00552 tmp = (1.0-omega)*e*kappa-f*(1.0-omega*kappa*h/uf+(ul-uf)*sqr(omega*kappa*h)/(uf*ul*ul));
00553 }
00554 else
00555 {
00556 dtmp = f*(ul-uf)*kappa*h/(b*sqrt(a+omega*b*kappa*h))-(e+f*h/b)*kappa;
00557 tmp = (1.0-omega)*e*kappa-f*(1.0+(2.0*uf*(ul-uf))/(b*b)-(2.0*(ul-uf)*sqrt(a+b*omega*h*kappa))/(b*b)+omega*h*kappa/b);
00558 }
00559 omega += - tmp/dtmp;
00560 if (omega < 0.0)
00561 {
00562 print_err("during iteration, omega reaches value < 0", __FILE__, __LINE__, __func__);
00563 abort();
00564 }
00565 if (fabs(tmp) <= ft*err)
00566 break;
00567 }
00568 if (fabs(tmp) > ft*err)
00569 print_err("iteration of omega doesn't converge with required error", __FILE__, __LINE__, __func__);
00570
00571 if (omega > 1.0)
00572 {
00573 if ((omega-1.0) > err)
00574 print_err("iteration of returns value > 1", __FILE__, __LINE__, __func__);
00575 else
00576 omega = 1.0;
00577 }
00578 if (omega < 0.0)
00579 {
00580 if (fabs(omega) > err)
00581 print_err("iteration of omega returns value < 0", __FILE__, __LINE__, __func__);
00582 else
00583 omega = 0.0;
00584 }
00585
00586
00587 return omega;
00588 }
00589
00590
00591
00592
00593
00594
00595
00596
00597
00598
00599
00600
00601
00602
00603
00604
00605
00606
00607 void ortodam::princ_dam(long ipp, vector &peps, double aat, double aac, vector &pdamt, vector &pdamc)
00608 {
00609 long i;
00610 double tmp, ltmp;
00611 double e, nu, af, auf, aul;
00612 double kappa;
00613
00614 switch(damevf)
00615 {
00616 case brittle:
00617 e = Mm->give_actual_ym(ipp);
00618 nu = Mm->give_actual_nu(ipp);
00619 for(i=0; i<3; i++)
00620 {
00621 kappa = dam_eq_strain(nu, i, peps);
00622
00623 if (peps(i) >= 0.0)
00624 {
00625
00626 af = Mm->give_actual_ft(ipp);
00627
00628 auf = uft*(af/ft);
00629 pdamt(i) = brittle_damage(ipp, kappa, e, af, auf, pdamt(i));
00630 }
00631
00632 if (peps(i) < 0.0)
00633 {
00634
00635 af = Mm->give_actual_fc(ipp);
00636
00637 auf = ufc*(af/fc);
00638 pdamc(i) = brittle_damage(ipp, kappa, e, af, auf, pdamc(i));
00639 }
00640 }
00641 break;
00642 case quasi_brittle:
00643 nu = Mm->give_actual_nu(ipp);
00644 for(i=0; i<3; i++)
00645 {
00646 kappa = dam_eq_strain(nu, i, peps);
00647
00648 if (peps(i) >= 0.0)
00649 {
00650 if (kappa > y0t)
00651 {
00652 ltmp = log(kappa - y0t);
00653 pdamt(i) = aat*exp(ltmp*bt);
00654 tmp = 1.0/(1.0+pdamt(i));
00655 pdamt(i) = 1.0 - tmp;
00656 }
00657 if (pdamt(i) < 0.0)
00658 {
00659 fprintf(stdout, "\nDamage < 0.0 (ipp=%ld)\n", ipp);
00660 abort();
00661 }
00662 }
00663
00664 if (peps(i) < 0.0)
00665 {
00666 if (kappa > y0c)
00667 {
00668 ltmp = log(kappa - y0c);
00669 pdamc(i) = aac*exp(ltmp*bc);
00670 tmp = 1.0/(1.0+pdamc(i));
00671 pdamc(i) = 1.0-tmp;
00672 }
00673 if (pdamc(i) < 0.0)
00674 {
00675 fprintf(stdout, "\nDamage < 0.0 (ipp=%ld)\n", ipp);
00676 abort();
00677 }
00678 }
00679 }
00680 break;
00681 case quadbezier:
00682 e = Mm->give_actual_ym(ipp);
00683 nu = Mm->give_actual_nu(ipp);
00684 for(i=0; i<3; i++)
00685 {
00686 kappa = dam_eq_strain(nu, i, peps);
00687
00688 if (kappa >= 0.0)
00689 {
00690
00691 af = Mm->give_actual_ft(ipp);
00692
00693 auf = uft*(af/ft);
00694
00695 aul = ult*(af/ft);
00696 pdamt(i) = qbezier_damage(ipp, fabs(kappa), e, af, auf, aul, pdamt(i),pqt);
00697 }
00698
00699 if (kappa < 0.0)
00700 {
00701
00702 af = Mm->give_actual_fc(ipp);
00703
00704 auf = ufc*(af/fc);
00705
00706 aul = ulc*(af/fc);
00707 pdamc(i) = qbezier_damage(ipp, fabs(kappa), e, af, auf, aul, pdamc(i),pqc);
00708 }
00709 }
00710 break;
00711 default:
00712 print_err(" unknown type of damage evolution function is required", __FILE__, __LINE__,__func__);
00713 }
00714 }
00715
00716
00717
00718
00719
00720
00721
00722
00723
00724
00725
00726
00727
00728
00729
00730
00731
00732
00733
00734
00735 void ortodam::give_actual_param_a(long ipp, long ido, double &aat, double &aac)
00736 {
00737 long eid;
00738 double e, h, tmp;
00739
00740 aat = Mm->ip[ipp].eqother[ido];
00741 aac = Mm->ip[ipp].eqother[ido+1];
00742 if ((aac == 0.0) || (aat == 0.0))
00743 {
00744 if (cde == corr_off)
00745
00746 {
00747 aat = at;
00748 aac = ac;
00749 Mm->ip[ipp].eqother[ido] = aat;
00750 Mm->ip[ipp].eqother[ido+1] = aac;
00751 return;
00752 }
00753 eid = Mm->elip[ipp];
00754 switch (Mt->give_dimension(eid))
00755 {
00756 case 1:
00757 h = Mt->give_length(eid);
00758 break;
00759 case 2:
00760 h = sqrt(Mt->give_area(eid));
00761 break;
00762 case 3:
00763 h = pow(Mt->give_volume(eid), 1.0/3.0);
00764 break;
00765 default:
00766 print_err("cannot determine dimension of the element.", __FILE__, __LINE__, __func__);
00767 }
00768 e = Mm->give_actual_ym(ipp);
00769
00770 tmp = (gft/h - y0t/e)*bt*e*sin(M_PI/bt)/M_PI;
00771 if (tmp < 0.0)
00772 print_err("cannot express variable softening for tensile damage.", __FILE__, __LINE__, __func__);
00773 aat = pow(tmp, -bt);
00774
00775 tmp = (gfc/h - y0c/e)*bc*e*sin(M_PI/bc)/M_PI;
00776 if (tmp < 0.0)
00777 print_err("cannot express variable softening for compressive damage.", __FILE__, __LINE__, __func__);
00778 aac = pow(tmp, -bc);
00779 Mm->ip[ipp].eqother[ido] = aat;
00780 Mm->ip[ipp].eqother[ido+1] = aac;
00781 }
00782 }
00783
00784
00785
00786
00787
00788
00789
00790
00791
00792
00793
00794
00795
00796 double ortodam::give_actual_ft(long ipp)
00797 {
00798 return ft;
00799 }
00800
00801
00802
00803
00804
00805
00806
00807
00808
00809
00810
00811
00812
00813 double ortodam::give_actual_fc(long ipp)
00814 {
00815 return fc;
00816 }
00817
00818
00819
00820
00821
00822
00823
00824
00825
00826
00827
00828
00829
00830
00831 void ortodam::initvalues(long ipp, long ido)
00832 {
00833 double aat=0.0, aac=0.0;
00834 if (damevf == quasi_brittle)
00835 give_actual_param_a(ipp, ido, aac, aat);
00836 }
00837
00838
00839
00840
00841
00842
00843
00844
00845
00846
00847
00848
00849
00850
00851 void ortodam::matstiff (matrix &d,long ipp,long ido)
00852 {
00853
00854
00855 switch (Mp->nlman->stmat)
00856 {
00857 case initial_stiff:
00858 elmatstiff(d, ipp);
00859 break;
00860 case tangent_stiff:
00861 elmatstiff(d, ipp);
00862
00863 break;
00864 default:
00865 print_err("unknown type of stifness matrix is required.", __FILE__, __LINE__, __func__);
00866 }
00867 }
00868
00869
00870
00871
00872
00873
00874
00875
00876
00877
00878
00879
00880
00881 void ortodam::elmatstiff (matrix &d,long ipp)
00882 {
00883 double e, e0;
00884 long idem = Mm->ip[ipp].gemid();
00885 long idoem=Mm->givencompeqother(ipp,0);
00886 long tmp=Mm->givencompeqother(ipp,idem);
00887 idoem -= tmp;
00888
00889 Mm->elmatstiff (d,ipp);
00890 e = Mm->give_actual_ym(ipp);
00891 e0 = Mm->give_actual_ym(ipp,idem,idoem);
00892 cmulm(e/e0, d);
00893 }
00894
00895
00896
00897
00898
00899
00900
00901
00902
00903
00904
00905
00906
00907
00908
00909
00910 void ortodam::pelmatstiff (long ipp, matrix &d, double e, double nu)
00911 {
00912 double b;
00913
00914 switch (Mm->ip[ipp].ssst)
00915 {
00916 case planestress:
00917 b = e/(1.0-nu*nu);
00918 d[0][0] = b; d[0][1] = b*nu; d[0][2] = b*nu;
00919 d[1][0] = b*nu; d[1][1] = b; d[1][2] = b*nu;
00920 d[2][0] = b*nu; d[2][1] = b*nu; d[2][2] = b;
00921 break;
00922 case planestrain:
00923 case spacestress:
00924 b = e /((1.0 + nu)*(1.0 - 2.0*nu));
00925 d[0][0] = d[1][1] = d[2][2] = 1.0 - nu;
00926 d[0][1] = d[0][2] = d[1][0] = d[1][2] = d[2][0] = d[2][1] = nu;
00927 cmulm(b, d);
00928 break;
00929 default:
00930 print_err("Unknown type of stress/strain state is required", __FILE__, __LINE__, __func__);
00931 break;
00932 }
00933 }
00934
00935
00936
00937
00938
00939
00940
00941
00942
00943
00944
00945
00946
00947
00948
00949 void ortodam::nlstresses (long ipp, long im, long ido)
00950 {
00951 long i;
00952 long ncompstr=Mm->ip[ipp].ncompstr;
00953
00954
00955 double aat=0.0, aac=0.0;
00956
00957
00958 vector pyto(3), pyco(3);
00959
00960 vector pdto(3), pdco(3);
00961
00962 vector pdt(3), pdc(3);
00963
00964 vector epsn(ncompstr);
00965
00966 matrix epst(3,3), t(3,3);
00967
00968 matrix sigt(3,3), d(3,3);
00969 vector sigma(ncompstr);
00970
00971 double e, nu;
00972
00973
00974 vector peps(3), psig(3), lft(3), lfc(3),pepsp(3);
00975 double h;
00976 long eid;
00977
00978 if ((Mp->matmodel == nonlocal) && (Mm->ip[ipp].hmt & 2))
00979
00980 {
00981 if (Mp->nonlocphase == 1)
00982 {
00983
00984
00985 return;
00986 }
00987 else
00988 {
00989
00990 for (i=0;i<ncompstr;i++)
00991 epsn[i] = Mm->ip[ipp].nonloc[i];
00992 }
00993 }
00994 else
00995
00996 {
00997
00998 for (i=0;i<ncompstr;i++)
00999 epsn[i] = Mm->ip[ipp].strain[i];
01000 }
01001
01002
01003 if (damevf == quasi_brittle)
01004 give_actual_param_a(ipp, ido, aat, aac);
01005 e = Mm->give_actual_ym(ipp);
01006 nu = Mm->give_actual_nu(ipp);
01007 if (Mm->ip[ipp].ssst == planestress)
01008 Mm->ip[ipp].strain[3] = -nu / (1.0 - nu) * (Mm->ip[ipp].strain[0]+Mm->ip[ipp].strain[1]);
01009
01010
01011 for (i=0;i<ncompstr;i++)
01012 epsn[i] = Mm->ip[ipp].strain[i];
01013
01014
01015 for (i=0; i<3; i++)
01016 {
01017
01018 pdt[i] = pdto[i] = Mm->ip[ipp].eqother[ido+2+i];
01019 pdc[i] = pdco[i] = Mm->ip[ipp].eqother[ido+2+3+i];
01020 }
01021
01022
01023 vector_tensor(epsn, epst, Mm->ip[ipp].ssst, strain);
01024 princ_val (epst,peps,t,nijac,limit,Mp->zero,3, 1);
01025
01026
01027 if (fat == fatigue_on)
01028 {
01029 for (i=0; i<3; i++)
01030 {
01031 pepsp[i] = 0.5*betat*pdto[i]*pdto[i]/(1.0-pdto[i]);
01032 pepsp[i] -= 0.5*betac*pdco[i]*pdco[i]/(1.0-pdco[i]);
01033 }
01034 subv(peps,pepsp,peps);
01035 }
01036
01037 princ_dam(ipp, peps, aat, aac, pdt, pdc);
01038 for (i=0; i<3; i++)
01039 {
01040 if (pdt[i] < pdto[i])
01041 pdt[i] = pdto[i];
01042 if (pdc[i] < pdco[i])
01043 pdc[i] = pdco[i];
01044 }
01045
01046
01047 if (fat == fatigue_on)
01048 {
01049
01050
01051
01052
01053
01054
01055 }
01056
01057 pelmatstiff (ipp, d, e, nu);
01058 mxv(d, peps, psig);
01059 for (i=0; i<3; i++)
01060 {
01061 if (psig(i) > 0.0)
01062 psig(i) *= 1.0-pdt(i);
01063 if (psig(i) < 0.0)
01064 psig(i) *= 1.0-pdc(i);
01065 }
01066
01067
01068 fillm(0.0, sigt);
01069 for (i = 0; i < 3; i++)
01070 sigt[i][i] = psig[i];
01071 lgmatrixtransf(sigt, t);
01072 tensor_vector(sigma, sigt, Mm->ip[ipp].ssst, stress);
01073 if (Mm->ip[ipp].ssst == planestress)
01074 sigma[3] = 0.0;
01075
01076
01077 for (i=0;i<ncompstr;i++)
01078 Mm->ip[ipp].stress[i]=sigma[i];
01079 for (i=0; i<3; i++)
01080 {
01081
01082 Mm->ip[ipp].other[ido+2+i] = pdt[i];
01083 Mm->ip[ipp].other[ido+2+3+i] = pdc[i];
01084 }
01085
01086
01087 fillm(0.0, sigt);
01088 for (i = 0; i < 3; i++)
01089 sigt[i][i] = pdt[i];
01090 lgmatrixtransf(sigt, t);
01091 for (i = 0; i < 3; i++)
01092 Mm->ip[ipp].other[ido+8+i] = sigt[i][i];
01093 Mm->ip[ipp].other[ido+11] = sigt[1][2];
01094 Mm->ip[ipp].other[ido+12] = sigt[0][2];
01095 Mm->ip[ipp].other[ido+13] = sigt[0][1];
01096
01097 eid = Mm->elip[ipp];
01098 switch (Mt->give_dimension(eid))
01099 {
01100 case 1:
01101 h = Mt->give_length(eid);
01102 break;
01103 case 2:
01104 h = sqrt(Mt->give_area(eid));
01105 break;
01106 case 3:
01107 h = pow(Mt->give_volume(eid), 1.0/3.0);
01108 break;
01109 default:
01110 print_err("unknown dimension of element", __FILE__, __LINE__, __func__);
01111 }
01112
01113 for (i = 0; i < 3; i++)
01114 {
01115 if (peps(i) > 0.0)
01116 Mm->ip[ipp].other[ido+14+i] = pdt[i]*(peps[i]-psig[i]/e)*h;
01117 else
01118 Mm->ip[ipp].other[ido+14+i] = 0.0;
01119 }
01120
01121
01122 }
01123
01124
01125
01126
01127
01128
01129
01130
01131
01132
01133
01134
01135
01136 void ortodam::updateval (long ipp,long im,long ido)
01137 {
01138 long i;
01139 long n = Mm->givencompeqother(ipp,im);
01140
01141 for (i=0;i<n;i++)
01142 Mm->ip[ipp].eqother[ido+i]=Mm->ip[ipp].other[ido+i];
01143 }