1 : C MATER2TG SOURCE AM 12/02/28 00:28:35 7289 2 : SUBROUTINE mater2tg(cmat,ntot,lerror,idimpara3, 3 : * parahot3,npttot3,H66t,sigma6v,deps6v,lcp, 4 : * NOEL,NPT,Kinc) 5 : 6 : c cmat is the name (characters) of the material 7 : c ntot is the integration point (global number for the whole structure): counter 8 : c idimpara3 is the number of parameters for parahot3 9 : c parahot3 is the array that contains the material parameters and state variables 10 : c npttot3 is the number of integration points (number for the considered material 3D) 11 : c and is thus the number of columns in parahot3 12 : c H66t is the tangent matrix 13 : c sigma6v is the nominal stress vector 14 : c deps6v is the incremental strains vector 15 : 16 : c=========================================================================== 17 : 18 : c input 19 : c CMAT : name of the material 20 : c NTOT : # of the integration point ( global # for the entire structure ) 21 : c can be NTOT1 or NTOT2 or NTOT3 22 : 23 : c output 24 : c H66t : tangent stiffness matrix for 3D materials (THG add) 25 : 26 : c input and output 27 : c PARAHOT3 : the first columns contain the material properties at elevated temperature, 28 : c the last columns contain strains, stresses, .... 29 : 30 : c SIGMA6v : stress vector for 3D materials (THG add) 31 : c DEPS6V : incremental strain vector for 3D materials (THG add) 32 : 33 : c======================================================================= 34 : 35 : 36 : c This routine is run for every iteration, 37 : c for every point of integration. 38 : c 39 : C At the integration point NTOT is the material described by CMAT. 40 : c This material has, at the point NTOT, the caracteristics PARAHOT 41 : c at elevated temperature. 42 : 43 : c In multiaxial situation 44 : c We are going to calculate the stress vector sigm6v and the stiffness 45 : c matrix H66t as a function of the strain incr. vector deps6v, 46 : c calculated from the strain vector at the beginning of the time step. 47 : 48 : c i o 49 : c a : df/dsigma 50 : c df/dsigmax, df/dsigmay, df/dtauxy, df/dsigmaz 51 : c * Ai : Hardening parameter 52 : 53 : c THG add concrete model 3D 54 : 55 : IMPLICIT REAL*8 (A-B,D-H,O-Z) 56 : implicit integer (I-K,M,N) 57 : implicit logical (L) 58 : implicit character*10 (C) 59 : 60 : dimension parahot3(idimpara3,npttot3) 61 : dimension epstr6v(6), esigmac(6) 62 : dimension esigma6v(6),H66t(6,6),deps6v(6),sigma6v(6) 63 : dimension rt1(1,6),rc1(1,6),H66tsave(6,6),H66(6,6) 64 : 65 : i1 = 1 66 : i2 = 2 67 : i3 = 3 68 : i4 = 4 69 : i5 = 5 70 : i6 = 6 71 : i7 = 7 72 : i8 = 8 73 : i9 = 9 74 : i10=10 75 : i11=11 76 : i12=12 77 : i13=13 78 : i14=14 79 : i15=15 80 : i16=16 81 : i17=17 82 : i18=18 83 : i19=19 84 : i20=20 85 : i21=21 86 : i22=22 87 : i23=23 88 : i24=24 89 : i25=25 90 : i26=26 91 : i27=27 92 : i28=28 93 : i29=29 94 : i30=30 95 : i31=31 96 : i32=32 97 : i33=33 98 : i34=34 99 : i35=35 100 : i36=36 101 : i37=37 102 : i38=38 103 : i39=39 104 : 105 : r2 = 2. 106 : r3 = 3. 107 : 108 : c TEMPERATURE 109 : T = parahot3(idimpara3-39,ntot) 110 : TMAX = parahot3(idimpara3-38,ntot) 111 : 112 : 113 : c Compressive strength 114 : fc20 = parahot3(i3,ntot) 115 : rloc = rKfcsi(t,tmax) 116 : parahot3(i3,ntot) = rloc * fc20 117 : fc = parahot3(i3,ntot) 118 : 119 : c Tensile strength 120 : ft20 = parahot3(i4,ntot) 121 : rloc = rkftco(t,tmax) 122 : parahot3(i4,ntot) = rloc * ft20 123 : if (parahot3(i4,ntot).lt.(1.d5)) then 124 : c minimal value for ft = 0.10 MPa implemented for numerical stability 125 : c so if ft entered by the user < 0.10 MPa, the concrete is considered as 126 : c previously cracked with dt>0.96 and ft = 1.00 MPa 127 : parahot3(i4,ntot) = (1.d6) 128 : parahot3(idimpara3-41,ntot) = 129 : . max(0.96d0,parahot3(idimpara3-41,ntot)) 130 : endif 131 : ft = parahot3(i4,ntot) 132 : 133 : c Strain at peak stress 134 : c Strain at peak stress according to EC2: epsc1,ec2 135 : rloc = 1.0d0 136 : epsc1ec2 = EPSc1co(tmax,rloc) 137 : c Strain at peak stress according to ENV (minimal value): epsc1,env 138 : epsc1env = EPSc1co(tmax,-rloc) 139 : c Strain at peak stress for the instantaneous stress-strain relationship 140 : c epsc1,ETC = (2 * epsc1,min + epsc1,EC2)/3 141 : epsc1etc = (2.*epsc1env+epsc1ec2)/3. 142 : parahot3(i5,ntot) = epsc1etc 143 : 144 : epsc1ec220 = EPSc1co(20.d0,rloc) 145 : epsc1env20 = EPSc1co(20.d0,-rloc) 146 : epsc1etc20 = (2.*epsc1env20+epsc1ec220)/3. 147 : 148 : c Modulus at the origin 149 : parahot3(i1,ntot) = 2.0d0*parahot3(i3,ntot)/parahot3(i5,ntot) 150 : 151 : c Compressive limit of elasticity 152 : parahot3(i6,ntot) = (0.3d0)*parahot3(i3,ntot) 153 : fc0 = parahot3(i6,ntot) 154 : 155 : c Biaxial compression strength: considered as 1.20 x compressive strength 156 : parahot3(7,ntot) = (1.2d0)*parahot3(3,ntot) 157 : 158 : c thermal strain 159 : c first, save thermal strain at the precedent step 160 : parahot3(idimpara3-31,ntot) = parahot3(idimpara3,ntot) 161 : c then, calculate new thermal strain 162 : parahot3(idimpara3,ntot) = epsthsi(t,tmax) 163 : 164 : c Modulus used in tension 165 : parahot3(i12,ntot) = parahot3(i1,ntot) 166 : 167 : c Material parameter dc: damage in compression at peak stress 168 : dco = parahot3(i9,ntot) 169 : 170 : c Material parameter x,c: adimensional densitive crack energy in compression 171 : c This parameter is assumed to be constant with temperature 172 : rxc = parahot3(i10,ntot) 173 : 174 : c Model parameter kappa,c1: accumulated compressive plastic strain at peak stress 175 : denom = (2.-2.*dco)*((parahot3(i8,ntot))-1.) 176 : parahot3(i13,ntot) = -parahot3(i5,ntot)*(1.-2.*dco)/(denom) 177 : rkappa1 = parahot3(i13,ntot) 178 : 179 : c Model Parameter ac 180 : parahot3(i14,ntot) = ( -dlog(1-dco)/(rkappa1) ) 181 : 182 : c Model Parameter bc 183 : parahot3(i15,ntot) = 2*fc*rxc/((1-rxc)*((fc0*rkappa1)+ 184 : . ((fc-fc0)*rkappa1*dlog(r2)))) 185 : 186 : c Parameter at 187 : if (ft20.ge.1.d6) then 188 : parahot3(i16,ntot) = (7.0d0*ft20)/(12.0d0* 189 : . parahot3(i11,ntot)) 190 : else 191 : c for numerical reason at cannot be equal to 0 even if ft0 = 0 192 : parahot3(i16,ntot) = (7.0d0*1.0d6)/(12.0d0* 193 : . parahot3(i11,ntot)) 194 : endif 195 : 196 : c TRANSIENT CREEP STRAIN 197 : do iloc=1,6 198 : epstr6v(iloc) = 0.d0 199 : enddo 200 : c Phi function 201 : c First, phi(ti) is saved as phiprev 202 : phiprev = parahot3(17,ntot) 203 : c Then, phi(ti+1) is calculated and saved as phi 204 : c phi = 2*(epsc1,EC2 - epsc1,min) / (3*(fc/fck)) 205 : phi = 2.*(epsc1ec2-epsc1env)/(3.*(fc/fc20)) 206 : parahot3(17,ntot) = phi 207 : c Calculation of transient creep strain 208 : c If the temperature is bellow Tmax (T is decreasing or it is 209 : c not a first heating) the transient creep strain is unchanged 210 : if (T.lt.(tmax-1.0d-3)) then 211 : continue 212 : c Else, it is a first-time heating 213 : c so there is an increment of transient creep strain under compressive stress 214 : else 215 : c unless the concrete is in the descending branch of the material law 216 : c if (kappa,c < kappa,c1) (i.e. material in the ascending branch of the 217 : c compressive part of the material law) 218 : rkappac = parahot3(idimpara3-36,ntot) 219 : if (rkappac.le.rkappa1) then 220 : c epstr,i+1 = epstr,i + [phi(Ti+1) - phi(Ti)] * (1 - dc,i) * (sigma,eff-,i / fck) 221 : dc = parahot3(idimpara3-40,ntot) 222 : do iloc=i1,i6 223 : esigmac(iloc) = parahot3(29+iloc,ntot) 224 : epstr6v(iloc) = parahot3(17+iloc,ntot) + ((phi-phiprev)/ 225 : . fc20)*(1.0d0 - dc) * esigmac(iloc) 226 : end do 227 : c epstr cannot decrease in absolute value 228 : do iloc=i1,i6 229 : if (epstr6v(iloc).lt.parahot3(17+iloc,ntot)) then 230 : parahot3(17+iloc,ntot) = epstr6v(iloc) 231 : endif 232 : end do 233 : endif 234 : endif 235 : 236 : c CALCULATE THE INCREMENT IN MECHANICAL STRAIN DEPS 237 : c (the increment in thermal strain is subtracted) 238 : c Deps = Deps - (Epsth(n)-Epsth(n-1)) 239 : deps6v(1) = deps6v(1) - parahot3(idimpara3,ntot) 240 : . + parahot3(idimpara3-31,ntot) 241 : deps6v(2) = deps6v(2) - parahot3(idimpara3,ntot) 242 : . + parahot3(idimpara3-31,ntot) 243 : if (.not.lcp) deps6v(3) = deps6v(3) - parahot3(idimpara3,ntot) 244 : . + parahot3(idimpara3-31,ntot) 245 : 246 : lg = .false. 247 : lbroyden = .false. 248 : lfulldamage = .false. 249 : 250 : n = 1 251 : 23 continue 252 : do istep = 1,n 253 : c To improve convergence, the step can be subdivided 254 : 255 : 256 : c *********************************************************** 257 : c ********************* 258 : c ! 3D concrete model ! 259 : c ! PLASTIC BOX ! 260 : c ********************* 261 : 261 : call DPRAN3D(parahot3(i1,ntot),idimpara3,H66t,rt1,rc1, 262 : c ******* 263 : . esigma6v,deps6v,lg,ntot,lerror,i1,i2,i3,i4,i5, 264 : . i6,n,lnoconv,istep,lcp,lbroyden,H66,lfulldamage) 265 : c Plastic Box of the 3D plastic-damage concrete model by T. Gernay 266 : c Yield functions of Drucker Prager in compression and Rankine in tension 267 : c Calculate the effective stresses and the effective elastoplastic tangent 268 : c modulus Dt (implicit process) 269 : c *********************************************************** 270 : 271 : if (lnoconv) then 272 : n = n*2 273 : c the step is subdivided by factor 2 274 : c we come back to the state at the previous (converged) step: rkappac = rkappaic 275 : parahot3(idimpara3-35,ntot) = parahot3(idimpara3-37,ntot) 276 : parahot3(idimpara3-34,ntot) = parahot3(idimpara3-36,ntot) 277 : if (n.lt.32) go to 23 278 : c the step can be divided in maximum 16 "sub-steps" 279 : if (.not.lbroyden) then 280 : c if still no convergence when dividing the size of the step => try another algorithm 281 : lbroyden = .true. 282 : n = 1 283 : go to 23 284 : endif 285 : if (.not.lfulldamage) then 286 : c still no convergence => if completely damaged, "remove" the point (sigma = 0) 287 : dc = parahot3(idimpara3-40,ntot) 288 : dt = parahot3(idimpara3-41,ntot) 289 : if ((dt.gt.0.97).or.(dc.gt.0.97)) then 290 : lfulldamage = .true. 291 : n = 1 292 : go to 23 293 : endif 294 : endif 295 : c if we arrive here, all the strategies to improve convergence have been tested with 296 : c no success => lerror is true (send back an error from the MATERIAL subroutines) 297 : lerror = .true. 298 : return 299 : endif 300 : enddo 301 : 302 : c if ((NOEL.eq.18).and.(NPT.eq.23)) then 303 : * print *,'SORTIE DPRAN3D contraintes = ',esigma6v 304 : c endif 305 : 306 : c ************************************************************* 307 : c ********************* 308 : c ! 3D concrete model ! 309 : c ! DAMAGE BOX ! 310 : c ********************* 310 : call DAMTG3D(parahot3(i1,ntot),idimpara3,H66t,rt1,rc1, 311 : c ******* 312 : . esigma6v,lg,ntot,lerror,i1,i2,i3,i4,i5,i6,H66, 313 : . lfulldamage,lcp) 314 : c Damage Box of the 3D plastic-damage concrete model by T. Gernay 315 : c Calculate the nominal stresses and the nominal consistent algorithmic 316 : c tangent modulus Ct (explicit process) 317 : c ************************************************************* 318 : 319 : c if ((NOEL.eq.18).and.(NPT.eq.23)) then 320 : * print *,'SORTIE DAMTG3D contraintes = ', sigma6v 321 : c endif 322 : 323 : 324 : c The tangent matrix is symmetrized 325 : c because SAFIR uses a symmetric matrix for the resolution of the equilibrium 326 : c T 327 : c B = ( A + A ) / 2 328 : do iloc=2,6 329 : rloc = H66t(1,iloc) + H66t(iloc,1) 330 : H66t(1,iloc) = rloc/(2.d0) 331 : H66t(iloc,1) = H66t(1,iloc) 332 : enddo 333 : do iloc=3,6 334 : rloc = H66t(2,iloc) + H66t(iloc,2) 335 : H66t(2,iloc) = rloc/(2.d0) 336 : H66t(iloc,2) = H66t(2,iloc) 337 : enddo 338 : do iloc=4,6 339 : rloc = H66t(3,iloc) + H66t(iloc,3) 340 : H66t(3,iloc) = rloc/(2.d0) 341 : H66t(iloc,3) = H66t(3,iloc) 342 : enddo 343 : do iloc=5,6 344 : rloc = H66t(4,iloc) + H66t(iloc,4) 345 : H66t(4,iloc) = rloc/(2.d0) 346 : H66t(iloc,4) = H66t(4,iloc) 347 : enddo 348 : rloc = H66t(5,6) + H66t(6,5) 349 : H66t(5,6) = rloc/(2.d0) 350 : H66t(6,5) = H66t(5,6) 351 : 352 : c Save the nominal stresses 353 : do iloc=i1,i6 354 : sigma6v(iloc) = PARAHOT3(idimpara3-13+iloc,ntot) 355 : end do 356 : 357 : c if (lcp) then 358 : cc Plane stress model 359 : c parahot2(idimpara2-i3,ntot) = sigma6v(1) 360 : c parahot2(idimpara2-i2,ntot) = sigma6v(2) 361 : c parahot2(idimpara2-i1,ntot) = sigma6v(4) 362 : c PARAHOT2(idimpara2-i7,ntot)=parahot3(idimpara3-18,ntot) 363 : c PARAHOT2(idimpara2-i6,ntot)=parahot3(idimpara3-17,ntot) 364 : c PARAHOT2(idimpara2-i5,ntot)=parahot3(idimpara3-15,ntot) 365 : c endif 366 : END 367 : 368 :
© Cast3M 2003 - All rights reserved.
Disclaimer