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