OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
m21law.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!|| m21law ../engine/source/materials/mat/mat021/m21law.F
25!||--- called by ------------------------------------------------------
26!|| mmain ../engine/source/materials/mat_share/mmain.F90
27!||--- calls -----------------------------------------------------
28!|| finter ../engine/source/tools/curve/finter.F
29!||--- uses -----------------------------------------------------
30!|| constant_mod ../common_source/modules/constant_mod.F
31!||====================================================================
32 SUBROUTINE m21law (
33 1 PM ,OFF ,SIG ,EINT ,RHO ,
34 2 EPSQ ,EPXE ,VOL ,MAT ,SSP ,
35 3 DVOL ,VNEW ,D1 ,D2 ,D3 ,
36 4 D4 ,D5 ,D6 ,SOLD1 ,SOLD2 ,
37 5 SOLD3 ,SOLD4 ,SOLD5 ,SOLD6 ,TF ,
38 6 NPF ,SIGY ,DEFP ,IPM ,PNEW ,
39 7 PSH ,MU ,SEQ_OUTPUT,NEL ,NUMMAT ,
40 8 DPLA ,MU_BAK )
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
269 END
alpha
#define alpha
Definition eval.h:35
m21law
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)
Definition m21law.F:41
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