OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
szhour3_or.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!|| szhour3_or ../engine/source/elements/solid/solidez/szhour3_or.F
25!||--- called by ------------------------------------------------------
26!|| szforc3 ../engine/source/elements/solid/solidez/szforc3.F
27!||--- calls -----------------------------------------------------
28!|| gfhour_or ../engine/source/elements/solid/solidez/gfhour_or.F
29!|| mdama24 ../engine/source/elements/solid/solidez/mdama24.F
30!|| mmod_norm ../engine/source/elements/solid/solidez/mmod_norm.F
31!|| mmodul ../engine/source/elements/solid/solidez/mmodul.F
32!|| szetfac ../engine/source/elements/solid/solidez/szetfac.F
33!|| szstrainhg ../engine/source/elements/solid/solidez/szstrainhg.F
34!|| szsvm_or ../engine/source/elements/solid/solidez/szsvm_or.F
35!||--- uses -----------------------------------------------------
36!|| elbufdef_mod ../common_source/modules/mat_elem/elbufdef_mod.F90
37!|| matparam_def_mod ../common_source/modules/mat_elem/matparam_def_mod.F90
38!||====================================================================
39 SUBROUTINE szhour3_or(
40 1 ELBUF_STR,PM, RHO, OFF,
41 2 VX1, VX2, VX3, VX4,
42 3 VX5, VX6, VX7, VX8,
43 4 VY1, VY2, VY3, VY4,
44 5 VY5, VY6, VY7, VY8,
45 6 VZ1, VZ2, VZ3, VZ4,
46 7 VZ5, VZ6, VZ7, VZ8,
47 8 F11, F21, F31, F12,
48 9 F22, F32, F13, F23,
49 A F33, F14, F24, F34,
50 B F15, F25, F35, F16,
51 C F26, F36, F17, F27,
52 D F37, F18, F28, F38,
53 E PX1H1, PX1H2, PX1H3, PX1H4,
54 F PX2H1, PX2H2, PX2H3, PX2H4,
55 G PX3H1, PX3H2, PX3H3, PX3H4,
56 H PX4H1, PX4H2, PX4H3, PX4H4,
57 I VOL, MAT, CXX, VIS,
58 J VD2, DELTAX, PID, GEO,
59 K PARTSAV, IPARTS, DXX, DYY,
60 L DZZ, D4, D5, D6,
61 M FHOUR, JR0, JS0, JT0,
62 N EINT, VOL0, SIGY, SIG0,
63 O SIGOLD, ICP, DEFP, MATVIS,
64 P ET, D_MAX, NEL, GAMA,
65 Q UPARAM, STRHG, STRAIN, ISTRAIN,
66 R MTN, ISMSTR, JLAG, IINT ,
67 S MAT_PARAM )
68C-----------------------------------------------
69C M o d u l e s
70C-----------------------------------------------
71 USE elbufdef_mod
72 use matparam_def_mod
73C-----------------------------------------------
74C I m p l i c i t T y p e s
75C-----------------------------------------------
76#include "implicit_f.inc"
77C-----------------------------------------------
78C G l o b a l P a r a m e t e r s
79C-----------------------------------------------
80#include "mvsiz_p.inc"
81C-----------------------------------------------
82C C o m m o n B l o c k s
83C-----------------------------------------------
84#include "param_c.inc"
85#include "com04_c.inc"
86#include "com08_c.inc"
87#include "scr14_c.inc"
88C-----------------------------------------------
89C D u m m y A r g u m e n t s
90C-----------------------------------------------
91 INTEGER, INTENT(IN) :: MTN
92 INTEGER, INTENT(IN) :: ISMSTR
93 INTEGER, INTENT(IN) :: JLAG
94 INTEGER, INTENT(IN) :: IINT
95 INTEGER NEL,ISTRAIN
96 my_real
97 . PM(NPROPM,*),GEO(NPROPG,*), RHO(*),OFF(*),
98 . VX1(*),VX2(*),VX3(*),VX4(*),VX5(*),VX6(*),VX7(*),VX8(*),
99 . vy1(*),vy2(*),vy3(*),vy4(*),vy5(*),vy6(*),vy7(*),vy8(*),
100 . vz1(*),vz2(*),vz3(*),vz4(*),vz5(*),vz6(*),vz7(*),vz8(*),
101 . f11(*),f21(*),f31(*),f12(*),f22(*),f32(*),
102 . f13(*),f23(*),f33(*),f14(*),f24(*),f34(*),
103 . f15(*),f25(*),f35(*),f16(*),f26(*),f36(*),
104 . f17(*),f27(*),f37(*),f18(*),f28(*),f38(*),
105 . px1h1(*), px1h2(*), px1h3(*), px1h4(*),
106 . px2h1(*), px2h2(*), px2h3(*), px2h4(*),
107 . px3h1(*), px3h2(*), px3h3(*), px3h4(*),
108 . px4h1(*), px4h2(*), px4h3(*), px4h4(*),
109 . partsav(npsav,*),
110 . vol(*),cxx(*),vis(*),vd2(*),deltax(*),
111 . fhour(nel,3,4),jr0(*),js0(*),jt0(*) ,eint(*),
112 . dxx(*), dyy(*), dzz(*), d4(*), d5(*), d6(*) ,
113 . sigy(*) ,sig0(nel,6),vol0(*),sigold(nel,6),defp(*),et(*),
114 . d_max(*),strhg(nel,18),strain(nel,6)
115 my_real
116 . uparam(*),gama(mvsiz,6)
117 INTEGER MAT(*),PID(*),IPARTS(*),ICP,MATVIS
118 TYPE (ELBUF_STRUCT_), TARGET :: ELBUF_STR
119 type(matparam_struct_) , intent(in) :: mat_param
120C-----------------------------------------------
121C L o c a l V a r i a b l e s
122C-----------------------------------------------
123 INTEGER I, MX, J, IET, MT, IPLAST
124 my_real
125 . caq(mvsiz), fcl(mvsiz),deint(mvsiz),
126 . h11(mvsiz), h22(mvsiz), h33(mvsiz),
127 . h12(mvsiz), h13(mvsiz), h23(mvsiz),
128 . hgx1(mvsiz), hgx2(mvsiz), hgx3(mvsiz), hgx4(mvsiz),
129 . hgy1(mvsiz), hgy2(mvsiz), hgy3(mvsiz), hgy4(mvsiz),
130 . hgz1(mvsiz), hgz2(mvsiz), hgz3(mvsiz), hgz4(mvsiz),
131 . vx3478, vx2358, vx1467, vx1256,
132 . vy3478, vy2358, vy1467, vy1256,
133 . vz3478, vz2358, vz1467, vz1256,
134 . vx17, vy17, vz17,
135 . vx28, vy28, vz28,
136 . vx35, vy35, vz35,
137 . vx46, vy46, vz46,
138 . g_3dt(mvsiz),nu,gg(mvsiz),de,ds,dsig(6),
139 . sm1(mvsiz),sm2(mvsiz),smo1(mvsiz),smo2(mvsiz),smo,
140 . jr_1(mvsiz),js_1(mvsiz),jt_1(mvsiz),nfhour(mvsiz,3,4),
141 . dfhour(mvsiz,3,4),dt05,rho0,
142 . nus(mvsiz),nu2(mvsiz),nu4(mvsiz),e0(mvsiz),
143 . e_r,e_s,e_t,fac,fac1,fac2,coefh,hq13p,hq13n,hq24p,hq24n,ff
144 my_real
145 . cc(mvsiz,3,3),cg(mvsiz,3,3),g33(mvsiz,3,3),gm,gmin,dama_g(mvsiz,3)
146C-----------------------------------------------
147 iet =iint
148 coefh = zep9999
149 iplast = elbuf_str%GBUF%G_PLA
150C---- note: r->eta; s->zeta; t->ksi------------
151
152C +++ MAT Visco-----
153! DO I=1,NEL
154! SIG0V(I,1) = SIG0(I,1)
155! SIG0V(I,2) = SIG0(I,2)
156! SIG0V(I,3) = SIG0(I,3)
157! SIG0V(I,4) = SIG0(I,4)
158! SIG0V(I,5) = SIG0(I,5)
159! SIG0V(I,6) = SIG0(I,6)
160! ENDDO
161 mx = mat(1)
162 rho0=pm(1,mx)
163 nu=pm(21,mx)
164 CALL mmodul(1 ,nel ,pm ,mat ,mtn ,
165 . gama ,uparam ,cc ,cg ,g33 , mat_param )
166 DO i=1,nel
167 gm = third*(g33(i,1,1)+g33(i,2,2)+g33(i,3,3))
168 gg(i)=half*rho0*cxx(i)*cxx(i)*(one -two*nu)/(one-nu)
169C----- sometimes GM is too small
170 gmin = gg(i)*em04
171 gg(i)=min(gg(i),gm)
172 gg(i)=max(gg(i),gmin)
173 e0(i)=two*(one+nu)*gg(i)
174 ENDDO
175C +++ MAT Hydro---
176c IF (MTN==3 .OR.MTN==4 .OR.MTN==6) THEN
177c MX = MAT(1)
178c DO I=1,NEL
179c GG(I)=PM(22,MX)
180c GG(I)=MAX(GG(I),VIS(I))
181c ENDDO
182c ENDIF
183C
184 IF (iet > 1 .AND. matvis>0 ) THEN
185 CALL szetfac(1,nel,iet,mtn,et,gg )
186 ELSEIF (matvis==1.AND.ismstr<10) THEN
187 DO i=1,nel
188 ff=third*(dxx(i)+dyy(i)+dzz(i))
189 de =(dxx(i)-ff)*(dxx(i)-ff)+(dyy(i)-ff)*(dyy(i)-ff)+
190 . (dzz(i)-ff)*(dzz(i)-ff) + fourth*(d4(i)*d4(i)+
191 . d5(i)*d5(i)+d6(i)*d6(i))
192 de = de*dt1
193 dsig(1)=sig0(i,1)-sigold(i,1)
194 dsig(2)=sig0(i,2)-sigold(i,2)
195 dsig(3)=sig0(i,3)-sigold(i,3)
196 dsig(4)=sig0(i,4)-sigold(i,4)
197 dsig(5)=sig0(i,5)-sigold(i,5)
198 dsig(6)=sig0(i,6)-sigold(i,6)
199 ff= third*(dsig(1)+dsig(2)+dsig(3))
200 dsig(1)=dsig(1)-ff
201 dsig(2)=dsig(2)-ff
202 dsig(3)=dsig(3)-ff
203 ds =dsig(1)*dsig(1)+dsig(2)*dsig(2)+dsig(3)*dsig(3)+
204 . dsig(4)*dsig(4)+dsig(5)*dsig(5)+dsig(6)*dsig(6)
205 gg(i)=max(fiveem2*gg(i),sqrt(ds/max(de,em30)))
206 ENDDO
207 ENDIF
208C
209 IF(invstr>=35)THEN
210 mt = pid(1)
211 DO i=1,nel
212 caq(i)=fourth*off(i)*geo(13,mt)
213 ENDDO
214 ELSE
215 mx = mat(1)
216 DO i=1,nel
217 caq(i)=fourth*off(i)*pm(4,mx)
218 ENDDO
219 ENDIF
220 DO i=1,nel
221 g_3dt(i)=third*off(i)*gg(i)*dt1
222 ENDDO
223C---------normalized by GG-------
224 CALL mmod_norm(1,nel ,gg ,cc ,cg ,g33 )
225C
226C MATERIAUX NON FLUIDES CXX->SPP
227 IF (iet > 1 ) THEN
228C------now when IHKT=2, w.r.t. Isolid=1, fac=0.2*0.006666 ~0.001
229 DO i=1,nel
230 fcl(i)=onep1*caq(i)*rho(i)*vol(i)**third
231 fcl(i)=zep00666666667*fcl(i)*cxx(i)
232 ENDDO
233 ELSE
234 DO i=1,nel
235 fcl(i)=caq(i)*rho(i)*vol(i)**third
236 fcl(i)=zep00666666667*fcl(i)*cxx(i)
237 ENDDO
238 END IF
239C
240c IF(ICP==1)THEN
241C-------!!!curieux!!!!!!!!
242c DO I=1,NEL
243c NU1(I) =FOUR_OVER_3
244c NU2(I) =-TWO_THIRD
245c NU3(I) =ZEP444
246c NU4(I) =ZERO
247c ENDDO
248c ELSEIF(ICP==2)THEN
249c DO I=1,NEL
250c FAC1 = SIGY(I)/E0(I)+DEFP(I)
251c FAC2=ONE-DEFP(I)/FAC1
252c FF =(ONE +NU)/(ONE -TWO*NU)*FAC2
253c NU1(I) =TWO_THIRD*(TWO+FF)
254c NU2(I) =TWO_THIRD*(FF-ONE)
255c NU3(I) =ZEP222*(TWO+FF)
256c NU4(I) =ZERO
257c ENDDO
258c ELSE
259c DO I=1,NEL
260c NU1(I) =TWO/(ONE-NU)
261c NU2(I) =NU*NU1(I)
262c NU3(I) =TWO_THIRD*(ONE + NU)
263c NU4(I) =NU
264c ENDDO
265c ENDIF
266 IF(icp==1)THEN
267 DO i=1,nel
268 nus(i) =zep499
269 ENDDO
270 ELSEIF(icp==2.AND.iplast>0)THEN
271 DO i=1,nel
272 fac1 = sigy(i)/e0(i)+defp(i)
273 fac2 = defp(i)/fac1
274 nus(i)=nu+(half-nu)*fac2
275 ENDDO
276 ELSE
277 DO i=1,nel
278 nus(i) =nu
279 ENDDO
280 ENDIF
281 DO i=1,nel
282 nu2(i) =nus(i)/(one-nus(i))
283 nu4(i) =nus(i)
284 ENDDO
285 dt05 =half*dt1
286 DO i=1,nel
287 vx3478=vx3(i)-vx4(i)-vx7(i)+vx8(i)
288 vx2358=vx2(i)-vx3(i)-vx5(i)+vx8(i)
289 vx1467=vx1(i)-vx4(i)-vx6(i)+vx7(i)
290 vx1256=vx1(i)-vx2(i)-vx5(i)+vx6(i)
291C
292 vy3478=vy3(i)-vy4(i)-vy7(i)+vy8(i)
293 vy2358=vy2(i)-vy3(i)-vy5(i)+vy8(i)
294 vy1467=vy1(i)-vy4(i)-vy6(i)+vy7(i)
295 vy1256=vy1(i)-vy2(i)-vy5(i)+vy6(i)
296C
297 vz3478=vz3(i)-vz4(i)-vz7(i)+vz8(i)
298 vz2358=vz2(i)-vz3(i)-vz5(i)+vz8(i)
299 vz1467=vz1(i)-vz4(i)-vz6(i)+vz7(i)
300 vz1256=vz1(i)-vz2(i)-vz5(i)+vz6(i)
301!
302 hgx3(i)=(vx1467-vx2358)*one_over_8
303 hgx1(i)=(vx1467+vx2358)*one_over_8
304 hgx2(i)=(vx1256-vx3478)*one_over_8
305 hgx4(i)=-(vx1256+vx3478)*one_over_8
306C
307 hgy3(i)=(vy1467-vy2358)*one_over_8
308 hgy1(i)=(vy1467+vy2358)*one_over_8
309 hgy2(i)=(vy1256-vy3478)*one_over_8
310 hgy4(i)=-(vy1256+vy3478)*one_over_8
311C
312 hgz3(i)=(vz1467-vz2358)*one_over_8
313 hgz1(i)=(vz1467+vz2358)*one_over_8
314 hgz2(i)=(vz1256-vz3478)*one_over_8
315 hgz4(i)=-(vz1256+vz3478)*one_over_8
316 ENDDO
317 DO i=1,nel
318 vx17=vx1(i)-vx7(i)
319 vx28=vx2(i)-vx8(i)
320 vx35=vx3(i)-vx5(i)
321 vx46=vx4(i)-vx6(i)
322 vy17=vy1(i)-vy7(i)
323 vy28=vy2(i)-vy8(i)
324 vy35=vy3(i)-vy5(i)
325 vy46=vy4(i)-vy6(i)
326 vz17=vz1(i)-vz7(i)
327 vz28=vz2(i)-vz8(i)
328 vz35=vz3(i)-vz5(i)
329 vz46=vz4(i)-vz6(i)
330C alpha =1 ->eta zeta
331C 1 1 -1 -1 -1 -1 1 1
332 hgx1(i)= hgx1(i)
333 & -(px1h1(i)*vx17+px2h1(i)*vx28
334 & +px3h1(i)*vx35+px4h1(i)*vx46)
335 hgy1(i)= hgy1(i)
336 & -(px1h1(i)*vy17+px2h1(i)*vy28
337 & +px3h1(i)*vy35+px4h1(i)*vy46)
338 hgz1(i)= hgz1(i)
339 & -(px1h1(i)*vz17+px2h1(i)*vz28
340 & +px3h1(i)*vz35+px4h1(i)*vz46)
341C
342C alpha =2 ->zeta ksi
343C 1 -1 -1 1 -1 1 1 -1
344 hgx2(i)= hgx2(i)
345 & -(px1h2(i)*vx17+px2h2(i)*vx28
346 & +px3h2(i)*vx35+px4h2(i)*vx46)
347 hgy2(i)= hgy2(i)
348 & -(px1h2(i)*vy17+px2h2(i)*vy28
349 & +px3h2(i)*vy35+px4h2(i)*vy46)
350 hgz2(i)= hgz2(i)
351 & -(px1h2(i)*vz17+px2h2(i)*vz28
352 & +px3h2(i)*vz35+px4h2(i)*vz46)
353C alpha =3 ->ksi eta
354C 1 -1 1 -1 1 -1 1 -1
355 hgx3(i)= hgx3(i)
356 & -(px1h3(i)*vx17+px2h3(i)*vx28
357 & +px3h3(i)*vx35+px4h3(i)*vx46)
358 hgy3(i)= hgy3(i)
359 & -(px1h3(i)*vy17+px2h3(i)*vy28
360 & +px3h3(i)*vy35+px4h3(i)*vy46)
361 hgz3(i)= hgz3(i)
362 & -(px1h3(i)*vz17+px2h3(i)*vz28
363 & +px3h3(i)*vz35+px4h3(i)*vz46)
364C
365C alpha =4 ->ksi eta zeta
366C -1 1 -1 1 1 -1 1 -1
367 hgx4(i)= hgx4(i)
368 & -(px1h4(i)*vx17+px2h4(i)*vx28
369 & +px3h4(i)*vx35+px4h4(i)*vx46)
370 hgy4(i)= hgy4(i)
371 & -(px1h4(i)*vy17+px2h4(i)*vy28
372 & +px3h4(i)*vy35+px4h4(i)*vy46)
373 hgz4(i)= hgz4(i)
374 & -(px1h4(i)*vz17+px2h4(i)*vz28
375 & +px3h4(i)*vz35+px4h4(i)*vz46)
376 ENDDO
377!
378 DO i=1,nel
379 jr_1(i) = one/max(em20,jr0(i))
380 js_1(i) = one/max(em20,js0(i))
381 jt_1(i) = one/max(em20,jt0(i))
382 h11(i) = js0(i)*jt0(i)*jr_1(i)
383 h22(i) = jr0(i)*jt0(i)*js_1(i)
384 h33(i) = jr0(i)*js0(i)*jt_1(i)
385 h12(i) = jt0(i)
386 h13(i) = js0(i)
387 h23(i) = jr0(i)
388C SIGY(I) = 1.0E+30
389 ENDDO
390 DO i=1,nel
391 fhour(i,1,1) = fhour(i,1,1)*off(i)
392 fhour(i,1,2) = fhour(i,1,2)*off(i)
393 fhour(i,1,3) = fhour(i,1,3)*off(i)
394 fhour(i,1,4) = fhour(i,1,4)*off(i)
395 fhour(i,2,1) = fhour(i,2,1)*off(i)
396 fhour(i,2,2) = fhour(i,2,2)*off(i)
397 fhour(i,2,3) = fhour(i,2,3)*off(i)
398 fhour(i,2,4) = fhour(i,2,4)*off(i)
399 fhour(i,3,1) = fhour(i,3,1)*off(i)
400 fhour(i,3,2) = fhour(i,3,2)*off(i)
401 fhour(i,3,3) = fhour(i,3,3)*off(i)
402 fhour(i,3,4) = fhour(i,3,4)*off(i)
403 ENDDO
404 IF (iplast==1)
405 . CALL szsvm_or(
406 1 jr0, js0, jt0, cc,
407 2 cg, g33, fhour, sigy,
408 3 sigold, nu, smo1, smo2,
409 4 nel, iint)
410 !-----------For energy calculation------------
411 IF(jlag==1)THEN
412 CALL gfhour_or(1,nel,
413 . fhour,jr0,js0,jt0,fcl,
414 . hgx1, hgx2, hgx3, hgx4,
415 . hgy1, hgy2, hgy3, hgy4,
416 . hgz1, hgz2, hgz3, hgz4,
417 . h11 , h22 , h33 ,
418 . h12 , h13 , h23 ,
419 . jr_1,js_1 , jt_1, nu4,nu2 ,
420 . cc ,cg ,g33 ,nfhour,nel)
421C
422 DO i=1,nel
423 deint(i)=
424 . nfhour(i,3,1)*hgz1(i) + nfhour(i,3,2)*hgz2(i) +
425 . nfhour(i,3,3)*hgz3(i) + nfhour(i,3,4)*hgz4(i) +
426 . nfhour(i,1,1)*hgx1(i) + nfhour(i,1,2)*hgx2(i) +
427 . nfhour(i,1,3)*hgx3(i) + nfhour(i,1,4)*hgx4(i) +
428 . nfhour(i,2,1)*hgy1(i) + nfhour(i,2,2)*hgy2(i) +
429 . nfhour(i,2,3)*hgy3(i) + nfhour(i,2,4)*hgy4(i)
430 eint(i)= eint(i)+dt05*deint(i)/max(em20,vol0(i))
431 ENDDO
432 ENDIF
433C
434 IF (iet > 1 .AND. mtn == 24 ) THEN
435 CALL mdama24(elbuf_str,1,nel ,pm ,mat ,dama_g )
436 DO j=1,3
437 DO i=1,nel
438 fac1=one- dama_g(i,j)
439 fhour(i,j,1:4) = fhour(i,j,1:4)*fac1
440 ENDDO
441 ENDDO
442 END IF !(IET > 1) THEN
443C -------elstic increament----
444 DO i=1,nel
445 e_r =g_3dt(i)*jr_1(i)
446 e_s =g_3dt(i)*js_1(i)
447 e_t =g_3dt(i)*jt_1(i)
448 dfhour(i,1,1) = e_r*hgx1(i)
449 dfhour(i,1,2) = e_r*hgx2(i)
450 dfhour(i,1,3) = e_r*hgx3(i)
451 dfhour(i,1,4) = e_r*hgx4(i)
452 dfhour(i,2,1) = e_s*hgy1(i)
453 dfhour(i,2,2) = e_s*hgy2(i)
454 dfhour(i,2,3) = e_s*hgy3(i)
455 dfhour(i,2,4) = e_s*hgy4(i)
456 dfhour(i,3,1) = e_t*hgz1(i)
457 dfhour(i,3,2) = e_t*hgz2(i)
458 dfhour(i,3,3) = e_t*hgz3(i)
459 dfhour(i,3,4) = e_t*hgz4(i)
460C
461 fhour(i,1,1) = fhour(i,1,1) + dfhour(i,1,1)
462 fhour(i,1,2) = fhour(i,1,2) + dfhour(i,1,2)
463 fhour(i,1,3) = fhour(i,1,3) + dfhour(i,1,3)
464 fhour(i,1,4) = fhour(i,1,4) + dfhour(i,1,4)
465 fhour(i,2,1) = fhour(i,2,1) + dfhour(i,2,1)
466 fhour(i,2,2) = fhour(i,2,2) + dfhour(i,2,2)
467 fhour(i,2,3) = fhour(i,2,3) + dfhour(i,2,3)
468 fhour(i,2,4) = fhour(i,2,4) + dfhour(i,2,4)
469 fhour(i,3,1) = fhour(i,3,1) + dfhour(i,3,1)
470 fhour(i,3,2) = fhour(i,3,2) + dfhour(i,3,2)
471 fhour(i,3,3) = fhour(i,3,3) + dfhour(i,3,3)
472 fhour(i,3,4) = fhour(i,3,4) + dfhour(i,3,4)
473 ENDDO
474 IF (iplast==1)
475 . CALL szsvm_or(
476 1 jr0, js0, jt0, cc,
477 2 cg, g33, fhour, sigy,
478 3 sig0, nu, sm1, sm2,
479 4 nel, iint)
480 !-----------ELASTIC-PLASTIC yield criterion------------
481 IF (iplast==1) THEN
482 DO i=1,nel
483 IF (sm1(i)>sigy(i).AND.deint(i)>0) THEN
484 smo = zep9*smo1(i)+em01*smo2(i)
485 fac1 = sigy(i)-smo
486 fac2 = sm1(i)-smo
487 IF (fac2<=em20) THEN
488 fac=zero
489 ELSE
490 fac = one - max(em20,fac1/fac2)
491 ENDIF
492 IF (sm2(i)<sigy(i)) THEN
493 fac1 =(sm1(i)-sigy(i))/max((sm1(i)-sm2(i)),em20)
494 fac1 =half + sqrt(fac1)
495 fac = min(fac1,one)*fac
496 ENDIF
497 fhour(i,1,1) = fhour(i,1,1) - fac*dfhour(i,1,1)
498 fhour(i,1,2) = fhour(i,1,2) - fac*dfhour(i,1,2)
499 fhour(i,1,3) = fhour(i,1,3) - fac*dfhour(i,1,3)
500 fhour(i,1,4) = fhour(i,1,4) - fac*dfhour(i,1,4)
501 fhour(i,2,1) = fhour(i,2,1) - fac*dfhour(i,2,1)
502 fhour(i,2,2) = fhour(i,2,2) - fac*dfhour(i,2,2)
503 fhour(i,2,3) = fhour(i,2,3) - fac*dfhour(i,2,3)
504 fhour(i,2,4) = fhour(i,2,4) - fac*dfhour(i,2,4)
505 fhour(i,3,1) = fhour(i,3,1) - fac*dfhour(i,3,1)
506 fhour(i,3,2) = fhour(i,3,2) - fac*dfhour(i,3,2)
507 fhour(i,3,3) = fhour(i,3,3) - fac*dfhour(i,3,3)
508 fhour(i,3,4) = fhour(i,3,4) - fac*dfhour(i,3,4)
509 ENDIF
510 ENDDO
511 ENDIF
512 CALL gfhour_or(1,nel,
513 . fhour,jr0,js0,jt0,fcl,
514 . hgx1, hgx2, hgx3, hgx4,
515 . hgy1, hgy2, hgy3, hgy4,
516 . hgz1, hgz2, hgz3, hgz4,
517 . h11 , h22 , h33 ,
518 . h12 , h13 , h23 ,
519 . jr_1,js_1 , jt_1, nu4,nu2 ,
520 . cc ,cg ,g33 ,nfhour,nel)
521 DO i=1,nel
522 hq13p = (nfhour(i,1,1)+nfhour(i,1,3))*one_over_8
523 hq13n = (nfhour(i,1,1)-nfhour(i,1,3))*one_over_8
524 hq24p = (nfhour(i,1,2)+nfhour(i,1,4))*one_over_8
525 hq24n = (nfhour(i,1,2)-nfhour(i,1,4))*one_over_8
526 ff =-px1h1(i)*nfhour(i,1,1)-px1h2(i)*nfhour(i,1,2)
527 . -px1h3(i)*nfhour(i,1,3)-px1h4(i)*nfhour(i,1,4)
528 f11(i) =-(hq13p+hq24n+ff)
529 f17(i) =-(hq13p+hq24p-ff)
530 ff =-px2h1(i)*nfhour(i,1,1)-px2h2(i)*nfhour(i,1,2)
531 . -px2h3(i)*nfhour(i,1,3)-px2h4(i)*nfhour(i,1,4)
532 f12(i) =-(hq13n-hq24n+ff)
533 f18(i) =-(hq13n-hq24p-ff)
534 ff =-px3h1(i)*nfhour(i,1,1)-px3h2(i)*nfhour(i,1,2)
535 . -px3h3(i)*nfhour(i,1,3)-px3h4(i)*nfhour(i,1,4)
536 f13(i) =-(-hq13n-hq24p+ff)
537 f15(i) =-(-hq13n-hq24n-ff)
538 ff =-px4h1(i)*nfhour(i,1,1)-px4h2(i)*nfhour(i,1,2)
539 . -px4h3(i)*nfhour(i,1,3)-px4h4(i)*nfhour(i,1,4)
540 f14(i) =-(-hq13p+hq24p+ff)
541 f16(i) =-(-hq13p+hq24n-ff)
542 ENDDO
543 DO i=1,nel
544 hq13p = (nfhour(i,2,1)+nfhour(i,2,3))*one_over_8
545 hq13n = (nfhour(i,2,1)-nfhour(i,2,3))*one_over_8
546 hq24p = (nfhour(i,2,2)+nfhour(i,2,4))*one_over_8
547 hq24n = (nfhour(i,2,2)-nfhour(i,2,4))*one_over_8
548 ff =-px1h1(i)*nfhour(i,2,1)-px1h2(i)*nfhour(i,2,2)
549 . -px1h3(i)*nfhour(i,2,3)-px1h4(i)*nfhour(i,2,4)
550 f21(i) =-(hq13p+hq24n+ff)
551 f27(i) =-(hq13p+hq24p-ff)
552 ff =-px2h1(i)*nfhour(i,2,1)-px2h2(i)*nfhour(i,2,2)
553 . -px2h3(i)*nfhour(i,2,3)-px2h4(i)*nfhour(i,2,4)
554 f22(i) =-(hq13n-hq24n+ff)
555 f28(i) =-(hq13n-hq24p-ff)
556 ff =-px3h1(i)*nfhour(i,2,1)-px3h2(i)*nfhour(i,2,2)
557 . -px3h3(i)*nfhour(i,2,3)-px3h4(i)*nfhour(i,2,4)
558 f23(i) =-(-hq13n-hq24p+ff)
559 f25(i) =-(-hq13n-hq24n-ff)
560 ff =-px4h1(i)*nfhour(i,2,1)-px4h2(i)*nfhour(i,2,2)
561 . -px4h3(i)*nfhour(i,2,3)-px4h4(i)*nfhour(i,2,4)
562 f24(i) =-(-hq13p+hq24p+ff)
563 f26(i) =-(-hq13p+hq24n-ff)
564 ENDDO
565 DO i=1,nel
566 hq13p = (nfhour(i,3,1)+nfhour(i,3,3))*one_over_8
567 hq13n = (nfhour(i,3,1)-nfhour(i,3,3))*one_over_8
568 hq24p = (nfhour(i,3,2)+nfhour(i,3,4))*one_over_8
569 hq24n = (nfhour(i,3,2)-nfhour(i,3,4))*one_over_8
570 ff =-px1h1(i)*nfhour(i,3,1)-px1h2(i)*nfhour(i,3,2)
571 . -px1h3(i)*nfhour(i,3,3)-px1h4(i)*nfhour(i,3,4)
572 f31(i) =-(hq13p+hq24n+ff)
573 f37(i) =-(hq13p+hq24p-ff)
574 ff =-px2h1(i)*nfhour(i,3,1)-px2h2(i)*nfhour(i,3,2)
575 . -px2h3(i)*nfhour(i,3,3)-px2h4(i)*nfhour(i,3,4)
576 f32(i) =-(hq13n-hq24n+ff)
577 f38(i) =-(hq13n-hq24p-ff)
578 ff =-px3h1(i)*nfhour(i,3,1)-px3h2(i)*nfhour(i,3,2)
579 . -px3h3(i)*nfhour(i,3,3)-px3h4(i)*nfhour(i,3,4)
580 f33(i) =-(-hq13n-hq24p+ff)
581 f35(i) =-(-hq13n-hq24n-ff)
582 ff =-px4h1(i)*nfhour(i,3,1)-px4h2(i)*nfhour(i,3,2)
583 . -px4h3(i)*nfhour(i,3,3)-px4h4(i)*nfhour(i,3,4)
584 f34(i) =-(-hq13p+hq24p+ff)
585 f36(i) =-(-hq13p+hq24n-ff)
586 ENDDO
587 !----hourglass energy is included in internal energy------
588 IF(jlag==1)THEN
589 DO i=1,nel
590 eint(i)= eint(i)+dt05*(
591 . nfhour(i,3,1)*hgz1(i) + nfhour(i,3,2)*hgz2(i) +
592 . nfhour(i,3,3)*hgz3(i) + nfhour(i,3,4)*hgz4(i) +
593 . nfhour(i,1,1)*hgx1(i) + nfhour(i,1,2)*hgx2(i) +
594 . nfhour(i,1,3)*hgx3(i) + nfhour(i,1,4)*hgx4(i) +
595 . nfhour(i,2,1)*hgy1(i) + nfhour(i,2,2)*hgy2(i) +
596 . nfhour(i,2,3)*hgy3(i) + nfhour(i,2,4)*hgy4(i) )
597 . /max(em20,vol0(i))
598 ENDDO
599 ENDIF
600 IF(istrain>0 .AND.
601 . ((anim_n(iad_gps+400+1) == 1) .OR. (anim_n(iad_gps+400+2) == 1) .OR.
602 . (anim_n(iad_gps+400+3) == 1) .OR. (anim_n(iad_gps+400+4) == 1) .OR.
603 . (anim_n(iad_gps+400+5) == 1) .OR. (anim_n(iad_gps+400+6) == 1)) )THEN
604 CALL szstrainhg(
605 1 jr_1, js_1, jt_1, strhg,
606 2 nel, hgx1, hgx2, hgx3,
607 3 hgx4, hgy1, hgy2, hgy3,
608 4 hgy4, hgz1, hgz2, hgz3,
609 5 hgz4, nu4, nu2)
610c CALL PRSTRAIN(STRAIN,STRHG,CC ,CG ,G33, GG,NEL )
611 ENDIF
612C
613 RETURN
614 END
cg
subroutine cg(dim, mat, rhs, sol, max_iter, tol)
Definition conjugate_gradient.F:29
gfhour_or
subroutine gfhour_or(lft, llt, fhour, jr0, js0, jt0, fcl, hgx1, hgx2, hgx3, hgx4, hgy1, hgy2, hgy3, hgy4, hgz1, hgz2, hgz3, hgz4, h11, h22, h33, h12, h13, h23, jr_1, js_1, jt_1, nu, nu2, cc, cg, g33, nfhour, nel)
Definition gfhour_or.F:37
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
mdama24
subroutine mdama24(elbuf_str, jft, jlt, pm, mat, dama_g)
Definition mdama24.F:38
mmod_norm
subroutine mmod_norm(jft, jlt, gg, cc, cg, g33)
Definition mmod_norm.F:29
mmodul
subroutine mmodul(jft, jlt, pm, mat, mtn, gama, uparam, cc, cg, g33, mat_param)
Definition mmodul.F:37
szetfac
subroutine szetfac(lft, llt, ikt, mtn, et, g)
Definition szetfac.F:31
szhour3_or
subroutine szhour3_or(elbuf_str, pm, rho, off, vx1, vx2, vx3, vx4, vx5, vx6, vx7, vx8, vy1, vy2, vy3, vy4, vy5, vy6, vy7, vy8, vz1, vz2, vz3, vz4, vz5, vz6, vz7, vz8, f11, f21, f31, f12, f22, f32, f13, f23, f33, f14, f24, f34, f15, f25, f35, f16, f26, f36, f17, f27, f37, f18, f28, f38, px1h1, px1h2, px1h3, px1h4, px2h1, px2h2, px2h3, px2h4, px3h1, px3h2, px3h3, px3h4, px4h1, px4h2, px4h3, px4h4, vol, mat, cxx, vis, vd2, deltax, pid, geo, partsav, iparts, dxx, dyy, dzz, d4, d5, d6, fhour, jr0, js0, jt0, eint, vol0, sigy, sig0, sigold, icp, defp, matvis, et, d_max, nel, gama, uparam, strhg, strain, istrain, mtn, ismstr, jlag, iint, mat_param)
Definition szhour3_or.F:68
szstrainhg
subroutine szstrainhg(jr_1, js_1, jt_1, strhg, nel, hgx1, hgx2, hgx3, hgx4, hgy1, hgy2, hgy3, hgy4, hgz1, hgz2, hgz3, hgz4, nu, nu1)
Definition szstrainhg.F:37
szsvm_or
subroutine szsvm_or(jr0, js0, jt0, cc, cg, g33, fhour, sigy, sig0, nu, svm1, svm2, nel, iint)
Definition szsvm_or.F:33