OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
delamination.F
Go to the documentation of this file.
1Copyright> OpenRadioss
2Copyright> Copyright (C) 1986-2025 Altair Engineering Inc.
3Copyright>
4Copyright> This program is free software: you can redistribute it and/or modify
5Copyright> it under the terms of the GNU Affero General Public License as published by
6Copyright> the Free Software Foundation, either version 3 of the License, or
7Copyright> (at your option) any later version.
8Copyright>
9Copyright> This program is distributed in the hope that it will be useful,
10Copyright> but WITHOUT ANY WARRANTY; without even the implied warranty of
11Copyright> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12Copyright> GNU Affero General Public License for more details.
13Copyright>
14Copyright> You should have received a copy of the GNU Affero General Public License
15Copyright> along with this program. If not, see <https://www.gnu.org/licenses/>.
16Copyright>
17Copyright>
18Copyright> Commercial Alternative: Altair Radioss Software
19Copyright>
20Copyright> As an alternative to this open-source version, Altair also offers Altair Radioss
21Copyright> software under a commercial license. Contact Altair to discuss further if the
22Copyright> commercial version may interest you: https://www.altair.com/radioss/.
23!||====================================================================
24!|| delamination ../engine/source/properties/composite_options/stack/delamination.F
25!||--- called by ------------------------------------------------------
26!|| cmain3 ../engine/source/materials/mat_share/cmain3.F
27!||--- calls -----------------------------------------------------
28!|| delm01law ../engine/source/properties/composite_options/stack/delm01law.F
29!|| delm02law ../engine/source/properties/composite_options/stack/delm02law.F
30!|| delm24law ../engine/source/properties/composite_options/stack/delm24law.F
31!||--- uses -----------------------------------------------------
32!|| elbufdef_mod ../common_source/modules/mat_elem/elbufdef_mod.F90
33!|| mat_elem_mod ../common_source/modules/mat_elem/mat_elem_mod.F90
34!||====================================================================
35 SUBROUTINE delamination(ELBUF_STR,MAT_PARAM,
36 1 JFT ,JLT ,IR ,IS ,NPT ,
37 2 MAT_IPLY ,IPM ,PM ,BUFMAT ,NPF ,
38 3 TF ,DT1C ,NGL ,OFF ,TH_IPLY ,
39 4 DEL_PLY ,SIG ,OFFI ,A11 ,FOR ,
40 5 MOM ,PLY_F ,THK0 ,SHF ,EXZ ,
41 6 EYZ ,AREA ,PID ,GEO ,SSP ,
42 7 POSLY ,THKLY ,KXX ,KYY ,KXY ,
43 . DEXZ ,DEYZ ,EINT ,GSTR ,NEL ,
44 . NUMMAT )
45C----------------------------------------------
46C M o d u l e s
47C-----------------------------------------------
48 USE mat_elem_mod
49 USE elbufdef_mod
50C-----------------------------------------------
51C I m p l i c i t T y p e s
52C-----------------------------------------------
53#include "implicit_f.inc"
54#include "comlock.inc"
55C-----------------------------------------------
56C G l o b a l P a r a m e t e r s
57C-----------------------------------------------
58#include "mvsiz_p.inc"
59C-----------------------------------------------
60C C o m m o n B l o c k s
61C-----------------------------------------------
62#include "param_c.inc"
63#include "com08_c.inc"
64C-----------------------------------------------
65C D u m m y A r g u m e n t s
66C-----------------------------------------------
67 INTEGER ,INTENT(IN) :: NUMMAT
68 INTEGER JFT,JLT,IR,IS,NPT,NEL
69 INTEGER MAT_IPLY(MVSIZ,*),IPM(NPROPMI,*),NGL(*),NPF(*),PID(*)
70C REAL
71 my_real
72 . PM(NPROPM,*), TH_IPLY(MVSIZ,NPT),
73 . DEL_PLY(MVSIZ,3,NPT),OFF(*),TF(*),DT1C(*),
74 . BUFMAT(*), SIG(MVSIZ,3,NPT),OFFI(MVSIZ,*),A11(MVSIZ,NPT),
75 . THK0(*),FOR(NEL,5),PLY_F(MVSIZ,5,NPT),
76 . exz(*),eyz(*),shf(*),area(*),mom(nel,3),ssp(*),geo(npropg,*),
77 . posly(*),thkly(*), kxx(*), kyy(*), kxy(*), dexz(*), deyz(*),
78 . eint(jlt,2), gstr(nel,8)
79 TYPE (ELBUF_STRUCT_), TARGET :: ELBUF_STR
80 TYPE (MATPARAM_STRUCT_) ,DIMENSION(NUMMAT) ,INTENT(IN) :: MAT_PARAM
81C-----------------------------------------------
82C L o c a l V a r i a b l e s
83C-----------------------------------------------
84 INTEGER I,J,K,JJ,IG,IC,II,IFL,IPLY,NVARF,NFAIL,JVAR,NVARF_MAX,
85 . IPMAT_IPLY,JINF,JSUP,
86 . ILAW,IRUPT,IMAT,JJINF,JJSUP,NBDEL,NF,NL,KK(5)
87 INTEGER MAT(MVSIZ),IFAIL(MVSIZ),MATF(MVSIZ),
88 . nbiply_del(mvsiz),niply_del(mvsiz,npt),jly,
89 . nindx,indx(mvsiz)
90C REAL
91 my_real
92 . syz(mvsiz),sxz(mvsiz),szz(mvsiz) , du_iply(3,mvsiz),
93 . syz0(mvsiz),sxz0(mvsiz),szz0(mvsiz),
94 . epsyz_ip(mvsiz),epsxz_ip(mvsiz),epszz_ip(mvsiz),reduc(mvsiz,npt),
95 . degmb(mvsiz),degfx(mvsiz)
96 my_real
97 . nu, g, e33,f4,f5,vol,g0,dd,fac,dtinv,sspi,visc,
98 . wmc,thky,scale,th(npt),
99 . dezz_ip(mvsiz),deyz_ip(mvsiz),dexz_ip(mvsiz),thi(mvsiz)
100c
101 my_real,
102 . DIMENSION(:) ,POINTER :: uvarf
103 TYPE(buf_intlay_) ,POINTER :: INTLAY
104 TYPE(BUF_INTLOC_) ,POINTER :: ILBUF
105 TYPE(buf_fail_) ,POINTER :: FBUF
106 TYPE(L_BUFEL_) ,POINTER :: LBUF
107C======================================================================|
108C EPS POINT EQUIVALENT (au sens energetique)
109C-----------------------------------------------------------
110C interface
111C-----------------------------------------------------------
112 sig(jft:jlt,1:3,1:elbuf_str%NINTLAY)=zero
113!
114 DO i=1,5
115 kk(i) = nel*(i-1)
116 ENDDO
117!
118C
119 nindx=0
120 DO i=jft,jlt
121 thi(i)= zero
122 IF(off(i)/=zero)THEN
123 nindx=nindx+1
124 indx(nindx)=i
125 END IF
126 ENDDO
127
128 DO k=1,nindx
129 i=indx(k)
130 nbiply_del(i) = 0
131 DO iply = 1,elbuf_str%NINTLAY
132 niply_del(i,iply) = 0
133 thi(i) = thi(i) + th_iply(i,iply)
134 ENDDO
135 ENDDO
136
137 DO iply = 1,elbuf_str%NINTLAY
138 jinf = iply
139 jsup = iply + 1
140 intlay => elbuf_str%INTLAY(iply)
141 ilbuf => elbuf_str%INTLAY(iply)%ILBUF(ir,is)
142 fbuf => elbuf_str%INTLAY(iply)%FAIL(ir,is)
143 nfail = intlay%NFAIL
144 imat = intlay%IMAT
145 ilaw = intlay%ILAW
146c
147 DO k=1,nindx
148 i=indx(k)
149C
150C Actuallement on utilise qu'une interface elastique isotrope. (Law01)
151C
152C Only with law59 (K33, K13, K23)
153c
154 e33 = pm(20, mat_iply(i, iply))*off(i)
155 g = pm(22, mat_iply(i, iply))*off(i)
156 g0 = g*th_iply(i,iply)*off(i)
157!! FAC = TH_IPLY(I,IPLY)/THK0(I)
158 fac= th_iply(i,iply)/thi(i)
159C
160 a11(i,iply) = e33 + g
161 IF(ilaw /= 59) THEN
162 g = two*g/th_iply(i,iply)*off(i)
163 e33 = e33/th_iply(i,iply)*off(i)
164 a11(i,iply) = e33 + g
165 g0 = pm(22, mat_iply(i, iply))*off(i)
166 ENDIF
167 sxz0(i) = g*(del_ply(i,1,jsup) - del_ply(i,1,jinf))
168 syz0(i) = g*(del_ply(i,2,jsup) - del_ply(i,2,jinf))
169C
170C warns G*GSTR <=> G shell or G interply ? ::
171 sig(i,2,iply)= g*(del_ply(i,1,jsup) - del_ply(i,1,jinf))
172 . + g0*gstr(i,5)
173 sig(i,1,iply)= g*(del_ply(i,2,jsup) - del_ply(i,2,jinf))
174 . + g0*gstr(i,4)
175 sig(i,3,iply)= e33*(del_ply(i,3,jsup) - del_ply(i,3,jinf))
176C
177 mat(i) = mat_iply(i,iply)
178 syz(i) = sig(i,1,iply) ! syz
179 sxz(i) = sig(i,2,iply) ! sxz
180 szz(i) = sig(i,3,iply) ! szz
181 szz0(i) = szz(i)
182C
183C compute ply sliding
184C
185 du_iply(1,i) = (del_ply(i,1,jsup) - del_ply(i,1,jinf))
186 du_iply(2,i) = (del_ply(i,2,jsup) - del_ply(i,2,jinf))
187 du_iply(3,i) = (del_ply(i,3,jsup) - del_ply(i,3,jinf))
188C
189C compute strain
190 epszz_ip(i)= du_iply(3,i)/th_iply(i,iply)
191 epsyz_ip(i)= du_iply(2,i)/th_iply(i,iply)
192 epsxz_ip(i)= du_iply(1,i)/th_iply(i,iply)
193C
194 ii = 3*(i-1)
195 dezz_ip(i) = epszz_ip(i) - ilbuf%EPS(ii + 1)
196 deyz_ip(i) = epsyz_ip(i) - ilbuf%EPS(ii + 2)
197 dexz_ip(i) = epsxz_ip(i) - ilbuf%EPS(ii + 3)
198C
199 ilbuf%EPS(ii + 1) = epszz_ip(i)
200 ilbuf%EPS(ii + 2) = epsyz_ip(i)
201 ilbuf%EPS(ii + 3) = epsxz_ip(i)
202 ENDDO
203C----------------------- -----------------------
204C delamination damage model
205C-----------------------------------------------
206 DO ifl = 1, nfail
207 uvarf => fbuf%FLOC(ifl)%VAR
208 nvarf = fbuf%FLOC(ifl)%NVAR
209 irupt = fbuf%FLOC(ifl)%ILAWF
210C
211 IF (irupt == 18) THEN
212C-- ladeveze delamination model
213 CALL delm01law(mat_param(imat)%FAIL(ifl),
214 1 jlt ,nvarf ,tt ,dt1c ,
215 3 ngl ,iply ,
216 4 off ,syz0 ,sxz0 ,szz ,uvarf ,
217 5 offi(1,iply),reduc(1,iply),intlay%COUNT,
218 . syz, sxz)
219C--- power law
220 ELSEIF (irupt == 19) THEN
221 CALL delm02law(mat_param(imat)%FAIL(ifl),
222 1 jlt ,nvarf ,tt ,dt1c ,
223 3 ngl ,iply ,
224 4 off ,syz0 ,sxz0 ,szz ,du_iply,
225 5 uvarf ,offi(1,iply),reduc(1,iply),
226 . intlay%COUNT)
227C--- total strain per direction
228 ELSEIF(irupt == 24) THEN
229 CALL delm24law(mat_param(imat)%FAIL(ifl),
230 1 jlt ,nvarf ,tt ,dt1c ,
231 3 ngl ,iply ,
232 4 off ,syz0 ,sxz0 ,szz ,epsyz_ip,
233 5 epsxz_ip,epszz_ip,uvarf ,offi(1,iply),reduc(1,iply),
234 6 intlay%COUNT)
235 ENDIF
236c
237 ENDDO ! IFL = 1, NFAIL
238c----------------------------
239C
240 DO k=1,nindx
241 i=indx(k)
242C
243 IF (int(intlay%COUNT(i)) == 4) THEN
244 fac = th_iply(i,iply)/thk0(i)
245 IF(reduc(i,iply) < 0 ) reduc(i,iply) = zero
246 sig(i,1,iply) = syz0(i) *reduc(i,iply)
247 sig(i,2,iply) = sxz0(i) *reduc(i,iply)
248 sig(i,3,iply) = szz(i) *reduc(i,iply)
249 IF(szz0(i) < zero)sig(i,3,iply) = szz0(i)
250C
251 nbiply_del(i) = nbiply_del(i) + 1
252 nbdel = nbiply_del(i)
253 niply_del(i,nbdel) = iply
254C
255 offi(i,iply) = reduc(i,iply)
256C
257c DEGMB(I) = - FOR(I,4)*DEXZ(I) - FOR(I,5)*DEYZ(I)
258c FAC = TH_IPLY(I,IPLY)/THK0(I)
259c F4=FOR(I,4)-FAC*G0*EXZ(I)*SHF(I)*OFF(I)
260c F5=FOR(I,5)-FAC*G0*EYZ(I)*SHF(I)*OFF(I)
261c IF(F4*FOR(I,4)<=ZERO)THEN
262c FOR(I,4)=ZERO
263c ELSE
264c FOR(I,4)=F4
265c END IF
266c IF(F5*FOR(I,5)<=ZERO)THEN
267c FOR(I,5)=ZERO
268c ELSE
269c FOR(I,5)=F5
270c END IF
271c DEGMB(I) = DEGMB(I) + FOR(I,4)*DEXZ(I)+FOR(I,5)*DEYZ(I)
272c EINT(I,1)= EINT(I,1) + DEGMB(I)*HALF*THK0(I)*AREA(I)
273C
274 ELSEIF(offi(i,iply) < one) THEN
275C
276 sig(i,1,iply) = syz0(i) *offi(i,iply)
277 sig(i,2,iply) = sxz0(i) *offi(i,iply)
278 sig(i,3,iply) = szz(i) *offi(i,iply)
279 IF(szz0(i) < zero)sig(i,3,iply) = szz0(i)
280C
281 ELSE
282 sig(i,1,iply) = syz0(i)
283 sig(i,2,iply) = sxz0(i)
284 sig(i,3,iply) = szz(i)
285 ENDIF
286C viscosity
287 visc =geo(20,pid(i))
288!! VISC =GEO(16,PID(I))
289!! IF(VISC == ZERO) VISC = FIVEEM2
290 fac = visc + sqrt(visc**2 + one)
291 dtinv =dt1c(i)/max(dt1c(i)**2,em20)
292 sspi=sqrt(a11(i,iply)*th_iply(i,iply)/pm(1,mat_iply(i,iply)))
293 visc =onep414*visc*pm(1,mat_iply(i,iply))*sspi
294 visc =visc*sqrt(area(i))*dtinv
295C
296 sig(i,1,iply)= sig(i,1,iply) + visc*deyz_ip(i)*offi(i,iply)
297 sig(i,2,iply)= sig(i,2,iply) + visc*dexz_ip(i)*offi(i,iply)
298 IF(sig(i,3,iply) < zero)THEN
299 sig(i,3,iply)= sig(i,3,iply) + visc*dezz_ip(i)
300 ELSE
301 sig(i,3,iply)= sig(i,3,iply) + visc*dezz_ip(i)*offi(i,iply)
302 END IF
303C
304C stress interply (output only)
305!! II = 3*(I-1)
306 ilbuf%SIG(kk(1) + i) = sig(i,3,iply)
307 ilbuf%SIG(kk(2) + i) = sig(i,1,iply)
308 ilbuf%SIG(kk(3) + i) = sig(i,2,iply)
309C
310C stability condition
311 a11(i,iply) = a11(i,iply) * fac*fac
312
313C stress
314!! ILBUF%SIG(KK(1) + I) = SIG(3,I,IPLY)
315!! ILBUF%SIG(KK(2) + I) = SIG(1,I,IPLY)
316!! ILBUF%SIG(KK(3) + I) = SIG(2,I,IPLY)
317!! ILBUF%EPS(II + 1) = EPSZZ_IP(I)
318!! ILBUF%EPS(II + 2) = EPSYZ_IP(I)
319!! ILBUF%EPS(II + 3) = EPSXZ_IP(I)
320c
321C Energie
322 IF (ir == 1 .and. is == 1) intlay%EINT(i) = zero
323 vol = area(i)*th_iply(i,iply)
324 intlay%EINT(i) = intlay%EINT(i) + half*(
325 . sig(i,3,iply)*epszz_ip(i) +
326 . sig(i,1,iply)*epsyz_ip(i)+ sig(i,2,iply)*epsxz_ip(i))*vol
327c
328 ENDDO ! I=JFT,JLT
329 ENDDO ! inter ply
330C
331C reduction of mom for shell
332 DO k=1,nindx
333 i=indx(k)
334 nbdel = nbiply_del(i)
335 IF(nbdel > 0 ) THEN
336 DO jj =1,nbdel + 1
337 th(jj) = zero
338 ENDDO
339 nf = 1
340 DO jj =1,nbdel
341 nl = niply_del(i,jj)
342 DO j =nf,nl
343 jly = (j-1)*jlt + i
344 th(jj) = th(jj) + thkly(jly)
345 ENDDO
346 nf = nl
347 ENDDO
348 jj = nbdel + 1
349 nl = elbuf_str%NINTLAY + 1
350 DO j = nf,nl
351 jly = (j-1)*jlt + i
352 th(jj) = th(jj) + thkly(jly)
353 ENDDO
354 scale =zero
355 DO jj =1,nbdel + 1
356 scale = scale + th(jj)**3
357 ENDDO
358 scale = max(scale,em01)
359C
360 degfx(i) = - mom(i,1)*kxx(i) - mom(i,2)*kyy(i)
361 + - mom(i,3)*kxy(i)
362 mom(i,1) = scale*mom(i,1)
363 mom(i,2) = scale*mom(i,2)
364 mom(i,3) = scale*mom(i,3)
365C
366 degfx(i) = degfx(i) + mom(i,1)*kxx(i)+mom(i,2)*kyy(i)
367 + + mom(i,3)*kxy(i)
368 eint(i,2)= eint(i,2) + degfx(i)*half*thk0(i)*thk0(i)*area(i)
369 ENDIF
370 ENDDO ! I=JFT,JLT
371C -----------------------------------------------------------------
372 RETURN
373 END
delamination
subroutine delamination(elbuf_str, mat_param, jft, jlt, ir, is, npt, mat_iply, ipm, pm, bufmat, npf, tf, dt1c, ngl, off, th_iply, del_ply, sig, offi, a11, for, mom, ply_f, thk0, shf, exz, eyz, area, pid, geo, ssp, posly, thkly, kxx, kyy, kxy, dexz, deyz, eint, gstr, nel, nummat)
Definition delamination.F:45
delm01law
subroutine delm01law(fail, nel, nuvar, time, timestep, ngl, iply, off, signyz0, signxz0, signzz, uvar, offi, reduc, count, signyz, signxz)
Definition delm01law.F:36
delm02law
subroutine delm02law(fail, nel, nuvar, time, timestep, ngl, iply, off, signyz, signxz, signzz, du, uvar, offi, reduc, count)
Definition delm02law.F:36
delm24law
subroutine delm24law(fail, nel, nuvar, time, timestep, ngl, iply, off, signyz, signxz, signzz, epsyz, epsxz, epszz, uvar, offi, reduc, count)
Definition delm24law.F:37
area
subroutine area(d1, x, x2, y, y2, eint, stif0)
Definition hm_read_prop32.F:567
max
#define max(a, b)
Definition macros.h:21