OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
m2lawp.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!|| m2lawp ../engine/source/materials/mat/mat002/m2lawp.F
25!||--- called by ------------------------------------------------------
26!|| main_beam3 ../engine/source/elements/beam/main_beam3.F
27!||--- uses -----------------------------------------------------
28!|| matparam_def_mod ../common_source/modules/mat_elem/matparam_def_mod.f90
29!||====================================================================
30 SUBROUTINE m2lawp(mat_param ,
31 . FOR ,MOM ,EINT ,GEO ,
32 . OFF ,PLA ,EXX ,EXY ,EXZ ,
33 . KXX ,KYY ,KZZ ,AL ,FA1 ,
34 . FA2 ,FA3 ,MA1 ,MA2 ,MA3 ,
35 . NEL ,PID ,NGL ,
36 . NUVAR ,UVAR ,SIGY)
37!-----------------------------------------------
38! m o d u l e s
39!-----------------------------------------------
40 use matparam_def_mod
41C-----------------------------------------------
42C I m p l i c i t T y p e s
43C-----------------------------------------------
44#include "implicit_f.inc"
45#include "comlock.inc"
46C-----------------------------------------------
47C C o m m o n B l o c k s
48C-----------------------------------------------
49#include "units_c.inc"
50#include "param_c.inc"
51#include "scr17_c.inc"
52#include "com08_c.inc"
53#include "impl1_c.inc"
54C-----------------------------------------------
55C D u m m y A r g u m e n t s
56C-----------------------------------------------
57 INTEGER ,INTENT(IN) :: NEL,NUVAR
58 INTEGER ,INTENT(IN) :: PID(NEL),NGL(NEL)
59 my_real :: OFF(*), PLA(*),AL(NEL),
60 . FOR(NEL,3), MOM(NEL,3), EINT(NEL,2), GEO(NPROPG,*),
61 . EXX(NEL),EXY(NEL),EXZ(NEL),KXX(NEL),KYY(NEL),
62 . KZZ(NEL),FA1(NEL),FA2(NEL),FA3(NEL),
63 . ma1(nel),ma2(nel),ma3(nel),a1(nel)
64 my_real ,DIMENSION(NEL,NUVAR), INTENT(INOUT) :: uvar
65 my_real ,DIMENSION(NEL) ,INTENT(INOUT) :: sigy
66 type (matparam_struct_) ,intent(in) :: mat_param
67C-----------------------------------------------
68C L o c a l V a r i a b l e s
69C-----------------------------------------------
70 INTEGER :: I,J,NINDX,IFORM,ICC,VP
71 INTEGER :: INDX(NEL)
72 my_real :: epdr,cc,plap1,asrate
73 my_real :: degfx(nel),esp(nel),g(nel),
74 . yldtmp(nel),ym(nel),b1(nel), b2(nel), b3(nel), degmb(nel),
75 . dmpm(nel), dmpf(nel), shf(nel),
76 . f1(nel), m1(nel), m2(nel), m3(nel),
77 . degsh(nel), yeq(nel), yld(nel), dwpla(nel),
78 . epmx(nel), dwelm(nel), dwelf(nel),z3(nel),z4(nel),
79 . ca(nel), cb(nel), cn(nel), ymax(nel),
80 . rr(nel), epsp(nel),
81 . sh(nel), yma2(nel), sh10(nel), sh20(nel), sh0(nel),
82 . sh1(nel), sh2(nel),plap(nel),dpla(nel),
83 . tmp1(nel), tmp2(nel), tmp3(nel)
84C-----------------------------------------------
85 IF (impl_s == 0 .OR. idyna > 0) THEN
86 DO i=1,nel
87 dmpm(i)=geo(16,pid(i))*al(i)
88 dmpf(i)=geo(17,pid(i))*al(i)
89 ENDDO
90 ELSE
91 DO i=1,nel
92 dmpm(i)=zero
93 dmpf(i)=zero
94 ENDDO
95 ENDIF
96C
97 iform = mat_param%iparam(1)
98 icc = mat_param%iparam(2)
99 vp = mat_param%iparam(3)
100 cc = mat_param%uparam(6)
101 epdr = mat_param%uparam(7)
102 if (vp == 1) then
103 epdr = max(epdr, em20)
104 else
105 epdr = max(epdr*dt1, em20)
106 endif
107 asrate = mat_param%uparam(9)
108 asrate = min(one, asrate*dt1)
109 DO i=1,nel
110 g(i) = mat_param%shear
111 ym(i) = mat_param%young
112C--------------------------
113 ca(i) = mat_param%uparam(1)
114 cb(i) = mat_param%uparam(2)
115 cn(i) = mat_param%uparam(3)
116 epmx(i)= mat_param%uparam(4)
117 ymax(i)= mat_param%uparam(5)
118C-----------------------------
119 a1(i) =geo(1,pid(i))
120 b1(i) =geo(2,pid(i))
121 b2(i) =geo(18,pid(i))
122 b3(i) =geo(4,pid(i))
123 shf(i) =geo(37,pid(i))
124 z3(i) = mat_param%uparam(10)
125 z4(i) = mat_param%uparam(11)
126 ENDDO
127C-----------------------------
128C DAMPING terms removed to pforce3
129C-----------------------------
130 DO i=1,nel
131 esp(i) = (eint(i,1)+eint(i,2))/al(i)/a1(i)
132 ENDDO
133C
134 DO i=1,nel
135 degmb(i) = for(i,1)*exx(i)
136 degsh(i) = for(i,2)*exy(i)+for(i,3)*exz(i)
137 degfx(i) = mom(i,1)*kxx(i)+mom(i,2)*kyy(i)+mom(i,3)*kzz(i)
138 ENDDO
139C
140C TRANSVERSE SHEAR CALCULATED WITH K1=12EI/L**2 K2=5/6GA
141C
142 DO i=1,nel
143 sh(i)=five_over_6*g(i)*a1(i)
144 yma2(i)=twelve*ym(i)/al(i)**2
145 sh10(i)=yma2(i)*b1(i)
146 sh20(i)=yma2(i)*b2(i)
147 sh0(i)=(one - shf(i))*sh(i)
148 sh1(i)=sh0(i)*sh10(i)/(sh(i)+sh10(i)) + shf(i)*sh10(i)
149 sh2(i)=sh0(i)*sh20(i)/(sh(i)+sh20(i)) + shf(i)*sh20(i)
150C
151 f1(i) =for(i,1)+ exx(i)*a1(i)*ym(i)
152 for(i,2)=for(i,2)+ exy(i)*sh2(i)
153 for(i,3)=for(i,3)+ exz(i)*sh1(i)
154 m1(i) =mom(i,1)+ kxx(i)*g(i) *b3(i)
155 m2(i) =mom(i,2)+ kyy(i)*ym(i)*b1(i)
156 m3(i) =mom(i,3)+ kzz(i)*ym(i)*b2(i)
157 ENDDO
158C-------------
159C CRITERE
160C-------------
161 DO i=1,nel
162 yeq(i)= f1(i)*f1(i) + three * a1(i) *
163 + ( m1(i)*m1(i) / max(b3(i),em20)
164 + + m2(i)*m2(i) / max(b1(i),em20)
165 + + m3(i)*m3(i) / max(b2(i),em20) )
166 yeq(i)= sqrt(yeq(i))/a1(i)
167 ENDDO
168C-------------
169C STRAIN RATE (JOHNSON-COOK, ZERILLI-ARMSTRONG)
170C-------------
171 IF (cc /= zero) THEN
172 DO i=1,nel
173 IF(vp == 1)THEN
174 plap(i) = uvar(i,1)
175 plap(i) = max(plap(i),epdr)
176 epsp(i) = log(plap(i)/epdr)
177 ELSE
178 epsp(i)=abs(degmb(i)+degfx(i))/(yeq(i)+ em20)/a1(i)
179 tmp2(i)=abs(degmb(i)+degfx(i))
180 tmp3(i)=epsp(i)
181 epsp(i)= max(epsp(i),epdr)
182 epsp(i)= log(epsp(i)/epdr)
183 ENDIF
184 ENDDO
185 IF (iform == 1) THEN ! zerilli
186 DO i=1,nel
187 epsp(i)=(one + cc * epsp(i))
188 tmp1(i)=epsp(i)
189 epsp(i)= cc*exp((-z3(i)+z4(i) * epsp(i))*esp(i))
190 IF(icc==1)ymax(i)= ymax(i) + epsp(i)
191 ca(i) = ca(i) + epsp(i)
192 epsp(i)=one
193 ENDDO
194 ELSE ! j-c
195 DO i=1,nel
196 epsp(i)=(one + cc * epsp(i))
197 IF (icc == 1) ymax(i) = ymax(i) * epsp(i)
198 ENDDO
199 ENDIF
200 ELSE
201 DO i=1,nel
202 epsp(i)=one
203 ENDDO
204 ENDIF
205C-----------------------------------
206C YIELD
207C-----------------------------------
208 DO i=1,nel
209 yld(i)= ca(i) + cb(i) * exp(cn(i) * log(pla(i)+em30))
210 yldtmp(i)=yld(i)
211 yld(i)= min(yld(i)*epsp(i),ymax(i))
212 rr(i) = min(one,yld(i)/(yeq(i)+ em20))
213 sigy(i) = yld(i)
214 ENDDO
215C
216 DO i=1,nel
217 f1(i) = f1(i)*rr(i)
218 dwelm(i) =(f1(i)+for(i,1))*(f1(i)-for(i,1))/ym(i)/a1(i)
219 degmb(i) = degmb(i) + f1(i)*exx(i)
220 ENDDO
221C
222 DO i=1,nel
223 m1(i) = m1(i)*rr(i)
224 m2(i) = m2(i)*rr(i)
225 m3(i) = m3(i)*rr(i)
226 dwelf(i) =(m1(i)+mom(i,1))*(m1(i)-mom(i,1))/ g(i)/b3(i)+
227 . (m2(i)+mom(i,2))*(m2(i)-mom(i,2))/ym(i)/b1(i)+
228 . (m3(i)+mom(i,3))*(m3(i)-mom(i,3))/ym(i)/b2(i)
229 degfx(i) = degfx(i)+ m1(i)*kxx(i)+m2(i)*kyy(i)+m3(i)*kzz(i)
230 ENDDO
231C
232 DO i=1,nel
233 dwpla(i)= degmb(i)+degfx(i)-dwelm(i)-dwelf(i)
234 ENDDO
235C-----------------------
236C EPS PLASTIQUE
237C-----------------------
238 DO i=1,nel
239 tmp1(i)=epsp(i)*dwpla(i)/yld(i)
240c PLA(I)=PLA(I)+OFF(I)*MAX(ZERO,0.5*EPSP(I)*DWPLA(I)/YLD(I)/A1(I))
241 dpla(i) = max(zero,half*tmp1(i)/a1(i))
242 pla(i) = pla(i)+off(i) * dpla(i)
243 ENDDO
244 DO i=1,nel
245 IF (vp == 1) THEN
246 plap1 = dpla(i) / max(em20,dt1)
247 uvar(i,1) = asrate * plap1 + (one - asrate) * plap(i)
248 ENDIF
249 ENDDO
250
251C--------------------------------
252C ductile failure test
253C--------------------------------
254 DO i=1,nel
255 IF (off(i) < em01) off(i) = zero
256 IF (off(i) < one ) off(i) = off(i)*four_over_5
257 ENDDO
258C
259 nindx = 0
260 DO i=1,nel
261 IF (off(i) < one) cycle
262 IF (pla(i) < epmx(i)) cycle
263 off(i)=four_over_5
264 idel7nok = 1
265 nindx=nindx+1
266 indx(nindx)=i
267 ENDDO
268C
269 IF (nindx > 0 .AND. imconv == 1) THEN
270 DO j=1,nindx
271#include "lockon.inc"
272 WRITE(iout, 1000) ngl(indx(j))
273 WRITE(istdo,1100) ngl(indx(j)),tt
274#include "lockoff.inc"
275 ENDDO
276 ENDIF
277C
278 DO i=1,nel
279 fa1(i) = f1(i)
280 fa2(i) = for(i,2)
281 fa3(i) = for(i,3)
282 ma1(i) = m1(i)
283 ma2(i) = m2(i)
284 ma3(i) = m3(i)
285 ENDDO
286C
287 1000 FORMAT(1x,'-- RUPTURE OF BEAM ELEMENT NUMBER ',i10)
288 1100 FORMAT(1x,'-- RUPTURE OF BEAM ELEMENT :',i10,' AT TIME :',g11.4)
289C-----------------------------------------------
290 RETURN
291 END
ymax
subroutine ymax(idn, fac, npc, pld, stiffmin, stiffmax, stiffini, stiffavg)
Definition law100_upd.F:274
m2lawp
subroutine m2lawp(mat_param, for, mom, eint, geo, off, pla, exx, exy, exz, kxx, kyy, kzz, al, fa1, fa2, fa3, ma1, ma2, ma3, nel, pid, ngl, nuvar, uvar, sigy)
Definition m2lawp.F:37
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21