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