00001 #include "scaldamcc.h"
00002 #include "matrix.h"
00003 #include "vector.h"
00004 #include "tensor.h"
00005 #include "elastisomat.h"
00006 #include "global.h"
00007 #include "intpoints.h"
00008 #include "vecttens.h"
00009 #include <stdio.h>
00010 #include <stdlib.h>
00011 #include <math.h>
00012
00013
00014 #define nijac 20
00015 #define limit 1.0e-8
00016
00017
00018
00019
00020
00021 scaldamcc::scaldamcc (void)
00022 {
00023 k0 = at = bt = ac = bc = beta = k = 0.0;
00024 }
00025
00026
00027
00028
00029 scaldamcc::~scaldamcc (void)
00030 {
00031
00032 }
00033
00034
00035
00036
00037
00038
00039
00040
00041 void scaldamcc::read (XFILE *in)
00042 {
00043 xfscanf (in,"%lf %lf %lf %lf %lf %lf %m",&k0,&at, &bt, &ac, &bc, &beta, ¶mf_type_kwdset, (int *)(&ft));
00044 if (ft == vonmises)
00045 xfscanf(in, "%lf", &k);
00046 }
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056 void scaldamcc::elmatstiff (matrix &d,long ipp)
00057 {
00058 double e, e0;
00059 long idem = Mm->ip[ipp].gemid();
00060 long idoem=Mm->givencompeqother(ipp,0);
00061 long tmp=Mm->givencompeqother(ipp,idem);
00062 idoem -= tmp;
00063
00064 Mm->elmatstiff (d,ipp);
00065 e = Mm->give_actual_ym(ipp);
00066 e0 = Mm->give_actual_ym(ipp,idem,idoem);
00067 cmulm(e/e0, d);
00068 }
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082 void scaldamcc::damfuncpar(long ipp, vector &eps, vector &kappa)
00083 {
00084 vector poseps;
00085 vector princeps;
00086 matrix pvect;
00087 matrix epst;
00088 matrix dev;
00089 double i1e, j2e, nu, e, tmp;
00090 long ncomp, idem;
00091
00092 switch (ft)
00093 {
00094 case vonmises:
00095 ncomp=Mm->ip[ipp].ncompstr;
00096 allocm(3, 3, epst);
00097 allocm(3, 3, dev);
00098 vector_tensor(eps, epst, Mm->ip[ipp].ssst, strain);
00099 i1e = first_invar(epst);
00100 deviator (epst,dev);
00101 j2e = second_invar (dev);
00102 idem = Mm->ip[ipp].gemid();
00103 if (Mm->ip[ipp].tm[idem] != elisomat)
00104 {
00105 fprintf(stderr, "\n\n Invalid type of elastic material is required");
00106 fprintf(stderr, "\n in function scaldamcc::damfuncpar, (%s, line %d)\n", __FILE__, __LINE__);
00107 }
00108 e = Mm->eliso[idem].e;
00109 nu = Mm->eliso[idem].nu;
00110 kappa[0] = (k-1.0)/(2.0*k)/(1.0-2.0*nu)*i1e;
00111 tmp = (k-1.0)*(k-1.0)/(1.0-2.0*nu)/(1.0-2.0*nu)*i1e*i1e -12.0*k/(1.0+nu)/(1.0+nu)*j2e;
00112 kappa[0] += 1.0/(2.0*k)*sqrt(tmp);
00113 break;
00114 case normazar:
00115 allocm(3, 3, epst);
00116 allocm(3, 3, pvect);
00117 allocv(3, princeps);
00118 allocv(3,poseps);
00119 fillv(0.0, princeps);
00120 vector_tensor(eps, epst, Mm->ip[ipp].ssst, strain);
00121 princ_val (epst,princeps,pvect,nijac,limit,Mp->zero,3,1);
00122 extractposv (princeps, poseps);
00123 scprd(poseps, poseps, tmp);
00124 kappa[0] = sqrt(tmp);
00125 break;
00126 default:
00127 {
00128 fprintf(stderr, "\n\n Unknown damage parameter function type");
00129 fprintf(stderr, "\n in function scaldamcc::damfuncpar, (%s, line %d)\n", __FILE__, __LINE__);
00130 }
00131 }
00132 return;
00133 }
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143 double scaldamcc::damfunction(long ipp, vector &kappa, vector &eps)
00144 {
00145 long ncomp=Mm->ip[ipp].ncompstr;
00146 matrix epst, sigt, pvect, d(ncomp, ncomp),tmp3(ncomp, ncomp),tmp;
00147 vector princeps;
00148 vector poseps;
00149 vector tenseps;
00150 vector compeps;
00151 vector sigm(ncomp), princsig, tmp1, tmp2;
00152 vector negeps;
00153 double omega = 0.0, dt, dc, alphat, alphac;
00154
00155 if (kappa[0] < k0)
00156
00157 return 0.0;
00158
00159 allocm(3, 3, epst);
00160 allocm(3, 3, sigt);
00161 allocm(3, 3, tmp);
00162 allocv(3, poseps);
00163 allocv(3, negeps);
00164 allocv(3, tenseps);
00165 allocv(3, compeps);
00166 allocv(3, princeps);
00167 allocv(3, princsig);
00168 allocv(3, tmp1);
00169 allocv(3, tmp2);
00170 allocm(3, 3, pvect);
00171
00172 dt = 1.0 - k0 * (1.0 - at) / kappa[0] - at/exp(bt*(kappa[0]-k0));
00173 dc = 1.0 - k0 * (1.0 - ac) / kappa[0] - ac/exp(bc*(kappa[0]-k0));
00174 vector_tensor(eps, epst, Mm->ip[ipp].ssst, strain);
00175 princ_val (epst,princeps,pvect,nijac,limit,Mp->zero,3,1);
00176
00177
00178
00179 Mm->elmatstiff (d,ipp);
00180 mxv(d,eps,sigm);
00181 vector_tensor(sigm,sigt, Mm->ip[ipp].ssst, stress);
00182 princ_val(sigt,princsig,pvect,nijac,limit,Mp->zero,3,1);
00183 extractposv(princsig,tmp1);
00184 extractnegv(princsig,tmp2);
00185 invm(d,tmp3,Mp->zero);
00186
00187 extractm(tmp,tmp3,0,3);
00188 mxv(tmp,tmp1,tenseps);
00189 mxv(tmp,tmp2,compeps);
00190
00191 extractposv (princeps, poseps);
00192 subv(princeps, poseps, negeps);
00193 scprd(tenseps, poseps, alphat);
00194 alphat = alphat / kappa[0] / kappa[0];
00195 alphat = pow(alphat, beta);
00196 scprd(compeps, poseps, alphac);
00197 alphac = alphac / kappa[0] / kappa[0];
00198 alphac = pow(alphac, beta);
00199
00200 omega = alphat*dt + alphac*dc;
00201
00202 if ((omega > 1.0) || (omega < 0.0)) {
00203 fprintf(stderr, "\n\nError - iteration of omega doesn't convergate to range <0,1> %g ",omega);
00204 fprintf(stderr, "\n in function damfunction(), (%s, line %d)\n", __FILE__, __LINE__);
00205 }
00206 if(omega>1.0)
00207 omega = 1.0 ;
00208 if(omega<0.0)
00209 omega=0.0;
00210
00211 return omega;
00212 }
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222 void scaldamcc::matstiff (matrix &d,long ipp,long ido)
00223 {
00224 double dp;
00225
00226 switch (Mp->nlman->stmat)
00227 {
00228 case initial_stiff:
00229 Mm->elmatstiff (d,ipp);
00230 break;
00231 case tangent_stiff:
00232 Mm->elmatstiff (d,ipp);
00233 dp=Mm->ip[ipp].other[ido+1];
00234 cmulm (1.0-dp,d);
00235 break;
00236 default:
00237 fprintf(stderr, "\n\nError - unknown type of stifness matrix");
00238 fprintf(stderr, "\n in function scaldamcc::matstiff (%s, line %d)\n", __FILE__, __LINE__);
00239 }
00240 }
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252 void scaldamcc::nlstresses (long ipp,long im,long ido)
00253 {
00254 long i,ncomp=Mm->ip[ipp].ncompstr;
00255 double dp;
00256 vector epsn(ncomp),sigma(ncomp), kappa(1), omegao(1);
00257 matrix d(ncomp, ncomp);
00258
00259
00260 for (i=0;i<ncomp;i++){
00261 epsn[i] = Mm->ip[ipp].strain[i];
00262 }
00263 kappa[0] = Mm->ip[ipp].other[ido+0];
00264 omegao[0] = Mm->ip[ipp].other[ido+1];
00265
00266 elmatstiff(d,ipp);
00267
00268 dp = Mm->scal_dam_sol (ipp, im, ido, epsn, kappa, omegao, sigma,d);
00269
00270
00271 for (i=0;i<ncomp;i++){
00272 Mm->ip[ipp].stress[i]=sigma[i];
00273 }
00274 Mm->ip[ipp].other[ido+2]=kappa[0];
00275 Mm->ip[ipp].other[ido+3]=dp;
00276 }
00277
00278 void scaldamcc::updateval (long ipp,long im,long ido)
00279 {
00280 long i,n = Mm->givencompeqother(ipp, im);
00281
00282 for (i=0;i<n;i++){
00283 Mm->ip[ipp].eqother[ido+i]=Mm->ip[ipp].other[ido+i];
00284 }
00285
00286 }
00287
00288
00289
00290
00291
00292 double scaldamcc::epsefunction ()
00293 {
00294 return k0 ;
00295 }