1 /*
2 This routine enables to implement kinematic hardening, bounding
3 surface plasticity models into Tochnog.
4 Source code is fully compatible with single element program Triax.
5
6 Uses forward euler scheme with error controll and automatic substepping.
7 Initial intersection with the yield surface found using nephton-raphson
8 iterations for scalar multiplier alpha.
9
10 options_element_dof were added to store dofs in integration points, rather
11 than only in nodes.
12
13 These changes into standard distribution of Tochnog are necessary:
14 -include this file in makefile
15 -add functions matrix_inverse_general, make_dev, matrix4_ab, matrix_a_contr_b
16 matrix_ab4 and cinv into math.cc and declare these functions in tochnog.h
17 -include matrix.h in tochnog.h
18 -declare plasti_incr in tochnog.h
19 -cut-off backward--euler iterations in stress.cc
20 -add new models into database.cc
21 -add new elastic models into stress.cc
22 It is possible to use constitutive models source code of single element
23 program Triax after minor modification
24
25 For options_element_dof these files had to be modified:
26 -tochnog.h, database.cc, create.cc, top.cc, materi.cc, initia.cc, elem.cc
27 Modifications are mostly labeled with comment '// added for options_element_dof'
28 search for 'element_dof' as well. As comparison with the professional version
29 discovered some bug in node_dof from free version, whereas the
30 options_element_dof implemented here gives the same results as the
31 professional version, default was set to options_element_dof -yes.
32
33 For nonlocal calculation try to grep -i for:
34 options_nonlocal_softvar
35 materi_plasti_softvar_local, materi_plasti_softvar_nonlocal,
36 element_nonlocal_weight, element_nonlocal_ipoint, nonlocal_element_info
37 find_local_softvar, nonlocal_first_set, find_nonlocal_weights,
38 hope it is enough
39
40 To print materi_plasti_incremental_substeps into output change input.cc,
41 database.cc, tochnog.h, initia.cc, general.cc, materi.cc, plasti_incr.cc. Search
42 for the word 'substeps'
43
44 groundflow_addtopressure introduced in order to model hydrostatic pore pressure distribution
45 above phreatic level. Added into tochnog.h, database.cc and groundfl.cc. Search for 'addtopressure'.
46
47 Messages from SuperLU solver were suppressed in get_perm_c.c and dgssv.c. Labelled by
48 comment 'commented out!' :-). (Just in my copy)
49
50 David Masin masin@natur.cuni.cz
51
52 This program is free softvare; you can redistribute it and/or modify
53 it under the terms of the GNU General Public License as published by
54 the Free softvare Foundation; either version 2 of the License, or
55 (at your option) any later version.
56
57 This program is distributed in the hope that it will be useful,
58 but WITHOUT ANY WARRANTY; without even the implied warranty of
59 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
60 GNU General Public License for more details.
61
62
63 You should have received a copy of the GNU General Public License
64 along with this program; if not, write to the Free softvare Foundation
65 59 Temple Place, Suite 330, Boston, MA, 02111-1307, USA
66 */
67
68 /**********************header***************************************************************/
69
70 #include "tochnog.h"
71
72 const double PI=3.14159265358979323846264338327950288;
73 const double EULER = 2.7182818284590452354; /* e */
74 const double onethird=0.333333333333333333333333333333333333333333333333333333333333333333333333333 /*:-)*/;
75 const double LIMITF=1.e-5;
76 const int MAX_INTERSEC_ITER=40;
77
78 struct Variable {
79 double sig[MDIM*MDIM], new_sig[MDIM*MDIM],
80 geo_sig[MDIM*MDIM], d_epp[MDIM*MDIM], epp[MDIM*MDIM], d_ept[MDIM*MDIM];
81 double delta_sig_trial[MDIM*MDIM];
82 double df_trial[MDIM*MDIM];
83 double dgdsig[MDIM*MDIM];
84 double kron_delta[MDIM][MDIM];
85 double lambda_mltp [MDIM*MDIM];
86 double matrix_Cep[MDIM][MDIM][MDIM][MDIM];
87 double matrix_C[MDIM][MDIM][MDIM][MDIM];
88 long int element, usematrix;
89 double f, g, softv_nonl, softv_l, ddum[1];
90 bool plasti, nonloc;
91 vector<double> elasti_data;
92 vector<double> plasti_data;
93 vector<double> hisv;
94 vector<double> other_f;
95
VariableVariable96 Variable() {
97 inicialize();
98 }
99 void inicialize();
100 void recount();
101 };
102
103 struct F_euler_errorc_data {
104 double rotated_old_sig[MDIM*MDIM], inc_ept[MDIM*MDIM], init_dgdsig[MDIM*MDIM];
105 vector<double> old_hisv; //input
106 double double_sig[MDIM*MDIM], double_epp[MDIM*MDIM], double_dgdsig[MDIM*MDIM];
107 vector<double> double_hisv; //output after one step
108 double single_sig[MDIM*MDIM], single_epp[MDIM*MDIM], single_dgdsig[MDIM*MDIM];
109 vector<double> single_hisv; //output after half step
110 double total_inc_epp[MDIM*MDIM], df_init[9];
111 double init_lambda[MDIM*MDIM], single_lambda[MDIM*MDIM], double_lambda[MDIM*MDIM];
112 bool largeerror;
113 double num_substeps;
114 double max_substeps;
115 double min_substeps;
116 double sigerror;
F_euler_errorc_dataF_euler_errorc_data117 F_euler_errorc_data() {
118 inicialize();
119 }
120 void inicialize();
121 };
122
123 enum Num_pl_model {
124 camclay_incremental,
125 tskh_ai_plasti,
126 tskh_std_plasti,
127 };
128
129 enum Integration_scheme {
130 f_euler_errorcont,
131 forw_eul_subst,
132 forward_eul,
133 };
134
135 struct Plasti_rule {
136 double df[9];
Plasti_rulePlasti_rule137 Plasti_rule() {};
give_fPlasti_rule138 virtual void give_f(double sig[MDIM*MDIM], double epp[9], double d_epp[9],
139 double &f, double &g, bool recount_hisv) {};
give_matrixPlasti_rule140 virtual void give_matrix(bool only_f, bool elpl_matrix)
141 {cout<<"Plasti doesn't work virtually ";};
142 virtual void make_matrix_Cep (double df[], double dg[], double H, bool elpl_matrix);
inic_hisvPlasti_rule143 virtual void inic_hisv(double hisv[], bool only_substep) {};
erase_rec_histPlasti_rule144 virtual void erase_rec_hist() {};
145 };
146
147 struct Count {
148 Plasti_rule* plasti;
CountCount149 Count(Plasti_rule* b) : plasti(b) {
150 };
151 virtual void count_elasti(bool trial);
152 virtual void count_plasti();
153 virtual bool controll_plasti();
154 virtual void calc_results();
155 };
156
157 struct Forward_euler : Count {
Forward_eulerForward_euler158 Forward_euler(Plasti_rule* bn) : Count(bn) {};
159 };
160
161 struct Forw_eul_substep : Count {
Forw_eul_substepForw_eul_substep162 Forw_eul_substep(Plasti_rule* bn) : Count(bn) {};
163 virtual void calc_results();
164 virtual void finish_substeps();
165 virtual void apply_errorc();
166 };
167
168 struct F_euler_errorc : Forw_eul_substep {
F_euler_errorcF_euler_errorc169 F_euler_errorc(Plasti_rule* bn) : Forw_eul_substep(bn) {};
170 void apply_errorc();
171 void calc_substeps();
172 };
173
174 enum Tskh_mode {
175 tskh_std,
176 tskh_genai,
177 };
178
179 struct Camclay_incremental : Plasti_rule {
Camclay_incrementalCamclay_incremental180 Camclay_incremental () : Plasti_rule() {};
181 void give_f(double sig[MDIM*MDIM], double epp[9], double d_epp[9],
182 double &f, double &g, bool recount_hisv);
183 void give_matrix(bool only_f, bool elpl_matrix);
184 void inic_hisv(double hisv[], bool only_substep);
185 };
186
187 struct Tskh_general_plasti : Plasti_rule {
188 Tskh_mode mode;
189 double m, kappa, lambda, T, S, psi, N, m_flow, r_m, r_mfl, m_cmp, m_flow_K0,
190 rolodeb, rolodefl, drolodenbdlodenb, drolodefldlodenb;
191 double sig_b[9], sig_a[9], I_0, sig_0[9], v;
Tskh_general_plastiTskh_general_plasti192 Tskh_general_plasti (Tskh_mode m) : Plasti_rule(), mode(m){};
193 void give_f(double sig[MDIM*MDIM], double epp[9], double d_epp[9],
194 double &f, double &g, bool recount_hisv);
195 void give_matrix(bool only_f, bool elpl_matrix);
196 void inic_hisv(double hisv[], bool only_substep);
197 void erase_rec_hist();
198 void calc_rolode(double lode, double beta, double &rolode, double &drolodedlode,
199 bool alsoderivative);
200 };
201
202 /**********************end header***********************************************/
203
204 vector<Plasti_rule*> vector_plasti;
205 vector<Count*> countv;
206 Integration_scheme ischem;
207
208 Camclay_incremental camclay_incr;
209 Tskh_general_plasti tskh_ai_pl(tskh_genai);
210 Tskh_general_plasti tskh_std_pl(tskh_std);
211
212 //accesible for everyone
213 Variable vari;
214 F_euler_errorc_data errorc;
215
inicialize()216 void Variable :: inicialize() {
217 array_set(sig, 0, MDIM*MDIM);
218 array_set(new_sig, 0, MDIM*MDIM);
219 array_set(d_epp, 0, MDIM*MDIM);
220 array_set(d_ept, 0, MDIM*MDIM);
221 array_set(epp , 0, MDIM*MDIM);
222 array_set(delta_sig_trial, 0, MDIM*MDIM);
223 array_set(df_trial, 0, MDIM*MDIM);
224 array_set(dgdsig, 0, MDIM*MDIM);
225 array_set(lambda_mltp, 0, MDIM*MDIM);
226 array_set(&matrix_C[0][0][0][0], 0, MDIM*MDIM*MDIM*MDIM);
227 array_set(&matrix_Cep[0][0][0][0], 0, MDIM*MDIM*MDIM*MDIM);
228 f=-10;
229 g=-1.e10;
230 softv_nonl=softv_l=ddum[0]=0;
231 plasti=nonloc=false;
232 element=0;
233 usematrix=-NO;
234 for(int i=0; i<MDIM; i++) {
235 for(int j=0; j<MDIM; j++) {
236 if(i==j) kron_delta[i][j]=1;
237 else kron_delta[i][j]=0;
238 }
239 }
240 recount();
241
242 vector_plasti.push_back(&camclay_incr);
243 vector_plasti.push_back(&tskh_ai_pl);
244 vector_plasti.push_back(&tskh_std_pl);
245 }
246
inicialize()247 void F_euler_errorc_data :: inicialize() {
248 array_set(rotated_old_sig, 0, MDIM*MDIM);
249 array_set(inc_ept, 0, MDIM*MDIM);
250 array_set(double_sig, 0, MDIM*MDIM);
251 array_set(double_epp, 0, MDIM*MDIM);
252 array_set(single_sig, 0, MDIM*MDIM);
253 array_set(single_epp, 0, MDIM*MDIM);
254 array_set(init_dgdsig, 0, MDIM*MDIM);
255 array_set(single_dgdsig, 0, MDIM*MDIM);
256 array_set(double_dgdsig, 0, MDIM*MDIM);
257 array_set(total_inc_epp, 0, MDIM*MDIM);
258 array_set(df_init, 0, MDIM*MDIM);
259 array_set(init_lambda, 0, MDIM*MDIM);
260 array_set(single_lambda, 0, MDIM*MDIM);
261 array_set(double_lambda, 0, MDIM*MDIM);
262 largeerror=true;
263 num_substeps=1;
264 max_substeps=5000000;
265 min_substeps=2;
266 sigerror=0.00001;
267 }
268
recount()269 void Variable::recount() {
270 array_multiply(sig, geo_sig, -1, MDIM*MDIM);
271 }
272
273
plasti_incr(double rotated_old_sig[],double inc_ept[],double old_epp[],double new_ept[],double C[3][3][3][3],double old_hisv[],long int plasti_type,double plasti_dt[],long lenght_pl,long int gr,long int element,double softvar_nonl,double & softvar_l,double new_sigv[],double inc_epp[],double new_hisv[],double & new_f,double & new_substeps,double Cep_cons[3][3][3][3])274 void plasti_incr( double rotated_old_sig[], double inc_ept[], double old_epp[],
275 double new_ept[], double C[3][3][3][3], double old_hisv[], long int plasti_type,
276 double plasti_dt[], long lenght_pl, long int gr, long int element,
277 double softvar_nonl, double &softvar_l,
278 double new_sigv[], double inc_epp[], double new_hisv[], double &new_f, double &new_substeps,
279 double Cep_cons[3][3][3][3]) {
280
281 /*******************************************************
282 This function is an iterface between single element program
283 Triax (D.Masin@city.ac.uk) and Tochnog.
284 *******************************************************/
285 Num_pl_model num_plasti=tskh_std_plasti;
286 double elasti_dt[DATA_ITEM_SIZE], ddum[1]={0};
287 for(int i=0; i<DATA_ITEM_SIZE; i++) elasti_dt[i]=0;
288 long lenght_el=0;
289 long int ldum=0, idum[1]={0};
290
291 if(plasti_type==GROUP_MATERI_PLASTI_TSKH) {
292 num_plasti=tskh_std_plasti;
293 if(!get_group_data( GROUP_MATERI_ELASTI_TSKH, gr, element, ddum,
294 elasti_dt, lenght_el, GET_IF_EXISTS )) {
295 cout<<"3-SKH model must be used with the 3-SKH elasticity"<<endl;
296 exit(1);
297 }
298 }
299 else if(plasti_type==GROUP_MATERI_PLASTI_AITSKH) {
300 num_plasti=tskh_ai_plasti;
301 if(!get_group_data( GROUP_MATERI_ELASTI_TSKH, gr, element, ddum,
302 elasti_dt, lenght_el, GET_IF_EXISTS )) {
303 cout<<"AI3-SKH model must be used with the 3-SKH elasticity"<<endl;
304 exit(1);
305 }
306 }
307 else if(plasti_type==GROUP_MATERI_PLASTI_CAMCLAY_INCREMENTAL) {
308 num_plasti=camclay_incremental;
309 }
310 else {
311 cout<<"New model was not added into plasti_incr"<<endl;
312 exit(1);
313 }
314 F_euler_errorc f_euler_errorc(vector_plasti[num_plasti]);
315 Forw_eul_substep forw_eul_substep(vector_plasti[num_plasti]);
316 Forward_euler forward_euler(vector_plasti[num_plasti]);
317 countv.push_back(&f_euler_errorc);
318 countv.push_back(&forw_eul_substep);
319 countv.push_back(&forward_euler);
320
321 //ischem=forward_eul; // not in use actually - without substeps
322 ischem=f_euler_errorcont; // default with error 0.00001
323
324 if(get_group_data( GROUP_MATERI_PLASTI_INCREMENTAL_FESUBSTEPS, gr, element, ddum,
325 &errorc.num_substeps, ldum, GET_IF_EXISTS )) {
326 ischem=forw_eul_subst;
327 }
328 else if(get_group_data( GROUP_MATERI_PLASTI_INCREMENTAL_FEERROR, gr, element, ddum,
329 &errorc.sigerror, ldum, GET_IF_EXISTS )) {
330 ischem=f_euler_errorcont;
331 }
332 if(find_local_softvar) {
333 ischem=forw_eul_subst;
334 errorc.num_substeps=1;
335 }
336
337 get_group_data( GROUP_MATERI_PLASTI_INCREMENTAL_MAXSUBSTEPS, gr, element, ddum,
338 &errorc.max_substeps, ldum, GET_IF_EXISTS );
339 get_group_data( GROUP_MATERI_PLASTI_INCREMENTAL_MINSUBSTEPS, gr, element, ddum,
340 &errorc.min_substeps, ldum, GET_IF_EXISTS );
341 if (db_active_index( GROUP_MATERI_PLASTI_INCREMENTAL_USEMATRIX, gr, VERSION_NORMAL) )
342 db( GROUP_MATERI_PLASTI_INCREMENTAL_USEMATRIX, gr, &vari.usematrix, ddum, ldum,
343 VERSION_NORMAL, GET );
344 if (scalar_dabs(options_nonlocal_softvar)>TINY) {
345 //vari.usematrix=-NO;
346 if(!find_local_softvar) vari.nonloc=true;
347 else vari.nonloc=false;
348 }
349
350 //Feeding vari with current values
351 for(int i=0; i<lenght_el; i++) vari.elasti_data.push_back(elasti_dt[i]);
352 for(int i=0; i<lenght_pl; i++) vari.plasti_data.push_back(plasti_dt[i]);
353 array_move(rotated_old_sig, vari.sig, MDIM*MDIM);
354 array_move(rotated_old_sig, vari.new_sig, MDIM*MDIM);
355 array_move(inc_ept, vari.d_ept, MDIM*MDIM);
356 array_move(&C[0][0][0][0], &vari.matrix_C[0][0][0][0], MDIM*MDIM*MDIM*MDIM);
357 countv[ischem]->plasti->inic_hisv(old_hisv, false);
358 vari.element=element;
359 vari.softv_nonl=softvar_nonl;
360 vari.softv_l=softvar_l;
361 vari.recount();
362
363 //Erase the influence of recent history
364 if ( db_active_index( GROUP_MATERI_PLASTI_INCREMENTAL_ERASERECENTHISTORY,
365 gr, VERSION_NORMAL ) ) {
366 long int eraserecenthistory=-NO;
367 double time[1]={1};
368 db( GROUP_MATERI_PLASTI_INCREMENTAL_ERASERECENTHISTORY, gr,
369 &eraserecenthistory, ddum, ldum, VERSION_NORMAL, GET );
370 db( TIME_CURRENT, 0, idum, time, ldum, VERSION_NORMAL, GET );
371 if(time[0]==0 && eraserecenthistory==-YES) countv[ischem]->plasti->erase_rec_hist();
372 }
373
374 //calculation
375
376 countv[ischem]->calc_results();
377
378 long int printf;
379 if (db_active_index( GROUP_MATERI_PLASTI_INCREMENTAL_PRINTF, gr, VERSION_NORMAL) ) {
380 db( GROUP_MATERI_PLASTI_INCREMENTAL_PRINTF, gr, &printf, ddum, ldum,
381 VERSION_NORMAL, GET );
382 if ( printf==element && !find_local_softvar)
383 cout<<vari.f<<" "<<vari.other_f[0]<<" "<<vari.other_f[1]<<endl;
384 }
385
386 //Feeding output
387 array_move(vari.new_sig, new_sigv, MDIM*MDIM);
388 array_move(vari.d_epp, inc_epp, MDIM*MDIM);
389
390 double d_sig[MDIM*MDIM];
391 array_set(d_sig, 0, MDIM*MDIM);
392 array_subtract(vari.new_sig, rotated_old_sig, d_sig, MDIM*MDIM);
393
394 for(int i=0; i<int(vari.hisv.size()); i++)
395 new_hisv[i]=vari.hisv[i];
396 new_f=vari.f;
397
398 if(new_f>0) new_substeps=errorc.num_substeps;
399 else new_substeps=1;
400
401 array_move(&vari.matrix_Cep[0][0][0][0], &Cep_cons[0][0][0][0], MDIM*MDIM*MDIM*MDIM);
402
403 //cleaning vari for the next step
404 vari.elasti_data.clear();
405 vari.plasti_data.clear();
406 vari.hisv.clear();
407 vari.other_f.clear();
408 countv.clear();
409 }
410
calc_results()411 void Count::calc_results() {
412 if(controll_plasti()==true) count_plasti();
413 else count_elasti(false);
414
415 //Calculates new history variables
416 if(vari.plasti) plasti->give_f(vari.sig, vari.epp, vari.d_epp, vari.ddum[0], vari.ddum[0], true);
417 }
418
controll_plasti()419 bool Count::controll_plasti() {
420 double dfdsig;
421 plasti->give_f(vari.sig, vari.epp, vari.d_epp, vari.f, vari.g, false);
422 count_elasti(true); //calculates trial stress rate
423 if (vari.f>=0) {
424 plasti->give_matrix(false, false); //calculates df/dsig
425 dfdsig=array_inproduct(vari.df_trial, vari.delta_sig_trial, MDIM*MDIM);
426 if( dfdsig > 0) vari.plasti=true;
427 else vari.plasti=false;
428 }
429 else vari.plasti = false;
430 return(vari.plasti);
431 }
432
count_elasti(bool trial)433 void Count::count_elasti(bool trial) {
434 double delta_epe[MDIM*MDIM], delta_sig[MDIM*MDIM];
435 array_set(delta_epe, 0, MDIM*MDIM);
436 array_set(delta_sig, 0, MDIM*MDIM);
437 array_move(vari.d_ept, delta_epe, MDIM*MDIM);
438 matrix_a4b(vari.matrix_C, delta_epe, delta_sig);
439
440 if(trial) array_move(delta_sig, vari.delta_sig_trial, MDIM*MDIM);
441 else array_add(vari.sig, delta_sig, vari.new_sig, MDIM*MDIM);
442
443 array_set(vari.d_epp, 0, MDIM*MDIM);
444 }
445
count_plasti()446 void Count::count_plasti() {
447 array_move(vari.sig, vari.new_sig, MDIM*MDIM);
448 double delta_ept[MDIM*MDIM];
449 array_set(delta_ept, 0, MDIM*MDIM);
450
451 double delta_sig[MDIM*MDIM];
452 array_set(delta_sig, 0, MDIM*MDIM);
453 double delta_epe[MDIM*MDIM];
454 array_set(delta_epe, 0, MDIM*MDIM);
455
456 double delta_epp[MDIM*MDIM];
457 array_set(delta_epp, 0, MDIM*MDIM);
458
459 array_move(vari.d_ept, delta_ept, MDIM*MDIM);
460
461 double used_lambda=array_inproduct(vari.lambda_mltp, vari.d_ept, MDIM*MDIM);
462 array_multiply(vari.dgdsig, delta_epp, used_lambda, MDIM*MDIM);
463 array_subtract(delta_ept, delta_epp, delta_epe, MDIM*MDIM);
464 matrix_a4b(vari.matrix_C, delta_epe, delta_sig);
465
466 array_add(vari.sig, delta_sig, vari.new_sig, MDIM*MDIM);
467 array_move(delta_epp, vari.d_epp, 9);
468 }
469
470 /******Forward Euler, const given number of substeps*************/
471
calc_results()472 void Forw_eul_substep::calc_results() {
473 array_set(errorc.total_inc_epp, 0, MDIM*MDIM);
474
475 controll_plasti();
476
477 if(!vari.plasti) {
478 double testf_after = -10;
479 count_elasti(false); //is stress outside after elastic step?
480 plasti->give_f(vari.new_sig, vari.epp, vari.d_epp, testf_after, vari.ddum[0], false);
481
482 //if yes -> calculates initial intersection with the yield surface
483 if(testf_after>0 && vari.f<0) {
484 double delta_sig_e[MDIM*MDIM];
485 double sig_init[MDIM*MDIM];
486 array_subtract(vari.new_sig, vari.sig, delta_sig_e, MDIM*MDIM);
487 array_move(vari.sig, sig_init, MDIM*MDIM);
488 double alpha=-vari.f/(testf_after-vari.f);
489 double dalpha=0, dftrdsige=0;
490 double num_iter=1;
491 double dsige_times_alpha[MDIM*MDIM];
492 double df_transposed[MDIM*MDIM];
493 array_set(dsige_times_alpha, 0, MDIM*MDIM);
494 array_set(df_transposed, 0, MDIM*MDIM);
495
496 //newton-raphson iterations for scalar multiplier alpha
497 while(1) {
498 array_multiply(delta_sig_e, dsige_times_alpha, alpha, MDIM*MDIM);
499 array_add(sig_init, dsige_times_alpha, vari.sig, MDIM*MDIM);
500 vari.recount();
501 plasti->give_f(vari.sig, vari.epp, vari.d_epp, vari.f, vari.ddum[0], false);
502 plasti->give_matrix(true, false); //calculates df/dsig_i
503
504 if(num_iter>MAX_INTERSEC_ITER || (vari.f>0 && vari.f<LIMITF)) break;
505
506 for(int i=0; i<MDIM; i++) { //transpose df/dsig_i
507 for(int j=0; j<MDIM; j++) {
508 df_transposed[3*i+j]=vari.df_trial[3*j+i];
509 }
510 }
511 dftrdsige=array_inproduct(df_transposed, delta_sig_e, MDIM*MDIM);
512 if(!dftrdsige==0) dalpha=-vari.f/dftrdsige;
513 else dalpha=0;
514 num_iter++;
515 alpha+=dalpha;
516 }
517 //new d_ept for el-pl time integration, vari.sig is now on yield surface
518 array_multiply(vari.d_ept, vari.d_ept, (1-alpha), MDIM*MDIM);
519 vari.plasti=true;
520 }
521 }
522 if(vari.plasti) {
523 apply_errorc();
524 array_set(errorc.total_inc_epp, 0, MDIM*MDIM);
525 array_multiply(vari.d_ept, vari.d_ept, (1/errorc.num_substeps), MDIM*MDIM);
526
527 for(int i=0; i<errorc.num_substeps; i++)
528 finish_substeps();
529
530 array_move(vari.new_sig, vari.sig, MDIM*MDIM);
531 }
532
533 if(vari.plasti && vari.usematrix==-YES) plasti->give_matrix(false, true);
534 else array_move(&vari.matrix_C[0][0][0][0], &vari.matrix_Cep[0][0][0][0], MDIM*MDIM*MDIM*MDIM);
535 array_move(errorc.total_inc_epp, vari.d_epp, MDIM*MDIM);
536 }
537
finish_substeps()538 void Forw_eul_substep::finish_substeps() {
539 vari.recount();
540 plasti->give_f(vari.sig, vari.epp, vari.d_epp, vari.f, vari.g, false); //before step->calc f
541 plasti->give_matrix(false, false);
542 count_plasti();
543 plasti->give_f(vari.sig, vari.epp, vari.d_epp, vari.ddum[0], vari.ddum[0], true); //after step-> calc hisv
544 array_move(vari.new_sig, vari.sig, MDIM*MDIM);
545 array_add(errorc.total_inc_epp, vari.d_epp, errorc.total_inc_epp, MDIM*MDIM);
546 }
547
apply_errorc()548 void Forw_eul_substep::apply_errorc() {
549 }
550
551 /******************Forward Euler, control of local error, automatic substepping*********************/
552
apply_errorc()553 void F_euler_errorc::apply_errorc() {
554 errorc.largeerror=true;
555 errorc.num_substeps=2;
556
557 //initial values
558 array_move(vari.sig, errorc.rotated_old_sig, MDIM*MDIM);
559 array_move(vari.d_ept, errorc.inc_ept, MDIM*MDIM);
560 for(int i=0; i<int(vari.hisv.size()); i++) errorc.old_hisv.push_back(vari.hisv[i]);
561 for(int i=0; i<int(vari.hisv.size()); i++) errorc.single_hisv.push_back(vari.hisv[i]);
562 for(int i=0; i<int(vari.hisv.size()); i++) errorc.double_hisv.push_back(vari.hisv[i]);
563 plasti->give_matrix(false, false);
564 array_move(vari.dgdsig, errorc.init_dgdsig, MDIM*MDIM);
565 array_move(vari.lambda_mltp, errorc.init_lambda, MDIM*MDIM);
566 array_move(plasti->df, errorc.df_init, MDIM*MDIM);
567
568 count_plasti();
569 plasti->give_f(vari.sig, vari.epp, vari.d_epp, vari.f, vari.ddum[0], true); //new_hisv
570
571 //saves results after one step
572 array_move(vari.new_sig, errorc.double_sig, MDIM*MDIM);
573 array_move(vari.d_epp, errorc.double_epp, MDIM*MDIM);
574 for(int i=0; i<int(vari.hisv.size()); i++) errorc.double_hisv[i]=vari.hisv[i];
575
576 while(errorc.largeerror) calc_substeps();
577
578 //refeed with init. values
579 array_move(errorc.rotated_old_sig, vari.sig, MDIM*MDIM);
580 array_move(errorc.inc_ept, vari.d_ept, MDIM*MDIM);
581 array_set(errorc.total_inc_epp, 0, MDIM*MDIM);
582 for(int i=0; i<int(vari.hisv.size()); i++) vari.hisv[i]=errorc.old_hisv[i];
583 plasti->inic_hisv(vari.ddum, true);
584
585 errorc.old_hisv.clear();
586 errorc.double_hisv.clear();
587 errorc.single_hisv.clear();
588 }
589
calc_substeps()590 void F_euler_errorc::calc_substeps() {
591
592 //refeed vari with init values
593 array_move(errorc.rotated_old_sig, vari.sig, MDIM*MDIM);
594 array_move(errorc.df_init, plasti->df, MDIM*MDIM);
595 array_move(errorc.init_dgdsig, vari.dgdsig, MDIM*MDIM);
596 array_move(errorc.init_lambda, vari.lambda_mltp, MDIM*MDIM);
597 array_multiply(vari.d_ept, vari.d_ept, 0.5, MDIM*MDIM); //half step size
598 for(int i=0; i<int(vari.hisv.size()); i++) vari.hisv[i]=errorc.old_hisv[i];
599 plasti->inic_hisv(vari.ddum, true);
600
601 //calc first substep
602 plasti->give_f(vari.sig, vari.epp, vari.d_epp, vari.f, vari.g, false);
603 count_plasti();
604 plasti->give_f(vari.sig, vari.epp, vari.d_epp, vari.ddum[0], vari.ddum[0], true);
605 array_move(vari.new_sig, vari.sig, MDIM*MDIM);
606
607 //saves results after substep
608 array_move(vari.new_sig, errorc.single_sig, MDIM*MDIM);
609 array_move(vari.d_epp, errorc.single_epp, MDIM*MDIM);
610 for(int i=0; i<int(vari.hisv.size()); i++) errorc.single_hisv[i]=vari.hisv[i];
611 vari.recount();
612 plasti->give_matrix(false, false);
613 array_move(vari.dgdsig, errorc.single_dgdsig, MDIM*MDIM);
614 array_move(vari.lambda_mltp, errorc.single_lambda, MDIM*MDIM);
615
616 //calc second substep
617 plasti->give_f(vari.sig, vari.epp, vari.d_epp, vari.f, vari.ddum[0], false);
618 count_plasti();
619 plasti->give_f(vari.sig, vari.epp, vari.d_epp, vari.ddum[0], vari.ddum[0], true);
620 array_move(vari.new_sig, vari.sig, MDIM*MDIM);
621 array_add(errorc.total_inc_epp, vari.d_epp, errorc.total_inc_epp, MDIM*MDIM);
622
623 //calc local trucation error
624 double errorsig[MDIM*MDIM];
625 array_set(errorsig, 0, MDIM*MDIM);
626 array_subtract(errorc.double_sig, vari.new_sig, errorsig, MDIM*MDIM);
627 double error = sqrt(array_inproduct(errorsig, errorsig, MDIM*MDIM));
628
629 if((error<errorc.sigerror|| errorc.num_substeps>errorc.max_substeps)&&
630 (errorc.num_substeps>(errorc.min_substeps-0.1))) errorc.largeerror=false;
631 else {
632 errorc.largeerror=true;
633
634 //half step is main step for next turn
635 array_move(errorc.single_sig, errorc.double_sig, MDIM*MDIM);
636 array_move(errorc.single_epp, errorc.double_epp, MDIM*MDIM);
637 for(int i=0; i<int(vari.hisv.size()); i++) errorc.double_hisv[i]=errorc.single_hisv[i];
638 array_move(errorc.single_dgdsig, errorc.double_dgdsig, MDIM*MDIM);
639 array_move(errorc.single_lambda, errorc.double_lambda, MDIM*MDIM);
640
641 array_set(errorc.total_inc_epp, 0, MDIM*MDIM);
642 errorc.num_substeps*=2;
643 }
644 }
645
646 /********************Elasto-plastic stiffness matrix***********************************/
647
make_matrix_Cep(double df[],double dg[],double H,bool elpl_matrix)648 void Plasti_rule :: make_matrix_Cep (double df[], double dg[], double H,
649 bool elpl_matrix) {
650
651 double tmp1[MDIM*MDIM], tmp2[MDIM*MDIM], tmp3[MDIM][MDIM][MDIM][MDIM], chi=1;
652 double Cp[MDIM][MDIM][MDIM][MDIM];
653
654 matrix_a4b(vari.matrix_C, dg, tmp1);
655 matrix_ab4(df, vari.matrix_C, tmp2);
656 matrix4_ab(tmp1, tmp2, tmp3);
657 matrix_a_contr_b(df, tmp1, chi);
658
659 chi=chi+H;
660 if(chi==0) {
661 cout<<" Problem in calculation of elasto--plastic stiffness matrix -> chi=0 ";
662 //exit_tn( -YES );
663 }
664 array_multiply(tmp2, vari.lambda_mltp, 1/chi, MDIM*MDIM);
665
666 if(elpl_matrix) {
667 for (int i=0; i<3; i++ ) {
668 for (int j=0; j<3; j++ ) {
669 for (int k=0; k<3; k++ ) {
670 for (int l=0; l<3; l++ ) {
671 Cp[i][j][k][l] = tmp3[i][j][k][l]/chi;
672 }
673 }
674 }
675 }
676
677 for (int i=0; i<3; i++ ) {
678 for (int j=0; j<3; j++ ) {
679 for (int k=0; k<3; k++ ) {
680 for (int l=0; l<3; l++ ) {
681 vari.matrix_Cep[i][j][k][l] = vari.matrix_C[i][j][k][l] - Cp[i][j][k][l];
682 }
683 }
684 }
685 }
686 }
687 array_move(dg, vari.dgdsig, MDIM*MDIM);
688 }
689
690 /***************************MODELS***********************************************/
691 /******** Cam-Clay -- incremental solution for comparison with iterative ********/
692
inic_hisv(double hisv[],bool only_substep)693 void Camclay_incremental::inic_hisv(double hisv[], bool only_substep) {
694 array_set(df, 0, 9);
695 vari.hisv.push_back(hisv[0]);
696 vari.hisv.push_back(hisv[1]);
697 }
698
give_f(double sig[MDIM * MDIM],double epp[9],double d_epp[9],double & f,double & g,bool recount_hisv)699 void Camclay_incremental :: give_f(double sig[MDIM*MDIM], double epp[9], double d_epp[9],
700 double &f, double &g, bool recount_hisv) {
701
702 double geo_sigd[9], d_eppd[9];
703 array_set(geo_sigd, 0, 9);
704 array_set(d_eppd, 0, 9);
705 array_multiply(sig, geo_sigd, -1, MDIM*MDIM);
706 array_multiply(d_epp, d_eppd, -1, MDIM*MDIM);
707 double m = vari.plasti_data[0];
708 double kappa = vari.plasti_data[1];
709 double lambda = vari.plasti_data[2];
710 double N = vari.plasti_data[3];
711 if ( scalar_dabs(lambda-kappa)==0. ) {
712 cout<<"Error in Cam--Clay: lambda==kappa";
713 exit(1);
714 }
715 double dev = vari.d_ept[0]+vari.d_ept[4]+vari.d_ept[8];
716 double p = (geo_sigd[0]+geo_sigd[4]+geo_sigd[8])/3;
717 double e = vari.hisv[0];
718 //double p0 = vari.hisv[1];
719 double p0 = exp((N-kappa*log(p)-(1+e))/(lambda-kappa));
720 double devp = (d_eppd[0]+d_eppd[4]+d_eppd[8]);
721 double de = dev*(1.+e);
722 double dp0 = devp * p0 * (1.+e)/(lambda-kappa);
723 e += de;
724 p0 += dp0;
725 double q = sqrt(
726 0.5*( scalar_square(geo_sigd[0]-geo_sigd[4])+
727 scalar_square(geo_sigd[4]-geo_sigd[8])+
728 scalar_square(geo_sigd[0]-geo_sigd[8]) ) +
729 3.*( scalar_square(geo_sigd[1]) +
730 scalar_square(geo_sigd[2]) +
731 scalar_square(geo_sigd[5]) ) );
732 f = q*q - m*m*(p*(p0-p));
733 g = q*q - m*m*(p*(p0-p));
734
735 if(recount_hisv) {
736 vari.hisv[0] += de;
737 if (vari.f>0) vari.hisv[1] = p0;
738 }
739 }
740
give_matrix(bool only_f,bool elpl_matrix)741 void Camclay_incremental::give_matrix(bool only_f, bool elpl_matrix) {
742 double m = vari.plasti_data[0];
743 double kappa = vari.plasti_data[1];
744 double lambda = vari.plasti_data[2];
745 double N = vari.plasti_data[3];
746 double e = vari.hisv[0];
747 double p = (vari.geo_sig[0]+vari.geo_sig[4]+vari.geo_sig[8])/3;
748 //double p0 = vari.hisv[1];
749 double p0 = exp((N-kappa*log(p)-(1+e))/(lambda-kappa));
750
751 double df[MDIM*MDIM], dg[MDIM*MDIM], sigd[MDIM*MDIM], zero[MDIM*MDIM], eppd[MDIM*MDIM];
752 array_set(zero, 0, MDIM*MDIM);
753 array_set(dg, 0, MDIM*MDIM);
754 array_set(df, 0, MDIM*MDIM);
755 array_set(sigd, 0, MDIM*MDIM);
756 array_set(eppd, 0, MDIM*MDIM);
757 array_move(vari.sig, sigd, MDIM*MDIM);
758 double f_right=0, f_left=0;
759
760 double signorm=array_size(vari.sig, MDIM*MDIM)/100000;
761 if(signorm==0) signorm=0.0000001;
762 for(int i=0; i<MDIM*MDIM; i++) {
763 sigd[i] += signorm;
764 give_f(sigd, zero, zero, f_right, f_right, false);
765 sigd[i] -= 2*signorm;
766 give_f(sigd, zero, zero, f_left, f_left, false);
767 df[i]=(f_right-f_left)/(2*signorm);
768 array_move(vari.sig, sigd, MDIM*MDIM);
769 }
770 array_move(df, dg, MDIM*MDIM);
771
772 double dfdep[MDIM*MDIM];
773 array_set(dfdep, 0, MDIM*MDIM);
774
775 dfdep[0]=dfdep[4]=dfdep[8]=(1+e)*m*m*p*p0/(lambda-kappa);
776
777 double H=0;
778 for(int i=0; i<MDIM*MDIM; i++) {
779 H -= dfdep[i]*dg[i];
780 }
781 array_move(df, vari.df_trial, MDIM*MDIM);
782 if(!only_f) make_matrix_Cep(df, dg, H, elpl_matrix);
783 }
784
785 /***************************3-SKH********************************************************/
786
erase_rec_hist()787 void Tskh_general_plasti::erase_rec_hist() {
788 vari.hisv[2]=vari.hisv[8]=-vari.sig[0];
789 vari.hisv[3]=vari.hisv[9]=-vari.sig[1];
790 vari.hisv[4]=vari.hisv[10]=-vari.sig[2];
791 vari.hisv[5]=vari.hisv[11]=-vari.sig[4];
792 vari.hisv[6]=vari.hisv[12]=-vari.sig[5];
793 vari.hisv[7]=vari.hisv[13]=-vari.sig[8];
794
795 inic_hisv(vari.ddum, true);
796 }
797
inic_hisv(double hisv[],bool only_substep)798 void Tskh_general_plasti::inic_hisv(double hisv[], bool only_substep) {
799
800 if(!only_substep) {
801 m_cmp = vari.plasti_data[0];
802 kappa = vari.plasti_data[1];
803 lambda = vari.plasti_data[2];
804 T = vari.plasti_data[3];
805 S = vari.plasti_data[4];
806 psi = vari.plasti_data[5];
807 N = vari.plasti_data[6];
808 m_flow = m = m_flow_K0 = m_cmp;
809
810 r_m=r_mfl=3/(3+m_cmp);
811
812 array_set(df, 0, 9);
813 array_set(sig_b, 0, 9);
814 array_set(sig_a, 0, 9);
815 array_set(sig_0, 0, 9);
816
817 for(int i=0; i<14; i++) vari.hisv.push_back(hisv[i]);
818
819 vari.other_f.push_back(-10);
820 vari.other_f.push_back(-10);
821 }
822
823 v=vari.hisv[0]+1; //e stored as a history variable
824 I_0=3*vari.hisv[1];
825 sig_0[0]=sig_0[4]=sig_0[8]=I_0/3;
826 sig_a[0]=vari.hisv[2];
827 sig_a[1]=sig_a[3]=vari.hisv[3];
828 sig_a[2]=sig_a[6]=vari.hisv[4];
829 sig_a[4]=vari.hisv[5];
830 sig_a[5]=sig_a[7]=vari.hisv[6];
831 sig_a[8]=vari.hisv[7];
832 sig_b[0]=vari.hisv[8];
833 sig_b[1]=sig_b[3]=vari.hisv[9];
834 sig_b[2]=sig_b[6]=vari.hisv[10];
835 sig_b[4]=vari.hisv[11];
836 sig_b[5]=sig_b[7]=vari.hisv[12];
837 sig_b[8]=vari.hisv[13];
838 }
839
give_f(double sig[MDIM * MDIM],double epp[9],double d_epp[9],double & f,double & g,bool recount_hisv)840 void Tskh_general_plasti :: give_f(double sig[MDIM*MDIM], double epp[9], double d_epp[9],
841 double &f, double &g, bool recount_hisv) {
842
843 double geo_sig[9], d_geo_epp[9];
844 array_set(geo_sig, 0, 9);
845 array_set(d_geo_epp, 0, 9);
846 array_multiply(sig, geo_sig, -1, MDIM*MDIM);
847 array_multiply(d_epp, d_geo_epp, -1, MDIM*MDIM);
848 double f_hist=vari.other_f[0];
849 double f_bound=vari.other_f[1];
850 double f_yield=vari.f;
851 double dsig_b[9], dsig_a[9];
852 if ( scalar_dabs(lambda-kappa)==0. ) {
853 cout<<"Problem in the 3-SKH model: lambda must differ from kappa";
854 exit(1);
855 }
856
857 double p = (geo_sig[0]+geo_sig[4]+geo_sig[8])/3;
858 if(p<=0) p=TINY;
859 double A = vari.elasti_data[0];
860 double n= vari.elasti_data[1];
861 double m_elasti = vari.elasti_data[2];
862 double R0=2*I_0 / (3*p);
863 double G = A*scalar_power(p, n)*scalar_power(R0, m_elasti);
864 double K = p/kappa;
865 double poisson = (3*K - 2*G)/(2*G + 6*K);
866 if(poisson>=0.5) {
867 cout<<"error in AI3-SKH, location 1, element "<<vari.element<<endl;
868 //exit_tn( -YES );
869 }
870
871 //m_cmp = vari.plasti_data[0];
872
873 double mfl_dnm=
874 (((1-2*poisson)*(6-m_cmp)*lambda-m_cmp*kappa*(1+poisson))*((6-m_cmp)*(6-m_cmp)-9));
875 if(mfl_dnm<=0) {
876 cout<<"error in AI3-SKH, location 2, element "<<vari.element<<endl;
877 //exit_tn( -YES );
878 }
879 m_flow_K0=3*(6-m_cmp)*sqrt((1-2*poisson)*(lambda-kappa)*m_cmp/mfl_dnm);
880 if(m_flow_K0<=0) {
881 cout<<"error in AI3-SKH, location 3, element "<<vari.element<<endl;
882 //exit_tn( -YES );
883 }
884
885 if(mode==tskh_std) {
886 r_m=1;
887 r_mfl=1;
888 m_flow_K0=m_cmp;
889 }
890
891 double sig_sigb[9];
892 double sig_siga[9];
893 array_set(sig_sigb, 0, 9);
894 array_set(sig_siga, 0, 9);
895 array_subtract(geo_sig, sig_b, sig_sigb, 9);
896 array_subtract(geo_sig, sig_a, sig_siga, 9);
897
898 if(mode==tskh_genai) {
899 r_mfl=m_cmp/m_flow_K0;
900 // In 2D y axis is vertical
901 // in 3D z axis is vertical
902 if(ndim==2) {
903 if(!((sig_sigb[4]>sig_sigb[8])&&(sig_sigb[4]>sig_sigb[0]))) {
904 r_mfl=1;
905 m_flow_K0=m_cmp;
906 }
907 }
908 else if(ndim==3) {
909 if(!((sig_sigb[8]>sig_sigb[4])&&(sig_sigb[8]>sig_sigb[0]))) {
910 r_mfl=1;
911 m_flow_K0=m_cmp;
912 }
913 }
914 else {
915 cout<<"Can not use AI3-SKH in 1D"<<endl;
916 exit_tn( -YES );
917 }
918 }
919 /*if(mode==tskh_genai) { //MN3-SKH
920 r_mfl=1;
921 m_flow_K0=m_cmp;
922 }*/
923
924 if(recount_hisv && vari.plasti) {
925
926 /************* isotropic hardening ******************/
927
928 double devp = (d_geo_epp[0]+d_geo_epp[4]+d_geo_epp[8]);
929 double dI0 = devp * I_0/(lambda-kappa);
930
931 if(I_0<=0) {
932 cout<<"error in AI3-SKH, location 4, element "<<vari.element<<endl;
933 //exit_tn( -YES );
934 }
935
936 array_set(dsig_b, 0, 9);
937 array_set(dsig_a, 0, 9);
938 for (int j=0; j<9; j++ ) {
939 dsig_b[j]=dI0*sig_b[j]/I_0;
940 dsig_a[j]=dI0*sig_a[j]/I_0;
941 }
942 for (int j=0; j<9; j++ ) {
943 sig_b[j]+=dsig_b[j];
944 sig_a[j]+=dsig_a[j];
945 }
946 I_0+=dI0;
947 sig_0[0]=sig_0[4]=sig_0[8]=I_0/3;
948
949 /****************** direction of translation ******************/
950
951 double beta[9]; //direction for history surface
952 array_set(beta,0,9);
953 for (int j=0; j<9; j++ ) {
954 beta[j]=(vari.geo_sig[j]-sig_a[j])/T+sig_0[j]-vari.geo_sig[j];
955 }
956 double gama[9]; //direction for yield surface
957 array_set(gama,0,9);
958 for (int j=0; j<9; j++ ) {
959 gama[j]=(vari.geo_sig[j]-sig_b[j])/S+sig_a[j]-vari.geo_sig[j];
960 }
961
962 /************* kinematic hardening ******************/
963
964 array_set(dsig_b, 0, 9);
965 array_set(dsig_a, 0, 9);
966 double d_geosig[9];
967 array_set(d_geosig, 0, 9);
968 array_subtract(vari.new_sig, vari.sig, d_geosig, MDIM*MDIM);
969 array_multiply(d_geosig, d_geosig, -1, 9);
970 double tmpup=0, tmpdown=0;
971
972 if (f_yield>=0 && f_hist<0 && f_bound<0) {
973 tmpup=0;
974 tmpdown=0;
975 double dfdI0=-T*T*S*S*I_0/9;
976 for (int j=0; j<9; j++ ) {
977 tmpup +=-df[j]*d_geosig[j]+df[j]*dI0*sig_b[j]/I_0;
978 tmpdown +=-df[j]*gama[j]; //-df <- back to geotechnical convention
979 }
980 if(tmpdown==0) {
981 cout<<"error in AI3-SKH, location 5, element "<<vari.element<<endl;
982 //exit_tn( -YES );
983 }
984
985 double Zs=(tmpup+dfdI0*dI0)/tmpdown;
986 for (int j=0; j<9; j++ ) {
987 dsig_b[j]=Zs*gama[j];
988 }
989 for (int j=0; j<9; j++ ) {
990 sig_b[j]+=dsig_b[j];
991 }
992 }
993 else if (f_hist>=0 && f_bound < 0 && f_yield>=0) {
994 tmpup=0;
995 tmpdown=0;
996 double dfhdI0=-T*T*I_0/9;
997 double dfhdsig[9];
998 array_set(dfhdsig, 0, 9);
999 for (int j=0; j<9; j++ ) {
1000 dfhdsig[j]=df[j]/S;
1001 }
1002 for (int j=0; j<9; j++ ) {
1003 tmpup +=-dfhdsig[j]*d_geosig[j]+dfhdsig[j]*dI0*sig_a[j]/I_0;
1004 tmpdown +=-dfhdsig[j]*beta[j]; //-df <- back to geotechnical convention
1005 }
1006 if(tmpdown==0) {
1007 cout<<"error in AI3-SKH, location 6, element "<<vari.element<<endl;
1008 //exit_tn( -YES );
1009 }
1010 double Ws=(tmpup+dfhdI0*dI0)/tmpdown;
1011 for (int j=0; j<9; j++ ) {
1012 dsig_a[j]=Ws*beta[j];
1013 }
1014 for (int j=0; j<9; j++ ) {
1015 sig_a[j]+=dsig_a[j];
1016 sig_b[j]=geo_sig[j]-(geo_sig[j]-sig_a[j])*S;
1017 }
1018 }
1019 else if(f_bound >= 0) {
1020 //A position of the yield and hist. surf. on bounding s., touching at conj. points
1021 for (int j=0; j<9; j++ ) {
1022 sig_b[j]=geo_sig[j]-(geo_sig[j]-sig_0[j])*T*S;
1023 sig_a[j]=geo_sig[j]-(geo_sig[j]-sig_0[j])*T;
1024 }
1025 }
1026 else if(f_yield>0) {
1027 cout<<endl<<"some problem"<<endl;
1028 //exit_tn( -YES );
1029 }
1030
1031 /*********************************************/
1032
1033 //double dev = vari.d_ept[0]+vari.d_ept[4]+vari.d_ept[8];
1034 //double dv = dev*v;
1035 //v += dv;
1036
1037 vari.hisv[1]=I_0/3;
1038 vari.hisv[2]=sig_a[0];
1039 vari.hisv[3]=sig_a[1];
1040 vari.hisv[4]=sig_a[2];
1041 vari.hisv[5]=sig_a[4];
1042 vari.hisv[6]=sig_a[5];
1043 vari.hisv[7]=sig_a[8];
1044 vari.hisv[8]=sig_b[0];
1045 vari.hisv[9]=sig_b[1];
1046 vari.hisv[10]=sig_b[2];
1047 vari.hisv[11]=sig_b[4];
1048 vari.hisv[12]=sig_b[5];
1049 vari.hisv[13]=sig_b[8];
1050 }
1051 double lodeb=0, lodea=0, I=0, J=0, lode=0;
1052 double rolode=0, rolodea=0;
1053 double tmparray[3][3]={{0,0,0},{0,0,0},{0,0,0}};
1054 calc_IJlode(geo_sig, I, J, lode, false, tmparray, tmparray, tmparray);
1055 double Jnorm_a=0, Jnorm_b=0, Inorm_a=0, Inorm_b;
1056 calc_IJlode(sig_sigb, Inorm_b, Jnorm_b, lodeb, false, tmparray, tmparray, tmparray);
1057 calc_IJlode(sig_siga, Inorm_a, Jnorm_a, lodea, false, tmparray, tmparray, tmparray);
1058 calc_rolode(lode, r_m, rolode, vari.ddum[0], false);
1059 calc_rolode(lodeb, r_m, rolodeb, drolodenbdlodenb, true);
1060 calc_rolode(lodea, r_m, rolodea, vari.ddum[0], false);
1061 calc_rolode(lodeb, r_mfl, rolodefl, drolodefldlodenb, true);
1062
1063 m=m_cmp*rolodeb;
1064 m_flow=m_flow_K0*rolodefl;
1065 double m_bound=m_cmp*rolode;
1066 double m_hist=m_cmp*rolodea;
1067 double m_yield=m;
1068 if((m_bound<=0)||(m_hist<=0)||(m_yield<=0)||(m_flow<=0)) {
1069 cout<<"error in AI3-SKH, location 7, element "<<vari.element<<endl;
1070 //exit_tn( -YES );
1071 }
1072
1073 f_yield=(3*Jnorm_b*Jnorm_b/(m_yield*m_yield)+scalar_power(Inorm_b/3, 2)-
1074 T*T*S*S*I_0*I_0/9)/2;
1075 f_hist=(3*Jnorm_a*Jnorm_a/(m_hist*m_hist)+scalar_power(Inorm_a/3, 2)-
1076 T*T*I_0*I_0/9)/2;
1077 f_bound=(3*J*J/(m_bound*m_bound)+scalar_power((I-I_0)/3, 2)-
1078 I_0*I_0/9)/2;
1079 g=(3*Jnorm_b*Jnorm_b/(m_flow*m_flow)+scalar_power(Inorm_b/3, 2)-
1080 T*T*S*S*I_0*I_0/9)/2;
1081 f=f_yield;
1082
1083 vari.other_f[0] = f_hist;
1084 vari.other_f[1] = f_bound;
1085
1086 double powv=N-lambda*log(2*I_0/3)+kappa*log(2*I_0/I);
1087 v=scalar_power(EULER, powv);
1088 vari.hisv[0]=v-1;
1089 }
1090
give_matrix(bool only_f,bool elpl_matrix)1091 void Tskh_general_plasti::give_matrix(bool only_f, bool elpl_matrix) {
1092
1093 double dg[MDIM*MDIM];
1094 array_set(dg, 0, MDIM*MDIM);
1095 array_set(df, 0, MDIM*MDIM);
1096 double lodeb=0;
1097 double sig_sigb[9];
1098 array_set(sig_sigb, 0, 9);
1099 array_subtract(vari.geo_sig, sig_b, sig_sigb, 9);
1100 double Jnorm_b=0, Inorm_b;
1101 double dJnbdsig[3][3];
1102 double dInbdsig[3][3];
1103 double dlodenbdsig[3][3];
1104 array_set(*dJnbdsig, 0, 9);
1105 array_set(*dInbdsig, 0, 9);
1106 array_set(*dlodenbdsig, 0, 9);
1107 calc_IJlode(sig_sigb, Inorm_b, Jnorm_b, lodeb, true, dInbdsig, dJnbdsig, dlodenbdsig);
1108 double dfdJnb=0, dfdInb=0, dfdlodenb=0, dgdJnb=0, dgdInb=0, dgdlodenb=0;
1109 if((m<=0)||(m_flow<=0)||(rolodeb<=0)||(rolodefl<=0)) {
1110 cout<<"error in AI3-SKH, location 8, element "<<vari.element<<endl;
1111 //exit_tn( -YES );
1112 }
1113 dfdInb=dgdInb=Inorm_b/9;
1114 dfdJnb=3*Jnorm_b/(m*m);
1115 dgdJnb=3*Jnorm_b/(m_flow*m_flow);
1116
1117 double dfdrolodenb=-3*Jnorm_b*Jnorm_b/(m_cmp*m_cmp*scalar_power(rolodeb,3));
1118 double dgdrolodenb=-3*Jnorm_b*Jnorm_b/(m_flow_K0*m_flow_K0*scalar_power(rolodefl,3));
1119 dfdlodenb = dfdrolodenb * drolodenbdlodenb;
1120 dgdlodenb = dgdrolodenb * drolodefldlodenb;
1121
1122 for (int i=0; i<3; i++ ) {
1123 for (int j=0; j<3; j++ ) {
1124 df[3*i+j]= dfdInb*dInbdsig[i][j]+dfdJnb*dJnbdsig[i][j]+
1125 dfdlodenb*dlodenbdsig[i][j];
1126 df[3*i+j]*=-1; //Triax requires solid mechanics convention
1127 dg[3*i+j]= dgdInb*dInbdsig[i][j]+dgdJnb*dJnbdsig[i][j]+
1128 dgdlodenb*dlodenbdsig[i][j];
1129 dg[3*i+j]*=-1; //Triax requires solid mechanics convention
1130 }
1131 }
1132
1133 double dfdepp[9];
1134 array_set(dfdepp, 0, 9);
1135 double dfdI0=-T*T*S*S*I_0/9;
1136 double dI0devp=I_0/(lambda-kappa);
1137 double dsigbdepv[3][3];
1138 array_set(*dsigbdepv, 0, 9);
1139 for (int i=0; i<3; i++ ) {
1140 for (int j=0; j<3; j++ ) {
1141 dsigbdepv[i][j]=sig_b[3*i+j]/(lambda-kappa);
1142 }
1143 }
1144 double tmpA[3][3];
1145 array_set(*tmpA, 0, 9);
1146 for (int i=0; i<3; i++ ) {
1147 for (int j=0; j<3; j++ ) {
1148 tmpA[i][j]=-dfdInb*dInbdsig[i][j]-dfdJnb*dJnbdsig[i][j]-dfdlodenb*dlodenbdsig[i][j];
1149 }
1150 }
1151 double tmpB=array_inproduct(*tmpA, *dsigbdepv, 9)+dfdI0*dI0devp;
1152 for (int i=0; i<3; i++ ) {
1153 for (int j=0; j<3; j++ ) {
1154 dfdepp[3*i+j]=tmpB*vari.kron_delta[i][j];
1155 dfdepp[3*i+j] *= -1; //Triax requires solid mechanics convention
1156 }
1157 }
1158
1159 double b1max=2*I_0*(1-T)/3;
1160 if(b1max==0) b1max=TINY;
1161 double beta[9];
1162 array_set(beta,0,9);
1163 for (int j=0; j<9; j++ ) {
1164 beta[j]=(vari.geo_sig[j]-sig_b[j])/(T*S)+sig_0[j]-
1165 (vari.geo_sig[j]-sig_b[j])/S - sig_a[j];
1166 }
1167 double b1=array_inproduct(beta, df, 9)*3/(S*T*I_0);
1168 if(vari.other_f[1]>0) b1=0;
1169 double H1=S*S*scalar_power(scalar_dabs(b1)/b1max, psi)*I_0*I_0*I_0/((lambda-kappa)*27);
1170 double b2max=2*T*I_0*(1-S)/3;
1171 if(b2max==0) b2max=TINY;
1172 double gama[9];
1173 array_set(gama,0,9);
1174 for (int j=0; j<9; j++ ) {
1175 gama[j]=(vari.geo_sig[j]-sig_b[j])/S+sig_a[j]-vari.geo_sig[j];
1176 }
1177 double b2=array_inproduct(gama, df, 9)*3/(T*S*I_0);
1178 if(vari.other_f[0]>0) b2=0;
1179 double H2=scalar_power(T*scalar_dabs(b2)/b2max, psi)*I_0*I_0*I_0/((lambda-kappa)*27);
1180
1181 double H0=0;
1182 for(int i=0; i<MDIM*MDIM; i++) H0-=dg[i]*dfdepp[i];
1183 double H = H0+H1+H2;
1184
1185 array_move(df, vari.df_trial, MDIM*MDIM);
1186 if (!only_f) make_matrix_Cep(df, dg, H, elpl_matrix);
1187 }
1188
calc_rolode(double lode,double beta,double & rolode,double & drolodedlode,bool alsoderivative)1189 void Tskh_general_plasti :: calc_rolode(double lode, double beta, double &rolode,
1190 double &drolodedlode, bool alsoderivative) {
1191 //eliptical lode dependence after william and warnkle
1192 lode*=-1; //w&w dependence defined for stress negative in compression
1193 double tmpA = 2*(1-beta*beta)*cos(PI/6-lode);
1194 double insqrt=(4*(1-beta*beta)*scalar_power(cos(PI/6-lode),2)+beta*(5*beta-4));
1195 if(insqrt<0) {
1196 cout<<"error in AI3-SKH, location 9, element "<<vari.element<<endl;
1197 //exit_tn( -YES );
1198 }
1199 double tmpB = (2*beta-1)*sqrt(insqrt);
1200 double tmpC = 4*(1-beta*beta)*scalar_power(cos(PI/6-lode),2)+(2*beta-1)*(2*beta-1);
1201 if(tmpC==0) {
1202 cout<<"error in AI3-SKH, location 10, element "<<vari.element<<endl;
1203 //exit_tn( -YES );
1204 }
1205
1206 rolode=(tmpA+tmpB)/tmpC;
1207 drolodedlode=0;
1208
1209 if(alsoderivative) {
1210 double insqrt=4*(1-beta*beta)*scalar_power(cos(PI/6-lode),2)+beta*(5*beta-4);
1211 if(insqrt<0) {
1212 cout<<"error in AI3-SKH, location 11, element "<<vari.element<<endl;
1213 //exit_tn( -YES );
1214 }
1215 double u = (2*beta-1)*sqrt(insqrt)+ 2*(1-beta*beta)*cos(PI/6-lode);
1216 double v = 4*(1-beta*beta)*scalar_power(cos(PI/6-lode),2)+(2*beta-1)*(2*beta-1);
1217 double udash_den=4*(1-beta*beta)*scalar_power(cos(PI/6-lode),2)+beta*(5*beta-4);
1218 if((udash_den<=0)||(v==0)) {
1219 cout<<"error in AI3-SKH, location 12, element "<<vari.element<<endl;
1220 //exit_tn( -YES );
1221 }
1222 double udash = 2*(1-beta*beta)*sin(PI/6-lode)+(2*(2*beta-1)*(1-beta*beta)*sin(PI/3-2*lode))/
1223 sqrt(udash_den);
1224 double vdash = 4*(1-beta*beta)*sin(PI/3-2*lode);
1225 drolodedlode = (udash*v-vdash*u)/(v*v);
1226 }
1227 drolodedlode*=-1; //bloody solid mechanics convention :-)
1228 }
1229