OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
law62_upd.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!|| law62_upd ../starter/source/materials/mat/mat062/law62_upd.F
25!||--- called by ------------------------------------------------------
26!|| updmat ../starter/source/materials/updmat.F
27!||--- uses -----------------------------------------------------
28!|| message_mod ../starter/share/message_module/message_mod.F
29!|| table_mod ../starter/share/modules1/table_mod.F
30!||====================================================================
31 SUBROUTINE law62_upd(IOUT,TITR ,MAT_ID,NUPARAM,UPARAM, PM ,
32 . GAMA_INF )
33C-----------------------------------------------
34C M o d u l e s
35C-----------------------------------------------
36 USE message_mod
37 USE table_mod
38 USE names_and_titles_mod , ONLY : nchartitle
39C-----------------------------------------------
40C I m p l i c i t T y p e s
41C-----------------------------------------------
42#include "implicit_f.inc"
43C-----------------------------------------------
44C C o m m o n B l o c k s
45C-----------------------------------------------
46#include "param_c.inc"
47C-----------------------------------------------
48C D u m m y A r g u m e n t s
49C-----------------------------------------------
50 CHARACTER(LEN=NCHARTITLE) :: TITR
51 INTEGER MAT_ID,IOUT,NUPARAM
52 my_real, DIMENSION(NUPARAM), INTENT(INOUT) :: uparam
53 my_real, INTENT(IN) :: gama_inf
54 my_real, DIMENSION(:), INTENT(INOUT) :: pm(npropm)
55C-----------------------------------------------
56C L o c a l V a r i a b l e s
57C-----------------------------------------------
58 INTEGER NDATA,I,K,INDX,ITER,NORDRE,J,NPRONY
59 my_real
60 . lam_min,lam_max,strain_min,strain_max,
61 . d11,d22,d12,lam1,lam2,lam3,lam12,invd1,invd2,a1,a2,a3,d33,d13,d23,
62 . aa,b,bb,dt2,dt3,t2,t3,invd3,rv,bulk,gs,nug,young
63 my_real , DIMENSION(:), ALLOCATABLE :: stress,stretch,mu,al,beta,
64 . muoveral,albeta,mubeta,nu
65 INTEGER , DIMENSION(:), ALLOCATABLE :: ITAB_ON_A,INDEX
66
67C=======================================================================
68! CHECK the material stability (Drucker proguer conditions)
69!====================================================================
70 !! NDATA = 5000 ! La=0.01,10, EM3!!
71 lam_min =em03
72 lam_max =four
73 ndata = nint((lam_max-lam_min)/em03)
74C
75 nordre = nint(uparam(2))
76 nprony = nint(uparam(3))
77 ALLOCATE ( mu(nordre),al(nordre),beta(nordre),
78 . muoveral(nordre),albeta(nordre),mubeta(nordre),nu(nordre) )
79 DO j= 1,nordre
80 mu(j) = uparam(6 + j)/gama_inf
81 al(j) = uparam(6 + nordre + j)
82 beta(j) = uparam(2*nordre + 2*nprony + 7 + j)
83 muoveral(j) = two*mu(j)/al(j)
84 albeta(j) = al(j)*beta(j)
85 mubeta(j) = mu(j)*beta(j)
86 nu(j) = one*beta(j)/(one + two*beta(j))
87 ENDDO
88 gs = zero
89 bulk = zero
90 IF(gama_inf /= one) THEN
91 DO j= 1,nordre
92 uparam(6 + j) = mu(j)
93 gs = gs + mu(j)
94 bulk = bulk + two*mu(j)*(third + beta(j))
95 ENDDO
96 nug = half*(three*bulk - two*gs)/(three*bulk+ gs)
97 uparam(1) = nug
98 uparam(5) = bulk
99C parameters
100 gs = gs*two
101 young = gs*(one + nug)
102 pm(20) = young
103 pm(21) = nug
104 pm(22) = gs !! TWO*G
105 pm(24) = young/(one - nug**2)
106 pm(32) = bulk
107 pm(100) = bulk
108C Formulation for solid elements time step computation.
109 pm(105) = gs/(bulk + two_third*gs)
110 ENDIF
111C
112 ALLOCATE (stretch(ndata))
113 ALLOCATE (stress(ndata))
114 ALLOCATE (itab_on_a(ndata))
115 ALLOCATE (index(ndata))
116 stretch(1)=lam_min
117 itab_on_a = 0
118 DO k= 2,ndata
119 stretch(k)=stretch(k-1) + em03
120 ENDDO
121 stress=zero
122C
123 IF(gama_inf /= one ) THEN
124 WRITE(iout,1000) mat_id
125 WRITE(iout,1100) nug,gs*half,bulk,nordre
126 WRITE(iout,1200) (mu(j),al(j),nu(j),j=1,nordre)
127 ENDIF
128 WRITE(iout,2000)mat_id
129
130 ! uniaxial Lam1= lam, lam2=lam3 (find lam2 (T2(lam2) = 0)
131 indx =0
132 DO i = 1, ndata
133 d11= zero
134 d22= zero
135 d33 = zero
136 d12 = zero
137 d13 = zero
138 d23 = zero
139 lam1 =stretch(i)
140 lam2 = lam1 ! starting point
141 lam3 = lam2
142 DO iter = 1,20
143 rv = lam1*lam2*lam3 ! lam3=lam2
144 t2= zero
145 dt2 = zero
146 DO j=1,nordre
147 aa = exp(al(j)*log(lam2))
148 bb = exp(-albeta(j)*log(rv))
149 t2 = t2 + muoveral(j)*(aa - bb)
150 dt2 = dt2 + two*mu(j)*aa + two*mubeta(j)*bb
151 dt2 = dt2 /lam2
152 ENDDO ! NORDRE
153 IF(dt2 /= zero) lam2 = lam2 - t2/dt2
154 lam2 = max(em03, lam2)
155 lam3 = lam2
156 ENDDO ! iter
157 rv= lam1*lam2*lam3
158 a1 = zero
159 a2 = zero
160 a3 = zero
161 b = zero
162 DO j=1,nordre
163 a1 = a1 + two*mu(j)*exp(al(j)*log(lam1))
164 a2 = a2 + two*mu(j)*exp(al(j)*log(lam2))
165 b = b + two*mubeta(j)*exp(-albeta(j)*log(rv))
166 ENDDO
167 d11 = a1 + b
168 d22 = a2 + b
169 d33 = d22
170 !D12 = D13= D21 = D23 = D31 = D32 = B
171 aa = b**2
172 invd1 = d11 + d22 + d33
173 invd2 = d11*d22 + d22*d33 + d11*d33 - three*aa
174 invd3 = d11*d22*d33 + two*aa*b - two*d22*aa - d11*aa
175 IF (invd1 > zero .AND. invd2 > zero .AND. invd3 > zero) THEN
176 indx = indx +1
177 itab_on_a(indx) = 1
178 index(indx) = i
179 ENDIF
180 ENDDO ! nDATA
181 IF(indx > 0 .AND. indx < ndata) THEN
182 WRITE(iout,2010)
183 i = index(1) - 1
184 IF(i > 1 .AND. i < ndata) THEN
185 strain_min = stretch(i) - one
186 WRITE(iout,2100)strain_min
187 ENDIF
188 i = index(indx) + 1
189 IF(i < ndata)THEN
190 strain_max = stretch(i) - one
191 WRITE(iout,2200)strain_max
192 ENDIF
193 ENDIF
194C Biaxial
195 indx =0
196 DO i = 1, ndata
197 d11 = zero
198 d22 = zero
199 d33 = zero
200 d12 = zero
201 d13 = zero
202 d23 = zero
203 lam1 = stretch(i)
204 lam2 = lam1
205 lam3 = lam1 ! staring point
206C ! bi-axial Lam1=lam2= lam, lam3 (find lam3 corresponding to (T3(lam3) = 0)
207 DO iter = 1,20
208 rv = lam1*lam2*lam3 !
209 t3 = zero
210 dt3 = zero
211 DO j=1,nordre
212 aa = exp(al(j)*log(lam3))
213 bb = exp(-albeta(j)*log(rv))
214 t3 = t3 + muoveral(j)*(aa - bb)
215 dt3 = dt3 + two*mu(j)* aa + two*mubeta(j)*bb
216 dt3 = dt3 /lam3
217 ENDDO ! NORDRE
218 IF(dt3 /= zero) lam3 = lam3 - t3/dt3
219 lam3 = max(em03, lam3)
220 ENDDO ! Iter
221 rv= lam1*lam2*lam3
222 a1 = zero
223 a2 = zero
224 a3 = zero
225 b = zero
226 DO j=1,nordre
227 a1 = a1 + two*mu(j)*exp(al(j)*log(lam1))
228 a3 = a3 + two*mu(j)*exp(al(j)*log(lam3))
229 b = b + two*mubeta(j)*exp(-albeta(j)*log(rv))
230 ENDDO
231 d11 = a1 + b
232 d22 = d11
233 d33 = a3 + b
234 !D12 = D13= D21 = D23 = D31 = D32 = B
235 aa = b**2
236 invd1 = d11 + d22 + d33
237 invd2 = d11*d22 + d22*d33 + d11*d33 - three*aa
238 invd3 = d11*d22*d33 + two*aa*b - two*d22*aa - d11*aa
239 IF (invd1 > zero .AND. invd2 > zero .AND. invd3 > zero) THEN
240 indx = indx +1
241 itab_on_a(indx) = 1
242 index(indx) = i
243 ENDIF
244 ENDDO ! NDATA
245 IF(indx > 0 .AND. indx < ndata) THEN
246 WRITE(iout,2020)
247 i = index(1) - 1
248 IF(i > 1 .AND. i < ndata) THEN
249 strain_min = stretch(i) - one
250 WRITE(iout,2100)strain_min
251 ENDIF
252 i = index(indx) + 1
253 IF(i < ndata)THEN
254 strain_max = stretch(i) - one
255 WRITE(iout,2200)strain_max
256 ENDIF
257 ENDIF
258C plane test
259 indx =0
260 DO i = 1, ndata
261 d11 = zero
262 d22 = zero
263 d33 = zero
264 d12 = zero
265 d13 = zero
266 d23 = zero
267 lam1 = stretch(i)
268 lam2 = one
269 lam3 = lam1 ! starting point
270C ! bi-axial Lam1=LAM, lam2= ONE, lam3 (find lam3 corresponding to (T3(lam3) = 0)
271 DO iter = 1,20
272 rv = lam1*lam2*lam3 !
273 t3 = zero
274 dt3 = zero
275 DO j=1,nordre
276 aa = exp(al(j)*log(lam3))
277 bb = exp(-albeta(j)*log(rv))
278 t3 = t3 + muoveral(j)*(aa - bb )
279 dt3 = dt3 + two*mu(j)* aa + two*mubeta(j)*bb
280 dt3 = dt3 /lam3
281 ENDDO ! NORDRE
282 IF(dt3 /= zero) lam3 = lam3 - t3/dt3
283 lam3 = max(em03, lam3)
284 ENDDO ! iter
285 rv= lam1*lam2*lam3
286 a1 = zero
287 a2 = zero
288 a3 = zero
289 b = zero
290 DO j=1,nordre
291 a1 = a1 + two*mu(j)*exp(al(j)*log(lam1))
292 a2 = a2 + two*mu(j)*exp(al(j)*log(lam2))
293 a3 = a3 + two*mu(j)*exp(al(j)*log(lam3))
294 b = b + two*mubeta(j)*exp(-albeta(j)*log(rv))
295 ENDDO
296 d11 = a1 + b
297 d22 = a2 + b
298 d33 = a3 + b
299 !D12 = D13= D21 = D23 = D31 = D32 = B
300 aa = b**2
301 invd1 = d11 + d22 + d33
302 invd2 = d11*d22 + d22*d33 + d11*d33 - three*aa
303 invd3 = d11*d22*d33 + two*aa*b - two*d22*aa - d11*aa
304 IF (invd1 > zero .AND. invd2 > zero .AND. invd3 > zero) THEN
305 indx = indx +1
306 itab_on_a(indx) = 1
307 index(indx) = i
308 ENDIF
309 ENDDO ! NDATA
310!!
311 IF(indx > 0 .AND. indx < ndata) THEN
312 WRITE(iout,2030)
313 i = index(1) - 1
314 IF(i > 1 .AND. i < ndata) THEN
315 strain_min = stretch(i) - one
316 WRITE(iout,2100)strain_min
317 ENDIF
318 i = index(indx) + 1
319 IF(i < ndata)THEN
320 strain_max = stretch(i) - one
321 WRITE(iout,2200)strain_max
322 ENDIF
323 ENDIF
324C
325 DEALLOCATE (stretch)
326 DEALLOCATE (stress)
327 DEALLOCATE (itab_on_a)
328 DEALLOCATE (index)
329 DEALLOCATE ( mu,al,beta, muoveral,albeta,mubeta)
330 RETURN
331 1000 FORMAT
332 & (//5x, 'MODIFIED MATERIAL RIGIDITY ' ,/,
333 & 5x, ' ---------------------------', /,
334 & 5x, 'MATERIAL LAW = LAW62 ',/,
335 & 5x, 'MATERIAL NUMBER =',i10,//)
336 1100 FORMAT
337 &(5x,'EQUIVALENT POISSON RATIO . . . . . . .=',e12.4/
338 &,5x,'INITIAL SHEAR MODULUS . . . . . . . . .=',e12.4/
339 & 5x,'INITIAL BULK MODULUS. . .. . . . . . .=',e12.4//
340 &,5x,'ORDER OF STRAIN ENERGY. . . . . . . . .=',i8)
341 1200 FORMAT(
342 & 7x,'MATERIAL PARAMETER (MU). . . . . . . . =',e12.4/
343 & 7x,'MATERIAL PARAMETER (ALPHA) . . . . . . =',e12.4/
344 & 7x,'MATERIAL PARAMETER (NU) . . . . . . . =',e12.4/)
345 2000 FORMAT
346 & (//5x, 'CHECK THE DRUCKER PRAGER STABILITY CONDITIONS ' ,/,
347 & 5x, ' -----------------------------------------------', /,
348 & 5x, 'MATERIAL LAW = LAW62 ',/,
349 & 5x, 'MATERIAL NUMBER =',i10,//)
350 2010 FORMAT(
351 & 7x,'TEST TYPE = UNIXIAL ')
352 2020 FORMAT(//,
353 & 7x,'TEST TYPE = BIAXIAL ')
354 2030 FORMAT(//,
355 & 7x,'TEST TYPE = PLANAR (SHEAR)')
356 2100 FORMAT(
357 & 8x,'COMPRESSION: UNSTABLE AT A NOMINAL STRAIN LESS THAN ',1pg20.13)
358 2200 FORMAT(
359 & 8x,'TENSION: UNSTABLE AT A NOMINAL STRAIN LARGER THAN ',1pg20.13)
360c-----------
361 END
my_real
#define my_real
Definition cppsort.cpp:32
law62_upd
subroutine law62_upd(iout, titr, mat_id, nuparam, uparam, pm, gama_inf)
Definition law62_upd.F:33
max
#define max(a, b)
Definition macros.h:21
message_mod
Definition message_mod.F:1249
names_and_titles_mod
Definition names_and_titles_mod.F:997
names_and_titles_mod::nchartitle
integer, parameter nchartitle
Definition names_and_titles_mod.F:1003
table_mod
Definition table_mod.F:112