Télécharger damtg3d.eso

Retour à la liste

Numérotation des lignes :

  1. C DAMTG3D SOURCE CB215821 16/04/21 21:16:14 8920
  2. SUBROUTINE DAMTG3D(parahot3,idimpara3,H66t,rt1,rc1,
  3. . esigma6v,lg,ntot,lerror,i1,i2,i3,i4,i5,i6,H66,
  4. . lfulldamage,lcp)
  5.  
  6. c SUBROUTINE added by THG
  7.  
  8. c parahot3 is the array that contains the material parameters and state variables
  9. c idimpara3 is the number of parameters for parahot3
  10. c H66t is the tangent matrix
  11. c rt1 and rc1 are the vectors built in DPRAN3D and used in DAMTG3D for calculation of tangent matrix
  12. c esigma6v is the effective stress vector
  13. c lg and lerror are logical
  14. c ntot is the integration point (global number for the whole structure): counter
  15.  
  16. c ==========================================================================================
  17. c ! DAMAGE BOX of the 3D plastic-damage concrete model by T. Gernay !
  18. c ! This subroutine calculates !
  19. c ! the nominal stress at the end of the time step !
  20. c ! the nominal consistent algorithmic tangent modulus Ct !
  21. c ! !
  22. c ! The effective stress esigma6v and the eff. tangent modulus H66t given by plastic box !
  23. C ! !
  24. c ==========================================================================================
  25.  
  26. IMPLICIT REAL*8 (A-B,D-H,O-Z)
  27. implicit integer (I-K,M,N)
  28. implicit logical (L)
  29. implicit character*10 (C)
  30.  
  31. dimension H66t(6,6)
  32. c Tangent matrix (output)
  33. dimension parahot3(idimpara3)
  34. c parahot3 : Material prop. and various info. at elevated temp. (input and output)
  35. c parahot3(idimpara3-34): equivalent plastic strain in compression at the end of the current step
  36. c parahot3(idimpara3-35): equivalent plastic strain in tension at the end of the current step
  37. c parahot3(idimpara3-6:idimpara3-1) = effective stress obtained by dpran3D
  38. dimension esigma6v(6)
  39. c Effective stress (input)
  40.  
  41. dimension rc1(1,6),rt1(1,6)
  42. dimension v1(3),v2(3),v3(3),rcossigmapr(3,3)
  43. dimension p1(3,3),p2(3,3),p3(3,3)
  44. dimension pt1(3,3,3,3),pt2(3,3,3,3),pt3(3,3,3,3)
  45. dimension projt(3,3,3,3),w4(3,3,3,3)
  46. dimension Qtens(3,3,3,3),Qcomp(3,3,3,3)
  47. dimension esigmat(6),esigmac(6)
  48. dimension sigma(6),projtm(6,6),projcm(6,6)
  49. dimension p12(3,3),p23(3,3),p13(3,3)
  50. dimension pp12(3,3,3,3),pp23(3,3,3,3),pp13(3,3,3,3)
  51. dimension W4m(6,6),H66temp(6,6)
  52. dimension sigrt(6,6),sigrc(6,6)
  53. dimension p121(3,3),p231(3,3),p131(3,3)
  54. dimension p122(3,3),p232(3,3),p132(3,3)
  55. dimension damtensor(6,6), damtensor2(6,6),rt1mat(3,3)
  56. dimension rc1mat(3,3),esigmatmat(3,3),esigmacmat(3,3)
  57. dimension sigrttens(3,3,3,3),sigrctens(3,3,3,3),H66(6,6)
  58. dimension Qtensm(6,6),Qcompm(6,6),rc1t(6)
  59. dimension U(6,6)
  60. * dimension id4(3,3,3,3)
  61.  
  62. c Matrix of unity I, or U
  63. c -----------------------
  64.  
  65. do jloc=1,6
  66. do iloc=1,6
  67. if (iloc.eq.jloc) then
  68. U(iloc,jloc) = 1.0d0
  69. else
  70. U(iloc,jloc) = 0.0d0
  71. endif
  72. enddo
  73. enddo
  74.  
  75. c Fourth identity tensor id4
  76. C do iloc=1,3
  77. C do jloc=1,3
  78. C do kloc=1,3
  79. C do mloc=1,3
  80. C if ((iloc.eq.jloc).and.(kloc.eq.mloc)) then
  81. C id4(iloc,jloc,kloc,mloc) = 1.0d0
  82. C else
  83. C id4(iloc,jloc,kloc,mloc) = 0.0d0
  84. C endif
  85. C end do
  86. C end do
  87. C end do
  88. C end do
  89.  
  90. i1 = 1
  91. i2 = 2
  92. i3 = 3
  93. i4 = 4
  94. i5 = 5
  95. i6 = 6
  96. i7 = 7
  97. i8 = 8
  98. i9 = 9
  99. i10=10
  100. i11=11
  101. i12=12
  102. i13=13
  103. i14=14
  104. i15=15
  105. i16=16
  106. i17=17
  107. i18=18
  108. i19=19
  109. i20=20
  110. i21=21
  111. i22=22
  112. i23=23
  113. i24=24
  114. i25=25
  115. i26=26
  116. i27=27
  117. i28=28
  118. i29=29
  119. i30=30
  120. i31=31
  121. i32=32
  122. i33=33
  123. i34=34
  124. i35=35
  125. i36=36
  126. i37=37
  127. i38=38
  128. i39=39
  129.  
  130. r0 = 0.
  131. r1 = 1.
  132.  
  133. c ================================================
  134. c ! !
  135. c ! EIGEN VALUES OF THE EFFECTIVE STRESSES !
  136. c ! ====================================== !
  137. c ================================================
  138.  
  139. c Effective stresses
  140. do iloc=i1,i6
  141. ESIGMA6V(iloc) = parahot3(idimpara3-i7+iloc)
  142. end do
  143.  
  144.  
  145. c Eigen values and eigen vector
  146. call sigmapr3ETC(esigma6v,S1,S2,S3,rcossigmapr,lcp)
  147. do iloc=i1,i3
  148. v1(iloc) = rcossigmapr(iloc,i1)
  149. v2(iloc) = rcossigmapr(iloc,i2)
  150. v3(iloc) = rcossigmapr(iloc,i3)
  151. end do
  152.  
  153. c ===============================================
  154. c ! !
  155. c ! UPDATE DAMAGE VARIABLES !
  156. c ! ======================= !
  157. c ===============================================
  158.  
  159. rkappac = parahot3(idimpara3-34)
  160. rkappat = parahot3(idimpara3-35)
  161.  
  162. c Tensile damage
  163. at = parahot3(16)
  164. dt = 1.0d0 - ((0.5d0)*exp(-at*rkappat)+(0.5d0)*
  165. . exp(-6.0d0*at*rkappat))
  166. rdti = parahot3(idimpara3-41)
  167. if (dt.lt.rdti) then
  168. dt = rdti
  169. endif
  170. parahot3(idimpara3-33) = dt
  171.  
  172.  
  173. c Compressive damage
  174. rkappaic = parahot3(idimpara3-36)
  175. rdeltakc = rkappac - rkappaic
  176. rdci = parahot3(idimpara3-40)
  177. rdamprev = LOG(1.0d0-rdci)
  178. dc = 1.0d0 - exp((rdamprev)-parahot3(i14)*rdeltakc)
  179. if (dc.lt.parahot3(idimpara3-40)) dc = parahot3(idimpara3-40)
  180. parahot3(idimpara3-32) = dc
  181.  
  182. c ===================================================================
  183. c ! !
  184. c ! PROJECTION OPERATORS OF THE EFFECTIVE STRESSES P+ AND P- !
  185. c ! ======================================================== !
  186. c ===================================================================
  187.  
  188. c p1(i3,i3) = v1 X v1 (tensorial product)
  189. call mulABT(v1,v1,p1,i3,i1,i3,i1)
  190. c [3;3]=[3;1]x[1;3]
  191. call mulABT(v2,v2,p2,i3,i1,i3,i1)
  192. c [3;3]=[3;1]x[1;3]
  193. call mulABT(v3,v3,p3,i3,i1,i3,i1)
  194. c [3;3]=[3;1]x[1;3]
  195.  
  196. c P+ calculated as the sum = H(sig,i) p,i X p,i
  197. c First we calculate the PT,i = p,i X p,i
  198. do mloc=i1,i3
  199. do kloc=i1,i3
  200. do jloc=i1,i3
  201. do iloc=i1,i3
  202. pt1(iloc,jloc,kloc,mloc) = p1(iloc,jloc)*p1(kloc,mloc)
  203. c [3;3;3;3] = [3;3] x [3;3] tensorial product
  204. pt2(iloc,jloc,kloc,mloc) = p2(iloc,jloc)*p2(kloc,mloc)
  205. pt3(iloc,jloc,kloc,mloc) = p3(iloc,jloc)*p3(kloc,mloc)
  206. end do
  207. end do
  208. end do
  209. end do
  210.  
  211. if (S1.gt.r0) then
  212. h1 = r1
  213. else
  214. h1 = r0
  215. endif
  216. if (S2.gt.r0) then
  217. h2 = r1
  218. else
  219. h2 = r0
  220. endif
  221. if (S3.gt.r0) then
  222. h3 = r1
  223. else
  224. h3 = r0
  225. endif
  226.  
  227. c Projection operators of the effective stresses P+ [3;3;3;3]
  228. do mloc=i1,i3
  229. do kloc=i1,i3
  230. do jloc=i1,i3
  231. do iloc=i1,i3
  232. projt(iloc,jloc,kloc,mloc) = h1*pt1(iloc,jloc,kloc,mloc)+
  233. . h2*pt2(iloc,jloc,kloc,mloc) + h3*pt3(iloc,jloc,kloc,mloc)
  234. end do
  235. end do
  236. end do
  237. end do
  238.  
  239. c Transform the [3;3;3;3] tensor into [6;6] matrix (Voigt notation)
  240. c projtm is the tensile projection matrix [6;6] of the effective stress
  241. do jloc=i1,i3
  242. do iloc=i1,i3
  243. projtm(iloc,jloc) = projt(iloc,iloc,jloc,jloc)
  244. end do
  245. end do
  246. do iloc=i1,i3
  247. projtm(iloc,i4) = 0.5d0*(projt(iloc,iloc,i1,i2) +
  248. . projt(iloc,iloc,i2,i1))
  249. projtm(iloc,i5) = 0.5d0*(projt(iloc,iloc,i1,i3) +
  250. . projt(iloc,iloc,i3,i1))
  251. c 5th component is eps,xz
  252. projtm(iloc,i6) = 0.5d0*(projt(iloc,iloc,i2,i3) +
  253. . projt(iloc,iloc,i3,i2))
  254. c 6th component is eps,yz
  255. end do
  256. do jloc=i1,i3
  257. projtm(i4,jloc) = 0.5d0*(projt(i1,i2,jloc,jloc) +
  258. . projt(i2,i1,jloc,jloc))
  259. projtm(i5,jloc) = 0.5d0*(projt(i1,i3,jloc,jloc) +
  260. . projt(i3,i1,jloc,jloc))
  261. projtm(i6,jloc) = 0.5d0*(projt(i2,i3,jloc,jloc) +
  262. . projt(i3,i2,jloc,jloc))
  263.  
  264. end do
  265. projtm(i4,i4) = 0.5d0*(projt(i1,i2,i1,i2) +
  266. . projt(i1,i2,i2,i1))
  267. projtm(i4,i5) = 0.5d0*(projt(i1,i2,i1,i3) +
  268. . projt(i1,i2,i3,i1))
  269. projtm(i4,i6) = 0.5d0*(projt(i1,i2,i2,i3) +
  270. . projt(i1,i2,i3,i2))
  271. projtm(i5,i4) = 0.5d0*(projt(i1,i3,i1,i2) +
  272. . projt(i1,i3,i2,i1))
  273. projtm(i5,i5) = 0.5d0*(projt(i1,i3,i1,i3) +
  274. . projt(i1,i3,i3,i1))
  275. projtm(i6,i5) = 0.5d0*(projt(i2,i3,i1,i3) +
  276. . projt(i2,i3,i3,i1))
  277. projtm(i6,i4) = 0.5d0*(projt(i2,i3,i1,i2) +
  278. . projt(i2,i3,i2,i1))
  279. projtm(i5,i6) = 0.5d0*(projt(i1,i3,i2,i3) +
  280. . projt(i1,i3,i3,i2))
  281. projtm(i6,i6) = 0.5d0*(projt(i2,i3,i2,i3) +
  282. . projt(i2,i3,i3,i2))
  283.  
  284. c Identity matrix [6,6] U used to calculate P- and Q-
  285. c for P-: not really necessary as sig- = sig - sig+
  286. c for Q-: necessary for calculation of consistent algorit. tangent
  287. c Projection operators of the effective stresses P-
  288. do jloc=i1,i6
  289. do iloc=i1,i6
  290. projcm(iloc,jloc) = U(iloc,jloc) - projtm(iloc,jloc)
  291. end do
  292. end do
  293.  
  294. c ===================================================================
  295. c ! !
  296. c ! POSITIVE AND NEGATIVE PARTS OF THE EFFECTIVE STRESS TENSOR !
  297. c ! ========================================================== !
  298. c ===================================================================
  299.  
  300. c projtm has to be transposed
  301. call mulATB(projtm,esigma6v,esigmat,i6,i6,i6,i1)
  302. c [6;1] = [6;6] x [6x1]
  303. do iloc=i1,i6
  304. esigmac(iloc) = esigma6v(iloc) - esigmat(iloc)
  305. end do
  306.  
  307. c Save the negative part of the effective stress tensor to be used in
  308. c the transient creep strain calculation
  309. do iloc=i1,i6
  310. parahot3(23+iloc) = esigmac(iloc)
  311. end do
  312.  
  313. c ========================================
  314. c ! !
  315. c ! NOMINAL STRESSES !
  316. c ! ================ !
  317. c ========================================
  318.  
  319. c SIGMA = (1-dt) * SIG,eff,tension + (1-dc) * SIG,eff,comp
  320.  
  321. c concrete completely damaged
  322. if ((lfulldamage).or.(dc.gt.0.98).or.(dt.gt.0.98)) then
  323. do iloc=1,6
  324. SIGMA(iloc) = 0.d0
  325. end do
  326. do iloc=i1,i6
  327. parahot3(idimpara3-i13+iloc) = sigma(iloc)
  328. end do
  329. do jloc=1,6
  330. do iloc=1,6
  331. H66t(iloc,jloc) = 0.d0
  332. end do
  333. end do
  334. return
  335. endif
  336.  
  337. dt = parahot3(idimpara3-33)
  338. do iloc = i1,i6
  339. SIGMA(iloc) = (1.0d0-dt)*esigmat(iloc) +
  340. . (1.0d0-dc)*esigmac(iloc)
  341. end do
  342.  
  343. do iloc=i1,i6
  344. parahot3(idimpara3-i13+iloc) = sigma(iloc)
  345. end do
  346.  
  347. c ***************************************************************
  348. c ========================================
  349. c ! !
  350. c ! TANGENT MODULUS !
  351. c ! =============== !
  352. c ========================================
  353.  
  354. c Calculation of (nominal) consistent algorithmic tangent modulus Ct.
  355. c In DPRan3D we have calculated the effective elastoplastic tangent
  356. c modulus Dt. Now, we have to include the damage contribution.
  357.  
  358. c APEX of Drucker-Prager or equivalent stress of Rankine = 0
  359. c First, verify that we are not at the apex of the surfaces
  360. c If we are, Dt = 0 and in this case we should have Ct = 0 too
  361. rloc = 0.d0
  362. do jloc=1,6
  363. do iloc=1,6
  364. rloc = rloc+ABS(H66t(iloc,jloc))
  365. enddo
  366. enddo
  367. if (rloc.lt.1.e+4) then
  368. do jloc=1,6
  369. do iloc=1,6
  370. H66t(iloc,jloc) = 0.d0
  371. end do
  372. end do
  373. return
  374. endif
  375. c End of test APEX
  376.  
  377. do iloc=i1,6
  378. rt1(i1,iloc) = ( (at/2.)*EXP(-at*rkappat)+(6.*at/2.)*
  379. . EXP(-6.*at*rkappat) ) * rt1(i1,iloc)
  380. rc1(i1,iloc) = ( parahot3(14)*EXP(-parahot3(14)*rkappac) )
  381. . * rc1(i1,iloc)
  382. end do
  383.  
  384. c Calculation of W4
  385. c Calculation of p12, p23, p13
  386. c p12(i3,i3) = (1/2) (v1 X v2 + v2 X v1) (tensorial product)
  387. call mulABT(v1,v2,p121,i3,i1,i3,i1)
  388. c [3;3]=[3;1]x[1;3]
  389. call mulABT(v2,v1,p122,i3,i1,i3,i1)
  390. call mulABT(v2,v3,p231,i3,i1,i3,i1)
  391. call mulABT(v3,v2,p232,i3,i1,i3,i1)
  392. call mulABT(v1,v3,p131,i3,i1,i3,i1)
  393. call mulABT(v3,v1,p132,i3,i1,i3,i1)
  394. do jloc=i1,i3
  395. do iloc=i1,i3
  396. p12(iloc,jloc) = (1.0d0/2.0d0) * (p121(iloc,jloc) +
  397. . p122(iloc,jloc))
  398. p23(iloc,jloc) = (1.0d0/2.0d0) * (p231(iloc,jloc) +
  399. . p232(iloc,jloc))
  400. p13(iloc,jloc) = (1.0d0/2.0d0) * (p131(iloc,jloc) +
  401. . p132(iloc,jloc))
  402. end do
  403. end do
  404.  
  405. c Calculation of PPijkl = pij X pij
  406. c [3;3;3;3]=[3;3] X [3;3]
  407. do mloc=i1,i3
  408. do kloc=i1,i3
  409. do jloc=i1,i3
  410. do iloc=i1,i3
  411. pp12(iloc,jloc,kloc,mloc)=p12(iloc,jloc)*p12(kloc,mloc)
  412. pp23(iloc,jloc,kloc,mloc)=p23(iloc,jloc)*p23(kloc,mloc)
  413. pp13(iloc,jloc,kloc,mloc)=p13(iloc,jloc)*p13(kloc,mloc)
  414. end do
  415. end do
  416. end do
  417. end do
  418.  
  419. c Calculation of Q+ ! [3;3;3;3]
  420. do mloc=i1,i3
  421. do kloc=i1,i3
  422. do jloc=i1,i3
  423. do iloc=i1,i3
  424. Qtens(iloc,jloc,kloc,mloc) = r0
  425. end do
  426. end do
  427. end do
  428. end do
  429. if ((S1-S2).ne.r0) then
  430. rloc = (h1*S1) - (h2*S2)
  431. do mloc=i1,i3
  432. do kloc=i1,i3
  433. do jloc=i1,i3
  434. do iloc=i1,i3
  435. Qtens(iloc,jloc,kloc,mloc)=Qtens(iloc,jloc,kloc,mloc)
  436. . + 2.0d0 * (rloc/(S1-S2))
  437. . * pp12(iloc,jloc,kloc,mloc)
  438. end do
  439. end do
  440. end do
  441. end do
  442. else
  443. continue
  444. endif
  445. if ((S2-S3).ne.r0) then
  446. rloc = (h2*S2) - (h3*S3)
  447. do mloc=i1,i3
  448. do kloc=i1,i3
  449. do jloc=i1,i3
  450. do iloc=i1,i3
  451. Qtens(iloc,jloc,kloc,mloc)=Qtens(iloc,jloc,kloc,mloc)
  452. . + 2.0d0 * (rloc/(S2-S3))
  453. . * pp23(iloc,jloc,kloc,mloc)
  454. end do
  455. end do
  456. end do
  457. end do
  458. else
  459. continue
  460. endif
  461. if ((S1-S3).ne.r0) then
  462. rloc = (h1*S1) - (h3*S3)
  463. do mloc=i1,i3
  464. do kloc=i1,i3
  465. do jloc=i1,i3
  466. do iloc=i1,i3
  467. Qtens(iloc,jloc,kloc,mloc)=Qtens(iloc,jloc,kloc,mloc)
  468. . + 2.0d0 * (rloc/(S1-S3))
  469. . * pp13(iloc,jloc,kloc,mloc)
  470. end do
  471. end do
  472. end do
  473. end do
  474. else
  475. continue
  476. endif
  477. do mloc=i1,i3
  478. do kloc=i1,i3
  479. do jloc=i1,i3
  480. do iloc=i1,i3
  481. Qtens(iloc,jloc,kloc,mloc) = projt(iloc,jloc,kloc,mloc)
  482. . + Qtens(iloc,jloc,kloc,mloc)
  483. c [3;3;3;3]
  484. end do
  485. end do
  486. end do
  487. end do
  488.  
  489. c Calculation of the matrix [6;6] Qtens
  490. c Transform the [3;3;3;3] Qtens tensor into [6;6] matrix (Voigt notation)
  491. do jloc=1,3
  492. do iloc=1,3
  493. Qtensm(iloc,jloc) = Qtens(iloc,iloc,jloc,jloc)
  494. end do
  495. end do
  496. do iloc=1,3
  497. Qtensm(iloc,4) = 0.5d0*(Qtens(iloc,iloc,1,2)
  498. . + Qtens(iloc,iloc,2,1))
  499. Qtensm(iloc,5) = 0.5d0*(Qtens(iloc,iloc,1,3)
  500. . + Qtens(iloc,iloc,3,1))
  501. Qtensm(iloc,6) = 0.5d0*(Qtens(iloc,iloc,2,3)
  502. . + Qtens(iloc,iloc,3,2))
  503. end do
  504. do jloc=1,3
  505. Qtensm(4,jloc) = 0.5d0*(Qtens(1,2,jloc,jloc)
  506. . + Qtens(2,1,jloc,jloc))
  507. Qtensm(5,jloc) = 0.5d0*(Qtens(1,3,jloc,jloc)
  508. . + Qtens(3,1,jloc,jloc))
  509. Qtensm(6,jloc) = 0.5d0*(Qtens(2,3,jloc,jloc)
  510. . + Qtens(3,2,jloc,jloc))
  511. end do
  512. Qtensm(4,4) = 0.5d0*(Qtens(1,2,1,2) + Qtens(1,2,2,1))
  513. Qtensm(4,5) = 0.5d0*(Qtens(1,2,1,3) + Qtens(1,2,3,1))
  514. Qtensm(4,6) = 0.5d0*(Qtens(1,2,2,3) + Qtens(1,2,3,2))
  515. Qtensm(5,4) = 0.5d0*(Qtens(1,3,1,2) + Qtens(1,3,2,1))
  516. Qtensm(5,5) = 0.5d0*(Qtens(1,3,1,3) + Qtens(1,3,3,1))
  517. Qtensm(5,6) = 0.5d0*(Qtens(1,3,2,3) + Qtens(1,3,3,2))
  518. Qtensm(6,4) = 0.5d0*(Qtens(2,3,1,2) + Qtens(2,3,2,1))
  519. Qtensm(6,5) = 0.5d0*(Qtens(2,3,1,3) + Qtens(2,3,3,1))
  520. Qtensm(6,6) = 0.5d0*(Qtens(2,3,2,3) + Qtens(2,3,3,2))
  521.  
  522. do jloc=1,6
  523. do iloc=1,6
  524. Qcompm(iloc,jloc) = U(iloc,jloc) - Qtensm(iloc,jloc)
  525. enddo
  526. enddo
  527. do iloc=1,3
  528. Qcompm(iloc,4) = 0.5d0*U(iloc,4) - Qtensm(iloc,4)
  529. Qcompm(iloc,5) = 0.5d0*U(iloc,5) - Qtensm(iloc,5)
  530. Qcompm(iloc,6) = 0.5d0*U(iloc,6) - Qtensm(iloc,6)
  531. enddo
  532. do jloc=1,3
  533. Qcompm(4,jloc) = 0.5d0*U(4,jloc) - Qtensm(4,jloc)
  534. Qcompm(5,jloc) = 0.5d0*U(5,jloc) - Qtensm(5,jloc)
  535. Qcompm(6,jloc) = 0.5d0*U(6,jloc) - Qtensm(6,jloc)
  536. enddo
  537. do jloc=4,6
  538. do iloc=4,6
  539. Qcompm(iloc,jloc) = 0.5d0*U(iloc,jloc) - Qtensm(iloc,jloc)
  540. enddo
  541. enddo
  542.  
  543. do jloc=i1,6
  544. do iloc=i1,6
  545. w4m(iloc,jloc) = dt*Qtensm(iloc,jloc) + dc*Qcompm(iloc,jloc)
  546. end do
  547. end do
  548.  
  549. c Calculation of (nominal) consistent algorithmic tangent modulus Ct.
  550. c [6;6]
  551. c Ct = (U - w4m) * Dt - sig,eff,+ * rt1 - sig,eff,- * rc1
  552. c In this subroutine Dt is called H66t (comes from DPRan3d)
  553. c and at the end, Ct has to be stored into H66t
  554.  
  555. c Calculation of the product (U - w4m) * Dt
  556. do jloc=i1,i6
  557. do iloc=i1,i6
  558. w4m(iloc,jloc) = U(iloc,jloc) - w4m(iloc,jloc)
  559. end do
  560. end do
  561. call mulAB(w4m,H66t,H66temp,i6,i6,i6,i6)
  562. c [6;6]
  563. do j=1,6
  564. do i=1,6
  565. H66t(i,j) = H66temp(i,j)
  566. enddo
  567. enddo
  568.  
  569. c esigmat and esigmac have been calculated above
  570. call mulAB(esigmat,rt1,sigrt,i6,i1,i1,i6)
  571. c [6;6] = [6;1] * [1;6]
  572. call mulAB(esigmac,rc1,sigrc,i6,i1,i1,i6)
  573.  
  574. do jloc=i1,i6
  575. do iloc=i1,i6
  576. H66t(iloc,jloc) = H66t(iloc,jloc)
  577. . - sigrt(iloc,jloc)
  578. . - sigrc(iloc,jloc)
  579. end do
  580. end do
  581.  
  582. c END OF TANGENT MODULUS
  583.  
  584. RETURN
  585. END
  586.  
  587.  
  588.  
  589.  
  590.  
  591.  
  592.  
  593.  
  594.  
  595.  
  596.  
  597.  
  598.  
  599.  
  600.  
  601.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales