00001 #include "hissplas.h"
00002 #include "global.h"
00003 #include "iotools.h"
00004 #include "matrix.h"
00005 #include "vector.h"
00006 #include "tensor.h"
00007 #include "vecttens.h"
00008 #include "intpoints.h"
00009 #include <stdlib.h>
00010 #include <math.h>
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020 hissplas::hissplas ()
00021 {
00022 beta = gamma_max = n = r = pa = alpha1 = eta1 = ksi_lim = k = 0.0;
00023 }
00024
00025
00026
00027
00028
00029
00030
00031
00032 hissplas::~hissplas ()
00033 {
00034 }
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048 void hissplas::read (XFILE *in)
00049 {
00050 xfscanf(in, "%k%lf %k%lf %k%lf% k%lf% k%lf% k%lf% k%lf %k%lf",
00051 "beta", &beta, "n", &n, "r", &r, "pa", &pa,
00052 "alpha1", &alpha1, "eta1", &eta1, "ksi_lim", &ksi_lim, "gamma_max", &gamma_max);
00053 sra.read (in);
00054 }
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075 double hissplas::yieldfunction (double i1, double j2, double cos3t, vector &q)
00076 {
00077 double f, fb, fc;
00078 double aux, alpha, gamma;
00079
00080
00081 alpha = q[0];
00082 gamma = q[1];
00083
00084 aux = (i1+r)/pa;
00085 fb = -alpha*pow(aux, n) + gamma*aux*aux;
00086
00087
00088
00089 fc = 1.0;
00090
00091 f = j2/(pa*pa)-fb*fc;
00092
00093 return f;
00094 }
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113 double hissplas::yieldfunction (matrix sigt, vector &q)
00114 {
00115 double i1, j2, j3, cos3t, f;
00116 matrix dev(3,3);
00117
00118
00119
00120
00121 i1 = first_invar(sigt);
00122 deviator(sigt, dev);
00123 j2 = second_invar(dev);
00124 j3 = third_invar(dev);
00125 cos3t = 3.0*sqrt(3.0)/2.0 * j3 /pow(j2, 3.0/2.0);
00126
00127 f = yieldfunction(i1, j2, cos3t, q);
00128
00129 return f;
00130 }
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145 double hissplas::g1function(double depsvp, double depsdp, double i1, double j2, double alpha, double gamma)
00146 {
00147 double ret, aux;
00148
00149
00150 ret = 2.0/(pa*pa) * sqrt(j2) * depsvp;
00151 aux = (i1+r)/pa;
00152 ret -= -(1.0/pa)*(-n*alpha*pow(aux, n-1.0) + 2.0*gamma*aux);
00153
00154 return ret;
00155 }
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174 double hissplas::dfdepsvp (long ipp, double i1, double alpha, double gamma, double ksi, matrix &depsp)
00175 {
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186 double aux, auxA, auxB, ret;
00187
00188 double e = Mm->give_actual_ym(ipp);
00189
00190 double nu = Mm->give_actual_nu(ipp);
00191
00192 double k = e/(3.0*(1-2*nu));
00193
00194 auxA = (1.0/pa)*(-n*alpha*pow((i1+r)/pa, n-1.0)+2.0*gamma*(i1+r)/pa);
00195 auxA *= -9.0*k;
00196
00197 auxB = -pow((i1+r)/pa, n);
00198
00199 aux = (1+pow(ksi,eta1));
00200 auxB *= -alpha1/(aux*aux)*(pow(ksi,eta1-1.0)*(eta1*ksi_lim + ksi*(1-eta1))+1.0);
00201 auxB *= first_invar(depsp)/tensornorm(depsp);
00202
00203 aux = -1.0;
00204 ret = aux*(auxA + auxB);
00205 return ret;
00206 }
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225 double hissplas::dfdepspd (long ipp, double i1, double j2, matrix &dev, double ksi, matrix &depsp)
00226 {
00227
00228
00229
00230
00231
00232
00233
00234 double aux, aux1, aux2, ret;
00235
00236 double e = Mm->give_actual_ym(ipp);
00237
00238 double nu = Mm->give_actual_nu(ipp);
00239
00240 double g = e/(2.0*(1.0+nu));
00241
00242
00243 aux1 = 1.0;
00244
00245 aux1 *= 2.0*sqrt(j2)/(pa*pa);
00246
00247 aux1 *= -g/2.0*tensor2_inner_prod(dev,dev)/j2;
00248
00249
00250 aux2 = -1.0;
00251
00252 aux2 *= -pow((i1+r)/pa, n);
00253
00254 aux = (1+pow(ksi,eta1));
00255 aux2 *= -alpha1/(aux*aux)*(pow(ksi,eta1-1.0)*(eta1*ksi_lim + ksi*(1-eta1))+1.0);
00256
00257
00258 aux2 *= tensor2_inner_prod(depsp, dev);
00259 aux2 /= 2.0*sqrt(j2)*tensornorm(depsp);
00260
00261 ret = aux1+aux2;
00262 return ret;
00263 }
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285 double hissplas::dg1depsvp (long ipp, double i1, double j2, double alpha, double gamma, matrix &dev, double depsdp, double ksi, matrix &depsp)
00286 {
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296 double ddfdsqrtJ2ddeps_vp;
00297 double aux, auxA, auxB, ret;
00298
00299 double e = Mm->give_actual_ym(ipp);
00300
00301 double nu = Mm->give_actual_nu(ipp);
00302
00303 double k = e/(3.0*(1-2*nu));
00304
00305
00306
00307 ddfdsqrtJ2ddeps_vp = 2.0/(pa*pa)*sqrt(j2);
00308
00309
00310
00311
00312
00313
00314 auxA = pow(i1+r, n-1.0);
00315
00316 aux = (1+pow(ksi,eta1));
00317 auxA *= -alpha1/(aux*aux)*(pow(ksi,eta1-1.0)*(eta1*ksi_lim + ksi*(1-eta1))+1.0);
00318
00319
00320 auxA *= first_invar(depsp)/tensornorm(depsp);
00321
00322
00323 auxA += -alpha*9.0*k*(n-1.0)*pow(i1+r, n-2.0);
00324
00325 auxA *= -n;
00326
00327
00328
00329 auxB = -18.0*gamma*k*pow(pa, n-2.0);
00330
00331 ret = ddfdsqrtJ2ddeps_vp-depsdp/pow(pa, n)*(auxA+auxB);
00332 return ret;
00333 }
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356 double hissplas::dg1depsdp (long ipp, double i1, double j2, double alpha, double gamma, matrix &dev, double depsvp, double depsdp, double ksi, matrix &depsp)
00357 {
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373 double aux, auxA, auxB, ret;
00374
00375 double e = Mm->give_actual_ym(ipp);
00376
00377 double nu = Mm->give_actual_nu(ipp);
00378
00379 double g = e/(2.0*(1.0+nu));
00380
00381 auxA = depsvp*2.0/(pa*pa);
00382 auxA *= -g/2.0*tensor2_inner_prod(dev, dev)/j2;
00383
00384 auxB = n/pow(pa, n);
00385 auxB *= pow(i1+r, n-1.0);
00386 aux = (1+pow(ksi,eta1));
00387 auxB *= -alpha1/(aux*aux)*(pow(ksi,eta1-1.0)*(eta1*ksi_lim + ksi*(1-eta1))+1.0);
00388 auxB *= tensor2_inner_prod(depsp, dev)/2.0/sqrt(j2);
00389 auxB /= tensornorm(depsp);
00390 auxB *= depsdp;
00391 auxB += -1.0/pa*(-n*alpha*pow((i1+r)/pa, n-1) + 2.0*gamma*(i1+r)/pa);
00392
00393 ret = auxA-auxB;
00394 return ret;
00395 }
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409 void hissplas::matstiff (matrix &d, long ipp, long ido)
00410 {
00411 if (Mp->nlman->stmat==initial_stiff)
00412 {
00413
00414 Mm->elmatstiff (d,ipp);
00415 }
00416 if (Mp->nlman->stmat==tangent_stiff)
00417 {
00418
00419 matrix ad(d.m,d.n);
00420
00421
00422 Mm->elmatstiff (d,ipp);
00423 }
00424 }
00425
00426
00427
00428
00429
00430
00431
00432
00433
00434
00435
00436
00437
00438
00439 void hissplas::nonloc_nlstresses (long ipp,long ido)
00440 {
00441 }
00442
00443
00444
00445
00446
00447
00448
00449
00450
00451
00452
00453
00454 void hissplas::nlstresses (long ipp,long ido)
00455 {
00456 long ncomp = Mm->ip[ipp].ncompstr;
00457 long eqstate;
00458 long i;
00459 vector epsn(ncomp);
00460 vector epsp(ncomp);
00461 vector depsp(ncomp);
00462 vector ddepsp(ncomp);
00463 vector epsa(ncomp);
00464 vector sig(ncomp);
00465 vector dsig(ncomp);
00466 vector qo(4);
00467 vector q(4);
00468 matrix epspt(3,3);
00469 matrix depspt(3,3);
00470 matrix ddepspt(3,3);
00471 matrix sigt(3,3);
00472 matrix dev(3,3);
00473 matrix d(ncomp,ncomp);
00474 matrix h(2,2);
00475 vector e(2);
00476 vector r(2);
00477 double f;
00478 double g1;
00479 double lambda;
00480 double norsig;
00481 double nordsig;
00482 double norepsp;
00483 double norddepsp;
00484 double pnordsig;
00485 double pnorddepsp;
00486 double i1;
00487 double j2, j3, cos3t;
00488 double depsvp;
00489 double depsdp;
00490 double &alpha = q[0];
00491 double &gamma = q[1];
00492 double &ksi = q[2];
00493 double &wpf = q[3];
00494
00495
00496
00497
00498 for (i=0; i<ncomp; i++)
00499 {
00500 epsn[i]=Mm->ip[ipp].strain[i];
00501 epsp[i]=Mm->ip[ipp].eqother[ido+i];
00502 }
00503 lambda = Mm->ip[ipp].eqother[ido+ncomp];
00504 qo[0] = Mm->ip[ipp].eqother[ido+ncomp+1];
00505 qo[1] = Mm->ip[ipp].eqother[ido+ncomp+2];
00506 qo[2] = Mm->ip[ipp].eqother[ido+ncomp+3];
00507 qo[3] = Mm->ip[ipp].eqother[ido+ncomp+4];
00508 copyv(qo, q);
00509
00510
00511
00512
00513 Mm->elmatstiff(d, ipp);
00514 subv(epsn, epsp, epsa);
00515 mxv(d,epsa,sig);
00516 vector_tensor (sig,sigt,Mm->ip[ipp].ssst,stress);
00517
00518
00519
00520
00521 i1 = first_invar(sigt);
00522 deviator(sigt, dev);
00523 j2 = second_invar(dev);
00524 j3 = third_invar(dev);
00525 cos3t = 3.0*sqrt(3.0)/2.0 * j3 /pow(j2, 3.0/2.0);
00526
00527
00528
00529
00530 eqstate = 1;
00531 f = yieldfunction(i1, j2, cos3t, q);
00532 if (f > sra.err)
00533 {
00534
00535
00536
00537 vector_tensor (epsp,epspt,Mm->ip[ipp].ssst,strain);
00538 depsvp = 0.0;
00539 depsdp = 0.0;
00540
00541 eqstate = 0;
00542
00543 if (alpha > 0.0)
00544 {
00545
00546
00547
00548 g1 = g1function(depsvp, depsdp, i1, j2, alpha, gamma);
00549 for (i=0; i<sra.ni; i++)
00550 {
00551
00552 h[0][0] = dg1depsvp(ipp, i1, j2, alpha, gamma, dev, depsdp, ksi, depspt);
00553 h[0][1] = dg1depsdp(ipp, i1, j2, alpha, gamma, dev, depsvp, depsdp, ksi, depspt);
00554 h[1][0] = dfdepsvp (ipp, i1, alpha, gamma, ksi, depspt);
00555 h[1][1] = dfdepspd (ipp, i1, j2, dev, ksi, depspt);
00556
00557
00558 e[0] = depsvp;
00559 e[1] = depsdp;
00560
00561 r[0] = 0.0 - g1;
00562 r[1] = 0.0 - f;
00563
00564
00565 gause(h,h,e,r,1.0e-10);
00566
00567
00568 norsig = normv(sig);
00569 norepsp = normv(epsp);
00570
00571
00572 copym(dev, ddepspt);
00573 cmulm(r[1]/(2.0*sqrt(j2)), ddepspt);
00574 ddepspt[0][0] = r[0];
00575 ddepspt[1][1] = r[0];
00576 ddepspt[2][2] = r[0];
00577 tensor_vector(ddepsp, ddepspt, Mm->ip[ipp].ssst, strain);
00578
00579
00580 depsvp += r[0];
00581 depsdp += r[1];
00582 depspt[0][0] = depspt[1][1] = depspt[2][2] = depsvp;
00583 cmulm(depsdp/(2.0*sqrt(j2)), dev);
00584 addm(depspt, dev, depspt);
00585 tensor_vector(depsp, depspt, Mm->ip[ipp].ssst, strain);
00586
00587
00588 addv(epsp, depsp, epsp);
00589
00590 mxv(d, ddepsp, dsig);
00591
00592 subv(epsn, epsp, epsa);
00593 mxv(d,epsa,sig);
00594 vector_tensor (sig,sigt,Mm->ip[ipp].ssst,stress);
00595
00596 i1 = first_invar(sigt);
00597 deviator(sigt, dev);
00598 j2 = second_invar(dev);
00599 j3 = third_invar(dev);
00600 cos3t = 3.0*sqrt(3.0)/2.0 * j3 /pow(j2, 3.0/2.0);
00601
00602 copyv(qo, q);
00603 updateq(ipp, depspt, sigt, depspt, q);
00604
00605
00606 f = yieldfunction(i1, j2, cos3t, q);
00607
00608 g1 = g1function(depsvp, depsdp, i1, j2, alpha, gamma);
00609
00610 norddepsp = normv(ddepsp);
00611
00612 nordsig = normv(dsig);
00613
00614
00615 pnordsig = nordsig/norsig;
00616 pnorddepsp = norddepsp/norepsp;
00617 if ((pnordsig < sra.err) && (pnorddepsp < sra.err) && (f < sra.err) && (fabs(g1) < sra.err))
00618 {
00619
00620
00621
00622
00623 lambda = depsdp/(2.0*sqrt(j2)) * pa * pa;
00624 eqstate = 1;
00625 break;
00626 }
00627 }
00628 }
00629 else
00630 {
00631
00632
00633
00634 }
00635 }
00636 if (!eqstate)
00637 {
00638 print_err("Stress return algorithm does not converge on element %ld (ipp %ld)",
00639 __FILE__, __LINE__, __func__, Mm->elip[ipp]+1, ipp+1);
00640 }
00641
00642
00643
00644
00645
00646 Mm->storestress(0,ipp,sig);
00647
00648 for (i=0; i<ncomp; i++)
00649 Mm->ip[ipp].other[ido+i] = epsp[i];
00650
00651
00652 Mm->ip[ipp].other[ido+ncomp+0] = lambda;
00653
00654 Mm->ip[ipp].other[ido+ncomp+1] = alpha;
00655 Mm->ip[ipp].other[ido+ncomp+2] = gamma;
00656 Mm->ip[ipp].other[ido+ncomp+3] = ksi;
00657 Mm->ip[ipp].other[ido+ncomp+4] = wpf;
00658 }
00659
00660
00661
00662
00663
00664
00665
00666
00667
00668
00669
00670
00671
00672
00673
00674
00675
00676
00677 void hissplas::updateq (long ipp, matrix &depspt, matrix &sigt, matrix &depst, vector &q)
00678 {
00679 double alpha, gamma, ksi, wpf;
00680
00681
00682 ksi = q[2];
00683
00684 ksi += tensornorm(depspt);
00685 if (ksi_lim >= ksi)
00686 alpha = alpha1*(ksi_lim - ksi)/(1+pow(ksi, eta1));
00687 else
00688 alpha = 0.0;
00689
00690
00691 if (q[0] == 0.0)
00692 {
00693 wpf = q[3];
00694
00695 wpf += tensor2_inner_prod(sigt, depst);
00696 gamma = sqrt(gamma_max*gamma_max-k*pow(wpf,3.0/2.0));
00697 }
00698 else
00699 gamma = gamma_max;
00700
00701
00702 q[0] = alpha;
00703 q[1] = gamma;
00704 q[2] = ksi;
00705 q[3] = wpf;
00706 }
00707
00708
00709
00710
00711
00712
00713
00714
00715
00716
00717
00718
00719
00720
00721 void hissplas::updateval (long ipp, long ido)
00722 {
00723 long i,n = Mm->ip[ipp].ncompstr;
00724
00725 for (i=0;i<n;i++){
00726 Mm->ip[ipp].eqother[ido+i]=Mm->ip[ipp].other[ido+i];
00727 }
00728 Mm->ip[ipp].eqother[ido+n] = Mm->ip[ipp].other[ido+n+0];
00729 Mm->ip[ipp].eqother[ido+n+1] = Mm->ip[ipp].other[ido+n+1];
00730 Mm->ip[ipp].eqother[ido+n+2] = Mm->ip[ipp].other[ido+n+2];
00731 Mm->ip[ipp].eqother[ido+n+3] = Mm->ip[ipp].other[ido+n+3];
00732 Mm->ip[ipp].eqother[ido+n+4] = Mm->ip[ipp].other[ido+n+4];
00733 }
00734
00735
00736
00737
00738
00739
00740
00741
00742
00743
00744
00745
00746
00747
00748
00749
00750 void hissplas::giveirrstrains (long ipp, long ido, vector &epsp)
00751 {
00752 long i;
00753 long ncomp = Mm->ip[ipp].ncompstr;
00754
00755 if (epsp.n != ncomp)
00756 {
00757 print_err("wrong number of components of plastic strain vector", __FILE__, __LINE__, __func__);
00758 abort();
00759 }
00760 for(i=0; i<epsp.n; i++)
00761 epsp[i] = Mm->ip[ipp].eqother[ido+i];
00762 }
00763
00764
00765
00766
00767
00768
00769
00770
00771
00772
00773
00774
00775
00776 double hissplas::give_consparam (long ipp,long ido)
00777 {
00778 long ncomp = Mm->ip[ipp].ncompstr;
00779 return Mm->ip[ipp].eqother[ido+ncomp];
00780 }
00781
00782
00783
00784
00785
00786
00787
00788
00789
00790
00791
00792
00793
00794
00795 void changeparam (atsel &atm,vector &val)
00796 {
00797 }
00798
00799