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

Go to the source code of this file.

Functions/Subroutines

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

Function/Subroutine Documentation

◆ m21law()

subroutine m21law ( dimension(npropm,nummat), intent(in) pm,
dimension(nel), intent(in) off,
dimension(nel,6), intent(inout) sig,
dimension(nel), intent(in) eint,
dimension(nel), intent(in) rho,
dimension(nel), intent(inout) epsq,
dimension(nel), intent(inout) epxe,
dimension(nel), intent(in) vol,
integer, dimension(nel), intent(in) mat,
dimension(nel), intent(inout) ssp,
dimension(nel), intent(in) dvol,
dimension(nel), intent(in) vnew,
dimension(nel), intent(in) d1,
dimension(nel), intent(in) d2,
dimension(nel), intent(in) d3,
dimension(nel), intent(in) d4,
dimension(nel), intent(in) d5,
dimension(nel), intent(in) d6,
dimension(nel), intent(in) sold1,
dimension(nel), intent(in) sold2,
dimension(nel), intent(in) sold3,
dimension(nel), intent(in) sold4,
dimension(nel), intent(in) sold5,
dimension(nel), intent(in) sold6,
dimension(*), intent(in) tf,
integer, dimension(*), intent(in) npf,
dimension(nel), intent(inout) sigy,
dimension(nel), intent(inout) defp,
integer, dimension(npropmi,nummat), intent(in) ipm,
dimension(nel), intent(inout) pnew,
dimension(nel), intent(inout) psh,
dimension(nel), intent(in) mu,
dimension(nel), intent(in) seq_output,
integer, intent(in) nel,
integer, intent(in) nummat,
dimension(nel), intent(inout) dpla,
dimension(nel), intent(inout) mu_bak )
Parameters
[in]nummatnumber of aterial laws (array size)

Definition at line 32 of file m21law.F.

41C-----------------------------------------------
42C M o d u l e s
43C-----------------------------------------------
44 USE constant_mod , ONLY : em14, zero, em20, three, third, two, half, one, onep333, ep20
45C-----------------------------------------------
46C I m p l i c i t T y p e s
47C-----------------------------------------------
48 implicit none
49C-----------------------------------------------
50C I n c l u d e d F i l e s
51C-----------------------------------------------
52#include "my_real.inc"
53#include "param_c.inc"
54#include "com08_c.inc"
55C-----------------------------------------------
56C D u m m y A r g u m e n t s
57C-----------------------------------------------
58 INTEGER,INTENT(IN) :: NPF(*),MAT(NEL),IPM(NPROPMI,NUMMAT),NEL
59 my_real,INTENT(IN) :: pm(npropm,nummat)
60 my_real,INTENT(IN) :: off(nel),eint(nel)
61 my_real,INTENT(IN) :: rho(nel), vol(nel), tf(*), seq_output(nel)
62 my_real,INTENT(IN) :: vnew(nel)
63 my_real,INTENT(IN) :: d1(nel), d2(nel), d3(nel), d4(nel), d5(nel), d6(nel)
64 my_real,INTENT(IN) :: dvol(nel), mu(nel)
65 my_real,INTENT(IN) :: sold1(nel), sold2(nel), sold3(nel)
66 my_real,INTENT(IN) :: sold4(nel), sold5(nel), sold6(nel)
67 my_real,INTENT(INOUT) :: mu_bak(nel) !< history for unloading path
68 my_real,INTENT(INOUT) :: psh(nel) !< pressure shift
69 my_real,INTENT(INOUT) :: ssp(nel) !< sound speed
70 my_real,INTENT(INOUT) :: pnew(nel) !< current pressure
71 my_real,INTENT(INOUT) :: sig(nel,6) !< stress tensor
72 my_real,INTENT(INOUT) :: epxe(nel),defp(nel),epsq(nel),dpla(nel)
73 my_real,INTENT(INOUT) :: sigy(nel) !< yield surface
74 INTEGER,INTENT(IN) :: NUMMAT !< number of aterial laws (array size)
75C-----------------------------------------------
76C L o c a l V a r i a b l e s
77C-----------------------------------------------
78 INTEGER I, MX, IFUNC, NPOINT
79 my_real :: dpdm(nel),dpdm0 !< total derivative dP/d(mu)
80 my_real :: t1(nel), t2(nel), t3(nel), t4(nel),t5(nel), t6(nel) !< deviatoric stress tensor
81 my_real :: pold(nel), p(nel), pne1, ptot, pmin !< pressure and derivated values
82 my_real :: aj2(nel) !< yield criteria
83 my_real :: g,gg,g43 !< shear modulus and derivated values
84 my_real :: bmin, bmax !< minimum and maximum unload modulus (tension : bmin , max compaction bmax)
85 my_real :: bulk(nel) !< unload modulus
86 my_real :: mumin, mumax !< volumetric strain for max compaction
87 my_real :: g0(nel),a0(nel),a1(nel),a2(nel),amax, yield2(nel) !< yield criterion parameter
88 my_real :: pstar !< yield criterion root
89 my_real :: facy !< function scale factor
90 my_real :: rho0(nel) !< initial density
91 my_real,EXTERNAL :: finter !< get /FUNCT value
92 my_real :: alpha
93 my_real :: p_
94 my_real :: svrt
95 my_real :: ratio(nel)
96C-----------------------------------------------
97C D e s c r i p t i o n
98C-----------------------------------------------
99C This subroutine uses Drucker-Prager criteria
100C F = J2 - A0 + A1*P + A2*P**2
101C to calculate deviatoric stresses.
102C If F > 1 then deviatoric tensor is projected on
103C Yield surface using scale factor RATIO(I).
104C Pressure is calculated from compaction EoS (using user function)
105C Energy integration is made in MEINT subroutine (warning : EoS not depending on internal energy)
106C
107C VARIABLE DEFINITIONS :
108C
109C G0 : YIELD ENVELOPE
110C AJ2 : 2ND INVARIANT FROM DEVIATORIC TENSOR
111C EPX : OLD MU
112C EPSEQ : VOLUMETRIC PLASTIC STRAIN
113C-----------------------------------------------
114C S o u r c e L i n e s
115C-----------------------------------------------
116 !----------------------------------------------------------------!
117 ! PARAMETER INIT. !
118 !----------------------------------------------------------------!
119 mx = mat(1)
120 g = dt1*pm(22,mx)
121 g43 = onep333*pm(22,mx)
122 gg = two*g
123 bmin = pm(32,mx)
124 bmax = pm(35,mx)
125 mumin = -ep20
126 mumax = pm(36,mx)
127 pmin = pm(37,mx)
128 psh(1:nel) = pm(43,mx)
129 pstar = pm(44,mx)
130 a0(1:nel) = pm(38,mx)
131 a1(1:nel) = pm(39,mx)
132 a2(1:nel) = pm(40,mx)
133 amax = pm(41,mx)
134 rho0(1:nel) = pm(1,mx)
135 facy = pm(42,mx)
136 ifunc = ipm(11,mx)
137
138 !----------------------------------------------------------------!
139 ! STATE INIT. !
140 !----------------------------------------------------------------!
141 DO i=1,nel
142 pold(i)=-third*(sig(i,1)+sig(i,2)+sig(i,3))
143 ENDDO
144
145 !----------------------------------------------------------------!
146 ! DEVIATORIC STRESS TENSOR !
147 !----------------------------------------------------------------!
148 DO i=1,nel
149 t1(i)=sig(i,1)+pold(i)
150 t2(i)=sig(i,2)+pold(i)
151 t3(i)=sig(i,3)+pold(i)
152 t4(i)=sig(i,4)
153 t5(i)=sig(i,5)
154 t6(i)=sig(i,6)
155 ENDDO
156
157 !----------------------------------------------------------------!
158 ! PRESSURE !
159 !----------------------------------------------------------------!
160 DO i=1,nel
161 p(i) = facy*finter(ifunc,mu(i),npf,tf,dpdm(i))
162 p_ = finter(ifunc,mu_bak(i),npf,tf,dpdm0)
163 !linear unload modulus
164 alpha = one
165 if(mumax > zero)then
166 alpha=mu_bak(i)/mumax
167 endif
168 bulk(i) = alpha*bmax+(one-alpha)*bmin
169 pne1 = p_-(mu_bak(i)-mu(i))*bulk(i)
170 if(mu_bak(i) > mumin) p(i) = min(pne1, p(i))
171 p(i) = max(p(i),pmin) *off(i)
172 if(mu(i) > mu_bak(i)) mu_bak(i) = min(mumax,mu(i))
173 ENDDO
174
175 !----------------------------------------------------------------!
176 ! sound speed !
177 !----------------------------------------------------------------!
178 do i=1,nel
179 dpdm(i) = max(bulk(i),dpdm(i))
180 dpdm(i)= g43 + max(bulk(i),dpdm(i))
181 ssp(i)=sqrt(abs(dpdm(i))/rho0(i))
182 enddo !next i
183
184 !----------------------------------------------------------------!
185 ! OUTPUT !
186 !----------------------------------------------------------------!
187 do i=1,nel
188 p(i)=max(pmin,p(i))*off(i)
189 pnew(i) = p(i)-psh(i)
190 enddo !next i
191
192 !----------------------------------------------------------------!
193 ! YIELD CRITERIA IF P > PMIN !
194 !----------------------------------------------------------------!
195 DO i=1,nel
196 IF(p(i) < pmin)THEN
197 a0(i)=zero
198 a1(i)=zero
199 a2(i)=zero
200 ENDIF
201 ENDDO
202
203 !----------------------------------------------------------------!
204 ! DEVIATORIC TENSOR - ELASTIC INCREMENT !
205 !----------------------------------------------------------------!
206 DO i=1,nel
207 svrt = third*(d1(i)+d2(i)+d3(i))
208 t1(i)=t1(i)+gg*(d1(i)-svrt)
209 t2(i)=t2(i)+gg*(d2(i)-svrt)
210 t3(i)=t3(i)+gg*(d3(i)-svrt)
211 t4(i)=t4(i)+g*d4(i)
212 t5(i)=t5(i)+g*d5(i)
213 t6(i)=t6(i)+g*d6(i)
214 ENDDO
215
216 !----------------------------------------------------------------!
217 ! YIELD SURFACE !
218 !----------------------------------------------------------------!
219 DO i=1,nel
220 aj2(i)=half*(t1(i)**2+t2(i)**2+t3(i)**2)+t4(i)**2+t5(i)**2+t6(i)**2
221 ptot = p(i) + psh(i)
222 g0(i) =a0(i)+a1(i)*ptot+a2(i)*ptot*ptot
223 g0(i)= min(amax,g0(i))
224 g0(i)= max(zero,g0(i))
225 IF(ptot <= pstar)g0(i)=zero
226 yield2(i)=aj2(i)-g0(i)
227 ENDDO
228
229 !----------------------------------------------------------------!
230 ! PROJECTION FACTOR ON YIELD SURFACE !
231 !----------------------------------------------------------------!
232 DO i=1,nel
233 ratio(i) = zero
234 IF(yield2(i)<=zero .AND. g0(i)>zero)THEN
235 ratio(i) = one
236 ELSE
237 ratio(i) = sqrt(g0(i)/(aj2(i)+ em14))
238 ENDIF
239 ENDDO
240
241 !----------------------------------------------------------------!
242 ! DEVIATORIC STRESS TENSOR !
243 !----------------------------------------------------------------!
244 DO i=1,nel
245 p(i)=p(i)*off(i)
246 pnew(i) = p(i)
247 sig(i,1)=ratio(i)*t1(i)*off(i)
248 sig(i,2)=ratio(i)*t2(i)*off(i)
249 sig(i,3)=ratio(i)*t3(i)*off(i)
250 sig(i,4)=ratio(i)*t4(i)*off(i)
251 sig(i,5)=ratio(i)*t5(i)*off(i)
252 sig(i,6)=ratio(i)*t6(i)*off(i)
253 dpla(i)=(one-ratio(i))*sqrt(aj2(i))*dt1 / max(em20,three*g)
254 ENDDO
255
256 !----------------------------------------------------------------!
257 ! OUTPUT / MISC !
258 !----------------------------------------------------------------!
259 DO i=1,nel
260 sigy(i) = g0(i) !YIELD SURFACE
261 defp(i) = defp(i) + dpla(i)
262 epxe(i) = defp(i)
263 epsq(i) = mu_bak(i) ! volumetric plastic strain
264 ENDDO
265
266
267 RETURN
268
my_real
#define my_real
Definition cppsort.cpp:32
alpha
#define alpha
Definition eval.h:35
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
constant_mod
Definition constant_mod.F:294