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 28 of file m10law.F.

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