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-1)%3));
00172 kappa /= (1-2*nu)*(1+nu);
00173 kappa = fabs(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
00615 switch(damevf)
00616 {
00617 case brittle:
00618 e = Mm->give_actual_ym(ipp);
00619 nu = Mm->give_actual_nu(ipp);
00620 for(i=0; i<3; i++)
00621 {
00622 kappa = dam_eq_strain(nu, i, peps);
00623
00624 if (peps(i) >= 0.0)
00625 {
00626
00627 af = Mm->give_actual_ft(ipp);
00628
00629 auf = uft*(af/ft);
00630 pdamt(i) = brittle_damage(ipp, kappa, e, af, auf, pdamt(i));
00631 }
00632
00633 if (peps(i) < 0.0)
00634 {
00635
00636 af = Mm->give_actual_fc(ipp);
00637
00638 auf = ufc*(af/fc);
00639 pdamc(i) = brittle_damage(ipp, kappa, e, af, auf, pdamc(i));
00640 }
00641 }
00642 break;
00643 case quasi_brittle:
00644 nu = Mm->give_actual_nu(ipp);
00645 for(i=0; i<3; i++)
00646 {
00647 kappa = dam_eq_strain(nu, i, peps);
00648
00649 if (peps(i) >= 0.0)
00650 {
00651 if (kappa > y0t)
00652 {
00653 ltmp = log(kappa - y0t);
00654 pdamt(i) = aat*exp(ltmp*bt);
00655 tmp = 1.0/(1.0+pdamt(i));
00656 pdamt(i) = 1.0 - tmp;
00657 }
00658 if (pdamt(i) < 0.0)
00659 {
00660 fprintf(stdout, "\nDamage < 0.0 (ipp=%ld)\n", ipp);
00661 abort();
00662 }
00663 }
00664
00665 if (peps(i) < 0.0)
00666 {
00667 if (kappa > y0c)
00668 {
00669 ltmp = log(kappa - y0c);
00670 pdamc(i) = aac*exp(ltmp*bc);
00671 tmp = 1.0/(1.0+pdamc(i));
00672 pdamc(i) = 1.0-tmp;
00673 }
00674 if (pdamc(i) < 0.0)
00675 {
00676 fprintf(stdout, "\nDamage < 0.0 (ipp=%ld)\n", ipp);
00677 abort();
00678 }
00679 }
00680 }
00681 break;
00682 case quadbezier:
00683 e = Mm->give_actual_ym(ipp);
00684 nu = Mm->give_actual_nu(ipp);
00685 for(i=0; i<3; i++)
00686 {
00687 kappa = dam_eq_strain(nu, i, peps);
00688
00689 if (peps(i) >= 0.0)
00690 {
00691
00692 af = Mm->give_actual_ft(ipp);
00693
00694 auf = uft*(af/ft);
00695
00696 aul = ult*(af/ft);
00697 pdamt(i) = qbezier_damage(ipp, kappa, e, af, auf, aul, pdamt(i),pqt);
00698 }
00699
00700 if (peps(i) < 0.0)
00701 {
00702
00703 af = Mm->give_actual_fc(ipp);
00704
00705 auf = ufc*(af/fc);
00706
00707 aul = ulc*(af/fc);
00708 pdamc(i) = qbezier_damage(ipp, kappa, e, af, auf, aul, pdamc(i),pqc);
00709 }
00710 }
00711 break;
00712 default:
00713 print_err(" unknown type of damage evolution function is required", __FILE__, __LINE__,__func__);
00714 }
00715 }
00716
00717
00718
00719
00720
00721
00722
00723
00724
00725
00726
00727
00728
00729
00730
00731
00732
00733
00734
00735
00736 void ortodam::give_actual_param_a(long ipp, long ido, double &aat, double &aac)
00737 {
00738 long eid;
00739 double e, h, tmp;
00740
00741 aat = Mm->ip[ipp].eqother[ido];
00742 aac = Mm->ip[ipp].eqother[ido+1];
00743 if ((aac == 0.0) || (aat == 0.0))
00744 {
00745 if (cde == corr_off)
00746
00747 {
00748 aat = at;
00749 aac = ac;
00750 Mm->ip[ipp].eqother[ido] = aat;
00751 Mm->ip[ipp].eqother[ido+1] = aac;
00752 return;
00753 }
00754 eid = Mm->elip[ipp];
00755 switch (Mt->give_dimension(eid))
00756 {
00757 case 1:
00758 h = Mt->give_length(eid);
00759 break;
00760 case 2:
00761 h = sqrt(Mt->give_area(eid));
00762 break;
00763 case 3:
00764 h = pow(Mt->give_volume(eid), 1.0/3.0);
00765 break;
00766 default:
00767 print_err("cannot determine dimension of the element.", __FILE__, __LINE__, __func__);
00768 }
00769 e = Mm->give_actual_ym(ipp);
00770
00771 tmp = (gft/h - y0t/e)*bt*e*sin(M_PI/bt)/M_PI;
00772 if (tmp < 0.0)
00773 print_err("cannot express variable softening for tensile damage.", __FILE__, __LINE__, __func__);
00774 aat = pow(tmp, -bt);
00775
00776 tmp = (gfc/h - y0c/e)*bc*e*sin(M_PI/bc)/M_PI;
00777 if (tmp < 0.0)
00778 print_err("cannot express variable softening for compressive damage.", __FILE__, __LINE__, __func__);
00779 aac = pow(tmp, -bc);
00780 Mm->ip[ipp].eqother[ido] = aat;
00781 Mm->ip[ipp].eqother[ido+1] = aac;
00782 }
00783 }
00784
00785
00786
00787
00788
00789
00790
00791
00792
00793
00794
00795
00796
00797 double ortodam::give_actual_ft(long ipp)
00798 {
00799 return ft;
00800 }
00801
00802
00803
00804
00805
00806
00807
00808
00809
00810
00811
00812
00813
00814 double ortodam::give_actual_fc(long ipp)
00815 {
00816 return fc;
00817 }
00818
00819
00820
00821
00822
00823
00824
00825
00826
00827
00828
00829
00830
00831
00832 void ortodam::initvalues(long ipp, long ido)
00833 {
00834 double aat, aac;
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, aac;
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
00979 give_actual_param_a(ipp, ido, aat, aac);
00980 e = Mm->give_actual_ym(ipp);
00981 nu = Mm->give_actual_nu(ipp);
00982 if (Mm->ip[ipp].ssst == planestress)
00983 Mm->ip[ipp].strain[3] = -nu / (1.0 - nu) * (Mm->ip[ipp].strain[0]+Mm->ip[ipp].strain[1]);
00984
00985
00986 for (i=0;i<ncompstr;i++)
00987 epsn[i] = Mm->ip[ipp].strain[i];
00988
00989
00990 for (i=0; i<3; i++)
00991 {
00992
00993 pdt[i] = pdto[i] = Mm->ip[ipp].eqother[ido+2+i];
00994 pdc[i] = pdco[i] = Mm->ip[ipp].eqother[ido+2+3+i];
00995 }
00996
00997
00998 vector_tensor(epsn, epst, Mm->ip[ipp].ssst, strain);
00999 princ_val (epst,peps,t,nijac,limit,Mp->zero,3, 1);
01000
01001
01002 if (fat == fatigue_on)
01003 {
01004 for (i=0; i<3; i++)
01005 {
01006 pepsp[i] = 0.5*betat*pdto[i]*pdto[i]/(1.0-pdto[i]);
01007 pepsp[i] -= 0.5*betac*pdco[i]*pdco[i]/(1.0-pdco[i]);
01008 }
01009 subv(peps,pepsp,peps);
01010 }
01011
01012 princ_dam(ipp, peps, aat, aac, pdt, pdc);
01013 for (i=0; i<3; i++)
01014 {
01015 if (pdt[i] < pdto[i])
01016 pdt[i] = pdto[i];
01017 if (pdc[i] < pdco[i])
01018 pdc[i] = pdco[i];
01019 }
01020
01021
01022 if (fat == fatigue_on)
01023 {
01024
01025
01026
01027
01028
01029
01030 }
01031
01032 pelmatstiff (ipp, d, e, nu);
01033 mxv(d, peps, psig);
01034 for (i=0; i<3; i++)
01035 {
01036 if (psig(i) > 0.0)
01037 psig(i) *= 1.0-pdt(i);
01038 if (psig(i) < 0.0)
01039 psig(i) *= 1.0-pdc(i);
01040 }
01041
01042
01043 fillm(0.0, sigt);
01044 for (i = 0; i < 3; i++)
01045 sigt[i][i] = psig[i];
01046 lgmatrixtransf(sigt, t);
01047 tensor_vector(sigma, sigt, Mm->ip[ipp].ssst, stress);
01048 if (Mm->ip[ipp].ssst == planestress)
01049 sigma[3] = 0.0;
01050
01051
01052 for (i=0;i<ncompstr;i++)
01053 Mm->ip[ipp].stress[i]=sigma[i];
01054 for (i=0; i<3; i++)
01055 {
01056
01057 Mm->ip[ipp].other[ido+2+i] = pdt[i];
01058 Mm->ip[ipp].other[ido+2+3+i] = pdc[i];
01059 }
01060
01061
01062 fillm(0.0, sigt);
01063 for (i = 0; i < 3; i++)
01064 sigt[i][i] = pdt[i];
01065 lgmatrixtransf(sigt, t);
01066 for (i = 0; i < 3; i++)
01067 Mm->ip[ipp].other[ido+8+i] = sigt[i][i];
01068 Mm->ip[ipp].other[ido+11] = sigt[1][2];
01069 Mm->ip[ipp].other[ido+12] = sigt[0][2];
01070 Mm->ip[ipp].other[ido+13] = sigt[0][1];
01071
01072 eid = Mm->elip[ipp];
01073 switch (Mt->give_dimension(eid))
01074 {
01075 case 1:
01076 h = Mt->give_length(eid);
01077 break;
01078 case 2:
01079 h = sqrt(Mt->give_area(eid));
01080 break;
01081 case 3:
01082 h = pow(Mt->give_volume(eid), 1.0/3.0);
01083 break;
01084 default:
01085 print_err("unknown dimension of element", __FILE__, __LINE__, __func__);
01086 }
01087
01088 for (i = 0; i < 3; i++)
01089 {
01090 if (peps(i) > 0.0)
01091 Mm->ip[ipp].other[ido+14+i] = pdt[i]*(peps[i]-psig[i]/e)*h;
01092 else
01093 Mm->ip[ipp].other[ido+14+i] = 0.0;
01094 }
01095
01096
01097 }
01098
01099
01100
01101
01102
01103
01104
01105
01106
01107
01108
01109
01110
01111 void ortodam::updateval (long ipp,long im,long ido)
01112 {
01113 long i;
01114 long n = Mm->givencompeqother(ipp,im);
01115
01116 for (i=0;i<n;i++)
01117 Mm->ip[ipp].eqother[ido+i]=Mm->ip[ipp].other[ido+i];
01118 }