OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
ruser36.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!|| ruser36 ../engine/source/elements/spring/ruser36.F
25!||--- called by ------------------------------------------------------
26!|| rforc3 ../engine/source/elements/spring/rforc3.F
27!||--- calls -----------------------------------------------------
28!|| get_u_func ../engine/source/user_interface/ufunc.F
29!|| get_u_mid ../engine/source/user_interface/upidmid.F
30!|| get_u_mnu ../engine/source/user_interface/upidmid.F
31!|| get_u_pid ../engine/source/user_interface/upidmid.F
32!|| get_u_pnu ../engine/source/user_interface/upidmid.F
33!||====================================================================
34 SUBROUTINE ruser36(
35 1 NEL ,IPROP ,UVAR ,NUVAR ,FR_WAVE,
36 2 FX ,FY ,FZ ,XMOM ,YMOM ,
37 3 ZMOM ,E ,OFF ,STIFM ,STIFR ,
38 4 VISCM ,VISCR ,MASS ,XINER ,DT ,
39 5 XL ,VX ,RY1 ,RZ1 ,RX ,
40 6 RY2 ,RZ2 )
41C-------------------------------------------------------------------------
42C I m p l i c i t T y p e s
43C-----------------------------------------------
44#include "implicit_f.inc"
45#include "mvsiz_p.inc"
46C----------------------------------------------------------
47C d u m m y a r g u m e n t s a n d f u n c t i o n
48C----------------------------------------------------------
49 INTEGER NEL,NUVAR,IPROP,
50 . GET_U_PNU,GET_U_PID,GET_U_MID,GET_U_MNU,
51 . KFUNC,KMAT,KPROP
52 MY_REAL
53 . UVAR(NUVAR,*),DT ,
54 . FX(*), FY(*), FZ(*), E(*), VX(*),MASS(*) ,XINER(*),
55 . ry1(*), rz1(*), off(*), xmom(*), ymom(*),
56 . zmom(*), rx(*), ry2(*), rz2(*),xl(*),
57 . stifm(*) ,stifr(*) , viscm(*) ,viscr(*) ,fr_wave(*) ,
58 . get_u_mat, get_u_geo, get_u_func
59 EXTERNAL get_u_mnu,get_u_pnu,get_u_mid,get_u_pid,
60 . get_u_mat,get_u_geo, get_u_func
61 parameter(kfunc=29)
62 parameter(kmat=31)
63 parameter(kprop=33)
64C-----------------------------------------------
65C L O C A L V A R I A B L E S
66C-----------------------------------------------
67 INTEGER I,KCST,N0,
68 . imat1,iprop1,iutyp1,
69 . imat2,iprop2,iutyp2
70 INTEGER J,N,NINDX,NMAX,INDEX(MVSIZ),IFUNC
71 MY_REAL
72 . FAC,RHO,AREA,IXX,IYY,IZZ,IMYZ,YOUNG,G,
73 . area1,ixx1,iyy1,izz1,rho1,young1,g1,
74 . area2,ixx2,iyy2,izz2,rho2,young2,g2,
75 . k11,k22,k26,k33,k35,k44,k55,k5b,k66,k6c
76 my_real r_y1,r_y2,r_z1,r_z2,r_y,r_z,r_x,facx,facy,facz,
77 . ay1,ay2,az1,az2,ay,az,by1,by2,bz1,bz2,by,bz,g3,
78 . cx1,cx2,cx,area_y(mvsiz),area_z(mvsiz),
79 . yld(mvsiz),dre,drg,dpla_i,dpla_j(mvsiz),dfe,dfg,
80 . err,f,df,yld_i,c1,fn,fm,dfn,dfm,y1,y2,
81 . pnx,pnx2,pny,pny2,pnz,pnz2,pmx,pmx2,pmy,pmy2,pmz,pmz2,
82 . nx2,ny2,nz2,mx2,my2,mz2,tempy,tempz,xl3,ktran,krot,
83 . nx(mvsiz),ny(mvsiz),nz(mvsiz),youngm(mvsiz),gm(mvsiz),
84 . mx(mvsiz),my(mvsiz),mz(mvsiz),svm(mvsiz),h(mvsiz),
85 . epsvxx(mvsiz),epsvxy(mvsiz),epsvxz(mvsiz),
86 . epscbx(mvsiz),epscby(mvsiz),epscbz(mvsiz),
87 . signxx(mvsiz),signxy(mvsiz),signxz(mvsiz),
88 . momnxx(mvsiz),momnyy(mvsiz),momnzz(mvsiz),fxap,
89 . gama(mvsiz),gamac(mvsiz),xl2(mvsiz),devol,fxav,
90 . pc1,dc1,pr1,ps1,dc2,pr2,ps2,dc,pr,ps,ueq,dtemp,
91 . sig01,sig02, sig0,hpla,h1,h2,m1,m2,m,dsigy,cc1,
92 . fac1,fac2
93 DATA nmax/10/
94C-----------------------------------------------
95 kcst = 0
96C
97 iprop1 = get_u_pnu(1,iprop,kprop)
98 area1 = get_u_geo(2,iprop1)
99 ixx1 = get_u_geo(3,iprop1)
100 iyy1 = get_u_geo(4,iprop1)
101 izz1 = get_u_geo(5,iprop1)
102 r_y1= get_u_geo(6,iprop1)
103 r_z1= get_u_geo(7,iprop1)
104 iprop2 = get_u_pnu(2,iprop,kprop)
105 area2 = get_u_geo(2,iprop2)
106 ixx2 = get_u_geo(3,iprop2)
107 iyy2 = get_u_geo(4,iprop2)
108 izz2 = get_u_geo(5,iprop2)
109 imat1 = get_u_pnu(1,iprop1,kmat)
110 ifunc = get_u_mnu(1,imat1,kfunc)
111 young1 = get_u_mat(7,imat1)
112 g1 = get_u_mat(6,imat1)
113 rho1 = get_u_mat(0,imat1)
114 ay1 = get_u_mat(8,imat1)
115 az1 = get_u_mat(9,imat1)
116 by1 = get_u_mat(10,imat1)
117 bz1 = get_u_mat(11,imat1)
118 cx1 = get_u_mat(12,imat1)
119 dc1 = get_u_mat(13,imat1)
120 pr1 = get_u_mat(14,imat1)
121 ps1 = get_u_mat(15,imat1)
122 sig01 = get_u_mat(16,imat1)
123 h1 = get_u_mat(17,imat1)
124 m1 = get_u_mat(18,imat1)
125 fac1 = get_u_mat(19,imat1)
126 imat2 = get_u_pnu(1,iprop2,kmat)
127 young2 = get_u_mat(7,imat2)
128 g2 = get_u_mat(6,imat2)
129 rho2 = get_u_mat(0,imat2)
130 ay2 = get_u_mat(8,imat2)
131 az2 = get_u_mat(9,imat2)
132 by2 = get_u_mat(10,imat2)
133 bz2 = get_u_mat(11,imat2)
134 cx2 = get_u_mat(12,imat2)
135 dc2 = get_u_mat(13,imat2)
136 pr2 = get_u_mat(14,imat2)
137 ps2 = get_u_mat(15,imat2)
138 sig02 = get_u_mat(16,imat2)
139 h2 = get_u_mat(17,imat2)
140 m2 = get_u_mat(18,imat2)
141 r_y2 = get_u_geo(6,iprop2)
142 r_z2 = get_u_geo(7,iprop2)
143 rho = half*(rho1+rho2)
144 area = half*(area1+area2)
145 fac = rho*area
146 ixx = half*(ixx1+ixx2)
147 iyy = half*(iyy1+iyy2)
148 izz = half*(izz1+izz2)
149 r_y = half*(r_y1+r_y2)
150 r_z = half*(r_z1+r_z2)
151 r_x = half*(r_y+r_z)
152 cx = half*(cx1+cx2)
153 ay = half*(ay1+ay2)
154 az = half*(az1+az2)
155 by = half*(by1+by2)
156 bz = half*(bz1+bz2)
157 imyz = max(iyy,izz)
158 young = half*(young1+young2)
159 g = half*(g1+g2)
160 facx = r_x/max(ixx,em20)
161 facy = r_y/max(iyy,em20)
162 facz = r_z/max(izz,em20)
163 tempy = g*area/max(twelve*young*iyy,em20)
164 tempz = g*area/max(twelve*young*izz,em20)
165 g3 = three*g
166 n0 = 11
167 dc = half*(dc1+dc2)
168 pr = half*(pr1+pr2)
169 ps = half*(ps1+ps2)
170 sig0 = half*(sig01+sig02)
171 hpla = half*(h1+h2)
172 m = half*(m1+m2)
173C
174C---- UVAR(N0,I) -> PLA , N0+1 -> SY0 , N0+2 -> EXX_P , N0+3 -> EYY_P
175C---- UVAR(N0+4,I) -> ENDOMMAGEMENT
176C
177 DO i=1,nel
178 youngm(i) = young*(1-uvar(n0+4,i))
179 gm(i) = g*(1-uvar(n0+4,i))
180C== ON CALCULE D'ABORD LES VITESSES DES DEFORMATIONS GENERALISEES ====
181 dtemp = one/max(xl(i),em20)
182 epsvxx(i) = vx(i)*dtemp
183 epsvxy(i) = -half*(rz1(i) + rz2(i))
184 epsvxz(i) = half*(ry1(i) + ry2(i))
185 epscbx(i) = r_x*rx(i)*dtemp
186 epscby(i) = r_y*(ry2(i) - ry1(i))*dtemp
187 epscbz(i) = r_z*(rz2(i) - rz1(i))*dtemp
188C== ON CALCULE CONTRAINTES GENERALISEES ====
189 xl2(i) = xl(i)*xl(i)
190 area_y(i) = area/(one +tempy*xl2(i))
191 area_z(i) = area/(one +tempz*xl2(i))
192 signxx(i) = fx(i)/max(area,em20)+youngm(i)*epsvxx(i)*dt
193 signxy(i) = fy(i)/max(area_y(i),em20)+gm(i)*epsvxy(i)*dt
194 signxz(i) = fz(i)/max(area_z(i),em20)+gm(i)*epsvxz(i)*dt
195 momnxx(i) = xmom(i)*facx+gm(i)*epscbx(i)*dt
196 momnyy(i) = ymom(i)*facy+youngm(i)*epscby(i)*dt
197 momnzz(i) = zmom(i)*facz+youngm(i)*epscbz(i)*dt
198 ENDDO
199C
200C--- UVAR(11,I)-> DEF PLASTIQUE EQUIVALENTE -
201C--- UVAR(N0,I)-> PLA, N0+1->SY0; N0+2->EXX_P; N0+3->EYY_P;
202C--- UVAR(N0+3,I)-> COURBURE PLASTIQUE Y,Z UVAR(N0+2,I)-> COURBURE PLASTIQUE X -
203C
204 DO i=1,nel
205 IF (ifunc /= 0) THEN
206 y1= fac1*get_u_func(ifunc,uvar(n0,i),y2)
207 y2= fac1*y2
208 yld(i) = max(y1,em20)
209 m=1
210 ELSE
211 y1 = sig0+hpla*(uvar(n0,i))**m
212 y2 = hpla
213 ENDIF
214 yld(i) = max(y1,em20)
215 h(i) = max(y2,zero)
216 uvar(n0+1,i)=min(yld(i),uvar(n0+1,i))
217C
218 pc1 = abs(six*youngm(i)*(uvar(n0+2,i)/uvar(n0+1,i)))
219 c1=one-fourth*exp(-pc1)
220 gamac(i)=cx*three_over_4/max(c1,em20)
221C
222 pc1 = six*youngm(i)*(uvar(n0+3,i)/uvar(n0+1,i))
223 cc1 = (three*pi/sixteen)**2
224 c1=one - (one - cc1)*exp(-pc1)
225 gama(i)=cc1/max(c1,em20)
226C
227 nx(i)=signxx(i)
228 ny(i)=ay*signxy(i)
229 nz(i)=az*signxz(i)
230 mx(i)=gamac(i)*momnxx(i)
231 my(i)=by*gama(i)*momnyy(i)
232 mz(i)=bz*gama(i)*momnzz(i)
233 svm(i)=nx(i)*nx(i)+three*(ny(i)*ny(i)+nz(i)*nz(i)+mx(i)*mx(i))+
234 1 my(i)*my(i)+mz(i)*mz(i)
235 svm(i)=sqrt(svm(i))
236 ENDDO
237 fac2 = g/max(em20,young)
238 DO i=1,nel
239C
240C TIME STEP
241C
242 dtemp = one/max(xl(i),em20)
243 ktran = max(area,fac2*area_y(i),fac2*area_z(i))*dtemp
244 krot = four *max(iyy*dtemp,izz*dtemp)
245 stifm(i) = youngm(i)*ktran
246 stifr(i) = max( gm(i)*ixx*dtemp,youngm(i)*krot)
247 viscm(i) = zero
248 viscr(i) = zero
249 mass(i) = xl(i)*fac
250 xiner(i) = xl(i)*rho*max(ixx,imyz+area*xl2(i)/twelve)
251 ENDDO
252C
253 nindx=0
254 DO i=1,nel
255 IF (svm(i) > yld(i) .AND. off(i) == one) THEN
256 nindx=nindx+1
257 index(nindx)=i
258 ENDIF
259 ENDDO
260C
261 IF (nindx /= 0) THEN
262 DO j=1,nindx
263 i=index(j)
264 dpla_j(i)=(svm(i)-yld(i))/(g3+h(i))
265 ENDDO
266C
267C ATTENTION, LA VERSION DU MODELE DE RESSORT (Q.ZENG/A.COMBESCURE/C.BENOIT)
268C NE CONSIDERE PAS D'ERREUR CUMULEE SUR PLUSIEURS ELEMENTS...
269C
270C --> EXISTE VERSION AVEC TEST MULTI-ELEMENTS (ERR_MULTI/RUSER29.F)
271C
272C
273C
274 n=0
275 err =2*em4
276 DO WHILE (err > em4 .AND. n < nmax)
277 DO j=1,nindx
278 i=index(j)
279 dpla_i=dpla_j(i)
280 IF (ifunc /= 0) THEN
281 yld_i =yld(i)+h(i)*dpla_i
282 ELSEIF (uvar(n0,i) > em20) THEN
283 yld_i =yld(i)+h(i)*m*dpla_i*(uvar(n0,i)+dpla_i)**(m-1)
284 ELSE
285 yld_i =yld(i)
286 ENDIF
287 dre =youngm(i)*dpla_i/yld_i
288 drg =g*dpla_i/yld_i
289 pnx =one/(one+dre)
290 pny =one/(one+three*drg*ay**2)
291 pnz =one/(one+three*drg*az**2)
292 pmx =one/(one+three*drg*gamac(i)**2)
293 pmy =one/(one+dre*(by*gama(i))**2)
294 pmz =one/(one+dre*(bz*gama(i))**2)
295 pnx2=pnx*pnx
296 pny2=pny*pny
297 pnz2=pnz*pnz
298 pmx2=pmx*pmx
299 pmy2=pmy*pmy
300 pmz2=pmz*pmz
301 nx2=nx(i)*nx(i)
302 ny2=ny(i)*ny(i)
303 nz2=nz(i)*nz(i)
304 mx2=mx(i)*mx(i)
305 my2=my(i)*my(i)
306 mz2=mz(i)*mz(i)
307 fn=nx2*pnx2+3.*(ny2*pny2+nz2*pnz2)
308 fm=3.*mx2*pmx2+my2*pmy2+mz2*pmz2
309 f = fn+fm-yld_i*yld_i
310 dfe = nx2*pnx2*pnx+my2*pmy2*pmy*(by*gama(i))**2
311 . +mz2*pmz2*pmz*(bz*gama(i))**2
312 dfg = ny2*pny2*pny*ay**2+nz2*pnz2*pnz*az**2
313 . +mx2*pmx2*pmx*gamac(i)**2
314 IF (ifunc /= 0) THEN
315 dsigy = -h(i)*yld_i*two
316 ELSEIF (uvar(n0,i)>em20) THEN
317 dsigy = -(h(i)*yld_i*m*(dpla_i+uvar(n0,i))**(m-1))*two
318 ELSE
319 dsigy = zero
320 ENDIF
321 df =(-dfe*(youngm(i)-dre*h(i))/yld_i
322 . -nine*dfg*(g-drg*h(i))/yld_i)*two
323 df=df+dsigy
324 err=abs(f/df)
325 dpla_j(i)=max(zero,dpla_i-f/df)
326 ENDDO ! DO J=1,NINDX
327 n=n+1
328 ENDDO ! DO WHILE (ERR > EM4 .AND. N < NMAX)
329C
330C CONTRAINTES PLASTIQUEMENT ADMISSIBLES
331C
332 DO j=1,nindx
333 i=index(j)
334 uvar(n0,i) = uvar(n0,i) + dpla_j(i)
335 dpla_i=dpla_j(i)
336C INCREMENT N+1
337 IF (ifunc /= 0) THEN
338 yld_i =yld(i)+h(i)*dpla_i
339 ELSEIF (uvar(n0,i) > em20) THEN
340 yld_i =yld(i)+h(i)*m*dpla_i*(uvar(n0,i))**(m-1)
341 ELSE
342 yld_i =yld(i)
343 ENDIF
344 c1 = dpla_i/yld_i
345 dre =youngm(i)*c1
346 drg =g*c1
347 pnx =one/(one+dre)
348 pny =one/(one+three*drg*ay**2)
349 pnz =one/(one+three*drg*az**2)
350 pmx =one/(one+three*drg*gamac(i)**2)
351 pmy =one/(one+dre*(by*gama(i))**2)
352 pmz =one/(one+dre*(bz*gama(i))**2)
353C
354 signxx(i) = signxx(i)*pnx
355 signxy(i) = signxy(i)*pny
356 signxz(i) = signxz(i)*pnz
357 momnxx(i) = momnxx(i)*pmx
358 momnyy(i) = momnyy(i)*pmy
359 momnzz(i) = momnzz(i)*pmz
360C
361 uvar(n0+2,i) = uvar(n0+2,i)+three*c1*gamac(i)**2*abs(momnxx(i))
362 ueq = c1*sqrt( signxx(i)**2
363 1 + gama(i)**4
364 1 *(by**4*momnyy(i)**2 + bz**4*momnzz(i)**2) )
365 uvar(n0+3,i) = uvar(n0+3,i) + ueq
366C
367C ENDOMAGEMENT
368C
369 IF (uvar(n0,i) > ps) THEN
370 devol = dc*dpla_i/(pr - ps)
371 ELSE
372 devol = zero
373 ENDIF
374C
375 uvar(n0+4,i) = uvar(n0+4,i)+devol
376C
377 IF (uvar(n0,i) >= pr .OR. uvar(n0+4,i) >= dc) THEN
378 off(i) = zero
379 uvar(n0+4,i) = dc
380 ENDIF
381 ENDDO ! DO J=1,NINDX
382 ENDIF ! IF (NINDX /= 0)
383C
384C TRANSLATE FORCES AND MOMENTS :
385C
386 DO i=1,nel
387 fx(i) = signxx(i)*area*off(i)
388 fy(i) = signxy(i)*area_y(i)*off(i)
389 fz(i) = signxz(i)*area_z(i)*off(i)
390 xmom(i) = momnxx(i)*off(i)/facx
391 ymom(i) = momnyy(i)*off(i)/facy
392 zmom(i) = momnzz(i)*off(i)/facz
393 ENDDO
394C---
395 1110 FORMAT(
396 & 5x,'N. . . . . . . =',i5/,
397 & 5x,'ERR. . . . . . =',e12.4/,
398 & 5x,'DPPLAS . . . . =',e12.4/,
399 & 5x,'F PLAS . . . . =',e12.4/,
400 & 5x,'DF PLAS. . . . =',e12.4///)
401 1120 FORMAT(
402 & 5x,'MODULE YOUNG . . . . . ..=',e12.4/,
403 & 5x,'PPLAS CUMULEE. . . . . . =',e12.4/,
404 & 5x,'PPLAS TORSION . . . . . =',e12.4/,
405 & 5x,'PPLAS FLEXION. . . . . . =',e12.4/,
406 & 5x,'ENDOMMAGEMENT . . . . . =',e12.4///)
407 1130 FORMAT(
408 & 5x,'PPLAS CUMULEE. . . . . . =',e12.4/,
409 & 5x,'DEVOL . . . . . =',e12.4/,
410 & 5x,'PS . . . . . ..=',e12.4/,
411 & 5x,'PR . . . . . =',e12.4///)
412 1140 FORMAT(
413 & 3x,e16.10,
414 & 3x,f14.4)
415C---
416 RETURN
417 END
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
ruser36
subroutine ruser36(nel, iprop, uvar, nuvar, fr_wave, fx, fy, fz, xmom, ymom, zmom, e, off, stifm, stifr, viscm, viscr, mass, xiner, dt, xl, vx, ry1, rz1, rx, ry2, rz2)
Definition ruser36.F:41