00001 #include <math.h>
00002 #include <stdio.h>
00003 #include <stdlib.h>
00004 #include "mohrcparab.h"
00005 #include "global.h"
00006 #include "matrix.h"
00007 #include "vector.h"
00008 #include "tensor.h"
00009 #include "intpoints.h"
00010
00011
00012 #define nijac 20
00013 #define limit 1.0e-8
00014
00015
00016
00017
00018
00019
00020 mohrcoulombpar::mohrcoulombpar (void)
00021 {
00022 phi=0.0; c=0.0; psi=0.0;
00023 alphaphi = 0.0; betaphi = 0.0; sigcphi = 0.0;
00024 alphapsi = 0.0; betapsi = 0.0; sigcpsi = 0.0;
00025 }
00026
00027
00028
00029
00030
00031
00032 mohrcoulombpar::~mohrcoulombpar (void)
00033 {
00034
00035 }
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050 void mohrcoulombpar::read (XFILE *in)
00051 {
00052 xfscanf (in,"%lf %lf %lf",&phi,&c,&psi);
00053 sra.read (in);
00054
00055 alphaphi = sin(phi)/sqrt(9+3*sin(phi)*sin(phi));
00056 betaphi = alphaphi / sqrt(3*(3*c*c - alphaphi*alphaphi));
00057 sigcphi = sqrt(3*(c*c - alphaphi*alphaphi/3.0));
00058
00059 alphapsi = sin(psi)/sqrt(9+3*sin(psi)*sin(psi));
00060 betapsi = alphapsi / sqrt(3*(3*c*c - alphapsi*alphapsi));
00061 sigcpsi = sqrt(3*(c*c - alphapsi*alphapsi/3.0));
00062 }
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075 double mohrcoulombpar::yieldfunction (matrix &sig)
00076 {
00077 double f, tmp;
00078 double j2, j1;
00079 matrix dev(3,3);
00080
00081 deviator(sig, dev);
00082 j1 = first_invar(dev);
00083 j2 = second_invar(dev);
00084 tmp = 3 * j2 + sqrt(3.0) * betaphi * sigcphi * j1;
00085 f = sqrt(tmp) - sigcphi;
00086
00087 return f;
00088 }
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103 void mohrcoulombpar::deryieldfsigma (matrix &sig,matrix &dfds)
00104 {
00105 double j2, j1, k1, k2;
00106 matrix dev(3,3);
00107
00108 deviator(sig, dev);
00109 j1 = first_invar(dev);
00110 j2 = second_invar(dev);
00111 k1 = 1.0 / (sqrt(j2 + sqrt(3.0)*betaphi*sigcphi*j1));
00112 k2 = sqrt(3.0) * betaphi * sigcphi;
00113 dfds[0][0] = k1 * (3.0 * dev[0][0] + k2);
00114 dfds[0][1] = k1 * (3.0 * dev[0][1]);
00115 dfds[0][2] = k1 * (3.0 * dev[0][2]);
00116 dfds[1][0] = k1 * (3.0 * dev[1][0]);
00117 dfds[1][1] = k1 * (3.0 * dev[1][1] + k2);
00118 dfds[1][2] = k1 * (3.0 * dev[1][2]);
00119 dfds[2][0] = k1 * (3.0 * dev[2][0]);
00120 dfds[2][1] = k1 * (3.0 * dev[2][1]);
00121 dfds[2][2] = k1 * (3.0 * dev[2][2] + k2);
00122 return;
00123 }
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138 void mohrcoulombpar::derplaspotsigma (matrix &sig,matrix &dfds)
00139 {
00140 double j2, j1, k1, k2;
00141 matrix dev(3,3);
00142
00143 deviator(sig, dev);
00144 j1 = first_invar(dev);
00145 j2 = second_invar(dev);
00146 k1 = 1.0 / (sqrt(j2 + sqrt(3.0)*betapsi*sigcpsi*j1));
00147 k2 = sqrt(3.0) * betapsi * sigcpsi;
00148 dfds[0][0] = k1 * (3.0 * dev[0][0] + k2);
00149 dfds[0][1] = k1 * (3.0 * dev[0][1]);
00150 dfds[0][2] = k1 * (3.0 * dev[0][2]);
00151 dfds[1][0] = k1 * (3.0 * dev[1][0]);
00152 dfds[1][1] = k1 * (3.0 * dev[1][1] + k2);
00153 dfds[1][2] = k1 * (3.0 * dev[1][2]);
00154 dfds[2][0] = k1 * (3.0 * dev[2][0]);
00155 dfds[2][1] = k1 * (3.0 * dev[2][1]);
00156 dfds[2][2] = k1 * (3.0 * dev[2][2] + k2);
00157 return;
00158 }
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172 void mohrcoulombpar::matstiff (matrix &d, long ipp, long ido)
00173 {
00174 if (Mp->nlman->stmat==0)
00175 {
00176
00177 Mm->elmatstiff (d,ipp);
00178 }
00179 if (Mp->nlman->stmat==1)
00180 {
00181
00182
00183 Mm->elmatstiff (d,ipp);
00184 }
00185 }
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201 void mohrcoulombpar::nlstresses (long ipp, long im, long ido)
00202
00203 {
00204 long i,ni, nc=Mm->ip[ipp].ncompstr;
00205 double gamma,err;
00206 vector epsn(nc),epsp(nc),q(0);
00207
00208
00209 for (i=0; i<nc; i++)
00210 {
00211 epsn[i]=Mm->ip[ipp].strain[i];
00212 epsp[i]=Mm->ip[ipp].eqother[ido+i];
00213 }
00214 gamma=Mm->ip[ipp].eqother[ido+nc];
00215
00216
00217 if (sra.give_tsra () == cp){
00218 ni=sra.give_ni ();
00219 err=sra.give_err ();
00220
00221 Mm->cutting_plane (ipp,im,ido,gamma,epsn,epsp,q,ni,err);
00222 }
00223 else{
00224 print_err("wrong type of stress return algorithm is required", __FILE__, __LINE__, __func__);
00225 abort ();
00226 }
00227
00228
00229 for (i=0;i<nc;i++){
00230 Mm->ip[ipp].other[ido+i]=epsp[i];
00231 }
00232 Mm->ip[ipp].other[ido+nc]=gamma;
00233 }
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251 void mohrcoulombpar::nonloc_nlstresses (long ipp, long im, long ido)
00252
00253 {
00254 long i,ni, nc=Mm->ip[ipp].ncompstr;
00255 double gamma,err;
00256 vector epsn(nc),epsp(nc),q(0);
00257
00258
00259 for (i=0; i<nc; i++)
00260 {
00261 epsn[i]=Mm->ip[ipp].strain[i];
00262 epsp[i]=Mm->ip[ipp].nonloc[i];
00263 }
00264 gamma=Mm->ip[ipp].eqother[ido+nc];
00265
00266
00267 if (sra.give_tsra () == cp){
00268 ni=sra.give_ni ();
00269 err=sra.give_err ();
00270
00271 Mm->cutting_plane (ipp,im,ido,gamma,epsn,epsp,q,ni,err);
00272 }
00273 else{
00274 print_err("wrong type of stress return algorithm is required", __FILE__, __LINE__, __func__);
00275 abort ();
00276 }
00277
00278
00279 for (i=0;i<nc;i++){
00280 Mm->ip[ipp].other[ido+i]=epsp[i];
00281 }
00282 Mm->ip[ipp].other[ido+nc]=gamma;
00283 }
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299 void mohrcoulombpar::updateval (long ipp, long im, long ido)
00300 {
00301 long i,n = Mm->givencompeqother(ipp, im);
00302
00303 for (i=0;i<n;i++){
00304 Mm->ip[ipp].eqother[ido+i]=Mm->ip[ipp].other[ido+i];
00305 }
00306 }
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321 void mohrcoulombpar::giveirrstrains (long ipp, long ido, vector &epsp)
00322 {
00323 long i;
00324 for (i=0;i<epsp.n;i++)
00325 epsp[i] = Mm->ip[ipp].eqother[ido+i];
00326 }
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341 double mohrcoulombpar::give_consparam (long ipp, long ido)
00342 {
00343 long ncompstr;
00344 double gamma;
00345
00346 ncompstr=Mm->ip[ipp].ncompstr;
00347 gamma = Mm->ip[ipp].eqother[ido+ncompstr];
00348
00349 return gamma;
00350 }
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364 void mohrcoulombpar::changeparam (atsel &atm,vector &val)
00365 {
00366 long i;
00367
00368 for (i=0;i<atm.num;i++){
00369 switch (atm.atrib[i]){
00370 case 0:{
00371 phi=val[0];
00372 break;
00373 }
00374 case 1:{
00375 c=val[1];
00376 break;
00377 }
00378 case 2:{
00379 psi=val[2];
00380 break;
00381 }
00382 default:
00383 print_err("wrong number of atribute", __FILE__, __LINE__, __func__);
00384 }
00385 }
00386
00387 alphaphi = sin(phi)/sqrt(9+3*sin(phi)*sin(phi));
00388 betaphi = alphaphi / sqrt(3*(3*c*c - alphaphi*alphaphi));
00389 sigcphi = sqrt(3*(c*c - alphaphi*alphaphi/3.0));
00390
00391 alphapsi = sin(psi)/sqrt(9+3*sin(psi)*sin(psi));
00392 betapsi = alphapsi / sqrt(3*(3*c*c - alphapsi*alphapsi));
00393 sigcpsi = sqrt(3*(c*c - alphapsi*alphapsi/3.0));
00394 }