OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
m10law.F File Reference
#include "implicit_f.inc"
#include "param_c.inc"
#include "com08_c.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine m10law (pm, off, sig, eint, rho, epsq, epxe, vol, mat, ssp, dvol, vnew, d1, d2, d3, d4, d5, d6, sold1, sold2, sold3, sold4, sold5, sold6, sigy, defp, pnew, psh, mu_new, seq_output, nel, dpdm, dpla, mu_bak)

Function/Subroutine Documentation

◆ m10law()

subroutine m10law ( pm,
dimension(nel), intent(in) off,
sig,
eint,
rho,
dimension(nel), intent(inout) epsq,
epxe,
vol,
integer, dimension(nel), intent(in) mat,
ssp,
dvol,
vnew,
d1,
d2,
d3,
d4,
d5,
d6,
sold1,
sold2,
sold3,
sold4,
sold5,
sold6,
sigy,
defp,
dimension(nel), intent(inout) pnew,
dimension(nel), intent(inout) psh,
mu_new,
seq_output,
integer, intent(in) nel,
dpdm,
dimension(nel), intent(inout) dpla,
dimension(nel), intent(inout) mu_bak )

Definition at line 30 of file m10law.F.

37C-----------------------------------------------
38C D e s c r i p t i o n
39C-----------------------------------------------
40Constitutive relations :
41C YIELD CRITERIA : Drucker-Prager Yield J2=A0+A1*P+A2*P**2
42C EOS : legacy input has compaction eos embedded.
43C
44C F = J2 - A0 + A1*P + A2*P**2
45C If F > 1 then deviatoric tensor is projected on F=1. Otherwise Elastic behavior if F<1
46C Yield surface using scale factor RATIO(I).
47C Pressure is cubic in compression linear in tension.
48C Energy integration is made in MEINT subroutine, but eos is not energy dependent...
49C
50C
51C VARIABLE DEFINITIONS :
52C
53C G0 : YIELD ENVELOPE
54C AJ2 : 2ND INVARIANT FROM DEVIATORIC TENSOR
55C-----------------------------------------------
56C M o d u l e s
57C-----------------------------------------------
58 use debug_mod
59C-----------------------------------------------
60C I m p l i c i t T y p e s
61C-----------------------------------------------
62#include "implicit_f.inc"
63C-----------------------------------------------
64C C o m m o n B l o c k s
65C-----------------------------------------------
66#include "param_c.inc"
67#include "com08_c.inc"
68C-----------------------------------------------
69C D u m m y A r g u m e n t s
70C-----------------------------------------------
71 INTEGER,INTENT(IN) :: NEL
72 INTEGER,INTENT(IN) :: MAT(NEL)
73 my_real pm(npropm,*), sig(nel,6), mu_old(nel), epxe(nel), eint(nel), rho(nel), vol(nel),seq_output(nel)
74 my_real,INTENT(INOUT) :: pnew(nel)
75 my_real,INTENT(INOUT) :: psh(nel)
76 my_real,INTENT(IN) :: off(nel)
77 my_real,INTENT(INOUT) :: dpla(nel)
78 my_real,INTENT(INOUT) :: epsq(nel)
79 my_real,INTENT(INOUT) :: mu_bak(nel)
80 my_real vnew(nel), ssp(nel), sigy(nel),defp(nel),
81 . d1(nel), d2(nel), d3(nel), d4(nel), d5(nel), d6(nel),
82 . dvol(nel), mu_new(nel),
83 . sold1(nel), sold2(nel), sold3(nel),
84 . sold4(nel), sold5(nel), sold6(nel),
85 . dpdm(nel)
86C-----------------------------------------------
87C L o c a l V a r i a b l e s
88C-----------------------------------------------
89 INTEGER I, MX
90 my_real t1(nel), t2(nel), t3(nel), t4(nel),
91 . t5(nel), t6(nel), pold(nel), p(nel), pne1(nel),
92 . g(nel), bulk(nel), a0(nel), a1(nel),
93 . a2(nel), amx(nel), aj2(nel), g0(nel), gg(nel),
94 . mu2(nel), svrt(nel), ratio(nel),
95 . yield2(nel), g43(nel),
96 . rho0(nel),ptot,pstar(nel),
97 . g43_1,c0_1,c1_1,c2_1,c3_1,
98 . bulk_1,bulk2_1,mu_max_1,psh_1,
99 . pstar_1,a0_1,a1_1,a2_1,amx_1,
100 . rho0_1,pfrac
101
102 integer nc
103 nc = nc_debug
104C-----------------------------------------------
105C D e s c r i p t i o n
106C-----------------------------------------------
107 !----------------------------------------------------------------!
108 ! PARAMETER INITIALIZATION !
109 !----------------------------------------------------------------!
110 mx=mat(1)
111 g43_1 = onep333*pm(22,mx)
112 pstar_1 = pm(44,mx)
113 a0_1 = pm(38,mx)
114 a1_1 = pm(39,mx)
115 a2_1 = pm(40,mx)
116 amx_1 = pm(41,mx)
117 rho0_1 = pm(1,mx)
118 psh_1 = pm(88,mx)
119 pfrac = pm(37,mx)
120 DO i=1,nel
121 g(i) = dt1*pm(22,mx)
122 g43(i) = g43_1
123 gg(i) = two*g(i)
124 psh(i) = psh_1
125 pstar(i) = pstar_1
126 a0(i) = a0_1
127 a1(i) = a1_1
128 a2(i) = a2_1
129 amx(i) = amx_1
130 rho0(i) = rho0_1
131 ENDDO !next I
132
133 !----------------------------------------------------------------!
134 ! STATE INIT. !
135 !----------------------------------------------------------------!
136 DO i=1,nel
137 pold(i)=-third*(sig(i,1)+sig(i,2)+sig(i,3))
138 svrt(i)= third*(d1(i)+d2(i)+d3(i))
139 mu2(i) = mu_new(i) * max(zero,mu_new(i))
140 ENDDO !next I
141
142 !----------------------------------------------------------------!
143 ! TEMPORARY DEVIATORIC STRESS TENSOR : T(1:6) !
144 !----------------------------------------------------------------!
145 DO i=1,nel
146 t1(i)=sig(i,1)+pold(i)
147 t2(i)=sig(i,2)+pold(i)
148 t3(i)=sig(i,3)+pold(i)
149 t4(i)=sig(i,4)
150 t5(i)=sig(i,5)
151 t6(i)=sig(i,6)
152 ENDDO !next I
153
154 !----------------------------------------------------------------!
155 ! SOUND SPEED !
156 !----------------------------------------------------------------!
157 DO i=1,nel
158 dpdm(i) = g43(i) + dpdm(i)
159 ssp(i) = sqrt(abs(dpdm(i))/rho0(i))
160 ENDDO !next I
161
162 !----------------------------------------------------------------!
163 ! DEVIATORIC TENSOR - ELASTIC INCREMENT !
164 !----------------------------------------------------------------!
165 DO i=1,nel
166 t1(i)=t1(i)+gg(i)*(d1(i)-svrt(i))
167 t2(i)=t2(i)+gg(i)*(d2(i)-svrt(i))
168 t3(i)=t3(i)+gg(i)*(d3(i)-svrt(i))
169 t4(i)=t4(i)+g(i)*d4(i)
170 t5(i)=t5(i)+g(i)*d5(i)
171 t6(i)=t6(i)+g(i)*d6(i)
172 ENDDO !next I
173
174 !----------------------------------------------------------------!
175 ! YIELD SURFACE !
176 !----------------------------------------------------------------!
177 DO i=1,nel
178 aj2(i)= half*(t1(i)**2+t2(i)**2+t3(i)**2)+t4(i)**2+t5(i)**2+t6(i)**2
179 ptot = pnew(i)+psh(i)
180 g0(i) = a0(i)+a1(i)*ptot+a2(i)*ptot*ptot
181 g0(i) = min(amx(i),g0(i))
182 g0(i) = max(zero,g0(i))
183 IF(pnew(i)<=pfrac)g0(i)=zero
184 IF(ptot <= pstar(i))g0(i)=zero
185 yield2(i)=aj2(i)-g0(i)
186 ENDDO !next I
187
188 !----------------------------------------------------------------!
189 ! PROJECTION FACTOR ON YIELD SURFACE !
190 !----------------------------------------------------------------!
191 DO i=1,nel
192 ratio(i)=zero
193 IF(yield2(i)<=zero .AND. g0(i)>zero)THEN
194 ratio(i)=one
195 ELSE
196 ratio(i)=sqrt(g0(i)/(aj2(i)+ em14))
197 ENDIF
198 ENDDO !next I
199
200 !----------------------------------------------------------------!
201 ! UPDATE DEVIATORIC STRESS TENSOR IN SIG(:,:) !
202 !----------------------------------------------------------------!
203 DO i=1,nel
204 sig(i,1)=ratio(i)*t1(i)*off(i)
205 sig(i,2)=ratio(i)*t2(i)*off(i)
206 sig(i,3)=ratio(i)*t3(i)*off(i)
207 sig(i,4)=ratio(i)*t4(i)*off(i)
208 sig(i,5)=ratio(i)*t5(i)*off(i)
209 sig(i,6)=ratio(i)*t6(i)*off(i)
210 dpla(i) = (one -ratio(i))*sqrt(three*abs(aj2(i)))*dt1 / max(em20,three*g(i)) !G <- G*DT1
211 ENDDO !next I
212
213 !----------------------------------------------------------------!
214 ! OUTPUT / MISC. !
215 !----------------------------------------------------------------!
216 DO i=1,nel
217 sigy(i) = g0(i) !YIELD SURFACE
218 epxe(i) = epxe(i) + dpla(i)
219 defp(i) = epxe(i)
220 epsq(i) = mu_bak(i) ! updated if compaction EoS defined
221 ENDDO !next I
222
223 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
debug_mod
Definition debug_mod.F:45
debug_mod::nc_debug
integer nc_debug
Engine Cycle number.
Definition debug_mod.F:49