00001 #include <math.h>
00002 #include <stdio.h>
00003 #include <stdlib.h>
00004 #include "drprag2.h"
00005 #include "global.h"
00006 #include "vecttens.h"
00007 #include "intpoints.h"
00008 #include "matrix.h"
00009 #include "vector.h"
00010 #include "tensor.h"
00011
00012
00013
00014 #define nijac 20
00015 #define limit 1.0e-8
00016
00017
00018
00019
00020
00021
00022
00023 drprag2::drprag2 (void)
00024 {
00025 phi=0.0; c=0.0; psi=0.0;
00026 theta = 0.0; clim = 0.0;
00027 }
00028
00029
00030
00031
00032
00033
00034 drprag2::~drprag2 (void)
00035 {
00036
00037 }
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050 void drprag2::read (XFILE *in)
00051 {
00052 xfscanf (in,"%k%lf%k%lf%k%lf%k%lf%k%lf","phi",&phi,"coh",&c,"psi",&psi,"theta",&theta,"clim",&clim);
00053 sra.read (in);
00054
00055 alpha=6.0*sin(phi)/(3.0-sin(phi));
00056
00057 alpha1=6.0*sin(psi)/(3.0-sin(psi));
00058
00059 beta=6.0*cos(phi)/(3.0-sin(phi));
00060
00061 }
00062
00063
00064
00065 void drprag2::giveIota(strastrestate ssst, vector &iota)
00066 {
00067 fillv(0.0, iota);
00068 switch (ssst)
00069 {
00070 case bar:
00071 iota[0]=1.0;
00072 break;
00073 case planestrain:
00074 iota[0]=iota[1]=iota[3]=1.0;
00075 break;
00076 case axisymm:
00077 case spacestress:
00078 iota[0]=iota[1]=iota[2]=1.0;
00079 break;
00080 default:
00081 print_err("stress/strain state %d is not supported in this model", __FILE__, __LINE__, __func__, ssst);
00082 }
00083 }
00084
00085
00086
00087 void drprag2::giveIdev(strastrestate ssst, matrix &idev)
00088 {
00089 fillm(0.0, idev);
00090 switch (ssst)
00091 {
00092 case bar:
00093 idev[0][0]=1.0;
00094 break;
00095 case planestrain:
00096 idev[0][0] = 2.0/3.0; idev[0][1] = -1.0/3.0; idev[0][3] = -1.0/3.0;
00097 idev[1][0] = -1.0/3.0; idev[1][1] = 2.0/3.0; idev[1][3] = -1.0/3.0;
00098 idev[2][2] = 0.5;
00099 idev[3][0] = -1.0/3.0; idev[3][1] = -1.0/3.0; idev[3][3] = 2.0/3.0;
00100 break;
00101 case axisymm:
00102 idev[0][0] = 2.0/3.0; idev[0][1] = -1.0/3.0; idev[0][2] = -1.0/3.0;
00103 idev[1][0] = -1.0/3.0; idev[1][1] = 2.0/3.0; idev[1][2] = -1.0/3.0;
00104 idev[2][0] = -1.0/3.0; idev[2][1] = -1.0/3.0; idev[2][2] = 2.0/3.0;
00105 idev[3][3] = 0.5;
00106 break;
00107 case spacestress:
00108 idev[0][0] = 2.0/3.0; idev[0][1] = -1.0/3.0; idev[0][2] = -1.0/3.0;
00109 idev[1][0] = -1.0/3.0; idev[1][1] = 2.0/3.0; idev[1][2] = -1.0/3.0;
00110 idev[2][0] = -1.0/3.0; idev[2][1] = -1.0/3.0; idev[2][2] = 2.0/3.0;
00111 idev[3][3] = 0.5;
00112 idev[4][4] = 0.5;
00113 idev[5][5] = 0.5;
00114 break;
00115 default:
00116 print_err("stress/strain state %d is not supported in this model", __FILE__, __LINE__, __func__, ssst);
00117 }
00118 }
00119
00120
00121
00122 void drprag2::getImportantValues(long ipp, long ido, double &K_con, double &G_con, double &rho_tr, double &p_tr)
00123 {
00124
00125 double E, nu;
00126
00127 double scalar_prod1, scalar_prod2;
00128
00129
00130 long n=Mm->ip[ipp].ncompstr;
00131 strastrestate ssst = Mm->ip[ipp].ssst;
00132 vector epsilon(n), epsilon_p(n), epsilon_etr(n), aux_epsilon_etr1(n);
00133 vector iota(n);
00134 matrix Idev(n,n);
00135 vector aux_epsilon_etr2(n);
00136
00137 E = Mm->give_actual_ym(ipp);
00138 nu = Mm->give_actual_nu(ipp);
00139 K_con = E/(3*(1-2*nu));
00140 G_con = E/(2*(1+nu));
00141
00142 giveIdev(ssst, Idev);
00143
00144
00145 giveirrstrains (ipp, ido, epsilon_p);
00146 Mm->givestrain(0, ipp, epsilon);
00147 subv(epsilon, epsilon_p, epsilon_etr);
00148
00149
00150 mxv(Idev, epsilon_etr, aux_epsilon_etr1);
00151 scprd(epsilon_etr, aux_epsilon_etr1, scalar_prod1);
00152 rho_tr = 2.0*G_con*sqrt(scalar_prod1);
00153
00154
00155 giveIota(ssst, iota);
00156 scprd(epsilon_etr, iota, scalar_prod2);
00157 p_tr = K_con*scalar_prod2;
00158
00159 }
00160
00161
00162 void drprag2::getImportantValues(long ipp, long ido, double &K_con, double &G_con, double &rho_tr, double &p_tr,
00163 vector &iota, vector &epsilon_etr, matrix &Idev)
00164 {
00165
00166 double E, nu;
00167
00168 double scalar_prod1, scalar_prod2;
00169
00170
00171 long n=Mm->ip[ipp].ncompstr;
00172 strastrestate ssst = Mm->ip[ipp].ssst;
00173 vector epsilon(n), epsilon_p(n), aux_epsilon_etr1(n);
00174
00175 E = Mm->give_initial_ym(ipp);
00176 nu = Mm->give_initial_nu(ipp);
00177 K_con = E/(3*(1-2*nu));
00178 G_con = E/(2*(1+nu));
00179
00180
00181 giveIdev(ssst, Idev);
00182
00183 giveirrstrains (ipp, ido, epsilon_p);
00184 Mm->givestrain(0, ipp, epsilon);
00185 subv(epsilon, epsilon_p, epsilon_etr);
00186
00187
00188 mxv(Idev, epsilon_etr, aux_epsilon_etr1);
00189 scprd(epsilon_etr, aux_epsilon_etr1, scalar_prod1);
00190 rho_tr = 2.0*G_con*sqrt(scalar_prod1);
00191
00192
00193
00194 giveIota(ssst, iota);
00195 scprd(epsilon_etr, iota, scalar_prod2);
00196 p_tr = K_con*scalar_prod2;
00197 }
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212 double drprag2::getQvalue(double gamma, long ipp, long ido)
00213 {
00214 double K_con, G_con, rho_tr, p_tr, q_value;
00215
00216
00217
00218 double H_value, Hder_value;
00219 double part_q1, part_q2, part_q3;
00220 getImportantValues(ipp, ido, K_con, G_con, rho_tr, p_tr);
00221
00222
00223
00224
00225 getValueOfHardening(ipp, ido, H_value, Hder_value);
00226
00227 part_q1 = sqrt(0.5)*max((rho_tr-gamma*G_con*sqrt(2.0)),0);
00228 part_q2 = alpha*(p_tr - gamma*K_con*alpha1);
00229 part_q3 = beta*(c+H_value);
00230
00231 q_value = part_q1 + part_q2 - part_q3;
00232
00233
00234 return q_value;
00235
00236 }
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251 void drprag2::getValueOfHardening(long ipp, long ido, double &H_value, double &Hder_value)
00252 {
00253 double eps_eqpl;
00254 long ncomp=Mm->ip[ipp].ncompstr;
00255 strastrestate ssst = Mm->ip[ipp].ssst;
00256 vector epsp(ncomp);
00257 matrix epst(ncomp,ncomp);
00258 double tgt = tan(theta);
00259 double cq;
00260
00261 giveirrstrains (ipp, ido, epsp);
00262
00263 vector_tensor (epsp, epst, ssst, strain);
00264
00265 cumulstrain(epst, eps_eqpl);
00266
00267
00268 cq = c + tgt*eps_eqpl;
00269 if ((cq > clim) && (c < clim))
00270 {
00271 cq = clim;
00272 H_value = cq - c;
00273 Hder_value = 0;
00274 }
00275 if ((cq < clim) && (c > clim))
00276 {
00277 cq = clim;
00278 H_value = cq - c;
00279 Hder_value = 0;
00280 }
00281
00282 H_value = cq - c;
00283 Hder_value = tgt * sqrt(1.0+2.0/9.0*alpha1*alpha1);
00284 }
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298 double drprag2::newton_one_element_case2(long ipp, long ido)
00299 {
00300 double K_con, G_con, rho_tr, p_tr, q_value;
00301 double H_value, Hder_value;
00302 double delta_lambda, delta_lambda_next, tau;
00303 double part_q2, part_q3, qDer_value;
00304 int counter_iter;
00305
00306 getImportantValues(ipp, ido, K_con, G_con, rho_tr, p_tr);
00307
00308 delta_lambda = rho_tr/(G_con*sqrt(2));
00309 delta_lambda_next = 0.0;
00310 tau = 1e-4;
00311 counter_iter = 0;
00312
00313 while(1) {
00314
00315 getValueOfHardening(ipp, ido, H_value, Hder_value);
00316 part_q2 = alpha*(p_tr - delta_lambda*K_con*alpha1);
00317 part_q3 = beta*(c+H_value);
00318 q_value = part_q2 - part_q3;
00319 qDer_value = -alpha*alpha1*K_con - beta*Hder_value;
00320
00321 delta_lambda_next = delta_lambda - q_value/qDer_value;
00322
00323 if ((delta_lambda_next-delta_lambda)/(delta_lambda_next+delta_lambda)<=tau || counter_iter >= 20) {
00324 break;
00325 }
00326 else {
00327 delta_lambda = delta_lambda_next;
00328 }
00329 counter_iter= counter_iter+1;
00330 }
00331
00332 return delta_lambda_next;
00333 }
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347 double drprag2::newton_one_element_case3(long ipp, long ido)
00348 {
00349 double K_con, G_con, rho_tr, p_tr, q_value;
00350 double H_value, Hder_value;
00351 double delta_lambda, delta_lambda_next, tau;
00352 double part_q1, part_q2, part_q3, qDer_value;
00353 int counter_iter;
00354
00355 getImportantValues(ipp, ido, K_con, G_con, rho_tr, p_tr);
00356
00357 delta_lambda = 0.0;
00358 delta_lambda_next = 0.0;
00359 tau = 1e-4;
00360 counter_iter = 0;
00361
00362 while(1) {
00363
00364 getValueOfHardening(ipp, ido, H_value, Hder_value);
00365 part_q1 = sqrt(0.5)*max((rho_tr-delta_lambda*G_con*sqrt(2.0)),0);
00366 part_q2 = alpha*(p_tr - delta_lambda*K_con*alpha1);
00367 part_q3 = beta*(c+H_value);
00368 q_value = part_q1 + part_q2 - part_q3;
00369 qDer_value = -G_con - alpha*alpha1*K_con - beta*Hder_value;
00370
00371 delta_lambda_next = delta_lambda - q_value/qDer_value;
00372
00373 if ((delta_lambda_next-delta_lambda)/(delta_lambda_next+delta_lambda)<=tau || counter_iter >= 20) {
00374 break;
00375 }
00376 else {
00377 delta_lambda = delta_lambda_next;
00378 }
00379 counter_iter= counter_iter+1;
00380 }
00381
00382 return delta_lambda_next;
00383 }
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397 void drprag2::matstiff (matrix &d, long ipp,long ido)
00398 {
00399 if (Mp->nlman->stmat==initial_stiff)
00400 {
00401
00402 Mm->elmatstiff (d,ipp);
00403 }
00404 if (Mp->nlman->stmat==tangent_stiff)
00405 {
00406
00407 matrix ad(d.m,d.n);
00408 Mm->elmatstiff (ad,ipp);
00409 tangentstiff (ad,d,ipp,ido);
00410 }
00411 }
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427 void drprag2::tangentstiff (matrix &d,matrix &td,long ipp,long ido)
00428 {
00429 long n=Mm->ip[ipp].ncompstr;
00430
00431 matrix Idev(n,n);
00432 double K_con, G_con, rho_tr, p_tr, q_value0, q_value1, delta_lambda;
00433 double H_value, H1, gamma;
00434 vector iota(n);
00435 vector epsilon_etr(n);
00436 double aux_con1, aux_con2;
00437 vector s_tr(n), n_tr(n), aux_s_tr(n);
00438 matrix aux_td2(td), aux_td3(td), aux_Ntr(td);
00439 vector vec_td31(n), vec_td32(n);
00440
00441 getImportantValues(ipp, ido, K_con, G_con, rho_tr, p_tr, iota, epsilon_etr, Idev);
00442
00443 gamma = 0.0;
00444 q_value0 = getQvalue(gamma, ipp, ido);
00445 gamma = rho_tr/(G_con*sqrt(2));
00446 q_value1 = getQvalue(gamma, ipp, ido);
00447
00448
00449 if (q_value0 < 1e-10) {
00450 copym (d,td);
00451 }
00452
00453 else if (q_value1 >= 1e-10) {
00454
00455
00456 delta_lambda = newton_one_element_case2(ipp, ido);
00457
00458 getValueOfHardening(ipp, ido, H_value, H1);
00459
00460
00461 aux_con1 = (beta*beta*K_con*H1)/(K_con*alpha*alpha1 + beta*beta*H1);
00462 vxv(iota, iota, td);
00463 cmulm(aux_con1, td);
00464 }
00465 else if ((q_value0 > 1e-10) && (q_value1 < 1e-10)) {
00466
00467
00468 delta_lambda = newton_one_element_case3(ipp, ido);
00469
00470 getValueOfHardening(ipp, ido, H_value, H1);
00471
00472 mxv(Idev, epsilon_etr, aux_s_tr);
00473 cmulv((2*G_con), aux_s_tr, s_tr);
00474 cmulv((1.0/rho_tr), s_tr, n_tr);
00475
00476
00477 aux_con1 = delta_lambda * (G_con*sqrt(2)/rho_tr);
00478 vxv(n_tr, n_tr, aux_Ntr);
00479 subm(Idev, aux_Ntr, aux_td2);
00480 cmulm(aux_con1, aux_td2);
00481
00482 aux_con2 = 1.0/(G_con + K_con*alpha*alpha1 + beta*beta*H1);
00483 fillv(0.0, vec_td31);
00484 addmultv (vec_td31, n_tr, (G_con*sqrt(2)));
00485 addmultv (vec_td31, iota, (K_con*alpha1));
00486 fillv(0.0, vec_td32);
00487 addmultv (vec_td32, n_tr, (G_con*sqrt(2)));
00488 addmultv (vec_td32, iota, (K_con*alpha));
00489 cmulv(aux_con2, vec_td32);
00490 vxv(vec_td31, vec_td32, aux_td3);
00491
00492
00493 subm(d, aux_td2, td);
00494 subm(td, aux_td3, td);
00495
00496 }
00497 else {
00498
00499 }
00500
00501
00502 }
00503
00504
00505
00506
00507
00508
00509
00510
00511
00512
00513
00514
00515
00516
00517
00518
00519
00520 void drprag2::nlstresses (long ipp, long im, long ido)
00521 {
00522 long i,n=Mm->ip[ipp].ncompstr;
00523 vector epsp(n),q(1),q_next(1);
00524
00525 matrix Idev(n,n);
00526 matrix De(n,n);
00527 double K_con, G_con, rho_tr, p_tr, q_value0, q_value1, delta_lambda;
00528 double gamma, aux_con1;
00529 vector iota(n), epsilon(n), epsilon_etr(n), epsilon_p_previous(n), epsilon_p_next(n), aux_vec(n);
00530 vector s_tr(n), n_tr(n), aux_s_tr(n), sigma_tr(n), T_lok(n);
00531
00532 q[0] = Mm->ip[ipp].eqother[ido+n+1];
00533
00534 getImportantValues(ipp, ido, K_con, G_con, rho_tr, p_tr, iota, epsilon_etr, Idev);
00535
00536 gamma = 0.0;
00537 q_value0 = getQvalue(gamma, ipp, ido);
00538 gamma = rho_tr/(G_con*sqrt(2));
00539 q_value1 = getQvalue(gamma, ipp, ido);
00540
00541 Mm->elmatstiff(De,ipp);
00542 mxv(De, epsilon_etr, sigma_tr);
00543
00544
00545 Mm->givestrain(0, ipp, epsilon);
00546 giveirrstrains (ipp, ido, epsilon_p_previous);
00547
00548
00549 if (q_value0 < 1e-10) {
00550 copyv (sigma_tr,T_lok);
00551 copyv (epsilon_p_previous, epsilon_p_next);
00552 copyv (q, q_next);
00553 }
00554
00555 else if (q_value1 >= 1e-10) {
00556
00557
00558 delta_lambda = newton_one_element_case2(ipp, ido);
00559
00560
00561 aux_con1 = (p_tr - delta_lambda*K_con*alpha1);
00562 cmulv(aux_con1, iota, T_lok);
00563
00564
00565 aux_con1 = (p_tr - delta_lambda*K_con*alpha1)/(3.0*K_con);
00566 fillv(0.0, aux_vec);
00567 addmultv (aux_vec, iota, aux_con1);
00568 subv(epsilon, aux_vec, epsilon_p_next);
00569
00570
00571 q_next[0] = q[0] + delta_lambda*beta;
00572
00573 }
00574 else if ((q_value0 > 1e-10) && (q_value1 < 1e-10)) {
00575
00576
00577 delta_lambda = newton_one_element_case3(ipp, ido);
00578
00579
00580 fillv(0.0, aux_vec);
00581 addmultv (aux_vec, n_tr, (delta_lambda*G_con*sqrt(2)));
00582 addmultv (aux_vec, iota, (delta_lambda*K_con*alpha1));
00583 subv(sigma_tr, aux_vec, T_lok);
00584
00585
00586 fillv(0.0, aux_vec);
00587 addmultv (aux_vec, n_tr, (delta_lambda*sqrt(0.5)));
00588 addmultv (aux_vec, iota, (delta_lambda*alpha1/3.0));
00589 addv(epsilon_p_previous, aux_vec, epsilon_p_next);
00590
00591
00592 q_next[0] = q[0] + delta_lambda*beta;
00593
00594 }
00595 else {
00596
00597 }
00598
00599 Mm->storestress(0, ipp, T_lok);
00600
00601
00602 for (i=0;i<n;i++){
00603 Mm->ip[ipp].other[ido+i]=epsilon_p_next[i];
00604 }
00605 Mm->ip[ipp].other[ido+n]+=delta_lambda;
00606 Mm->ip[ipp].other[ido+n+1]=q_next[0];
00607
00608 }
00609
00610
00611
00612
00613
00614
00615
00616
00617
00618
00619
00620
00621
00622
00623
00624
00625
00626 void drprag2::updateval (long ipp, long im, long ido)
00627 {
00628 long i,n = Mm->givencompeqother(ipp,im);
00629
00630 for (i=0;i<n;i++){
00631 Mm->ip[ipp].eqother[ido+i]=Mm->ip[ipp].other[ido+i];
00632 }
00633 }
00634
00635
00636
00637
00638
00639
00640
00641
00642
00643
00644
00645
00646
00647
00648 void drprag2::giveirrstrains (long ipp, long ido, vector &epsp)
00649 {
00650 long i;
00651 for (i=0;i<epsp.n;i++)
00652 epsp[i] = Mm->ip[ipp].eqother[ido+i];
00653 }
00654
00655
00656
00657
00658
00659
00660
00661
00662
00663
00664
00665
00666
00667
00668 double drprag2::give_consparam (long ipp, long ido)
00669 {
00670 long ncompstr;
00671 double gamma;
00672
00673 ncompstr=Mm->ip[ipp].ncompstr;
00674 gamma = Mm->ip[ipp].eqother[ido+ncompstr];
00675
00676 return gamma;
00677 }
00678
00679
00680
00681
00682
00683
00684
00685
00686
00687
00688
00689
00690
00691
00692
00693 void drprag2::changeparam (atsel &atm,vector &val)
00694 {
00695 long i;
00696
00697 for (i=0;i<atm.num;i++){
00698 switch (atm.atrib[i]){
00699 case 0:{
00700 phi=val[i];
00701 break;
00702 }
00703 case 1:{
00704 c=val[i];
00705 break;
00706 }
00707 case 2:{
00708 psi=val[i];
00709 break;
00710 }
00711 case 3:{
00712 theta=val[i];
00713 break;
00714 }
00715 case 4:{
00716 clim=val[i];
00717 break;
00718 }
00719 default:{
00720 fprintf (stderr,"\n\n wrong number of atribute in function changeparam (file %s, line %d).\n",__FILE__,__LINE__);
00721 }
00722 }
00723 }
00724
00725 alpha=6.0*sin(phi)/(3.0-sin(phi));
00726 alpha1=6.0*sin(psi)/(3.0-sin(psi));
00727
00728 beta=6.0*cos(phi)/(3.0-sin(phi));
00729 }