Télécharger endorf.eso

Retour à la liste

Numérotation des lignes :

endorf
  1. C ENDORF SOURCE PV090527 23/01/27 21:15:30 11574
  2. subroutine endorf(sigref,eplrf,syr,sur,epsrf,epu,hplr,
  3. # yor,wpr,lcr,damrt0,damrtf,damrt1,damrt2,damrc0,damrcf,refr1,
  4. # epsrt,epsrc,plr_iso,diff0,difff,Er,Dr,Hi,grad_sigr,
  5. # tauie,tyr,dami0,damif,epsmf,eplr0,sigr0,METHODNL,epsm0,epsr0,
  6. # depsa,Et0,Etf,Hs0,Hsf,sigrf,spref,epart,eparv,source_nl,
  7. # coeff_nlr,taui0,tauif,coeff_nli,eprkf,eprmf,
  8. # hypo,istep,err1,garder,endo_interface_explicit)
  9.  
  10. c tables de dimension fixe pour resolution des sytemes lineaires
  11. implicit real*8 (a-h,o-z)
  12. implicit integer (i-n)
  13.  
  14. c variables externes
  15. real*8 sigref,eplrf,syr,sur,epu,hplr,yor,wpr,lcr,coeff_nlr
  16. real*8 damrt1,damrt2,refr1,epsrt,epsrf,diff0,difff
  17. real*8 Er,Dr,Hi,lcgr,grad_sigr,dami0,damif,taui,tyr,damrtf
  18. real*8 eplr0,epsmf,damrcf,sigr0,depsa,epsm0,epsr0,eparv
  19. real*8 Et0,Etf,Hs0,Hsf,damrt0,damrc0,sigrf,source_nl,spref
  20. real*8 eprkf,eprmf,depspl,hypo,epsrc
  21. logical plr_iso,METHODNL,garder
  22. c variables locales
  23. real*8 damrt3,damrt4,epr,dmax,xref,epse,epart
  24. integer irest0
  25. real*8 precision1,epse_max,damrc3,damrt5,Hid,drf,dr0,Ert
  26. real*8 tauie,taui0,tauif,coeff_nli,eppr_et,coeffd
  27. c coeff de reduction d ecrouissage pour la phase 2
  28. real*8 deplrdt2,eprk
  29. integer err1,istep
  30. parameter (precision1=1.0d-5,dmax=1.d0-precision1)
  31. real*8 DPI,xx,depsrr
  32. PARAMETER (DPI=3.141592653589793D0)
  33. c gestion de la methode de calcul du module du renfort
  34. logical tangent,endo_interface_explicit
  35. parameter (tangent=.false.)
  36. c coeff correcteur de la longeur pour prendre en compte
  37. c l elargissement de la zone plasqtique endommagee
  38. c due a l intrpoation des contraintes non locales lorsque
  39. c l endo est local
  40. real*8 coeff_lc
  41. parameter (coeff_lc=1.0d0)
  42.  
  43. c*************** test des hypotheses de calculs prevues ***************
  44. c test calcul en module tangent
  45. if(tangent) then
  46. print*,'Module tangent non disponible'
  47. print*,'dans endorf'
  48. err1=1
  49. return
  50. end if
  51.  
  52. c test option de plasticite du renfort
  53. if (plr_iso) then
  54. print*,'Plasticite isotrope des renforts'
  55. print*,'incompatible avec endo des des renfort'
  56. print*,'dans endorf'
  57. err1=1
  58. return
  59. end if
  60.  
  61. c**************** initialisation ***************************************
  62. err1=0
  63.  
  64. c**************** endommagement du renfort *****************************
  65.  
  66. if (sur.ge.syr) then
  67. c endommagement de palier d ecrouissage
  68. if(garder) then
  69. c la solution actuelle est issue d'hypotheses coherentes
  70. if(sigref.gt.sur) then
  71. c pas d endo de compression
  72. damrc3=0.d0
  73. c cas de la traction
  74. damrt3=1.d0-sur/sigref
  75. else
  76. if(sigref.lt.-sur) then
  77. c pas d endo de traction
  78. damrt3=0.d0
  79. c cas de la compression
  80. damrc3=1.d0+sur/sigref
  81. else
  82. c pas d endo
  83. damrt3=0.d0
  84. damrc3=0.d0
  85. end if
  86. end if
  87. damrt1=max(damrt1,damrt3)
  88. damrcf=max(damrc0,damrc3)
  89. else
  90. c une hypothese incoherente detecteee
  91. c on reprend les endos initiaux
  92. damrcf=damrc0
  93. end if
  94. else
  95. print*,'Dans endorf il faur sur > syr'
  96. err1=1
  97. return
  98. end if
  99.  
  100. c increment de deformation plastique
  101. depspl=eplrf-eplr0
  102.  
  103. if((epu.ge.0.).and.(wpr.gt.0.)) then
  104. c les donnees pour la rupture sont coherentes
  105. c deformation caracteristique de rupture en traction
  106. c en fonction de la longueur caracteristique d integration lcr
  107. c - decroissance Gaussienne
  108. c eprk=epu+(wpr/(lcr*sur))*sqrt(2.d0/DPi)*(yor/(yor+hplr))
  109. c - decroissance exponentielle
  110. eprk=epu+(wpr/(coeff_lc*lcr*sur))*(yor/(yor+hplr))
  111. c - decroissance parabolique
  112. c eprk=epu+(3.d0*wpr/(2.d0*lcr*sur))*(yor/(yor+hplr))
  113. if (garder) then
  114. c les hypotheses d ecoulement etaient coherentes
  115. if(depspl.gt.0.d0) then
  116. c cumul des deformation positive pour endo tract
  117. epsrt=epsrt+depspl
  118. c remise a zero de la def de refermeture
  119. epsrc=0.d0
  120. else
  121. c actualisation de la refermeture
  122. epsrc=epsrc+depspl
  123. end if
  124. c evolution de l endo de traction
  125. if(epsrt.gt.epu) then
  126. c - endo Gaussien
  127. c xx=((epsrt-epu)/(eprk-epu))**2
  128. c damrt5=1.d0-exp(-0.5d0*xx)
  129. c - endo exponentielle
  130. xx=((epsrt-epu)/(eprk-epu))
  131. damrt5=1.d0-exp(-xx)
  132. c - endo parabolique
  133. c if(epsrt.lt.eprk) then
  134. c damrt5=((epsrt-epu)/(eprk-epu))**2
  135. c else
  136. c damrt5=dmax
  137. c end if
  138. damrt5=min(damrt5,dmax)
  139. else
  140. damrt5=damrt2
  141. end if
  142. c actualisation endo de traction localise
  143. damrt2=max(damrt5,damrt2)
  144. else
  145. c les hypotheses d ecoulement etaient incoherentes
  146. c on reprend le cas converge precedent
  147. continue
  148. end if
  149. else
  150. print*,'Pas de donnees pour la rupture du renfort'
  151. print*,'dans endorf'
  152. err1=1
  153. return
  154. end if
  155.  
  156. c endommagement resultant en traction
  157. damrt4=1.d0-(1.d0-damrt1)*(1.d0-damrt2)
  158.  
  159. c condition de croissance pour l endo de traction
  160. damrtf=max(damrt0,damrt4)
  161.  
  162. c endommagement maximum autorisé
  163. damrtf=min(damrtf,dmax)
  164.  
  165. c********** coeff de refermeture de fissure ***************************
  166.  
  167. c increment total de def plastique necessaire pour refermer
  168. c - cas de l endo Gaussien
  169. c depsrr=sqrt(0.2D1)*sqrt(-log((1.d0-damrt2))*(eprk-epu)**2)
  170. c - cas de l endo exponentiel
  171. depsrr=(epu-eprk)*log(1.d0-damrt2)
  172. c - cas de l endo parabolique
  173. c depsrr=sqrt(damrt2*((eprk-epu)**2))
  174. c test somme des increments neg - somme des incr positifs
  175. if(epsrc.gt.-depsrr) then
  176. c le cumul est en faveur de l ouverture
  177. c - endo Gaussien
  178. c xx=((depsrr+epsrc)/(eprk-epu))**2
  179. c refr1=exp(-0.5d0*xx)
  180. c - endo exponentiel
  181. xx=((depsrr+epsrc)/(eprk-epu))
  182. refr1=exp(-xx)
  183. c - endo parabolique
  184. c refr1=1.d0-((depsrr+epsrc)/(eprk-epu))**2
  185. c bornes de la fonction de refermeture
  186. refr1=min(1.d0,max(precision1,refr1))
  187. else
  188. c le cumul est une refermeture
  189. refr1=1.d0
  190. end if
  191.  
  192. c**************** endommagement de l interface ************************
  193.  
  194. if((METHODNL.and.(istep.gt.1).and.(garder).and.
  195. # (.not.endo_interface_explicit)).or.
  196. # ((METHODNL.and.(istep.eq.1).and.
  197. # (endo_interface_explicit))))then
  198. c les hypotheses etaient coherentes
  199. c approximation de la contrainte de cisaillement actuelle
  200. taui=-grad_sigr*Dr/4.d0
  201. c test de plasticification de l interface
  202. if((tyr.eq.0.).and.(taui.ne.0.)) then
  203. print*,'Dans endorf'
  204. print*,'Resistance de l interface non fournie'
  205. print*,'pour endorf appele par fldo3d'
  206. print*,'fournir TYRi non nul pour continuer...'
  207. err1=1
  208. return
  209. else
  210. c contrainte effective d interface
  211. tauie=taui/coeff_nli
  212. if(abs(tauie).gt.tyr) then
  213. damif=1.d0-tyr/abs(tauie)
  214. c condition de croissance de lendo de l interface
  215. damif=max(damif,dami0)
  216. damif=min(damif,dmax)
  217. else
  218. damif=dami0
  219. end if
  220. c raideur de l interface endommagee
  221. Hid=Hi*(1.d0-damif)
  222. end if
  223. else if (METHODNL)then
  224. if(((.not.endo_interface_explicit).and.
  225. # ((istep.eq.1).or.(.not.garder))).or.
  226. # (endo_interface_explicit.and.(istep.ne.1)))then
  227. c on reprend la valeur convergee
  228. damif=dami0
  229. Hid=Hi*(1.d0-dami0)
  230. tauie=taui0/coeff_nli
  231. else
  232. print*,'endorf: Cas non prevu / interface'
  233. print*,METHODNL,endo_interface_explicit
  234. print*,istep,garder
  235. err1=1
  236. return
  237. end if
  238. else if(.not.METHODNL) then
  239. c calcul local
  240. Hid=Hi
  241. dami0=0.d0
  242. damif=dami0
  243. taui=0.d0
  244. tauie=taui
  245. else
  246. c on initialise a la valeur convergee precedente
  247. tauie=taui0/coeff_nli
  248. damif=dami0
  249. Hid=Hi*(1.d0-damif)
  250. end if
  251. c module de l interface
  252. if(METHODNL) then
  253. c module secant
  254. Hsf=Hid
  255. else
  256. c module initial
  257. Hsf=Hi
  258. end if
  259. c contrainte dans l interface apres endommagement
  260. coeff_nli=(1.d0-damif)
  261. tauif=tauie*coeff_nli
  262.  
  263.  
  264. c********** coeff de diffusion pour l adherence ************************
  265.  
  266. c calcul en module secant
  267. if(.not.garder) then
  268. c on reprend la contrainte convergee pour calculer les coeffs
  269. c non lineaires necessaire au calcul de diffusion
  270. sigref=sigr0
  271. end if
  272. if (sigref.ge.0.d0) then
  273. c avec l endo de traction
  274. coeff_nlr=(1.0d0-damrtf)
  275. else
  276. c avec l endo de compression et la refermeture
  277. coeff_nlr=(1.0d0-damrcf)*refr1
  278. end if
  279.  
  280.  
  281. c module endommagee actuel pour le renfort
  282. Ert=Er*coeff_nlr
  283.  
  284. c coeff de diffusion
  285. if(METHODNL) then
  286. c .and.(istep.ne.1)) then
  287. c actualisation du coeff pour le pas suivant
  288. difff=(Ert*Dr)/(4.d0*Hid)
  289. c calcul suivant hypothese d ecoulement
  290. difff=difff*(hplr/(hplr+yor*abs(hypo)))
  291. c difff=4.0d-2
  292. c else if (METHODNL) then
  293. c print*,'Cas imprevu sur lnl dans endorf'
  294. c err1=1
  295. c return
  296. else if ((istep.eq.0).or.(.not.METHODNL)) then
  297. c calcul local
  298. difff=0.d0
  299. Etf=Et0
  300. Hsf=Hs0
  301. else
  302. print*,'Cas imprevu sur lnl dans endorf'
  303. err1=1
  304. return
  305. end if
  306. c print*,'endo renfort lnl',lnl,' istep',istep
  307.  
  308. c module actuel du renfort
  309. if(METHODNL) then
  310. Etf=Ert
  311. else
  312. Etf=Er
  313. end if
  314.  
  315. c*********** contraintes apparentes dans les renforts ******************
  316.  
  317. c calcul des contraintes apparentes dans les renforts
  318. sigrf=sigref*coeff_nlr
  319. c test coherence de la contrainte
  320. if(abs(sigrf).gt.sur*(1.d0+precision1)) then
  321. c incoherence detectee
  322. print*,'Contrainte anormale dans dans endorf'
  323. print*,'|sirf|>sur',sigrf,sur
  324. print*,'Coeff de non linearite',coeff_nlr
  325. print*,'valeur des endommagements'
  326. print*,'damrt1',damrt1
  327. print*,'damrt2',damrt2
  328. print*,'damrtf',damrtf
  329. print*,'damrcf',damrcf
  330. print*,'refermeture',refr1
  331. print*,'Def plastique actuelle',eplrf
  332. print*,'Def totale actuelle',epsrf
  333. if(damrt2.ge.dmax) then
  334. print*,'Def de refermeture actuelle',epsrt
  335. end if
  336. err1=1
  337. return
  338. end if
  339. if(difff.lt.0.) then
  340. print*,'diffusion negative dans endorf'
  341. err1=1
  342. return
  343. end if
  344.  
  345. c*********** source non locale *****************************************
  346.  
  347. if(METHODNL) then
  348. if (hypo.ne.0.) then
  349. c hypothese avec ecoulement plastique
  350. eppr_et=syr/hplr
  351. if (hypo.eq.1.) then
  352. source_nl=Hid*(epsmf-(eparv -eppr_et))
  353. else if(hypo.eq.-1.) then
  354. source_nl=Hid*(epsmf-(eparv+eppr_et))
  355. else
  356. print*,'valeur anormale de hypo ds endorf'
  357. err1=1
  358. return
  359. end if
  360. else
  361. c hypothese tir visco elastique
  362. source_nl=Hid*(epsmf-epart)
  363. end if
  364. c coeff multiplicateur des sources
  365. coeffd=4.d0*difff/Dr
  366. source_nl=(epsmf-epart)*Ert
  367. c print*,'endorf, source_nl=',source_nl
  368. c print*,'endorf, Difff=',difff
  369. else
  370. c calcul local: source inutile
  371. source_nl=0.d0
  372. end if
  373.  
  374. c***********************************************************************
  375.  
  376. return
  377. end
  378.  
  379.  
  380.  

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