00001 #include <math.h>
00002 #include <stdio.h>
00003 #include <stdlib.h>
00004 #include "camclay.h"
00005 #include "global.h"
00006 #include "globmat.h"
00007 #include "vecttens.h"
00008 #include "intpoints.h"
00009 #include "matrix.h"
00010 #include "vector.h"
00011 #include "tensor.h"
00012 #include "elastisomat.h"
00013
00014
00015
00016
00017
00018
00019 camclay::camclay (void)
00020 {
00021 m = 0.0;
00022 lambda = 0.0;
00023 kappa = 0.0;
00024 }
00025
00026
00027
00028
00029
00030
00031 camclay::~camclay (void)
00032 {
00033
00034 }
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045 void camclay::read (XFILE *in)
00046 {
00047 xfscanf (in, "%k%lf%k%lf%k%lf", "m", &m, "lambda", &lambda, "kappa", &kappa);
00048 sra.read (in);
00049 }
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060 void camclay::initval(long ipp, long ido)
00061 {
00062 double v_ini, i1s, v_pc0, v_lambda1;
00063 double v_kappa1, p1, pc_0;
00064 long i, ncompstr = Mm->ip[ipp].ncompstr;
00065 vector epse(ncompstr), sig(ncompstr);
00066 matrix d(ncompstr, ncompstr), sigt(3,3);
00067
00068 if (Mm->ip[ipp].eqother[ido+ncompstr+1] == 0.0)
00069 {
00070
00071 v_kappa1 = Mm->ic[ipp][0];
00072
00073 p1 = Mm->ic[ipp][1];
00074
00075 pc_0 = Mm->ic[ipp][2];
00076
00077
00078 for (i=0; i<ncompstr; i++)
00079 epse[i] = Mm->ic[ipp][3+i];
00080
00081
00082 Mm->elmatstiff (d,ipp);
00083
00084 mxv(d, epse, sig);
00085 vector_tensor (sig,sigt,Mm->ip[ipp].ssst,stress);
00086
00087 i1s = first_invar(sigt)/3.0;
00088
00089
00090 v_pc0 = v_kappa1 - kappa*log(pc_0/p1);
00091
00092 Mm->ip[ipp].eqother[ido+ncompstr+1] = pc_0;
00093
00094 Mm->ip[ipp].eqother[ido+ncompstr+2] = v_lambda1 = v_pc0 + lambda*log(pc_0/p1);
00095
00096 Mm->ip[ipp].eqother[ido+ncompstr+3] = v_ini = v_pc0 + kappa*log(pc_0/i1s);
00097 }
00098 return;
00099 }
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112 double camclay::yieldfunction (matrix &sig, vector &q)
00113 {
00114 double i1s,j2s,f;
00115
00116 matrix dev(3,3);
00117 deviator (sig,dev);
00118
00119 i1s = first_invar (sig)/3.0;
00120 j2s = second_invar (dev);
00121
00122 f = j2s/(m*m) + i1s * (i1s - q[0]);
00123 if (f > 0.0)
00124 return f;
00125
00126 return f;
00127 }
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142 void camclay::deryieldfsigma (matrix &sig, vector &q, matrix &dfds)
00143 {
00144 double volume,i1s;
00145 matrix dev(3,3);
00146
00147 i1s = first_invar (sig)/3.0;
00148
00149 deviator (sig,dfds);
00150
00151
00152
00153
00154
00155
00156
00157
00158 cmulm(1.0/(m*m), dfds);
00159 volume=1.0/3.0*(2.0*i1s - q[0]);
00160 dfds[0][0]+=volume;
00161 dfds[1][1]+=volume;
00162 dfds[2][2]+=volume;
00163 }
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175 void camclay::dderyieldfsigma (matrix &ddfds)
00176 {
00177 fillm(0.0, ddfds);
00178 ddfds[0][0] = ddfds[1][1] = ddfds[2][2] = 2.0/(3.0*m*m) + 2.0/9.0;
00179 ddfds[0][1] = ddfds[0][2] = ddfds[1][0] = ddfds[1][2] = -1.0/(3.0*m*m) + 2.0/9.0;
00180 ddfds[2][0] = ddfds[2][1] = ddfds[0][1];
00181 ddfds[3][3] = 1.0/(m*m);
00182 ddfds[4][4] = 1.0/(m*m);
00183 ddfds[5][5] = 1.0/(m*m);
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209 }
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221 void camclay::derpotsigma (matrix &sig, vector &q, matrix &dgds)
00222 {
00223 deryieldfsigma (sig, q, dgds);
00224 }
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238 void camclay::deryieldfq(matrix &sig, vector &dfq)
00239 {
00240 dfq[0] = -first_invar (sig)/3.0;
00241
00242
00243
00244
00245 return;
00246 }
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259 void camclay::deryieldfdqdq(matrix &dfq)
00260 {
00261 dfq[0][0] = 0.0;
00262 return;
00263 }
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277 void camclay::deryieldfdsdq(matrix &dfdsdqt)
00278 {
00279 dfdsdqt[0][0] = dfdsdqt[1][0] = dfdsdqt[2][0] = -1.0/3.0;
00280 return;
00281 }
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297 void camclay::dqdsigma(matrix &sigt, vector &qt, matrix &dqds)
00298 {
00299 double i1s;
00300
00301 i1s = first_invar (sigt)/3.0;
00302
00303 fillm(0.0, dqds);
00304 dqds[0][0] = dqds[1][1] = dqds[2][2] = qt[0]*kappa/i1s*1/3.0;
00305 return;
00306 }
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325 void camclay::dhardfdq(long ipp, long ido, double dgamma, vector &qt, matrix &dqdq)
00326 {
00327 double v_ini;
00328 long ncompstr = Mm->ip[ipp].ncompstr;
00329
00330 fillm(0.0, dqdq);
00331 v_ini = Mm->ip[ipp].eqother[ido+ncompstr+3];
00332 dqdq[0][0] = -dgamma * v_ini * qt[0]/(kappa-lambda);
00333 return;
00334 }
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351 void camclay::der_q_gamma(long ipp, long ido, matrix &sig, vector &qtr, vector &dqdg)
00352 {
00353 double v_ini, i1s;
00354 long ncompstr = Mm->ip[ipp].ncompstr;
00355
00356
00357 i1s = first_invar (sig)/3.0;
00358 v_ini = Mm->ip[ipp].eqother[ido+ncompstr+3];
00359
00360 dqdg[0] = qtr[0] * (v_ini)/(kappa-lambda) * (2.0*i1s-qtr[0]);
00361
00362
00363
00364 return;
00365 }
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379 double camclay::plasmodscalar(long ipp, long ido, matrix &sig, vector &qtr)
00380 {
00381 double ret;
00382 vector dfq(1);
00383 vector dqg(1);
00384
00385 deryieldfq(sig, dfq);
00386 der_q_gamma(ipp, ido, sig, qtr, dqg);
00387 scprd(dfq, dqg, ret);
00388
00389 return -ret;
00390 }
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406 void camclay::updateq(long ipp, long ido, vector &eps, vector &epsp, matrix &sig, vector &q)
00407 {
00408 matrix epst(3,3);
00409 vector depsp(epsp.n);
00410 double v_kappa1, v_lambda1, p1, depsvp, v_ini, pc_new;
00411 strastrestate ssst = Mm->ip[ipp].ssst;
00412 long ncompstr = Mm->ip[ipp].ncompstr;
00413 long i;
00414
00415
00416
00417 v_kappa1 = Mm->ic[ipp][0];
00418
00419 p1 = Mm->ic[ipp][1];
00420
00421 v_lambda1 = Mm->ip[ipp].eqother[ido+ncompstr+2];
00422
00423 v_ini = Mm->ip[ipp].eqother[ido+ncompstr+3];
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433 for (i=0; i<ncompstr; i++)
00434 depsp[i] = epsp[i] + Mm->ic[ipp][3+i];
00435 vector_tensor (depsp, epst, ssst, strain);
00436 depsvp = first_invar(epst);
00437
00438 pc_new = p1 * exp((v_kappa1 - v_lambda1 + v_ini*depsvp)/(kappa-lambda));
00439 if (pc_new < q[0])
00440 q[0] = pc_new;
00441
00442 return;
00443 }
00444
00445
00446
00447
00448
00449
00450
00451
00452
00453
00454
00455 double camclay::give_actual_ym(long ipp, long im, long ido)
00456 {
00457 long i;
00458 long idem = Mm->ip[ipp].gemid();
00459 long ncompo = Mm->givencompeqother(ipp, im);
00460 double e;
00461
00462
00463 double e_ini = Mm->give_initial_ym(ipp, idem, ido+ncompo);
00464
00465 double nu = give_actual_nu(ipp, im, ido);
00466
00467
00468 long ncomp=Mm->ip[ipp].ncompstr;
00469 strastrestate ssst = Mm->ip[ipp].ssst;
00470 vector aux(ncomp);
00471 matrix auxt(3,3);
00472
00473 Mm->givestrain(0, ipp, aux);
00474 vector_tensor(aux, auxt, ssst, strain);
00475 double epsv = first_invar(auxt);
00476
00477
00478 for (i=0; i<ncomp; i++)
00479 aux[i] = Mm->ic[ipp][3+i];
00480 vector_tensor(aux, auxt, ssst, strain);
00481 double epsv_ini = first_invar(auxt);
00482 if (epsv >= epsv_ini)
00483 return e_ini;
00484
00485
00486 double p_ini = e_ini/(3.0*(1.0-2.0*nu))*epsv_ini;
00487
00488 Mm->givestress(0, ipp, aux);
00489 vector_tensor (aux, auxt, ssst, stress);
00490 double p = first_invar(auxt)/3.0;
00491
00492 if (p >= p_ini)
00493 return e_ini;
00494
00495
00496 double v_ini = Mm->ip[ipp].eqother[ido+ncomp+3];
00497
00498 double v = v_ini*(1.0+epsv);
00499
00500 double pc_0 = Mm->ic[ipp][2];
00501 if (p >= pc_0)
00502 return e_ini;
00503
00504 e = 3.0*(1.0-nu)*v*(-p)/kappa;
00505
00506 return e;
00507 }
00508
00509
00510
00511
00512
00513
00514
00515
00516
00517
00518
00519 double camclay::give_actual_nu(long ipp, long im, long ido)
00520 {
00521 long idem = Mm->ip[ipp].gemid();
00522 long ncompo = Mm->givencompeqother(ipp, im);
00523 double nu = Mm->give_initial_nu(ipp, idem, ido+ncompo);
00524
00525 return nu;
00526 }
00527
00528
00529
00530
00531
00532
00533
00534
00535
00536
00537
00538
00539 void camclay::matstiff (matrix &d, long ipp, long im, long ido)
00540 {
00541 double zero=1.0e-20;
00542
00543 if (Mp->nlman->stmat==1)
00544 {
00545
00546 Mm->elmatstiff (d,ipp);
00547 }
00548 if (Mp->nlman->stmat > 1)
00549 {
00550
00551 long n = Mm->ip[ipp].ncompstr, i, j;
00552 double gamma = Mm->ip[ipp].other[ido+n];
00553
00554 matrix de(n,n), ddfdst(6,6), ddfds(n,n);
00555 vector sig(n), dfds(n), dgds(n);
00556 vector q(1);
00557
00558 matrix sigt(3,3), auxt(3,3);
00559 matrix auxm(n,n), tmp(n,n);
00560 vector auxv(n), tmpv(n);
00561 double aux1;
00562
00563
00564
00565
00566 Mm->elmatstiff (de,ipp);
00567
00568 if (Mp->stressstate == 0)
00569 nlstresses(ipp, im, ido);
00570
00571 Mm->givestress(0, ipp, auxv);
00572
00573 vector_tensor (auxv, sigt, Mm->ip[ipp].ssst, stress);
00574 q[0] = Mm->ip[ipp].eqother[ido+n+1];
00575 double gammap = Mm->ip[ipp].eqother[ido+n+5];
00576 double dgamma = gamma - gammap;
00577
00578 aux1 = yieldfunction(sigt, q);
00579 if ((aux1 < -1.0) || (dgamma == 0.0))
00580 {
00581 for (i=0; i<d.m; i++)
00582 {
00583 for(j=0; j<d.n; j++)
00584 d[i][j] = de[i][j];
00585 }
00586 return;
00587 }
00588 else
00589 {
00590 dderyieldfsigma (ddfdst);
00591
00592 tensor4_matrix (ddfds,ddfdst,Mm->ip[ipp].ssst);
00593 cmulm(gamma, de, tmp);
00594 mxm(tmp, ddfds, auxm);
00595 for (i = 0; i < n; i++)
00596 auxm[i][i] += 1.0;
00597 invm(auxm, tmp, zero);
00598 mxm(tmp, de, auxm);
00599
00600 }
00601
00602
00603 deryieldfsigma(sigt, q, auxt);
00604 tensor_vector(dfds, auxt, Mm->ip[ipp].ssst, stress);
00605 derpotsigma(sigt, q, auxt);
00606 tensor_vector(dgds, auxt, Mm->ip[ipp].ssst, stress);
00607
00608
00609
00610
00611
00612
00613
00614 mxv(auxm, dgds, auxv);
00615 vxm(dfds, auxm, tmpv);
00616 vxv(auxv, tmpv, tmp);
00617
00618
00619 mxv(auxm, dgds, auxv);
00620 scprd(dfds, auxv, aux1);
00621 aux1 += plasmodscalar(ipp, ido, sigt, q);
00622
00623 cmulm(1.0/aux1, tmp, tmp);
00624
00625
00626
00627
00628
00629 for (i=0; i<d.m; i++)
00630 {
00631 for(j=0; j<d.n; j++)
00632 d[i][j] = de[i][j] - tmp[i][j];
00633 }
00634 }
00635 }
00636
00637
00638
00639
00640
00641
00642
00643
00644
00645
00646
00647
00648
00649 void camclay::nlstresses (long ipp, long im, long ido)
00650
00651 {
00652 long i,ni,ncomp=Mm->ip[ipp].ncompstr;
00653 double gamma,err;
00654 vector epsn(ncomp),epsp(ncomp),q(1);
00655 vector epsa(ncomp),sig(ncomp),dfds(ncomp),dgds(ncomp);
00656 matrix d(ncomp,ncomp),sigt(3,3),dfdst(3,3),dgdst(3,3);
00657 matrix epst(3,3), epspt(3,3);
00658 matrix dev(3,3);
00659 double i1s, j2s, epsv, epsvp;
00660
00661
00662
00663
00664 for (i=0; i<ncomp; i++)
00665 {
00666 epsn[i] = Mm->ip[ipp].strain[i];
00667 epsp[i] = Mm->ip[ipp].eqother[ido+i];
00668
00669
00670
00671
00672 epsp[i] -= Mm->ic[ipp][3+i];
00673 }
00674
00675 gamma=Mm->ip[ipp].eqother[ido+ncomp];
00676 q[0] = Mm->ip[ipp].eqother[ido+ncomp+1];
00677
00678
00679 switch(sra.give_tsra ())
00680 {
00681 case cp:
00682 ni=sra.give_ni ();
00683 err=sra.give_err ();
00684 Mm->cutting_plane (ipp, im, ido, gamma, epsn, epsp, q, ni, err);
00685 break;
00686 case gsra:
00687 ni=sra.give_ni ();
00688 err=sra.give_err ();
00689 Mm->newton_stress_return (ipp, im, ido, gamma, epsn, epsp, q, ni, err);
00690 break;
00691 default:
00692 print_err("wrong type of stress return algorithm is required", __FILE__, __LINE__, __func__);
00693 abort ();
00694 }
00695
00696
00697
00698 for (i=0; i<ncomp; i++)
00699 {
00700 epsp[i] += Mm->ic[ipp][3+i];
00701 Mm->ip[ipp].other[ido+i]=epsp[i];
00702 }
00703
00704 Mm->ip[ipp].other[ido+ncomp]=gamma;
00705 Mm->ip[ipp].other[ido+ncomp+1]=q[0];
00706 Mm->givestress(0, ipp, sig);
00707 vector_tensor (sig,sigt,Mm->ip[ipp].ssst,strain);
00708 vector_tensor (epsn,epst,Mm->ip[ipp].ssst,strain);
00709 vector_tensor (epsp,epspt,Mm->ip[ipp].ssst,strain);
00710 deviator (sigt,dev);
00711 i1s = first_invar (sigt)/3.0;
00712 j2s = sqrt(second_invar (dev));
00713 epsv = first_invar (epst);
00714 epsvp = first_invar (epspt);
00715 Mm->ip[ipp].other[ido+ncomp+4] = i1s;
00716 Mm->ip[ipp].other[ido+ncomp+5] = give_actual_ym(ipp, im, ido);
00717 Mm->ip[ipp].other[ido+ncomp+6] = epsv;
00718 Mm->ip[ipp].other[ido+ncomp+7] = epsvp;
00719
00720 }
00721
00722
00723
00724
00725
00726
00727
00728
00729
00730
00731
00732 void camclay::updateval (long ipp,long ido)
00733 {
00734 long i,n = Mm->ip[ipp].ncompstr;
00735
00736 Mm->ip[ipp].eqother[ido+n+5] = Mm->ip[ipp].eqother[ido+n];
00737
00738 for (i=0;i<n;i++){
00739 Mm->ip[ipp].eqother[ido+i]=Mm->ip[ipp].other[ido+i];
00740 }
00741 Mm->ip[ipp].eqother[ido+n] = Mm->ip[ipp].other[ido+n];
00742 Mm->ip[ipp].eqother[ido+n+1] = Mm->ip[ipp].other[ido+n+1];
00743
00744
00745 Mm->ip[ipp].eqother[ido+n+4] = Mm->ip[ipp].other[ido+n+4];
00746 Mm->ip[ipp].eqother[ido+n+5] = Mm->ip[ipp].other[ido+n+5];
00747 Mm->ip[ipp].eqother[ido+n+6] = Mm->ip[ipp].other[ido+n+6];
00748 Mm->ip[ipp].eqother[ido+n+7] = Mm->ip[ipp].other[ido+n+7];
00749 }
00750
00751
00752
00753
00754
00755
00756
00757
00758
00759
00760
00761
00762 void camclay::giveirrstrains (long ipp, long ido, vector &epsp)
00763 {
00764 long i;
00765 for (i=0;i<epsp.n;i++)
00766 epsp[i] = Mm->ip[ipp].eqother[ido+i];
00767 }
00768
00769
00770
00771
00772
00773
00774
00775
00776
00777
00778
00779
00780
00781 double camclay::give_consparam (long ipp,long ido)
00782 {
00783 long ncompstr;
00784 double gamma;
00785
00786 ncompstr=Mm->ip[ipp].ncompstr;
00787 gamma = Mm->ip[ipp].eqother[ido+ncompstr];
00788
00789 return gamma;
00790 }
00791
00792
00793
00794
00795
00796
00797
00798
00799
00800
00801
00802
00803
00804
00805 double camclay::give_preconspress(long ipp, long ido)
00806 {
00807 long ncompstr;
00808 double pc;
00809
00810 ncompstr=Mm->ip[ipp].ncompstr;
00811 pc = Mm->ip[ipp].eqother[ido+ncompstr+1];
00812
00813 return pc;
00814 }
00815
00816
00817
00818
00819
00820
00821
00822
00823
00824
00825
00826
00827
00828
00829 double camclay::give_virgporosity(long ipp, long ido)
00830 {
00831 long ncompstr;
00832 double e_lambda1;
00833
00834 ncompstr=Mm->ip[ipp].ncompstr;
00835 e_lambda1 = Mm->ip[ipp].eqother[ido+ncompstr+2]-1.0;
00836
00837 return e_lambda1;
00838 }
00839
00840
00841
00842
00843
00844
00845
00846
00847
00848
00849
00850
00851
00852
00853 double camclay::give_iniporosity(long ipp, long ido)
00854 {
00855 long ncompstr;
00856 double e_ini;
00857
00858 ncompstr=Mm->ip[ipp].ncompstr;
00859
00860 e_ini = Mm->ip[ipp].eqother[ido+ncompstr+3]-1.0;
00861
00862 return e_ini;
00863 }
00864
00865
00866
00867
00868
00869
00870
00871
00872
00873
00874
00875
00876
00877 void camclay::changeparam (atsel &atm,vector &val)
00878 {
00879 long i;
00880
00881 for (i=0;i<atm.num;i++)
00882 {
00883 switch (atm.atrib[i])
00884 {
00885 case 0:
00886 m=val[i];
00887 break;
00888 case 1:
00889 lambda=val[i];
00890 break;
00891 case 2:
00892 kappa=val[i];
00893 break;
00894 default:
00895 print_err("wrong number of atribute", __FILE__, __LINE__, __func__);
00896 }
00897 }
00898 }