OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
law76_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!|| law76_upd ../starter/source/materials/mat/mat076/law76_upd.F
25!||--- called by ------------------------------------------------------
26!|| updmat ../starter/source/materials/updmat.F
27!||--- calls -----------------------------------------------------
28!|| ancmsg ../starter/source/output/message/message.F
29!|| func_comp ../starter/source/materials/mat/mat076/law76_upd.F
30!|| table2d_intersect ../starter/source/tools/curve/table2d_intersect.F
31!||--- uses -----------------------------------------------------
32!|| message_mod ../starter/share/message_module/message_mod.F
33!||====================================================================
34 SUBROUTINE law76_upd(IOUT ,TITR ,MAT_ID ,NUPARAM ,MATPARAM ,
35 . UPARAM ,NUMTABL ,ITABLE ,TABLE ,NFUNC ,
36 . IFUNC ,NPC ,PLD )
37C-----------------------------------------------
38C M o d u l e s
39C-----------------------------------------------
40 USE message_mod
41 USE matparam_def_mod
42 USE names_and_titles_mod , ONLY : nchartitle
43C-----------------------------------------------
44C I m p l i c i t T y p e s
45C-----------------------------------------------
46#include "implicit_f.inc"
47#include "com04_c.inc"
48C-----------------------------------------------
49C D u m m y A r g u m e n t s
50C-----------------------------------------------
51 INTEGER IOUT,MAT_ID,NUMTABL,NFUNC,NUPARAM
52 INTEGER ,DIMENSION(NUMTABL) :: ITABLE
53 INTEGER ,DIMENSION(NFUNC) :: IFUNC
54 INTEGER :: NPC(*)
55 my_real uparam(nuparam)
56 my_real pld(*)
57 CHARACTER(LEN=NCHARTITLE) :: TITR
58 TYPE(matparam_struct_) ,TARGET :: MATPARAM
59 TYPE(ttable), DIMENSION(NTABLE) ,INTENT(INOUT) ,TARGET :: TABLE
60C-----------------------------------------------
61C L o c a l V a r i a b l e s
62C-----------------------------------------------
63 INTEGER :: I,J,NDIM,NPT,NEPSP,FUNC_ID,FUNC_T,FUNC_C,FUNC_S,ICAS,ICONV,
64 . NPT_TRAC,NPT_COMP,NPT_SHEAR,NPTMAX,IFUN_NUP,IFX,IFY,STAT,LEN2
65 my_real :: xfac,epdt_min,epdt_max,epdc_min,epdc_max,epds_min,epds_max,
66 . nup,xint,yint
67 my_real ,DIMENSION(:) ,ALLOCATABLE :: x_comp,y_comp
68 TYPE(table_4d_), DIMENSION(:) ,POINTER :: TABLE_MAT
69C=======================================================================
70 func_t = itable(1)
71 func_c = itable(2)
72 func_s = itable(3)
73c
74 nup = uparam(9)
75 iconv = uparam(15)
76 icas = uparam(17)
77 xfac = uparam(18)
78c-----------------------------------------------------------------------------
79c Check yield stresses values (must be strictly positive)
80c-----------------------------------------------------------------------------
81 DO i = 1,numtabl
82 func_id = itable(i)
83 IF (func_id > 0) THEN
84 IF (table(func_id)%Y%VALUES(1) <= 0) THEN
85 CALL ancmsg(msgid=2063, msgtype=msgerror, anmode=aninfo_blind_1,
86 . i1=mat_id,
87 . c1=titr,
88 . i2=i)
89 ELSE IF (minval(table(func_id)%Y%VALUES) <= zero) THEN
90 ! Non strictly positive value was found
91 CALL ancmsg(msgid=2049, msgtype=msgwarning, anmode=aninfo_blind_1,
92 . i1=mat_id,
93 . c1=titr,
94 . i2=i)
95 ENDIF
96 ENDIF
97 ENDDO
98c---------------------------------------------------------------
99c check max and min strain rates by direction in yield tables
100c check if yield curves for different strain rates do not intersect
101c---------------------------------------------------------------
102 ndim = table(func_t)%NDIM
103 IF (ndim == 2) THEN
104 npt = SIZE(table(func_t)%X(1)%VALUES)
105 nepsp = SIZE(table(func_t)%X(2)%VALUES)
106 epdt_min = table(func_t)%X(2)%VALUES(1)*xfac
107 epdt_max = table(func_t)%X(2)%VALUES(nepsp)*xfac
108 uparam(19) = epdt_min
109 uparam(20) = epdt_max
110c
111 DO i = 2,nepsp
112 DO j = i+1,nepsp
113 CALL table2d_intersect(table(func_t) ,i ,j ,npt ,
114 . xfac ,xint ,yint)
115 IF (xint > zero .and. yint > zero) THEN
116 CALL ancmsg(msgid=3010, msgtype=msgwarning, anmode=aninfo,
117 . i1 = mat_id,
118 . i2 = table(func_t)%NOTABLE,
119 . c1 = titr ,
120 . r1 = table(func_t)%X(2)%VALUES(i)*xfac,
121 . r2 = table(func_t)%X(2)%VALUES(j)*xfac,
122 . r3 = xint,
123 . r4 = yint)
124 END IF
125 END DO
126 END DO
127 END IF
128c---
129 IF (func_c > 0) THEN
130 ndim = table(func_c)%NDIM
131 IF (ndim == 2) THEN
132 npt = SIZE(table(func_c)%X(1)%VALUES)
133 nepsp = SIZE(table(func_c)%X(2)%VALUES)
134 epdc_min = table(func_c)%X(2)%VALUES(1)*xfac
135 epdc_max = table(func_c)%X(2)%VALUES(nepsp)*xfac
136 uparam(21) = epdc_min
137 uparam(22) = epdc_max
138c
139 DO i = 2,nepsp
140 DO j = i+1,nepsp
141 CALL table2d_intersect(table(func_c) ,i ,j ,npt ,
142 . xfac ,xint ,yint)
143 IF (xint > zero .and. yint > zero) THEN
144 CALL ancmsg(msgid=3010, msgtype=msgwarning, anmode=aninfo,
145 . i1 = mat_id,
146 . i2 = table(func_t)%NOTABLE,
147 . c1 = titr ,
148 . r1 = table(func_t)%X(2)%VALUES(i)*xfac,
149 . r2 = table(func_t)%X(2)%VALUES(j)*xfac,
150 . r3 = xint,
151 . r4 = yint)
152 END IF
153 END DO
154 END DO
155 END IF
156 END IF
157c---
158 IF (func_s > 0) THEN
159 ndim = table(func_s)%NDIM
160 IF (ndim == 2) THEN
161 npt = SIZE(table(func_s)%X(1)%VALUES)
162 nepsp = SIZE(table(func_s)%X(2)%VALUES)
163 epds_min = table(func_s)%X(2)%VALUES(1)*xfac
164 epds_max = table(func_s)%X(2)%VALUES(nepsp)*xfac
165 uparam(23) = epds_min
166 uparam(24) = epds_max
167c
168 DO i = 2,nepsp
169 DO j = i+1,nepsp
170 CALL table2d_intersect(table(func_s) ,i ,j ,npt ,
171 . xfac ,xint ,yint)
172 IF (xint > zero .and. yint > zero) THEN
173 CALL ancmsg(msgid=3010, msgtype=msgwarning, anmode=aninfo,
174 . i1 = mat_id,
175 . i2 = table(func_t)%NOTABLE,
176 . c1 = titr ,
177 . r1 = table(func_t)%X(2)%VALUES(i)*xfac,
178 . r2 = table(func_t)%X(2)%VALUES(j)*xfac,
179 . r3 = xint,
180 . r4 = yint)
181 END IF
182 END DO
183 END DO
184 END IF
185 END IF
186c--------------------------------------------------------------------------
187c copy global function tables to local private storage for material law
188c--------------------------------------------------------------------------
189 ALLOCATE (matparam%TABLE(numtabl))
190 table_mat => matparam%TABLE(1:numtabl)
191 table_mat(1:numtabl)%NOTABLE = 0
192c
193c copy tension table
194c
195 IF (func_t > 0) THEN
196 table_mat(1)%NOTABLE = func_t
197 ndim = table(func_t)%NDIM
198 table_mat(1)%NDIM = ndim
199 ALLOCATE (table_mat(1)%X(ndim) ,stat=stat)
200c
201 DO i = 1,ndim
202 npt = SIZE(table(func_t)%X(i)%VALUES)
203 ALLOCATE (table_mat(1)%X(i)%VALUES(npt) ,stat=stat)
204 table_mat(1)%X(i)%VALUES(1:npt) = table(func_t)%X(i)%VALUES(1:npt)
205 END DO
206c
207 IF (ndim == 1) THEN
208 npt = SIZE(table(func_t)%X(1)%VALUES)
209 ALLOCATE (table_mat(1)%Y1D(npt) ,stat=stat)
210 table_mat(1)%Y1D(1:npt) = table(func_t)%Y%VALUES(1:npt)
211 ELSE IF (ndim == 2) THEN
212 npt = SIZE(table(func_t)%X(1)%VALUES)
213 len2 = SIZE(table(func_t)%X(2)%VALUES)
214 ALLOCATE (table_mat(1)%Y2D(npt,len2) ,stat=stat)
215 DO i=1,npt
216 DO j=1,len2
217 table_mat(1)%Y2D(i,j) = table(func_t)%Y%VALUES((j-1)*npt+i)
218 END DO
219 END DO
220 END IF
221 END IF
222c
223c copy compression table
224c
225 IF (func_c > 0) THEN
226
227 table_mat(2)%NOTABLE = func_c
228 ndim = table(func_c)%NDIM
229 table_mat(2)%NDIM = ndim
230 ALLOCATE (table_mat(2)%X(ndim) ,stat=stat)
231c
232 DO i = 1,ndim
233 npt = SIZE(table(func_c)%X(i)%VALUES)
234 ALLOCATE (table_mat(2)%X(i)%VALUES(npt) ,stat=stat)
235 table_mat(2)%X(i)%VALUES(1:npt) = table(func_c)%X(i)%VALUES(1:npt)
236 END DO
237c
238 IF (ndim == 1) THEN
239 npt = SIZE(table(func_c)%X(1)%VALUES)
240 ALLOCATE (table_mat(2)%Y1D(npt) ,stat=stat)
241 table_mat(2)%Y1D(1:npt) = table(func_c)%Y%VALUES(1:npt)
242 ELSE IF (ndim == 2) THEN
243 npt = SIZE(table(func_c)%X(1)%VALUES)
244 len2 = SIZE(table(func_c)%X(2)%VALUES)
245 ALLOCATE (table_mat(2)%Y2D(npt,len2) ,stat=stat)
246 DO i=1,npt
247 DO j=1,len2
248 table_mat(2)%Y2D(i,j) = table(func_c)%Y%VALUES((j-1)*npt+i)
249 END DO
250 END DO
251 END IF
252 END IF
253c
254c copy shear table
255c
256 IF (func_s > 0) THEN
257
258 table_mat(3)%NOTABLE = func_s
259 ndim = table(func_s)%NDIM
260 table_mat(3)%NDIM = ndim
261 ALLOCATE (table_mat(3)%X(ndim) ,stat=stat)
262c
263 DO i = 1,ndim
264 npt = SIZE(table(func_s)%X(i)%VALUES)
265 ALLOCATE (table_mat(3)%X(i)%VALUES(npt) ,stat=stat)
266 table_mat(3)%X(i)%VALUES(1:npt) = table(func_s)%X(i)%VALUES(1:npt)
267 END DO
268c
269 IF (ndim == 1) THEN
270 npt = SIZE(table(func_s)%X(1)%VALUES)
271 ALLOCATE (table_mat(3)%Y1D(npt) ,stat=stat)
272 table_mat(3)%Y1D(1:npt) = table(func_s)%Y%VALUES(1:npt)
273 ELSE IF (ndim == 2) THEN
274 npt = SIZE(table(func_s)%X(1)%VALUES)
275 len2 = SIZE(table(func_s)%X(2)%VALUES)
276 ALLOCATE (table_mat(3)%Y2D(npt,len2) ,stat=stat)
277 DO i=1,npt
278 DO j=1,len2
279 table_mat(3)%Y2D(i,j) = table(func_s)%Y%VALUES((j-1)*npt+i)
280 END DO
281 END DO
282 END IF
283 END IF
284c--------------------------------------------------------
285c Initialize plastic Poisson ratio if needed
286c--------------------------------------------------------
287 ifun_nup = ifunc(1)
288 IF (ifun_nup > 0) THEN
289 ifx = npc(ifun_nup)
290 ify = npc(ifun_nup + 1)
291 nup = pld(ify)
292 nup = max(zero, min(half, nup))
293 uparam(9) = nup
294 END IF
295c--------------------------------------------------------
296 IF (icas == 2) THEN
297 ! create new static compression table (1D)
298 ndim = 1
299 table_mat(2)%NOTABLE = 2
300 table_mat(2)%NDIM = ndim
301c
302 npt_trac = SIZE(table(func_t )%X(1)%VALUES)
303 npt_shear= SIZE(table(func_s )%X(1)%VALUES)
304 nptmax = npt_trac+npt_shear
305
306 ALLOCATE(x_comp(nptmax), stat=stat)
307 ALLOCATE(y_comp(nptmax), stat=stat)
308 x_comp(1:nptmax) = zero
309 y_comp(1:nptmax) = zero
310
311 CALL func_comp(table_mat,ntable ,nptmax ,npt_trac ,npt_shear ,
312 . npt_comp ,x_comp ,y_comp ,nup )
313
314
315 npt = npt_comp
316 ALLOCATE (table_mat(2)%X(ndim) ,stat=stat)
317 ALLOCATE (table_mat(2)%X(1)%VALUES(npt) ,stat=stat)
318 ALLOCATE (table_mat(2)%Y1D(npt) ,stat=stat)
319c
320 table_mat(2)%X(1)%VALUES(1:npt) = x_comp(1:npt)
321 table_mat(2)%Y1D(1:npt) = y_comp(1:npt)
322
323 icas =-1
324 iconv = 1
325
326 DEALLOCATE(x_comp,y_comp)
327
328 END IF
329c--------------------------------------------------------
330 uparam(15) = iconv
331 uparam(17) = icas
332c--------------------------------------------------------
333 RETURN
334 END
335!||====================================================================
336!|| func_comp ../starter/source/materials/mat/mat076/law76_upd.F
337!||--- called by ------------------------------------------------------
338!|| law76_upd ../starter/source/materials/mat/mat076/law76_upd.F
339!||--- uses -----------------------------------------------------
340!|| message_mod ../starter/share/message_module/message_mod.F
341!||====================================================================
342 SUBROUTINE func_comp(TABLE ,NTABLE ,NPTMAX ,NPT_TRAC ,NPT_SHEAR ,
343 . NPT_COMP,X_COMP ,Y_COMP ,NUP )
344C-----------------------------------------------
345C M o d u l e s
346C-----------------------------------------------
347 USE message_mod
348 USE matparam_def_mod
349C-----------------------------------------------
350C I m p l i c i t T y p e s
351C-----------------------------------------------
352#include "implicit_f.inc"
353C-----------------------------------------------
354C D u m m y A r g u m e n t s
355C-----------------------------------------------
356 INTEGER NTABLE,NPTMAX, NPT_TRAC,NPT_SHEAR,NPT_COMP
357 TYPE(TABLE_4D_), DIMENSION(NTABLE) :: TABLE
358 my_real :: NUP
359 my_real ,DIMENSION(NPTMAX) :: x_comp,y_comp
360C-----------------------------------------------
361C L o c a l V a r i a b l e s
362C-----------------------------------------------
363 INTEGER I,K,J,NT,NS,FUNC_TRAC,FUNC_SHEAR
364
365 my_real XI, XJ , NUM, DEN ,ALPHAMAX ,BETAMIN, SCALE_X_S,
366 . aa,bb,cc,delta,x1,x2,xa,xb,y1,y2,yc1,yc2,nupc,nup1,nup2
367 my_real x_t_shear(npt_shear)
368 my_real ,DIMENSION(NPTMAX) :: slope_trac,slope_shea,y_t, y_s,alpha,beta
369C=======================================================================
370 func_trac = 1
371 func_shear = 3
372 nt = table(func_trac)%NDIM
373 ns = table(func_shear)%NDIM
374 scale_x_s = sqrt(three)/(one+nup)
375 DO i = 1, npt_shear
376 x_t_shear(i) = scale_x_s* table(func_shear)%X(1)%VALUES(i)
377 ENDDO
378
379 i = 1
380 j = 1
381 k = 1
382 DO WHILE (i <= npt_trac .OR. j <= npt_shear)
383 IF (i <= npt_trac .AND. j <= npt_shear)THEN
384 xi = table(func_trac )%X(1)%VALUES(i)
385 xj = x_t_shear(j)
386 IF (xi < xj ) THEN
387 x_comp(k) = xi
388 y_t(k) = table(func_trac )%Y1D(i)
389 IF (j ==1) THEN
390 y_s(k) = table(func_shear)%Y1D(1) +
391 . (xi - x_t_shear(1) )*
392 . (table(func_shear)%Y1D(2) - table(func_shear)%Y1D(1))/
393 . (x_t_shear(2)- x_t_shear(1))
394 ELSE
395 y_s(k) = table(func_shear)%Y1D(j-1) +
396 . (xi - x_t_shear(j-1) )*
397 . (table(func_shear)%Y1D(j) - table(func_shear)%Y1D(j-1))/
398 . (x_t_shear(j)- x_t_shear(j-1))
399 ENDIF
400 i = i + 1
401 k = k + 1
402 ELSEIF (xj < xi ) THEN
403 x_comp(k) = xj
404 y_s(k) = table(func_shear )%Y1D(j)
405 IF (i ==1) THEN
406 y_t(k) = table(func_trac)%Y1D(1) +
407 . (xj - table(func_trac)%X(1)%VALUES(1) )*
408 . (table(func_trac)%Y1D(2) - table(func_trac)%Y1D(1))/
409 . (table(func_trac)%X(1)%VALUES(2)- table(func_trac)%X(1)%VALUES(1))
410 ELSE
411 y_t(k) = table(func_trac)%Y1D(i-1) +
412 . (xj - table(func_trac)%X(1)%VALUES(i-1) )*
413 . (table(func_trac)%Y1D(i) - table(func_trac)%Y1D(i-1))/
414 . (table(func_trac)%X(1)%VALUES(i)- table(func_trac)%X(1)%VALUES(i-1))
415 ENDIF
416
417 j = j + 1
418 k = k + 1
419 ELSEIF (xi == xj ) THEN
420 x_comp(k) = xi
421 y_t(k) = table(func_trac )%Y1D(i)
422 y_s(k) = table(func_shear )%Y1D(j)
423 i = i + 1
424 j = j + 1
425 k = k + 1
426 ENDIF
427 ELSEIF (i > npt_trac .AND. j <= npt_shear)THEN
428 xj=x_t_shear(j)
429 x_comp(k) = xj
430 y_s(k) = table(func_shear )%Y1D(j)
431 y_t(k) = table(func_trac)%Y1D(i-2) +
432 . (xj - table(func_trac)%X(1)%VALUES(i-2) )*
433 . (table(func_trac)%Y1D(i-1) - table(func_trac)%Y1D(i-2))/
434 . (table(func_trac)%X(1)%VALUES(i-1)- table(func_trac)%X(1)%VALUES(i-2))
435 j = j + 1
436 k = k + 1
437 ELSEIF (i <= npt_trac .AND. j > npt_shear)THEN
438 xi=table(func_trac )%X(1)%VALUES(i)
439 x_comp(k) = xi
440 y_t(k) = table(func_trac )%Y1D(i)
441 y_s(k) = table(func_shear)%Y1D(j-2) +
442 . (xi - x_t_shear(j-2) )*
443 . (table(func_shear)%Y1D(j-1) - table(func_shear)%Y1D(j-2))/
444 . (x_t_shear(j-1)- x_t_shear(j-2))
445 i = i + 1
446 k = k + 1
447 ELSE
448 EXIT
449 ENDIF
450 END DO
451 npt_comp = k - 1
452
453 alphamax = one
454 DO k= 2, npt_comp
455 slope_trac(k) = (y_t(k)-y_t(k-1)) / (x_comp(k)-x_comp(k-1))
456 slope_shea(k) = (y_s(k)-y_s(k-1)) / (x_comp(k)-x_comp(k-1))
457 IF( slope_trac(k)>zero .AND. slope_shea(k)>zero)THEN
458 alpha(k) = sqrt(three)*half *(slope_trac(k)/slope_shea(k)) * (y_s(k)/y_t(k))**2
459 alphamax = max(alphamax,alpha(k))
460 ENDIF
461 END DO
462 DO k= 1, npt_comp
463 num = sqrt(three) *alphamax *y_t(k)*y_s(k)
464 den = two * alphamax * y_t(k) - sqrt(three) * y_s(k)
465 y_comp(k) = num / max(em20, den)
466 END DO
467c-----------
468 !!
469c-----------
470 RETURN
471 END
472cc
473
my_real
#define my_real
Definition cppsort.cpp:32
alpha
#define alpha
Definition eval.h:35
func_comp
subroutine func_comp(table, ntable, nptmax, npt_trac, npt_shear, npt_comp, x_comp, y_comp, nup)
Definition law76_upd.F:344
law76_upd
subroutine law76_upd(iout, titr, mat_id, nuparam, matparam, uparam, numtabl, itable, table, nfunc, ifunc, npc, pld)
Definition law76_upd.F:37
min
#define min(a, b)
Definition macros.h:20
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
ancmsg
subroutine ancmsg(msgid, msgtype, anmode, i1, i2, i3, i4, i5, i6, i7, i8, i9, i10, i11, i12, i13, i14, i15, i16, i17, i18, i19, i20, r1, r2, r3, r4, r5, r6, r7, r8, r9, c1, c2, c3, c4, c5, c6, c7, c8, c9, prmode)
Definition message.F:889
table2d_intersect
subroutine table2d_intersect(table, i1, i2, npt, xfac, xint, yint)
Definition table2d_intersect.F:32