OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps158c.F File Reference
#include "implicit_f.inc"
#include "com04_c.inc"
#include "tabsiz_c.inc"
#include "vectorize.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine sigeps158c (nel, nuparam, nuvar, nfunc, ifunc, npf, tf, time, timestep, uparam, area, thkly, soundsp, viscmax, uvar, depsxx, depsyy, depsxy, depsyz, depszx, epsxx, epsyy, epsxy, epsyz, epszx, sigoxx, sigoyy, sigoxy, sigoyz, sigozx, signxx, signyy, signxy, signyz, signzx, sigvxx, sigvyy, sigvxy, tan_phi, offgg, rho0, etse, shf, aldt, nsensor, sensor_tab, niparam, iparam)

Function/Subroutine Documentation

◆ sigeps158c()

subroutine sigeps158c ( integer, intent(in) nel,
integer, intent(in) nuparam,
integer, intent(in) nuvar,
integer, intent(in) nfunc,
integer, dimension(nfunc), intent(in) ifunc,
integer, dimension(snpc), intent(in) npf,
intent(in) tf,
intent(in) time,
intent(in) timestep,
intent(in) uparam,
intent(in) area,
intent(in) thkly,
intent(out) soundsp,
intent(out) viscmax,
intent(inout) uvar,
intent(in) depsxx,
intent(in) depsyy,
intent(in) depsxy,
intent(in) depsyz,
intent(in) depszx,
intent(in) epsxx,
intent(in) epsyy,
intent(in) epsxy,
intent(in) epsyz,
intent(in) epszx,
intent(in) sigoxx,
intent(in) sigoyy,
intent(in) sigoxy,
intent(in) sigoyz,
intent(in) sigozx,
intent(out) signxx,
intent(out) signyy,
intent(out) signxy,
intent(out) signyz,
intent(out) signzx,
intent(out) sigvxx,
intent(out) sigvyy,
intent(out) sigvxy,
intent(out) tan_phi,
intent(in) offgg,
intent(in) rho0,
intent(out) etse,
intent(in) shf,
intent(in) aldt,
integer, intent(in) nsensor,
type (sensor_str_), dimension(nsensor), intent(in) sensor_tab,
integer, intent(in) niparam,
integer, dimension(niparam), intent(in) iparam )

Definition at line 32 of file sigeps158c.F.

43C-----------------------------------------------
44C M o d u l e s
45C-----------------------------------------------
46 USE sensor_mod
47C-----------------------------------------------
48C I m p l i c i t T y p e s
49C-----------------------------------------------
50#include "implicit_f.inc"
51C---------+---------+---+---+--------------------------------------------
52C VAR | SIZE |TYP| RW| DEFINITION
53C---------+---------+---+---+--------------------------------------------
54C NEL | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL
55C NUPARAM | 1 | I | R | SIZE OF THE USER PARAMETER ARRAY
56C NUVAR | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES
57C---------+---------+---+---+--------------------------------------------
58C NFUNC | 1 | I | R | NUMBER FUNCTION USED FOR THIS USER LAW
59C IFUNC | NFUNC | I | R | FUNCTION INDEX
60C NPF | * | I | R | FUNCTION ARRAY
61C TF | * | F | R | FUNCTION ARRAY
62C---------+---------+---+---+--------------------------------------------
63C TIMESTEP| 1 | F | R | CURRENT TIME STEP
64C UPARAM | NUPARAM | F | R | USER MATERIAL PARAMETER ARRAY
65C RHO0 | NEL | F | R | INITIAL DENSITY
66C AREA | NEL | F | R | AREA
67C THKLY | NEL | F | R | LAYER THICKNESS
68C ... | | | |
69C DEPSXX | NEL | F | R | STRAIN INCREMENT XX
70C DEPSYY | NEL | F | R | STRAIN INCREMENT YY
71C ... | | | |
72C EPSXX | NEL | F | R | STRAIN XX TRUE
73C EPSYY | NEL | F | R | STRAIN YY TRUE
74C ... | | | |
75C SIGOXX | NEL | F | R | OLD ELASTO PLASTIC STRESS XX
76C SIGOYY | NEL | F | R | OLD ELASTO PLASTIC STRESS YY
77C ... | | | |
78C---------+---------+---+---+--------------------------------------------
79C SIGNXX | NEL | F | W | NEW ELASTO PLASTIC STRESS XX
80C SIGNYY | NEL | F | W | NEW ELASTO PLASTIC STRESS YY
81C ... | | | |
82C SIGVXX | NEL | F | W | VISCOUS STRESS XX
83C SIGVYY | NEL | F | W | VISCOUS STRESS YY
84C ... | | | |
85C SOUNDSP | NEL | F | W | SOUND SPEED (NEEDED FOR TIME STEP)
86C VISCMAX | NEL | F | W | MAXIMUM DAMPING MODULUS(NEEDED FOR TIME STEP)
87C---------+---------+---+---+--------------------------------------------
88C UVAR |NEL*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
89C OFFG | NEL | F |R/W| DELETED ELEMENT FLAG (=1. ON, =0. OFF)
90C---------+---------+---+---+--------------------------------------------
91C C o m m o n B l o c k s
92C-----------------------------------------------
93#include "com04_c.inc"
94#include "tabsiz_c.inc"
95C-----------------------------------------------
96C I N P U T A r g u m e n t s
97C-----------------------------------------------
98 INTEGER ,INTENT(IN) :: NEL
99 INTEGER ,INTENT(IN) :: NFUNC
100 INTEGER ,INTENT(IN) :: NUPARAM
101 INTEGER ,INTENT(IN) :: NIPARAM
102 INTEGER ,INTENT(IN) :: NUVAR
103 INTEGER ,INTENT(IN) :: NSENSOR
104 my_real ,INTENT(IN) :: time,timestep
105 INTEGER ,DIMENSION(NFUNC) ,INTENT(IN) :: IFUNC
106 INTEGER ,DIMENSION(SNPC) ,INTENT(IN) :: NPF
107 my_real ,DIMENSION(STF) ,INTENT(IN) :: tf
108 INTEGER ,DIMENSION(NIPARAM) ,INTENT(IN) :: IPARAM
109 my_real ,DIMENSION(NUPARAM) ,INTENT(IN) :: uparam
110 my_real ,DIMENSION(NEL) ,INTENT(IN) :: area,rho0,thkly,offgg,
111 . depsxx,depsyy,depsxy,depsyz,depszx,epsxx,epsyy,epsxy,epsyz,epszx,
112 . sigoxx,sigoyy,sigoxy,sigoyz,sigozx,shf,aldt
113 TYPE (SENSOR_STR_) ,DIMENSION(NSENSOR) ,INTENT(IN) :: SENSOR_TAB
114C-----------------------------------------------
115C O U T P U T A r g u m e n t s
116C-----------------------------------------------
117 my_real ,DIMENSION(NEL) ,INTENT(OUT) :: tan_phi,etse,
118 . signxx,signyy,signxy,signyz,signzx,
119 . sigvxx,sigvyy,sigvxy,soundsp,viscmax
120C-----------------------------------------------
121C I N P U T O U T P U T A r g u m e n t s
122C-----------------------------------------------
123 my_real ,DIMENSION(NEL,NUVAR) ,INTENT(INOUT) :: uvar
124C-----------------------------------------------
125C L o c a l V a r i a b l e s
126C-----------------------------------------------
127 INTEGER :: I,ITER,NITER,ISENS,IFUNCC,IFUNCT,IFUNCS,IFUNF1,IFUNF2
128 INTEGER ,DIMENSION(NEL) :: IADF1,IADF2,IAD1,IAD2,IAD3,ILEN1,ILEN2,ILEN3,
129 . ILENF1,ILENF2,IPOS,IPOS1,IPOS2
130 my_real :: kmax,gmax,gfrot,gsh,stiff,kflex,kflex1,kflex2,
131 . mass,dyc,dyt,dc0,dt0,h0,hc0,ht0,visce,viscg,lc0,lt0,hdc,hdt,
132 . tfrot,sigg,ec2,et2,sinn,tan2,damp,v1,v2,dtinv,etc,ett,trace,
133 . zerostress,dsig,lmin,tstart,phi1,phi2,j11,j12,j21,j22,det
134 my_real ,DIMENSION(3) :: yfac
135 my_real ,DIMENSION(NEL) :: sigc,sigt,ec,et,lc,lt,fn,angle,kg,gxy,
136 . dtang,tfold,cvisc,cvist,flexc,flext,flexf,hc,ht,dc,dt,dcc,dtt,
137 . fc,ft,fpc,fpt,dflxc,dflxt,epsf,epsfc,epsft,xc,xt,beta,ph01,ph02
138C-----------------------------------------------
139C S t a t e V a r i a b l e s (U V A R)
140C-----------------------------------------------
141c UVAR(1) = SIGNXX + SIGVXX ! total stress in 1st direction (elastic + viscous)
142c UVAR(2) = SIGNYY + SIGVYY ! total stress in 2nd direction (elastic + viscous)
143c UVAR(3) = SIGNXY + SIGVXY ! total shear stress (elastic + friction)
144c UVAR(4) = EC ! total element elongation in 1st direction
145c UVAR(5) = ET ! total element elongation in 2nd direction
146c UVAR(6) = DEPSXY ! tan(alpha) of shear angle
147c UVAR(7) = YC ! total flex element elongation in 1st direction
148c UVAR(8) = YT ! total flex element elongation in 2nd direction
149c UVAR(9) = SIGVXY ! Friction stress (shear)
150c UVAR(10) = ! not used
151c UVAR(11) = SIGNXX ! elastic stress in 1st direction (used for zerostress relaxation)
152c UVAR(12) = SIGNYY ! elastic stress in 2nd direction (used for zerostress relaxation)
153c UVAR(13) = SIGNXY ! elastic stress in shear (used for zerostress relaxation)
154c UVAR(14) = ALDT ! initial characteristic length of the element
155c UVAR(15) = DC ! total fiber elongation in 1st direction
156c UVAR(16) = DT ! total fiber elongation in 2nd direction
157C======================================================================|
158C--- Initialisations
159C----------------------------------------------------------------------
160 isens = iparam(1)
161!
162 lc0 = one
163 lt0 = one
164 dc0 = uparam(1)
165 dt0 = uparam(2)
166 hc0 = uparam(3)
167 ht0 = uparam(4)
168 kflex = uparam(5)
169 kflex1 = uparam(6) ! stiffness factor for log function in dir1
170 kflex2 = uparam(7) ! stiffness factor for log function in dir1
171 zerostress = uparam(8)
172 kmax = uparam(10)
173 gmax = uparam(11)
174 yfac(1) = uparam(12)
175 yfac(2) = uparam(13)
176 yfac(3) = uparam(14)
177
178 ifuncc = ifunc(1)
179 ifunct = ifunc(2)
180 ifuncs = ifunc(3)
181 ifunf1 = ifunc(4)
182 ifunf2 = ifunc(5)
183
184 gfrot = gmax
185 gsh = shf(1)
186 visce = em02 ! fiber visc damping coefficient
187 viscg = em02 ! shear damping coefficient
188c
189 IF (zerostress > zero .and. isens > 0) THEN
190 tstart = sensor_tab(isens)%TSTART
191 ELSE
192 tstart = zero
193 ENDIF
194 dtinv = timestep/max(timestep**2,em20)
195c
196 h0 = hc0 + ht0
197 niter = 10
198C-----------------------------------------------------------
199#include "vectorize.inc"
200 DO i = 1,nel
201 tfold(i) = uvar(i,9)
202 mass = rho0(i)*area(i)*thkly(i)*fourth ! 2m = rho*V/4
203 cvisc(i) = sqrt(mass*kflex1)*dtinv*third
204 cvist(i) = sqrt(mass*kflex2)*dtinv*third
205 ENDDO
206
207c--- strain integration -> total engineering strains
208#include "vectorize.inc"
209 DO i = 1,nel
210 etc = uvar(i,4) + depsxx(i)
211 ett = uvar(i,5) + depsyy(i)
212 uvar(i,4) = etc
213 uvar(i,5) = ett
214 ec(i) = exp(etc) - one ! eng strain dir 1
215 et(i) = exp(ett) - one ! eng strain dir 2
216 lc(i) = lc0 * (one + ec(i)) ! element length dir 1
217 lt(i) = lt0 * (one + et(i)) ! element length dir 2
218 ENDDO
219
220 DO i = 1,nel
221 beta(i) = one ! increment scale factor for Newton iterations
222 ph01(i) = ep10
223 ph02(i) = ep10
224 ENDDO
225c
226#include "vectorize.inc"
227 DO i = 1,nel
228 ipos(i) = 1
229 iad1(i) = npf(ifuncc) / 2 + 1
230 iad2(i) = npf(ifunct) / 2 + 1
231 iad3(i) = npf(ifuncs) / 2 + 1
232 ilen1(i) = npf(ifuncc+1) / 2 - iad1(i) - ipos(i)
233 ilen2(i) = npf(ifunct+1) / 2 - iad2(i) - ipos(i)
234 ilen3(i) = npf(ifuncs+1) / 2 - iad3(i) - ipos(i)
235 END DO
236 IF (ifunf1 > 0) THEN
237 DO i = 1,nel
238 iadf1(i) = npf(ifunf1) / 2 + 1
239 ilenf1(i) = npf(ifunf1+1) / 2 - iadf1(i) - ipos(i)
240 END DO
241 END IF
242 IF (ifunf2 > 0) THEN
243 DO i = 1,nel
244 iadf2(i) = npf(ifunf2) / 2 + 1
245 ilenf2(i) = npf(ifunf2+1) / 2 - iadf2(i) - ipos(i)
246 END DO
247 END IF
248c
249c--- resolve nonlinear fiber equilibrium equations using Newton iterations
250c
251 epsfc(1:nel) = uvar(1:nel,7) ! eng strain flex dir 1
252 epsft(1:nel) = uvar(1:nel,8) ! eng strain flex dir 2
253
254 DO iter = 1,niter
255
256 epsf(1:nel) =(hc0 * epsfc(1:nel) + ht0 * epsft(1:nel)) / h0 ! eng strain coupling flex spring
257 flexf(1:nel) = kflex * epsf(1:nel)
258 DO i = 1,nel
259 xc(i) = epsfc(i)
260 xt(i) = epsft(i)
261 hc(i) = hc0 * (epsfc(i) + one) ! length of flex spring dir 1
262 ht(i) = ht0 * (epsft(i) + one) ! length of flex spring dir 2
263 dc(i) = sqrt(lc(i)**2 + hc(i)**2) ! fiber length in dir 1
264 dt(i) = sqrt(lt(i)**2 + ht(i)**2) ! fiber length in dir 2
265 dcc(i) = dc(i) - dc0
266 dtt(i) = dt(i) - dt0
267 END DO
268 ipos1(1:nel) = 1
269 ipos2(1:nel) = 1
270
271 CALL vinter2(tf,iad1,ipos1,ilen1,nel,dcc,fpc,fc) ! fiber force dir 1
272 CALL vinter2(tf,iad2,ipos2,ilen2,nel,dtt,fpt,ft) ! fiber force dir 2
273c
274 ! flex force dir 1
275 IF (ifunf1 > 0) THEN
276 ipos(1:nel) = 1
277 CALL vinter2(tf,iadf1,ipos,ilenf1,nel,epsfc,dflxc,flexc)
278 flexc(1:nel) = flexc(1:nel) * kflex1
279 dflxc(1:nel) = dflxc(1:nel) * kflex1
280 ELSE
281 DO i = 1,nel
282 flexc(i) = kflex1 * log(epsfc(i) + one)
283 dflxc(i) = kflex1 / (epsfc(i) + one)
284 END DO
285 END IF
286c flex force dir 2
287 IF (ifunf2 > 0) THEN
288 ipos(1:nel) = 1
289 CALL vinter2(tf,iadf2,ipos,ilenf2,nel,epsft,dflxt,flext)
290 flext(1:nel) = flext(1:nel) * kflex2
291 dflxt(1:nel) = dflxt(1:nel) * kflex2
292 ELSE
293 DO i = 1,nel
294 flext(i) = kflex2 * log(epsft(i) + one)
295 dflxt(i) = kflex2 / (epsft(i) + one)
296 END DO
297 END IF
298c
299#include "vectorize.inc"
300 DO i = 1,nel
301 dyc = epsfc(i) - uvar(i,7)
302 dyt = epsft(i) - uvar(i,8)
303 hdc = hc(i) / dc(i)
304 hdt = ht(i) / dt(i)
305 phi1 = fc(i) * hdc + flexc(i) + flexf(i) + cvisc(i)*dyc
306 phi2 = ft(i) * hdt + flext(i) + flexf(i) + cvist(i)*dyt
307 j12 = kflex * ht0 / h0
308 j21 = kflex * hc0 / h0
309 j11 = j12 + fpc(i)*hdc*hc0 + dflxc(i) + cvisc(i)
310 j22 = j21 + fpt(i)*hdt*ht0 + dflxt(i) + cvist(i)
311 det = j11 * j22 - j12 * j21
312
313 epsfc(i) = epsfc(i) - beta(i) * (j22 * phi1 - j12 * phi2) / det
314 epsft(i) = epsft(i) + beta(i) * (j12 * phi1 - j11 * phi2) / det
315
316 epsfc(i) = max(epsfc(i), em04 - one)
317 epsft(i) = max(epsft(i), em04 - one)
318
319 IF (abs(phi1) > ph01(i) .and. abs(phi2) > ph02(i)) THEN
320 epsfc(i) = xc(i)
321 epsft(i) = xt(i)
322 beta(i) = beta(i) * half
323 beta(i) = max(beta(i), em02)
324 END IF
325 ph01(i) = abs(phi1)
326 ph02(i) = abs(phi2)
327 END DO
328c
329 END DO ! ITER
330c
331#include "vectorize.inc"
332 DO i = 1,nel
333 sigc(i) = fc(i) * lc(i) / dc(i)
334 sigt(i) = ft(i) * lt(i) / dt(i)
335 uvar(i,7) = epsfc(i)
336 uvar(i,8) = epsft(i)
337 uvar(i,1) = sigc(i)
338 uvar(i,2) = sigt(i)
339 uvar(i,15) = dc(i)
340 uvar(i,16) = dt(i)
341 fn(i) = flexf(i) ! normal force for friction in shear
342 ENDDO !i,nel
343c
344C------------------------------------------------------------------
345C Trace = Eps1 + Eps2 (Trace of true principal strain tensor)
346C Trace = ec_true + ec2_true
347C ec2_eng = exp(Tr) / (ec_eng + 1) - 1
348C rfac = 2*Nc / (ec2_eng+1)
349C ec2_eng+1 > 0
350C------------------------------------------------------------------
351#include "vectorize.inc"
352 DO i = 1,nel
353 trace = exp(epsxx(i) + epsyy(i))! =exp(tr)
354 ec2 = max(trace / (ec(i) + one), em6)
355 et2 = max(trace / (et(i) + one), em6)
356C--- true stress membrane
357 signxx(i) = sigc(i) / ec2
358 signyy(i) = sigt(i) / et2
359 ENDDO !i,nel
360c------------------------------------------------------------------
361c--- SHEAR
362c------------------------------------------------------------------
363#include "vectorize.inc"
364 DO i = 1,nel
365 tan_phi(i)= depsxy(i)
366 angle(i) = atan(tan_phi(i))*hundred80/pi
367 ipos(i) = 1
368 ENDDO
369c
370 CALL vinter2(tf,iad3,ipos,ilen3,nel,angle,gxy,signxy) ! shear stress =f(angle)
371c
372#include "vectorize.inc"
373 DO i = 1,nel
374 tan2 = tan_phi(i)**2
375 kg(i) = gxy(i) * tan2 * yfac(3)
376c--- fiber visc damping
377 damp = sqrt(rho0(i)*area(i)*thkly(i)*half)
378 v1 = visce*damp*sqrt(kmax)
379 v2 = visce*damp*sqrt(kmax)
380 sigvxx(i) = dtinv*(depsxx(i))*v1
381 sigvyy(i) = dtinv*(depsyy(i))*v2
382c--- friction in shear
383 IF (fn(i) > zero) THEN
384 tfrot = two_third*viscg*fn(i)*(hc0+ht0)/(lc(i)+lt(i))
385 dtang(i) = depsxy(i) - tan_phi(i)
386 sigg = tfold(i) + gfrot*dtang(i)
387 IF (abs(sigg) > tfrot) THEN
388 sigvxy(i) = sign(tfrot,sigg)
389 ELSE
390 sigvxy(i) = sigg
391 ENDIF
392 ENDIF
393C---
394 sinn = tan_phi(i) / sqrt(one + tan2)
395 stiff = kmax*(one+sinn) + gmax
396 lmin = min(dc(i)/dc0,dt(i)/dt0)*uvar(i,14)
397 soundsp(i) = sqrt(stiff/(rho0(i)))*aldt(i)/lmin
398 viscmax(i) = max(v1,v2)
399 etse(i) = one
400 ENDDO
401C---
402#include "vectorize.inc"
403 DO i = 1,nel
404 uvar(i,3) = signxy(i) + sigvxy(i)
405 tan_phi(i)= depsxy(i)
406 uvar(i,6) = depsxy(i)
407 uvar(i,9) = sigvxy(i)
408c
409 signyz(i) = sigoyz(i) + gsh * depsyz(i)
410 signzx(i) = sigozx(i) + gsh * depszx(i)
411 ENDDO
412C-----------------------------------------------------------
413C REF-STATE ZEROSTRESS OPTION
414C-----------------------------------------------------------
415 IF (zerostress > zero)THEN
416 IF (time <= tstart) THEN
417#include "vectorize.inc"
418 DO i = 1,nel
419 uvar(i,11) = signxx(i)
420 uvar(i,12) = signyy(i)
421 uvar(i,13) = signxy(i)
422 signxx(i) = zero
423 signyy(i) = zero
424 signxy(i) = zero
425 ENDDO
426 ELSE
427#include "vectorize.inc"
428 DO i = 1,nel
429 dsig = signxx(i) - sigoxx(i) - uvar(i,11)
430 IF((uvar(i,11) > zero).AND.(dsig < zero))THEN
431 uvar(i,11) = max(zero,uvar(i,11)+zerostress*dsig)
432 ELSEIF((uvar(i,11) < zero).AND.(dsig > zero))THEN
433 uvar(i,11) = min(zero,uvar(i,11)+zerostress*dsig)
434 ENDIF
435 dsig = signyy(i) - sigoyy(i) - uvar(i,12)
436 IF((uvar(i,12) > zero).AND.(dsig < zero))THEN
437 uvar(i,12) = max(zero,uvar(i,12)+zerostress*dsig)
438 ELSEIF((uvar(i,12) < zero).AND.(dsig > zero))THEN
439 uvar(i,12) = min(zero,uvar(i,12)+zerostress*dsig)
440 ENDIF
441 dsig = signxy(i) - sigoxy(i) - uvar(i,13)
442 IF((uvar(i,13) > zero).AND.(dsig < zero))THEN
443 uvar(i,13) = max(zero,uvar(i,13)+zerostress*dsig)
444 ELSEIF((uvar(i,13) < zero).AND.(dsig > zero))THEN
445 uvar(i,13) = min(zero,uvar(i,13)+zerostress*dsig)
446 ENDIF
447 signxx(i) = signxx(i) - uvar(i,11)
448 signyy(i) = signyy(i) - uvar(i,12)
449 signxy(i) = signxy(i) - uvar(i,13)
450 ENDDO
451 ENDIF
452 ENDIF
453c
454 DO i = 1,nel
455 soundsp(i) = sqrt(kmax*two/(rho0(i)))
456 ENDDO
457C-----------
458 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
damping
subroutine damping(nodft, nodlt, v, vr, a, ar, damp, ms, in, igrnod, dim, itask, weight, tagslv_rby, wfext)
Definition damping.F:882
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
vinter2
subroutine vinter2(tf, iad, ipos, ilen, nel0, x, dydx, y)
Definition vinter.F:144