OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
main_beam18.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!|| main_beam18 ../engine/source/elements/beam/main_beam18.F
25!||--- called by ------------------------------------------------------
26!|| pforc3 ../engine/source/elements/beam/pforc3.F
27!||--- calls -----------------------------------------------------
28!|| fail_beam18 ../engine/source/elements/beam/fail_beam18.F
29!|| m2lawpi ../engine/source/materials/mat/mat002/m2lawpi.F
30!|| mulaw_ib ../engine/source/elements/beam/mulaw_ib.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 main_beam18(ELBUF_STR,
36 1 NEL ,NPT ,MTN ,IMAT ,
37 2 PID ,NGL ,PM ,IPM ,
38 3 GEO ,OFF ,FOR ,MOM ,
39 4 EINT ,AL ,EPSD ,BUFMAT ,NPF ,
40 5 TF ,EXX ,EXY ,EXZ ,KXX ,
41 6 KYY ,KZZ ,F1 ,F2 ,F3 ,
42 7 M1 ,M2 ,M3 ,JTHE ,TEMPEL ,
43 8 IFAIL ,SBUFMAT ,SNPC ,STF ,NUMMAT ,
44 9 NUMGEO ,IOUT ,ISTDO ,NPROPMI ,NPROPM ,
45 A NPROPG ,TIME ,DTIME ,IDEL7NOK,ISIGI ,
46 B IMCONV ,ISMSTR ,MAT_PARAM ,NTABLE,TABLE )
47C-----------------------------------------------
48C M o d u l e s
49C-----------------------------------------------
50 USE mat_elem_mod
51 USE elbufdef_mod
52C-----------------------------------------------
53C I m p l i c i t T y p e s
54C-----------------------------------------------
55#include "implicit_f.inc"
56#include "comlock.inc"
57C-----------------------------------------------
58C D u m m y A r g u m e n t s
59C-----------------------------------------------
60 INTEGER ,INTENT(IN) :: IMAT
61 INTEGER ,INTENT(IN) :: NEL,MTN,NPT,JTHE,IFAIL
62 INTEGER ,INTENT(IN) :: SBUFMAT
63 INTEGER ,INTENT(IN) :: SNPC
64 INTEGER ,INTENT(IN) :: STF
65 INTEGER ,INTENT(IN) :: NUMMAT
66 INTEGER ,INTENT(IN) :: NUMGEO
67 INTEGER ,INTENT(IN) :: NPROPMI
68 INTEGER ,INTENT(IN) :: NPROPM
69 INTEGER ,INTENT(IN) :: NPROPG
70 INTEGER ,INTENT(IN) :: IOUT
71 INTEGER ,INTENT(IN) :: ISTDO
72 INTEGER ,INTENT(IN) :: ISIGI
73 INTEGER ,INTENT(IN) :: IMCONV
74 INTEGER ,INTENT(IN) :: ISMSTR
75 INTEGER ,INTENT(IN) :: NTABLE
76 TYPE(ttable), DIMENSION(NTABLE), INTENT(INOUT) :: TABLE ! TABLE DATA
77 INTEGER ,DIMENSION(NEL) ,INTENT(IN) :: PID
78 INTEGER ,DIMENSION(NEL) ,INTENT(IN) :: NGL
79 INTEGER ,DIMENSION(SNPC) ,INTENT(IN) :: NPF
80 INTEGER IPM(NPROPMI,NUMMAT)
81 INTEGER ,INTENT(INOUT) :: IDEL7NOK
82 my_real ,INTENT(IN) :: time
83 my_real ,INTENT(IN) :: dtime
84 my_real ,INTENT(IN) :: pm(npropm,nummat)
85 my_real ,INTENT(IN) :: geo(npropg,numgeo)
86 my_real ,DIMENSION(NEL) ,INTENT(IN) :: al
87 my_real ,DIMENSION(NEL) ,INTENT(OUT) :: epsd
88 my_real ,DIMENSION(SBUFMAT) ,INTENT(IN) :: bufmat
89 my_real ,DIMENSION(NEL) ,INTENT(INOUT) :: tempel
90 my_real ,DIMENSION(NEL) ,INTENT(INOUT) :: off
91 my_real ,DIMENSION(NEL,2),INTENT(INOUT) :: eint
92 my_real ,DIMENSION(NEL,3),INTENT(INOUT) :: for,mom
93 my_real ,DIMENSION(NEL) ,INTENT(INOUT) :: exx,exy,exz
94 my_real ,DIMENSION(NEL) ,INTENT(INOUT) :: kxx,kyy,kzz
95 my_real ,DIMENSION(NEL) ,INTENT(INOUT) :: f1,f2,f3
96 my_real ,DIMENSION(NEL) ,INTENT(INOUT) :: m1,m2,m3
97 my_real ,DIMENSION(STF) ,INTENT(IN) :: tf
98C
99 TYPE (ELBUF_STRUCT_) ,INTENT(INOUT) :: ELBUF_STR
100 TYPE (MATPARAM_STRUCT_) ,INTENT(IN) :: MAT_PARAM
101C-----------------------------------------------
102C L o c a l V a r i a b l e s
103C-----------------------------------------------
104 INTEGER :: I,J,IPT,IPID,IPLA,ISRATE,IPY,IPZ,IPA,II(3)
105 my_real :: asrate,epsdi,fact,area,
106 . damage_loc,dfxx,dfxy,dfxz,
107 . signxx,signxy,signxz,sigoxx,sigoxy,sigoxz
108 my_real ,DIMENSION(NEL ):: degmb,degsh,degfx,ypt,zpt,apt,off_old
109 my_real ,DIMENSION(:,:) ,ALLOCATABLE :: dpla,sigy,vm,vm0
110C=======================================================================
111 ipy = 200
112 ipz = 300
113 ipa = 400
114 ipid = pid(1)
115 area = geo(1 ,ipid)
116 DO i=1,3
117 ii(i) = nel*(i-1)
118 ENDDO
119C-------------------
120C STRAIN RATE
121C-------------------
122 israte = ipm(3,imat)
123 asrate = min(one, pm(9,imat)*dtime)
124c calculate total strain rate on neutral fiber position
125 DO i = 1,nel
126 epsdi = sqrt(exx(i)**2 + half*(exy(i)**2 + exz(i)**2))
127 IF (israte > 0) THEN
128 epsd(i)= asrate*epsdi + (one - asrate)*epsd(i)
129 ELSE
130 epsd(i)= epsdi
131 ENDIF
132 ENDDO
133c
134 ! Save old OFF value
135 off_old(1:nel) = off(1:nel)
136c-------------------
137c STRAIN
138c-------------------
139 DO i = 1,nel
140 exx(i) = exx(i)*dtime
141 exy(i) = exy(i)*dtime
142 exz(i) = exz(i)*dtime
143 kxx(i) = kxx(i)*dtime
144 kyy(i) = kyy(i)*dtime
145 kzz(i) = kzz(i)*dtime
146 ENDDO
147C
148 DO i = 1,nel
149 degmb(i) = for(i,1)*exx(i)
150 degsh(i) = for(i,2)*exy(i) + for(i,3)*exz(i)
151 degfx(i) = mom(i,1)*kxx(i) + mom(i,2)*kyy(i) + mom(i,3)*kzz(i)
152 ENDDO
153C
154 IF (isigi == 0 .OR. (isigi /= 0 .AND. time /= zero)) THEN
155 DO i = 1,nel
156 for(i,1) = zero
157 for(i,2) = zero
158 for(i,3) = zero
159 mom(i,1) = zero
160 mom(i,2) = zero
161 mom(i,3) = zero
162 ENDDO
163 ENDIF
164c
165 ! Plasticity specific tables
166 ipla = elbuf_str%BUFLY(1)%L_PLA
167 ALLOCATE (dpla(nel*ipla,npt*ipla))
168 ALLOCATE (sigy(nel*ipla,npt*ipla))
169 ALLOCATE (vm0(nel*ipla,npt*ipla))
170 ALLOCATE (vm(nel*ipla,npt*ipla))
171 IF (ipla > 0) THEN
172 DO j=1,npt
173 DO i = 1,nel
174 dpla(i,j) = elbuf_str%BUFLY(1)%LBUF(1,1,j)%PLA(i)
175 sigy(i,j) = 1.e30
176 sigoxx = elbuf_str%BUFLY(1)%LBUF(1,1,j)%SIG(ii(1)+i)
177 sigoxy = elbuf_str%BUFLY(1)%LBUF(1,1,j)%SIG(ii(2)+i)
178 sigoxz = elbuf_str%BUFLY(1)%LBUF(1,1,j)%SIG(ii(3)+i)
179 vm0(i,j) = sqrt(sigoxx**2 + three*(sigoxy**2 + sigoxz**2))
180 ENDDO
181 ENDDO
182 ENDIF
183c---------------------------
184c material models
185c---------------------------
186 IF (mtn == 2) THEN
187 CALL m2lawpi(elbuf_str,mat_param,
188 1 nel ,npt ,geo ,eint ,off ,
189 3 pid ,epsd ,exx ,exy ,exz ,
190 4 kxx ,kyy ,kzz ,al ,asrate ,
191 5 dtime ,jthe ,tempel ,sigy )
192c
193 ELSE
194c
195 CALL mulaw_ib(elbuf_str,
196 1 nel ,npt ,mtn ,imat ,
197 2 pid ,ngl ,ipm ,
198 3 geo ,off ,
199 4 epsd ,bufmat ,npf ,tf ,
200 5 exx ,exy ,exz ,kxx ,
201 6 kyy ,kzz ,jthe ,tempel ,sigy )
202 ENDIF
203!
204c---------------------------
205c Computation of Plastic Work
206c---------------------------
207 IF (ipla > 0) THEN
208 DO j=1,npt
209 DO i = 1,nel
210 dpla(i,j) = elbuf_str%BUFLY(1)%LBUF(1,1,j)%PLA(i) - dpla(i,j)
211 signxx = elbuf_str%BUFLY(1)%LBUF(1,1,j)%SIG(ii(1)+i)
212 signxy = elbuf_str%BUFLY(1)%LBUF(1,1,j)%SIG(ii(2)+i)
213 signxz = elbuf_str%BUFLY(1)%LBUF(1,1,j)%SIG(ii(3)+i)
214 vm(i,j) = sqrt(signxx**2 + three*(signxy**2 + signxz**2))
215 IF (elbuf_str%GBUF%G_WPLA > 0) elbuf_str%GBUF%WPLA(i) = elbuf_str%GBUF%WPLA(i) +
216 . half*(vm(i,j)+vm0(i,j))*dpla(i,j)*al(i)*area/npt
217 ENDDO
218 ENDDO
219 ENDIF
220c---------------------------
221c failure models
222c---------------------------
223 IF (ifail > 0) THEN
224c
225 CALL fail_beam18(
226 . elbuf_str,mat_param%FAIL(1),mat_param,numgeo ,
227 . npropg ,snpc ,stf ,
228 . nel ,npt ,ipid ,jthe ,
229 . tempel ,ngl ,geo ,
230 . off ,epsd ,npf ,tf ,
231 . dpla ,eint ,time ,iout ,istdo ,
232 . al ,ismstr ,exx ,exy ,exz ,
233 . kxx ,kyy ,kzz ,dtime ,
234 . ntable ,table ,sigy )
235 END IF
236C------------------------------------
237c resultant force and moment
238c------------------------------------
239 IF (elbuf_str%BUFLY(1)%L_DMGSCL > 0 ) THEN
240 DO ipt = 1,npt
241 DO i=1,nel
242 ypt(i) = geo(ipy+ipt,pid(i))
243 zpt(i) = geo(ipz+ipt,pid(i))
244 apt(i) = geo(ipa+ipt,pid(i))
245 damage_loc = elbuf_str%BUFLY(1)%LBUF(1,1,ipt)%DMGSCL(i)
246 dfxx = apt(i)*elbuf_str%BUFLY(1)%LBUF(1,1,ipt)%SIG(ii(1)+i)*damage_loc
247 dfxy = apt(i)*elbuf_str%BUFLY(1)%LBUF(1,1,ipt)%SIG(ii(2)+i)*damage_loc
248 dfxz = apt(i)*elbuf_str%BUFLY(1)%LBUF(1,1,ipt)%SIG(ii(3)+i)*damage_loc
249 for(i,1) = for(i,1) + dfxx
250 for(i,2) = for(i,2) + dfxy
251 for(i,3) = for(i,3) + dfxz
252 mom(i,1) = mom(i,1) + dfxy*zpt(i) - dfxz*ypt(i)
253 mom(i,2) = mom(i,2) + dfxx*zpt(i)
254 mom(i,3) = mom(i,3) - dfxx*ypt(i)
255 ENDDO
256 ENDDO
257 ELSE
258C-----------------------
259C FORCES AND MOMENTS
260C-----------------------
261 DO ipt = 1,npt
262 DO i=1,nel
263 ypt(i) = geo(ipy+ipt,pid(i))
264 zpt(i) = geo(ipz+ipt,pid(i))
265 apt(i) = geo(ipa+ipt,pid(i))
266 dfxx = apt(i)*elbuf_str%BUFLY(1)%LBUF(1,1,ipt)%SIG(ii(1)+i)
267 dfxy = apt(i)*elbuf_str%BUFLY(1)%LBUF(1,1,ipt)%SIG(ii(2)+i)
268 dfxz = apt(i)*elbuf_str%BUFLY(1)%LBUF(1,1,ipt)%SIG(ii(3)+i)
269 for(i,1) = for(i,1) + dfxx
270 for(i,2) = for(i,2) + dfxy
271 for(i,3) = for(i,3) + dfxz
272 mom(i,1) = mom(i,1) + dfxy*zpt(i) - dfxz*ypt(i)
273 mom(i,2) = mom(i,2) + dfxx*zpt(i)
274 mom(i,3) = mom(i,3) - dfxx*ypt(i)
275 ENDDO
276 ENDDO
277 ENDIF
278
279
280C-------------------------------------
281 DO i = 1,nel
282 for(i,1) = for(i,1)*off(i)
283 for(i,2) = for(i,2)*off(i)
284 for(i,3) = for(i,3)*off(i)
285 mom(i,1) = mom(i,1)*off(i)
286 mom(i,2) = mom(i,2)*off(i)
287 mom(i,3) = mom(i,3)*off(i)
288 f1(i) = for(i,1)
289 f2(i) = for(i,2)
290 f3(i) = for(i,3)
291 m1(i) = mom(i,1)
292 m2(i) = mom(i,2)
293 m3(i) = mom(i,3)
294 ENDDO
295c---------------------------
296c element damping is removed outside with a common call
297C---------------------------
298C Internal energy
299C---------------------------
300 DO i = 1,nel
301 degmb(i) = degmb(i) + for(i,1)*exx(i)
302 degsh(i) = degsh(i) + for(i,2)*exy(i) + for(i,3)*exz(i)
303 degfx(i) = degfx(i)
304 . + mom(i,1)*kxx(i) + mom(i,2)*kyy(i) + mom(i,3)*kzz(i)
305 fact = half*off(i)*al(i)
306 eint(i,1) = eint(i,1) + fact*(degmb(i)+degsh(i))
307 eint(i,2) = eint(i,2) + fact* degfx(i)
308 ENDDO
309c----------------------------
310c Check element erosion (IMCONV = 1 => convergence of implicit)
311c----------------------------
312 DO i = 1,nel
313 IF (off(i) == four_over_5 .AND. imconv == 1) THEN
314#include "lockon.inc"
315 WRITE(iout, 1000) ngl(i)
316 WRITE(istdo,1100) ngl(i),time
317#include "lockoff.inc"
318 ENDIF
319 ENDDO
320C
321 DO i = 1,nel
322 IF (off(i) < em01) off(i) = zero
323 IF (off(i) < one ) off(i) = off(i)*four_over_5
324 ENDDO
325c
326c--------------------------------------------------------
327c SHOOTING NODES ALGORITHM ACTIVATION
328c--------------------------------------------------------
329 DO i = 1,nel
330 IF ((off_old(i) > zero) .AND. (off(i) == zero)) THEN
331 idel7nok = 1
332 ENDIF
333 ENDDO
334c
335 DEALLOCATE (dpla)
336C------------------------------------------
337 1000 FORMAT(1x,'-- RUPTURE OF BEAM ELEMENT NUMBER ',i10)
338 1100 FORMAT(1x,'-- RUPTURE OF BEAM ELEMENT :',i10,' AT TIME :',g11.4)
339C------------------------------------------
340 RETURN
341 END SUBROUTINE main_beam18
my_real
#define my_real
Definition cppsort.cpp:32
fail_beam18
subroutine fail_beam18(elbuf_str, fail, mat_param, numgeo, npropg, snpc, stf, nel, npt, iprop, jthe, tempel, ngl, geo, off, epsd, npf, tf, dpla, eint, time, iout, istdo, al, ismstr, exx, exy, exz, kxx, kyy, kzz, dtime, ntable, table, sigy)
Definition fail_beam18.F:54
area
subroutine area(d1, x, x2, y, y2, eint, stif0)
Definition hm_read_prop32.F:567
m2lawpi
subroutine m2lawpi(elbuf_str, mat_param, nel, npt, geo, eint, off, pid, epsp, exx, exy, exz, kxx, kyy, kzz, al, asrate, timestep, jthe, temp, sigy)
Definition m2lawpi.F:36
min
#define min(a, b)
Definition macros.h:20
main_beam18
subroutine main_beam18(elbuf_str, nel, npt, mtn, imat, pid, ngl, pm, ipm, geo, off, for, mom, eint, al, epsd, bufmat, npf, tf, exx, exy, exz, kxx, kyy, kzz, f1, f2, f3, m1, m2, m3, jthe, tempel, ifail, sbufmat, snpc, stf, nummat, numgeo, iout, istdo, npropmi, npropm, npropg, time, dtime, idel7nok, isigi, imconv, ismstr, mat_param, ntable, table)
Definition main_beam18.F:47
mulaw_ib
subroutine mulaw_ib(elbuf_str, nel, npt, mtn, imat, pid, ngl, ipm, geo, off, epsp, bufmat, npf, tf, exx, exy, exz, kxx, kyy, kzz, jthe, tempel, sigy)
Definition mulaw_ib.F:42
for
for(i8=*sizetab-1;i8 >=0;i8--)
Definition mumps_common.c:91