00001 #include "glasgowdam.h"
00002 #include "matrix.h"
00003 #include "vector.h"
00004 #include "elastisomat.h"
00005 #include "global.h"
00006 #include "intpoints.h"
00007 #include "tensor.h"
00008 #include "vecttens.h"
00009 #include <math.h>
00010
00011 #define nijac 20
00012 #define limit 1.0e-8
00013
00014
00015
00016
00017
00018
00019 glasgowdam::glasgowdam (void)
00020 {
00021 st=0.0;
00022 k=0.0;
00023 gf0=0.0;
00024 lc=0.0;
00025 reftemp = 0.0;
00026 ft = paramf_type(0);
00027 }
00028
00029
00030
00031
00032
00033
00034 glasgowdam::~glasgowdam (void)
00035 {
00036
00037 }
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047 void glasgowdam::read (XFILE *in)
00048 {
00049 xfscanf (in,"%lf %lf %lf %m", &st, &gf0, &lc, ¶mf_type_kwdset, (int *)(&ft));
00050 xfscanf(in, "%lf %lf", &k, &reftemp);
00051 }
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064 void glasgowdam::damfuncpar(long ipp, vector &eps, vector &kappa)
00065 {
00066 vector poseps;
00067 vector princeps;
00068 matrix pvect;
00069 matrix epst;
00070 matrix dev;
00071 double i1e, j2e, nu, e, tmp, a, b;
00072 long ncomp, idem, im;
00073
00074 switch (ft)
00075 {
00076 case vonmises :
00077
00078 ncomp=Mm->ip[ipp].ncompstr;
00079 allocm(3, 3, epst);
00080 allocm(3, 3, dev);
00081 vector_tensor(eps, epst, Mm->ip[ipp].ssst, strain);
00082 i1e = first_invar(epst);
00083 deviator (epst,dev);
00084 j2e = second_invar (dev);
00085 idem = Mm->ip[ipp].gemid();
00086 if (Mm->ip[ipp].tm[idem] != elisomat)
00087 {
00088 fprintf(stderr, "\n\n Invalid type of elastic material is required");
00089 fprintf(stderr, "\n in function glasgowdam::damfuncpar, (file %s, line %d)\n", __FILE__, __LINE__);
00090 }
00091 im = Mm->ip[ipp].idm[idem];
00092 e = Mm->eliso[im].e;
00093 nu = Mm->eliso[im].nu;
00094 a = (k-1.0)/(2.0*k)/(1.0-2.0*nu);
00095 b = 3.0/k/(1+nu)/(1.0+nu);
00096 tmp = a*a*i1e*i1e + b*j2e;
00097 kappa[0] = a*i1e+sqrt(tmp);
00098
00099
00100
00101
00102
00103 break;
00104 case normazar :
00105 allocm(3, 3, epst);
00106 allocm(3, 3, pvect);
00107 allocv(3, princeps);
00108 allocv(3,poseps);
00109 fillv(0.0, princeps);
00110 vector_tensor(eps, epst, Mm->ip[ipp].ssst, strain);
00111 princ_val (epst,princeps,pvect,nijac,limit,Mp->zero,3,1);
00112 extractposv (princeps, poseps);
00113 scprd(poseps, poseps, tmp);
00114 kappa[0] = sqrt(tmp);
00115 break;
00116 default :
00117 {
00118 fprintf(stderr, "\n\n Unknown damage parameter function type");
00119 fprintf(stderr, "\n in function glasgowdam::damfuncpar, (file %s, line %d)\n", __FILE__, __LINE__);
00120 }
00121 }
00122 return;
00123 }
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138 double glasgowdam::damfunction(long ipp,double tempr,double chi,vector &kappa)
00139 {
00140 double omega = 0.0, e0, kappa0, tn, gamma, gf,th;
00141 long idem = Mm->ip[ipp].gemid();
00142
00143 if (Mm->ip[ipp].tm[idem] != elisomat)
00144 {
00145 fprintf(stderr, "\n\nError - ip %ld has wrong type of elastic material combined with", ipp);
00146 fprintf(stderr, "\n scalar damge. The elastic isotropic material is requried");
00147 fprintf(stderr, "\n Function glasgowdam::damfunction (file %s, line %d)\n", __FILE__, __LINE__);
00148 return 0.0;
00149 }
00150
00151
00152
00153 e0 = Mm->eliso[Mm->ip[ipp].idm[idem]].e;
00154
00155
00156 tn = (tempr-20.0)/100.0;
00157
00158 th = (kappa[1]-20.0)/100.0;
00159 if (tn>th){
00160 th=tn;
00161 }
00162
00163 if ((th < 0.0) || (th > 7.9)){
00164 fprintf(stderr, "\n\nWarning - largest value of normalised temperature\n");
00165 fprintf(stderr, " is out of suggested range 0.0-7.9 (ipp %ld)\n",ipp);
00166 }
00167
00168
00169
00170 kappa0 = st/e0*(1.0 - 0.016*th*th);
00171
00172
00173 gf = gf0*(1.0+0.39*th-0.07*th*th);
00174
00175
00176 gamma = (1.0-chi)*st*lc/gf;
00177
00178 if (kappa[0] < kappa0)
00179
00180 return 0.0;
00181
00182
00183 omega = 1.0-(kappa0/kappa[0])*exp(-gamma*(kappa[0]-kappa0));
00184
00185 return omega;
00186 }
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196 void glasgowdam::compute_thermdilat (double t,double dt,matrix &eps)
00197 {
00198 double auxt,alpha;
00199
00200 auxt = (t-20.0)/100.0;
00201 if (auxt>=0.0 && auxt<=6.0){
00202 alpha = 6.0e-5/(7.0-auxt);
00203 }
00204 else{ alpha=0.0; }
00205
00206 fillm (0.0,eps);
00207
00208 eps[0][0]=alpha*dt;
00209 eps[1][1]=alpha*dt;
00210 eps[2][2]=alpha*dt;
00211
00212 }
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225 double glasgowdam::thermdamfunction (long ipp,double tempr,vector &kappa)
00226 {
00227 double chi = 0.0;
00228 double tkappa0 = 20.0;
00229
00230 if (tempr > kappa[0])
00231 kappa[0] = tempr;
00232 if (kappa[0] > tkappa0)
00233 chi = 2.0e-3*(kappa[0] - tkappa0)*(1.0-5.0e-4*(kappa[0]-tkappa0));
00234 else
00235 chi = 0.0;
00236
00237 return chi;
00238 }
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249 void glasgowdam::matstiff (matrix &d,long ipp,long ido)
00250 {
00251 double omega, chi, dp;
00252 vector kappa(1);
00253
00254 switch (Mp->nlman->stmat)
00255 {
00256 case initial_stiff :
00257 Mm->elmatstiff (d,ipp);
00258 break;
00259 case tangent_stiff :
00260 Mm->elmatstiff (d,ipp);
00261 omega =Mm->ip[ipp].eqother[ido+1];
00262 chi = 0.0;
00263 kappa[0] = 20.0;
00264 chi = thermdamfunction(ipp, reftemp, kappa);
00265 dp = (1.0 - chi)*(1.0-omega);
00266 if (dp < 0.000001)
00267 dp = 0.000001;
00268 cmulm (dp,d);
00269 break;
00270 default :
00271 fprintf(stderr, "\n\nError - unknown type of stifness matrix");
00272 fprintf(stderr, "\n in function glasgowdam::matstiff (%s, line %d)\n", __FILE__, __LINE__);
00273 }
00274 }
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287 void glasgowdam::nlstresses (long ipp, long im, long ido)
00288 {
00289 long i,ncomp=Mm->ip[ipp].ncompstr;
00290 vector epsn(ncomp),sigma(ncomp), kappa(2), updkappa(2), epsft(ncomp), eps(ncomp), tkappa(1);
00291 matrix d(ncomp, ncomp), epst(3,3);
00292 double omega, chi, dp;
00293
00294
00295 for (i=0;i<ncomp;i++){
00296 epsn[i] = Mm->ip[ipp].strain[i];
00297 }
00298 kappa[0] = Mm->ip[ipp].eqother[ido+0];
00299 updkappa[1] = kappa[1] = reftemp;
00300
00301
00302 if (Mp->matmodel == local)
00303
00304 {
00305 damfuncpar(ipp, epsn, updkappa);
00306 if (updkappa[0] > kappa[0])
00307 {
00308 kappa[0] = updkappa[0];
00309 omega = damfunction(ipp, reftemp, 0.0, kappa);
00310 }
00311 else
00312 omega = Mm->ip[ipp].eqother[ido+1];
00313
00314 copyv(epsn, eps);
00315
00316
00317
00318
00319
00320
00321
00322 tkappa[0] = 20.0;
00323
00324 chi = thermdamfunction(ipp, reftemp, tkappa);
00325
00326
00327 dp = (1.0 - omega)*(1.0-chi);
00328 Mm->elmatstiff (d,ipp);
00329 mxv(d, eps, sigma);
00330 cmulv(dp, sigma);
00331 }
00332
00333
00334 for (i=0;i<ncomp;i++){
00335 Mm->ip[ipp].stress[i]=sigma[i];
00336 }
00337 Mm->ip[ipp].other[ido+0]=kappa[0];
00338 Mm->ip[ipp].other[ido+1]=omega;
00339 }
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351 void glasgowdam::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 }