OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps117.F File Reference
#include "implicit_f.inc"
#include "units_c.inc"
#include "comlock.inc"
#include "sms_c.inc"
#include "lockon.inc"
#include "lockoff.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine sigeps117 (nel, nuparam, nuvar, jsms, time, timestep, uparam, uvar, area, off, offl, epszz, epsyz, epszx, depszz, depsyz, depszx, signzz, signyz, signzx, stifm, dmels, dmg, ipg, nfail, ngl, nfunc, ifunc, npf, tf)

Function/Subroutine Documentation

◆ sigeps117()

subroutine sigeps117 ( integer, intent(in) nel,
integer, intent(in) nuparam,
integer, intent(in) nuvar,
integer, intent(in) jsms,
intent(in) time,
intent(in) timestep,
intent(in) uparam,
intent(inout) uvar,
dimension(nel) area,
dimension(nel) off,
dimension(nel) offl,
dimension(nel) epszz,
dimension(nel) epsyz,
dimension(nel) epszx,
dimension(nel) depszz,
dimension(nel) depsyz,
dimension(nel) depszx,
dimension(nel) signzz,
dimension(nel) signyz,
dimension(nel) signzx,
intent(inout) stifm,
dimension(nel) dmels,
intent(out) dmg,
integer, intent(in) ipg,
integer, intent(out) nfail,
integer, dimension(nel), intent(in) ngl,
integer nfunc,
integer, dimension(nfunc) ifunc,
integer, dimension(*) npf,
tf )

Definition at line 30 of file sigeps117.F.

36C-----------------------------------------------
37C I m p l i c i t T y p e s
38C-----------------------------------------------
39#include "implicit_f.inc"
40C----------------------------------------------------------
41#include "units_c.inc"
42#include "comlock.inc"
43#include "sms_c.inc"
44C----------------------------------------------------------
45C D u m m y A r g u m e n t s
46C----------------------------------------------------------
47 INTEGER ,INTENT(IN) :: NEL,NUPARAM,NUVAR,JSMS,IPG
48 INTEGER ,INTENT(OUT) :: NFAIL
49 INTEGER ,DIMENSION(NEL) ,INTENT(IN) :: NGL
50 my_real ,INTENT(IN) :: time,timestep
51 my_real ,DIMENSION(NEL) :: off,offl,area,dmels,
52 . epszz,epsyz,epszx,depszz,depsyz,depszx,signzz,signyz,signzx
53 my_real ,DIMENSION(NEL) ,INTENT(OUT) :: dmg
54 my_real ,DIMENSION(NEL) ,INTENT(INOUT) :: stifm
55 my_real ,DIMENSION(NUPARAM) ,INTENT(IN) :: uparam
56 my_real ,DIMENSION(NEL,NUVAR) ,INTENT(INOUT) :: uvar
57C-----------------------------------------------
58 INTEGER NPF(*), NFUNC, IFUNC(NFUNC)
59 my_real finter ,tf(*)
60 EXTERNAL finter
61C-----------------------------------------------
62C L o c a l V a r i a b l e s
63C-----------------------------------------------
64 INTEGER :: I,J,K,II, I_REL, IRUPT,NINDXF
65 INTEGER ,DIMENSION(NEL) :: INDXF
66 my_real :: rho0, dm, dam, stf, et2,dtb,
67 . e_elas_n,dtinv , t_n,t_s,
68 . exp_g , exp_bk, delta0s_cst, delta0n_cst,
69 . gic , giic , e_elas_s, und_cst, utd_cst,gama
70c
71 my_real, DIMENSION(NEL) :: epsm,delta0m, beta,
72 . etha,deltafmax,epst,fac1, fac2,fac3,epsmmax ,epsn ,
73 . delta0n,delta0s,und,utd,dydx1,dydx2,length,tmax_n,tmax_s
74c-----------------------------------------------
75c UVAR(1) =
76c UVAR(2) =
77c UVAR(3) =
78c UVAR(4) =
79c UVAR(5) =
80c UVAR(6) =
81c UVAR(7) =
82c UVAR(8) =
83C=======================================================================
84C INPUT PARAMETERS INITIALIZATION
85C-----------------------------------------------
86 e_elas_n = uparam(1)
87 e_elas_s = uparam(2)
88 gama = uparam(3)
89 t_n = uparam(4)
90 t_s = uparam(5)
91 irupt = int(uparam(6) )
92 gic = uparam(7)
93 giic = uparam(8)
94 nfail = int(uparam(9))
95 exp_g = uparam(10)
96 exp_bk = uparam(11)
97 delta0n_cst = uparam(13) ! DISP puremode N
98 delta0s_cst = uparam(14) ! DISP puremode S
99 und_cst = uparam(15) ! ultimate displacement in normal direction (UND)
100 utd_cst = uparam(16) ! ultimate displacement in tangential direction (UTD)
101
102c
103 dtinv = one / (max(timestep, em20))
104 stf = e_elas_n + e_elas_s
105 nindxf = 0
106c-------------------------
107 IF (time == zero) THEN
108 DO i=1,nel
109 epsmmax(i)=zero
110 ENDDO
111 ELSE
112 DO i=1,nel
113 epsmmax(i) = uvar(i,4)
114 ENDDO
115 ENDIF
116
117 DO i=1,nel
118 epsn(i) = max(epszz(i) , zero)
119 epst(i) = sqrt(epsyz(i)**2 + epszx(i)**2)
120 epsm(i) = sqrt( epsn(i)**2 + epst(i)**2)
121 epsmmax(i) = max(epsm(i),epsmmax(i))
122 ENDDO !I=1,NEL
123 !---------------------------------------------------------
124 ! Compute damage initiation and max displacement
125 !---------------------------------------------------------
126 IF ( ifunc(1) /= 0.OR. ifunc(2) /= 0) THEN !depending on mesh size
127 length(1:nel) = sqrt(area(1:nel))
128 ENDIF
129 IF ( ifunc(1) /= 0) THEN
130 DO i=1,nel
131 tmax_n(i) = t_n * finter(ifunc(1),length(i),npf,tf,dydx1(i))
132 delta0n(i)= tmax_n(i) /e_elas_n
133 und(i) = two*gic /tmax_n(i)
134 ENDDO !I=1,NEL
135 ELSE
136 DO i=1,nel
137 delta0n(i) = delta0n_cst
138 und(i) = und_cst
139 ENDDO
140 ENDIF
141 IF ( ifunc(2) /= 0) THEN
142 DO i=1,nel
143 tmax_s(i) = t_s * finter(ifunc(2),length(i),npf,tf,dydx2(i))
144 delta0s(i)= tmax_s(i) /e_elas_s
145 utd(i) = two*giic/tmax_s(i)
146 ENDDO !I=1,NEL
147 ELSE
148 DO i=1,nel
149 delta0s(i) = delta0s_cst
150 utd(i) = utd_cst
151 ENDDO
152 ENDIF
153 !---------------------------------------------------------
154 ! Update failure initiation strains in mixed mode
155 !---------------------------------------------------------
156 DO i=1,nel
157 IF (epst(i) == zero) THEN
158 delta0m(i) = delta0n(i)
159 ELSE IF (epsn(i) == zero) THEN
160 delta0m(i) = delta0s(i)
161 ELSE ! mixed mode
162 beta(i) = abs(epst(i) / epsn(i))
163 delta0m(i)= delta0s(i)* delta0n(i)*sqrt((one + beta(i)**2)/
164 . ((delta0s(i)**2) + (beta(i)* delta0n(i))**2))
165 END IF
166 ENDDO !I=1,NEL
167!------------------------------------------------
168! damage : ultimate displacement
169!------------------------------------------------
170 ! IRUPT =2 = BK method -----------------
171 IF (irupt == 2) THEN
172 DO i=1,nel
173 IF (epst(i) == zero) THEN
174 deltafmax(i)= und(i)
175 ELSE IF (epsn(i) == zero) THEN
176 deltafmax(i)= utd(i)
177 ELSE ! mixed mode
178 fac1(i) = (e_elas_n**gama)/(one + beta(i)**2)
179 fac2(i) = (e_elas_s**gama)*(beta(i)**2)/(one + beta(i)**2)
180 fac3(i) = (fac1(i) + fac2(i) )**(one/gama)
181 !ultimate mixed mode displacement
182 deltafmax(i)= two/(delta0m(i)* fac3(i) )*
183 . (gic + (giic - gic)*((e_elas_s*beta(i)**2)/
184 . (e_elas_n + e_elas_s*beta(i)**2))**abs(exp_bk))
185 END IF
186 ENDDO
187 ELSE
188 ! IRUPT =1 = power law method-- default---------------
189 DO i=1,nel
190 IF (epst(i) == zero) THEN
191 deltafmax(i)= und(i)
192 ELSE IF (epsn(i) == zero) THEN
193 deltafmax(i)= utd(i)
194 ELSE ! mixed mode
195 fac1(i) = two*((one+beta(i)**2))/delta0m(i)
196 !ultimate mixed mode displacement
197 deltafmax(i)= fac1(i)*((e_elas_n/gic)**exp_g +
198 . (e_elas_s*(beta(i)**2)/giic)**exp_g)**(-one/exp_g)
199 END IF
200 ENDDO
201 ENDIF
202!---------------------------------------------------------------
203! damage evolution
204!---------------------------------------------------------------
205 DO i=1,nel
206 dm = epsmmax(i) - delta0m(i)
207 IF (dm > zero.AND.epsmmax(i) /= zero ) THEN
208 dam = (deltafmax(i) / epsmmax(i))*
209 . (epsmmax(i) - delta0m(i))/
210 . max((deltafmax(i) - delta0m(i)), em20)
211
212 dmg(i) = max(dmg(i), dam)!
213 dmg(i) = min(dmg(i), one)
214 IF (offl(i) == one .AND. epsmmax(i) > deltafmax(i) ) THEN
215 nindxf = nindxf+1
216 indxf(nindxf) = i
217 offl(i) = zero
218 END IF
219 END IF
220 ENDDO
221!---------------------------------------------------------------
222! Stress update
223!---------------------------------------------------------------
224 DO i=1,nel
225 IF (epszz(i) < zero ) THEN
226 signzz(i) = e_elas_n*epszz(i)
227 ELSE
228 signzz(i) = (one-dmg(i))*e_elas_n*epszz(i)
229 ENDIF
230 ! sigma shear Y
231 signyz(i) = (one-dmg(i))*e_elas_s*(epsyz(i))
232 ! sigma shear X
233 signzx(i) = (one-dmg(i))*e_elas_s*(epszx(i))
234 ENDDO
235
236 DO i=1,nel
237 uvar(i,1) = epszz(i)
238 uvar(i,2) = epszx(i)
239 uvar(i,3) = epsyz(i)
240 uvar(i,4)= epsmmax(i)
241
242 uvar(i,6)= signzz(i)
243 uvar(i,7)= signzx(i)
244 uvar(i,8)= signyz(i)
245 uvar(i,9)= delta0m(i)
246 uvar(i,10)= beta(i)
247 uvar(i,12) = deltafmax(i)
248 ENDDO
249c-----------------------------------------------------
250c-----------------------------------------------------
251c omega = sqrt(2k/2*dmels), dt=2/omega, 2*dmels=dt**2 * 2k / 4
252 IF (idtmins==2 .AND. jsms/=0) THEN
253 dtb = (dtmins/dtfacs)**2
254 DO i=1,nel
255 dmels(i)=max(dmels(i),half*dtb*stf*area(i)*off(i))
256 ENDDO
257 END IF
258 stifm(1:nel) = stifm(1:nel) + stf*area(1:nel)*off(1:nel)
259c-----------------------------------------------------
260 IF (nindxf > 0) THEN
261 DO ii=1,nindxf
262 i = indxf(ii)
263#include "lockon.inc"
264 WRITE(iout ,1000) ngl(i),ipg,epsm(i)
265 WRITE(istdo,1100) ngl(i),ipg,epsm(i),time
266#include "lockoff.inc"
267 END DO
268 END IF
269c-----------------------------------------------------
270 1000 FORMAT(5x,'FAILURE COHESIVE ELEMENT ',i10,
271 . ' INTEGRATION POINT',i2,', MIXED MODE STRAIN=',1pe16.9)
272 1100 FORMAT(5x,'FAILURE COHESIVE ELEMENT ',i10,
273 . ' INTEGRATION POINT',i2,', MIXED MODE STRAIN=',1pe16.9,
274 . ' AT TIME ',1pe16.9)
275c-----------
276 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
area
subroutine area(d1, x, x2, y, y2, eint, stif0)
Definition hm_read_prop32.F:567
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21