OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
law69_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!|| law69_upd ../starter/source/materials/mat/mat069/law69_upd.F
25!||--- called by ------------------------------------------------------
26!|| updmat ../starter/source/materials/updmat.F
27!||--- calls -----------------------------------------------------
28!|| ancmsg ../starter/source/output/message/message.F
29!|| arret ../starter/source/system/arret.F
30!|| law69_nlsqf ../starter/source/materials/tools/nlsqf.F
31!|| law69_nlsqf_auto ../starter/source/materials/tools/law69_nlsqf_auto.F
32!||--- uses -----------------------------------------------------
33!|| message_mod ../starter/share/message_module/message_mod.F
34!|| table_mod ../starter/share/modules1/table_mod.F
35!||====================================================================
36 SUBROUTINE law69_upd(IOUT,TITR ,MAT_ID,UPARAM,NFUNC,NFUNCT,
37 . IFUNC, FUNC_ID,NPC ,PLD ,PM,IPM, GAMA_INF)
38C-----------------------------------------------
39C M o d u l e s
40C-----------------------------------------------
41 USE message_mod
42 USE table_mod
43 USE names_and_titles_mod , ONLY : nchartitle
44C-----------------------------------------------
45C I m p l i c i t T y p e s
46C-----------------------------------------------
47#include "implicit_f.inc"
48C-----------------------------------------------
49C C o m m o n B l o c k s
50C-----------------------------------------------
51#include "param_c.inc"
52C-----------------------------------------------
53C D u m m y A r g u m e n t s
54C-----------------------------------------------
55 CHARACTER(LEN=NCHARTITLE) :: TITR
56 INTEGER MAT_ID,IOUT
57 INTEGER ,INTENT(IN) :: NFUNC
58 INTEGER ,INTENT(IN) :: NFUNCT
59 INTEGER NPC(*), FUNC_ID(NFUNCT) ,IPM(NPROPMI)
60 my_real
61 . uparam(*),pld(*),pm(npropm)
62 my_real , INTENT(IN) :: gama_inf
63 INTEGER, DIMENSION(NFUNC):: IFUNC
64C-----------------------------------------------
65C L o c a l V a r i a b l e s
66C-----------------------------------------------
67
68c
69 INTEGER N_NETWORK,N,K,ITEST,ICHECK,II,JJ,NSTART,IC1,IC2,NOGD,NDATA,NMULA,
70 . tab,nmual,iflag,is,i,indx,lawid, itag,is_encrypted
71 my_real
72 . mu(5),al(5),lam_min,lam_max,strain_min,strain_max,
73 . d11,d22,d12,lam1,lam2,lam3,lam12,invd1,invd2
74 my_real
75 . e,nu,gs,rbulk, d,young,
76 . errtol,ave_slope,g,mu_max,mu_min,dx,lam,beta,
77 . amula(2),mual(10),dy
78 my_real , DIMENSION(:), ALLOCATABLE :: stress,stretch
79 INTEGER , DIMENSION(:), ALLOCATABLE :: ITAB_ON_A,INDEX
80!
81C====================================================================
82C=======================================================================
83! IDENTIFICATION
84!====================================================================
85
86 IF(ifunc(2) == 0 ) RETURN
87C
88 nmual = 2*uparam(18)
89 iflag = nint(uparam(13))
90 icheck = nint(uparam(19))
91 nstart = nint(uparam(20))
92 is_encrypted = nint(uparam(21))
93 errtol = fiveem3
94C
95 ic1 = npc(ifunc(2))
96 ic2 = npc(ifunc(2)+1)
97C
98 nogd=(ic2-ic1)/2
99 ndata=nogd
100C
101 ALLOCATE (stretch(nogd))
102 ALLOCATE (stress(nogd))
103C
104 ave_slope = zero
105 jj=0
106 stretch=zero
107 stress=zero
108 g=zero
109 rbulk=zero
110 gs=zero
111 lam_max= zero
112 lam_min= zero
113 itag = 0
114 is = ifunc(2)
115 DO ii = ic1,ic2-2,2
116 jj=jj+1
117 stretch(jj) = pld(ii) + one
118 stress(jj) = pld(ii+1)
119 IF(pld(ii) < zero .AND. stress(jj) < zero ) itag = 1
120 IF(pld(ii) <= -one) THEN
121 CALL ancmsg(msgid=1175,
122 . msgtype=msgerror,
123 . anmode=aninfo,
124 . i1=mat_id,
125 . c1=titr,
126 . i2=func_id(is)) ! Id_function
127 CALL arret(2)
128 ENDIF
129 ENDDO
130C check if the curve is monotonic
131 DO jj =1,nogd - 1
132 dx = stretch(jj+1) - stretch(jj)
133 dy = stress(jj + 1) - stress(jj)
134 IF(dx <= zero .OR. dy <= zero) THEN
135 CALL ancmsg(msgid=1176,
136 . msgtype=msgerror,
137 . anmode=aninfo,
138 . i1=mat_id,
139 . c1=titr,
140 . i2=func_id(is))
141 ENDIF
142 ENDDO
143 ave_slope = zero
144
145 IF (iflag == 2 .AND. nmual /= 4) THEN
146 CALL ancmsg(msgid=126,
147 . msgtype=msgerror,
148 . anmode=aninfo_blind_1,
149 . i1=mat_id,
150 . c1=titr,
151 . i2=func_id(is))
152 ENDIF
153 DO ii=1,5
154 mual(2*ii-1)=five
155 mual(2*ii)=two
156 ENDDO
157!
158 WRITE(iout,277)
159 WRITE(iout,270)trim(titr),mat_id
160 IF(is_encrypted > 0) WRITE(iout,'(5X,A,//)')'CONFIDENTIAL DATA'
161!
162 IF (iflag == 1 .OR. iflag == 2) THEN
163 lawid = iflag
164 CALL law69_nlsqf(lawid,stretch,stress,10,nmual,nogd,mual,
165 $ icheck, nstart, errtol,mat_id,titr,is_encrypted)
166 ELSEIF(iflag == -1) THEN
167 lawid = 2
168 IF(itag == 0) icheck = abs(icheck)
169 CALL law69_nlsqf_auto(lawid,stretch,stress,10,nmual,nogd,mual,
170 $ icheck, nstart, errtol,mat_id,titr,is_encrypted)
171 uparam(18) = nmual / 2
172 ENDIF
173 DEALLOCATE (stretch)
174 DEALLOCATE (stress)
175 IF(gama_inf < one) THEN
176 DO ii=1,5
177 uparam(ii)=mual(2*ii-1)/gama_inf
178 uparam(ii+5)=mual(2*ii)
179 gs=gs + mual(2*ii-1)*mual(2*ii)
180 ENDDO
181 IF(is_encrypted == 0) WRITE(iout,100) mat_id
182 ELSE
183 DO ii=1,5
184 uparam(ii)=mual(2*ii-1)
185 uparam(ii+5)=mual(2*ii)
186 gs=gs+mual(2*ii-1)*mual(2*ii)
187 ENDDO
188 ENDIF
189 g=gs/two
190 nu = uparam(14)
191 rbulk=gs*(one+uparam(14))
192 & /max(em30,three*(one-two*uparam(14)))
193 uparam(11)=rbulk
194 uparam(17)=gs
195C parameters
196 young = gs*(one + nu)
197 pm(20) = young
198 pm(21) = nu
199 pm(22) = gs
200 pm(24) = young/(one - nu**2)
201 pm(32) = gs
202 pm(100) = rbulk
203C Formulation for solid elements time step computation.
204 ipm(252)= 2
205 pm(105) = gs/(rbulk + two_third*g)
206
207C
208 IF(is_encrypted == 0) THEN
209 IF (lawid /= 2) THEN
210 WRITE(iout,278)
211 & uparam(1),uparam(2),uparam(3),uparam(4),
212 & uparam(5),uparam(6),uparam(7),uparam(8),
213 & uparam(9),uparam(10),g,rbulk
214 ELSE
215 WRITE(iout,279)
216 & uparam(1),uparam(2),uparam(6),uparam(7),
217 & half*uparam(1),-half*uparam(2),g,rbulk
218 ENDIF
219 ENDIF
220c----------------
221c end of optimization loop
222c----------------
223C=======================================================================
224! CHECK the material stability (Drucker proguer conditions)
225!====================================================================
226 ndata = 10000 ! La=0.1,10, EM3!!
227 lam_min= em3
228 lam_max= ten
229C
230 DO i=1,5
231 mu(i) = uparam(i)
232 al(i) = uparam(5+i)
233 ENDDO
234C
235 ALLOCATE (stretch(ndata))
236 ALLOCATE (stress(ndata))
237 ALLOCATE (itab_on_a(ndata))
238 ALLOCATE (index(ndata))
239 stretch(1)=lam_min
240 itab_on_a = 0
241 DO k= 2,ndata
242 stretch(k)=stretch(k-1) + em3
243 ENDDO
244 stress=zero
245C
246 WRITE(iout,1000)mat_id
247C Tension/compression
248 indx =0
249 DO i = 1, ndata
250 d11= zero
251 d22= zero
252 d12= zero
253 lam1 =stretch(i)
254 lam2 = one/sqrt(lam1)
255 lam3 = lam2
256C
257 DO k=1,5
258 lam12 = (lam1*lam2)**(-al(k))
259 d11=d11+ al(k)*mu(k) * (lam1**al(k) + lam12 )
260 d22=d22+ al(k)*mu(k) * (lam2**al(k) + lam12 )
261 d12=d12+ al(k)*mu(k) * lam12
262 ENDDO
263
264 invd1 = d11 + d22
265 invd2 = d11*d22 - d12**2
266 IF (invd1 > 0 .AND. invd2 > 0) THEN
267 indx = indx +1
268 itab_on_a(indx) = 1
269 index(indx) = i
270 ENDIF
271 ENDDO
272 IF(indx > 0 .AND. indx < ndata) THEN
273 WRITE(iout,1010)
274 i = index(1) - 1
275 IF(i > 1) THEN
276 strain_min = stretch(i) - one
277 WRITE(iout,1100)strain_min
278 ENDIF
279 i = index(indx) + 1
280 IF(i <= ndata)THEN
281 strain_max = stretch(i) - one
282 WRITE(iout,1200)strain_max
283 ENDIF
284 ENDIF
285C Biaxial
286 indx =0
287 DO i = 1, ndata
288 d11= zero
289 d22= zero
290 d12= zero
291 lam1 =stretch(i)
292 lam2 = lam1
293 lam3 = one/lam1/lam1
294C
295 DO k=1,5
296 lam12 = (lam1*lam2)**(-al(k))
297 d11=d11+ al(k)*mu(k) * (lam1**al(k) + lam12 )
298 d22=d22+ al(k)*mu(k) * (lam2**al(k) + lam12 )
299 d12=d12+ al(k)*mu(k) * lam12
300 ENDDO
301
302 invd1 = d11 + d22
303 invd2 = d11*d22 - d12**2
304 IF (invd1 > 0 .AND. invd2 > 0) THEN
305 indx = indx + 1
306 itab_on_a(indx) = 1
307 index(indx) = i
308 ENDIF
309 ENDDO
310 IF(indx > 0 .AND. indx < ndata) THEN
311 WRITE(iout,1020)
312 i = index(1) - 1
313 IF(i > 1) THEN
314 strain_min = stretch(i) - one
315 WRITE(iout,1100)strain_min
316 ENDIF
317 i = index(indx) + 1
318 IF(i <= ndata)THEN
319 strain_max = stretch(i) - one
320 WRITE(iout,1200)strain_max
321 ENDIF
322 ENDIF
323C shear test
324 indx =0
325 DO i = 1, ndata
326 d11= zero
327 d22= zero
328 d12= zero
329 lam1 =stretch(i)
330 lam2 = one
331 lam3 = one/lam1
332C
333 DO k=1,5
334 lam12 = (lam1*lam2)**(-al(k))
335 d11=d11+ al(k)*mu(k) * (lam1**al(k) + lam12 )
336 d22=d22+ al(k)*mu(k) * (lam2**al(k) + lam12 )
337 d12=d12+ al(k)*mu(k) * lam12
338 ENDDO
339
340 invd1 = d11 + d22
341 invd2 = d11*d22 - d12**2
342 IF (invd1 > 0 .AND. invd2 > 0) THEN
343 indx = indx +1
344 itab_on_a(indx) = 1
345 index(indx) = i
346 ENDIF
347 ENDDO
348 IF(indx > 0 .AND. indx < ndata) THEN
349 WRITE(iout,1030)
350 i = index(1) - 1
351 IF(i > 1) THEN
352 strain_min = stretch(i) - one
353 WRITE(iout,1100)strain_min
354 ENDIF
355 i = index(indx) + 1
356 IF(i <= ndata)THEN
357 strain_max = stretch(i) - one
358 WRITE(iout,1200)strain_max
359 ENDIF
360 ENDIF
361 DEALLOCATE (stretch)
362 DEALLOCATE (stress)
363 DEALLOCATE (itab_on_a)
364 DEALLOCATE (index)
365C --------------------------------------------------------
366 RETURN
367 100 FORMAT
368 & (//5x, 'MODIFIED FITTED MATERIAL RIGIDITY ' ,/,
369 & 5x, ' ---------------------------', /,
370 & 5x, 'MATERIAL LAW = LAW69 ',/,
371 & 5x, 'MATERIAL NUMBER =',i10,//)
372 277 FORMAT
373 & (//5x, 'FITTED PARAMETERS FOR HYPERELASTIC_MATERIAL LAW69' ,/,
374 & 5x, ' ----------------------------------------')
375 270 FORMAT(
376 & 5x,a,/,
377 & 5x, 'material number =',I10,//)
378 278 FORMAT(
379 & 5X, 'ogden law parameters:',/
380 & ,5X, 'mu1 . . . . . . . . . . . . . . . . . .=',1PG20.13/
381 & ,5X, 'mu2 . . . . . . . . . . . . . . . . . .=',1PG20.13/
382 & ,5X, 'mu3 . . . . . . . . . . . . . . . . . .=',1PG20.13/
383 & ,5X, 'mu4 . . . . . . . . . . . . . . . . . .=',1PG20.13/
384 & ,5X, 'mu5 . . . . . . . . . . . . . . . . . .=',1PG20.13//
385 & ,5X, 'al1 . . . . . . . . . . . . . . . . . .=',1PG20.13/
386 & ,5X, 'al2 . . . . . . . . . . . . . . . . . .=',1PG20.13/
387 & ,5X, 'al3 . . . . . . . . . . . . . . . . . .=',1PG20.13/
388 & ,5X, 'al4 . . . . . . . . . . . . . . . . . .=',1PG20.13/
389 & ,5X, 'al5 . . . . . . . . . . . . . . . . . .=',1PG20.13//
390 & ,5X, 'initial shear modulus . . . . . . . . .=',1PG20.13/
391 & ,5X, 'bulk modulus. . . . . . . . . . . . . .=',1PG20.13//)
392 279 FORMAT(
393 & 5X, 'mooney-rivlin law parameters:',/
394 & ,5X, 'mu1 . . . . . . . . . . . . . . . . . .=',1PG20.13/
395 & ,5X, 'mu2 . . . . . . . . . . . . . . . . . .=',1PG20.13/
396 & ,5X, 'al1 . . . . . . . . . . . . . . . . . .=',1PG20.13/
397 & ,5X, 'al2 . . . . . . . . . . . . . . . . . .=',1PG20.13//
398 & ,5X, 'c10 . . . . . . . . . . . . . . . . . .=',1PG20.13/
399 & ,5X, 'c01 . . . . . . . . . . . . . . . . . .=',1PG20.13/
400 & ,5X, 'initial shear modulus . . . . . . . . .=',1PG20.13/
401 & ,5X, 'bulk modulus. . . . . . . . . . . . . .=',1PG20.13//)
402 1000 FORMAT
403 & (//5X, 'check the drucker prager stability conditions ' ,/,
404 & 5X, ' -----------------------------------------------', /,
405 & 5X, 'material law = ogden(law69) ',/,
406 & 5X, 'material number =',I10,//)
407 1010 FORMAT(
408 & 7X,'test TYPE = unixial ')
409 1020 FORMAT(//,
410 & 7X,'test TYPE = biaxial ')
411 1030 FORMAT(//,
412 & 7X,'test TYPE = planar (shear)')
413 1100 FORMAT(
414 & 8X,'compression: unstable at a nominal strain less than ',1pg20.13)
415 1200 FORMAT(
416 & 8x,'tension: unstable at a nominal strain larger than ',1PG20.13//)
417c-----------
418 END
my_real
#define my_real
Definition cppsort.cpp:32
the
end diagonal values have been computed in the(sparse) matrix id.SOL
law69_nlsqf_auto
subroutine law69_nlsqf_auto(lawid, stretch, y, ma, nmual, npt, mual, icheck, nstart, errtol, id, titr, is_encrypted)
Definition law69_nlsqf_auto.F:39
law69_upd
subroutine law69_upd(iout, titr, mat_id, uparam, nfunc, nfunct, ifunc, func_id, npc, pld, pm, ipm, gama_inf)
Definition law69_upd.F:38
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
law69_nlsqf
subroutine law69_nlsqf(lawid, stretch, y, ma, nmual, npt, mual, icheck, nstart, errtol, id, titr, is_encrypted)
Definition nlsqf.F:565
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
arret
subroutine arret(nn)
Definition arret.F:87