00001 #include <math.h>
00002 #include <stdio.h>
00003 #include <stdlib.h>
00004 #include "doubdp.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 doubdp::doubdp (void)
00017 {
00018 fc = 0.0; ft = 0.0; fb = 0.0; gtf = 0.0;
00019 }
00020
00021
00022
00023
00024
00025 doubdp::~doubdp (void)
00026 {
00027
00028 }
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039 void doubdp::read (XFILE *in)
00040 {
00041 xfscanf (in,"%lf %lf %lf %lf",&fc,&ft,&fb,>f);
00042 sra.read (in);
00043
00044 alphat = 1.73205/3*(fc-ft)/(fc+ft);
00045 taut = fc*(1.73205-3*alphat)/3;
00046
00047 alphac = 1.73205/3*(fb-fc)/(2*fb+fc);
00048 tauc = fc*(1.73205-3*alphac)/3;
00049
00050 icr = (taut-tauc)/(alphat-alphac);
00051 jcr = taut - alphat*icr;
00052
00053 ivrch = taut/alphat;
00054 jvrch = 0.0;
00055 }
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065 void doubdp::params (long ipp, vector &q)
00066 {
00067 long eid;
00068 double h;
00069
00070 eid = Mm->elip[ipp];
00071 switch (Mt->give_dimension(eid))
00072 {
00073 case 1:
00074 h = Mt->give_length(eid);
00075 break;
00076 case 2:
00077 h = sqrt(Mt->give_area(eid));
00078 break;
00079 case 3:
00080 h = pow(Mt->give_volume(eid), 1.0/3.0);
00081 break;
00082 default:
00083 print_err("unknown dimension of element", __FILE__, __LINE__, __func__);
00084 }
00085
00086 gamult = gtf*2/ft/h;
00087 q(1) = ft*(1-q(0)/gamult);
00088
00089 alphat = 1.73205/3*(fc-q(1))/(fc+q(1));
00090 taut = fc*(1.73205-3*alphat)/3;
00091
00092 alphac = 1.73205/3*(fb-fc)/(2*fb+fc);
00093 tauc = fc*(1.73205-3*alphac)/3;
00094
00095 icr = (taut-tauc)/(alphat-alphac);
00096 jcr = taut - alphat*icr;
00097
00098 ivrch = taut/alphat;
00099 jvrch = 0.0;
00100 }
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112 double doubdp::yieldfunction (matrix &sig, vector &q)
00113 {
00114 double i1, j2, f1, f2, n1, n2, n3;
00115 double tcr=0.0, tvrch=0.0;
00116 double a, b, c;
00117 double l, m, n;
00118 matrix dev(3,3);
00119
00120 i1 = first_invar (sig);
00121 deviator (sig,dev);
00122 j2 = sqrt(second_invar (dev));
00123
00124 f1 = alphat*i1+j2-taut;
00125 f2 = alphac*i1+j2-tauc;
00126 n1 = alphac*j2-i1-alphac*jcr+icr;
00127 n2 = alphat*j2-i1-alphat*jcr+icr;
00128 n3 = alphat*j2-i1-alphat*jvrch+ivrch;
00129
00130 a = i1-icr;
00131 b = j2-jcr;
00132 c = -a*icr-b*jcr;
00133 l = i1-ivrch;
00134 m = j2-jvrch;
00135 n = -l*ivrch-m*jvrch;
00136
00137 if ((f1<=0.0)&&(f2<=0.0)){v = 1; return f1;}
00138
00139 if ((f1>0.0) &&(f2<=0.0))
00140 {
00141 if (n3<0)
00142 {
00143 v = 5;
00144 tvrch = l*i1 + m*j2 + n;
00145 return tvrch;
00146 }
00147 else {v = 2; return f1;}
00148 }
00149
00150 if ((f1<=0.0)&&(f2>0.0)) {v = 3; return f2;}
00151
00152 if ((f1>0.0) &&(f2>0.0))
00153 {
00154 if (n1>=0) {v = 3; return f2;}
00155 if ((n2<=0)&&(n3>=0)) {v = 2; return f1;}
00156 if ((n1<0) &&(n2>0))
00157 {
00158 v = 4;
00159 tcr = a*i1 + b*j2 + c;
00160 return tcr;
00161 }
00162 if (n3<0)
00163 {
00164 v = 5;
00165 tvrch = l*i1 + m*j2 + n;
00166 return tvrch;
00167 }
00168 }
00169 return -1;
00170 }
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186 void doubdp::deryieldfsigma (matrix &sig, matrix &dfds, vector &q)
00187 {
00188 double i1, j2;
00189 double a, b, d;
00190 double l, m;
00191 long i,j;
00192 matrix dev(3,3);
00193
00194 for (i=0;i<3;i++)
00195 for (j=0;j<3;j++)
00196 dfds[i][j]=0.0;
00197
00198 i1 = first_invar(sig);
00199 deviator(sig,dev);
00200 j2 = sqrt(second_invar(dev));
00201
00202 a = i1-icr;
00203 b = j2-jcr;
00204 l = i1-ivrch;
00205 m = j2-jvrch;
00206
00207 dfds[0][0] += (sig[0][0]-sig[1][1])/3;
00208 dfds[0][1] += 2*sig[0][1];
00209 dfds[1][0] += 2*sig[0][1];
00210 dfds[1][1] += (sig[1][1]-sig[2][2])/3;
00211 dfds[1][2] += 2*sig[1][2];
00212 dfds[2][1] += 2*sig[1][2];
00213 dfds[2][0] += 2*sig[2][0];
00214 dfds[0][2] += 2*sig[2][0];
00215 dfds[2][2] += (sig[2][2]-sig[0][0])/3;
00216
00217 d = 1/(2*j2);
00218 cmulm(d, dfds);
00219
00220
00221 if (v==2)
00222 {
00223 dfds(0,0) += alphat;
00224 dfds(1,1) += alphat;
00225 dfds(2,2) += alphat;
00226 }
00227 if (v==3)
00228 {
00229 dfds(0,0) += alphac;
00230 dfds(1,1) += alphac;
00231 dfds(2,2) += alphac;
00232 }
00233 if (v==4)
00234 {
00235 cmulm(b, dfds);
00236 dfds(0,0) += a;
00237 dfds(1,1) += a;
00238 dfds(2,2) += a;
00239 }
00240 if (v==5)
00241 {
00242 cmulm(m, dfds);
00243 dfds(0,0) += l;
00244 dfds(1,1) += l;
00245 dfds(2,2) += l;
00246 }
00247 }
00248
00249
00250
00251
00252
00253 void doubdp::matstiff (matrix &d, long ipp,long ido)
00254 {
00255 if (Mp->nlman->stmat==initial_stiff)
00256 {
00257
00258 Mm->elmatstiff (d,ipp);
00259 }
00260 if (Mp->nlman->stmat==tangent_stiff)
00261 {
00262
00263
00264 Mm->elmatstiff (d,ipp);
00265
00266 }
00267 }
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284 void doubdp::nlstresses (long ipp, long im, long ido)
00285 {
00286 long i,j=1,ni,n=Mm->ip[ipp].ncompstr;
00287 double gamma,err,nu;
00288 vector epsn(n),epsp(n),q(2),sig(n),pom(n);
00289 matrix d(n,n), sigt(3,3);
00290
00291
00292 for (i=0;i<n;i++)
00293 {
00294 epsn[i]=Mm->ip[ipp].strain[i];
00295 epsp[i]=Mm->ip[ipp].eqother[ido+i];
00296 }
00297 gamma=Mm->ip[ipp].eqother [ido+n];
00298 q[0] = Mm->ip[ipp].eqother[ido+n+1];
00299 ni=sra.give_ni ();
00300 err=sra.give_err ();
00301
00302
00303 if (Mm->ip[ipp].ssst == platek)
00304 Mm->ip[ipp].ssst = planestress;
00305
00306 if ((Mm->ip[ipp].ssst == planestress))
00307 {
00308 for (i=0;i<n;i++)
00309 pom[i] = epsn[i];
00310 epsn[3] = pom[2];
00311
00312
00313
00314
00315
00316
00317 nu = Mm->give_actual_nu(ipp);
00318 epsn[2] = -nu/(1-nu)*(epsn[0]-epsp[0]+epsn[1]-epsp[1])+epsp[2];
00319 for (j=0;j<ni;j++)
00320 {
00321 Mm->cutting_plane (ipp,im,ido,gamma,epsn,epsp,q,ni,err);
00322 if (abs(Mm->ip[ipp].stress[3])<err)
00323 break;
00324 Mm->ip[ipp].ssst = axisymm;
00325 Mm->elmatstiff (d,ipp);
00326 Mm->ip[ipp].ssst = planestress;
00327 epsn[2]+= -Mm->ip[ipp].stress[3]/d(2,2);
00328 }
00329 }
00330 else
00331 Mm->cutting_plane (ipp,im,ido,gamma,epsn,epsp,q,ni,err);
00332
00333
00334
00335
00336
00337
00338
00339 for (i=0;i<n;i++){
00340 Mm->ip[ipp].other[ido+i]=epsp[i];
00341 }
00342 Mm->ip[ipp].other[ido+n]=gamma;
00343 Mm->ip[ipp].other[ido+n+1]=q[0];
00344 }
00345
00346
00347
00348
00349
00350
00351 void doubdp::updateval (long ipp, long im, long ido)
00352 {
00353 long i,n = Mm->givencompeqother(ipp,im);
00354
00355 for (i=0;i<n;i++)
00356 Mm->ip[ipp].eqother[ido+i]=Mm->ip[ipp].other[ido+i];
00357 }
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368 double doubdp::plasmodscalar(matrix &sig, vector &q)
00369 {
00370 double ret;
00371 double dfdq;
00372 double dqdg;
00373 double i1 = first_invar(sig);
00374
00375 dfdq = -(2*sqrt(3)/3)*fc/((fc+q(1))*(fc+q(1)))*(i1+fc);
00376 dqdg = -ft/gamult;
00377 ret = dfdq * dqdg;
00378
00379 ret = 0;
00380
00381 return -ret;
00382 }
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396 void doubdp::updateq(long ipp, double dgamma, vector &q)
00397 {
00398
00399
00400 }