OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
m49law.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!|| m49law ../engine/source/materials/mat/mat049/m49law.f
25!||--- called by ------------------------------------------------------
26!|| mmain ../engine/source/materials/mat_share/mmain.F90
27!||====================================================================
28 SUBROUTINE m49law(MAT ,PM ,OFF ,SIG ,EPXE ,
29 2 THETA ,EPD ,CXX ,DF ,D1 ,
30 3 D2 ,D3 ,D4 ,D5 ,D6 ,
31 4 RHO0 ,DPDM ,SIGY ,DEFP ,DPLA ,
32 5 ESPE ,NEL ,JLAG ,JTHE ,
33 6 FHEAT ,VOL)
34C-----------------------------------------------
35C I m p l i c i t T y p e s
36C-----------------------------------------------
37#include "implicit_f.inc"
38C-----------------------------------------------
39C G l o b a l P a r a m e t e r s
40C-----------------------------------------------
41#include "mvsiz_p.inc"
42C-----------------------------------------------
43C C o m m o n B l o c k s
44C-----------------------------------------------
45#include "com08_c.inc"
46#include "param_c.inc"
47C-----------------------------------------------
48C D u m m y A r g u m e n t s
49C-----------------------------------------------
50 INTEGER MAT(MVSIZ),NEL
51 INTEGER, INTENT(IN) :: JLAG !< flag for lagrangian framework
52 INTEGER, INTENT(IN) :: JTHE !< flag for heat equation
53 my_real PM(NPROPM,*),OFF(MVSIZ) ,SIG(NEL,6),EPXE(MVSIZ),
54 . THETA(MVSIZ),EPD(MVSIZ) ,CXX(MVSIZ) ,DF(MVSIZ) , D1(MVSIZ) ,
55 . d2(mvsiz) ,d3(mvsiz) ,d4(mvsiz) ,d5(mvsiz) , d6(mvsiz) ,
56 . rho0(mvsiz) ,dpdm(mvsiz),sigy(mvsiz) ,defp(mvsiz), dpla(mvsiz),
57 . espe(mvsiz)
58 my_real, DIMENSION(NEL) ,INTENT(INOUT) :: fheat ! Heat due to plastic work for Heat Equation with lagrangian framework
59 my_real, DIMENSION(NEL) ,INTENT(IN) :: vol ! Element Volume
60C-----------------------------------------------
61C L o c a l V a r i a b l e s
62C-----------------------------------------------
63 INTEGER I, MX,ITB, K, J, IFRAC(MVSIZ),IQUAD(MVSIZ), IQUADI, IFLAGR(MVSIZ), NBE
64 my_real g(mvsiz) ,g0 ,g1(mvsiz) ,g2(mvsiz) ,sig0,
65 . cb ,cn ,cb1 ,cb2 ,ch ,
66 . cf ,epmx ,sigmx ,sph ,t0 ,
67 . p(mvsiz) ,dav(mvsiz) ,tmelt ,emelt(mvsiz),qh(mvsiz) ,
68 . qd(mvsiz) ,aj2(mvsiz) ,scale(mvsiz) ,yld(mvsiz)
69 my_real qa, qb, qc, qe, rhocp
70C-----------------------------------------------
71C B o d y
72C-----------------------------------------------
73 mx = mat(1)
74 g0 = pm(22,mx)
75 sig0 = pm(38,mx)
76 cb = pm(39,mx)
77 cn = pm(40,mx)
78 epmx = pm(41,mx)
79 sigmx= pm(42,mx)
80 cb1 = pm(43,mx)
81 cb2 = pm(44,mx)
82 ch = pm(45,mx)
83 sph = pm(69,mx)
84 cf = pm(77,mx)
85 t0 = pm(78,mx)
86 tmelt= pm(46,mx)
87 rhocp= pm(69,mx)
88C
89 DO i=1,nel
90 p(i) =-(sig(i,1)+sig(i,2)+sig(i,3))*third
91 dav(i)=-(d1(i)+d2(i)+d3(i))*third
92 emelt(i) = sph*tmelt
93 ENDDO
94C----------------------------
95C SHEAR MODULUS & YIELD
96C----------------------------
97 DO i=1,nel
98 qa = p(i)*df(i)**third
99 qb = one - ch*(theta(i)-t0)
100 IF(emelt(i)<=zero .OR. cf<=zero) THEN
101 qc = one
102 ELSEIF(espe(i)>=emelt(i)) THEN
103 qc = zero
104 ELSE
105 qc = exp(cf * espe(i) / (espe(i)-emelt(i)))
106 ENDIF
107 g(i) = g0 * (cb1*qa + qb) * qc
108 qd(i) = (cb2*qa + qb) * qc
109 IF(epxe(i)<=zero) THEN
110 qe = sig0
111 ELSEIF(epxe(i)>epmx) THEN
112 qe = sig0 * ((one + cb*epmx)**cn)
113 ELSE
114 qe = sig0 * ((one + cb*epxe(i))**cn)
115 ENDIF
116 yld(i) = min(sigmx, qe) * qd(i)
117 ENDDO
118C--------------------------------
119C DEVIATORIC STRESS ESTIMATE
120C--------------------------------
121 DO i=1,nel
122 g1(i) = dt1*g(i)*off(i)
123 g2(i) = two*g1(i)*off(i)
124 sig(i,1)=sig(i,1)+p(i)+g2(i)*(d1(i)+dav(i))
125 sig(i,2)=sig(i,2)+p(i)+g2(i)*(d2(i)+dav(i))
126 sig(i,3)=sig(i,3)+p(i)+g2(i)*(d3(i)+dav(i))
127 sig(i,4)=sig(i,4)+g1(i)*d4(i)
128 sig(i,5)=sig(i,5)+g1(i)*d5(i)
129 sig(i,6)=sig(i,6)+g1(i)*d6(i)
130 ENDDO
131C--- dP/dRHO
132 DO i=1,nel
133 dpdm(i) = dpdm(i) + four_over_3*g(i)
134 cxx(i)=sqrt(abs(dpdm(i))/rho0(i))
135 ENDDO
136C
137 DO i=1,nel
138 epd(i)=off(i)*max(abs(d1(i)), abs(d2(i)), abs(d3(i)), half*abs(d4(i)),half*abs(d5(i)),half*abs(d6(i)))
139 aj2(i)=half*(sig(i,1)**2+sig(i,2)**2+sig(i,3)**2)+sig(i,4)**2+sig(i,5)**2+sig(i,6)**2
140 aj2(i)=sqrt(three*aj2(i))
141 ENDDO
142C---
143 DO i=1,nel
144C--- Melting
145 IF(theta(i)>=tmelt) THEN
146 qh(i) =zero
147 g(i) =zero
148 yld(i)=zero
149 aj2(i)=zero
150 scale(i)=zero
151 ELSE
152C--- Hardening rule
153 IF(cn>=1) THEN
154 qh(i)= qd(i)*sig0*cb*cn * ((one + cb*epxe(i))**(cn-one))
155 ELSEIF(epxe(i)>zero)THEN
156 qh(i)= qd(i)*sig0*cb*cn / ((one + cb*epxe(i))**(one-cn))
157 ELSE
158 qh(i)=zero
159 ENDIF
160C--- Radial return
161 IF(aj2(i)<=yld(i)) THEN
162 scale(i)=one
163 ELSEIF(aj2(i)/=zero) THEN
164 scale(i)=yld(i)/aj2(i)
165 ELSE
166 scale(i)=zero
167 ENDIF
168 ENDIF
169 ENDDO
170C--- Radial return
171 DO i=1,nel
172C--- plastic strain increment.
173 dpla(i)=(one -scale(i))*aj2(i)/max((three*g(i)+qh(i)),em15)
174C--- actual yield stress.
175 yld(i) =yld(i)+dpla(i)*qh(i)
176 sig(i,1)=scale(i)*sig(i,1)
177 sig(i,2)=scale(i)*sig(i,2)
178 sig(i,3)=scale(i)*sig(i,3)
179 sig(i,4)=scale(i)*sig(i,4)
180 sig(i,5)=scale(i)*sig(i,5)
181 sig(i,6)=scale(i)*sig(i,6)
182 epxe(i) =epxe(i)+dpla(i)
183 epxe(i) =epxe(i)*off(i)
184 ENDDO
185C
186 DO i=1,nel
187 sigy(i)=yld(i)
188 defp(i)=epxe(i)
189 ENDDO
190
191C----------------------------------------------
192C TEMPERATURE (Heating due to plastic work)
193C----------------------------------------------
194 IF (jthe /= 0 .AND. jlag /= 0) THEN
195 DO i=1,nel
196 fheat(i) = fheat(i) + sigy(i)*dpla(i)*vol(i)
197 ENDDO
198 ELSEIF(rhocp > zero)THEN
199 DO i=1,nel
200 theta(i) = theta(i) + sigy(i)*dpla(i) / rhocp
201 ! temperature and internal energy must be incremented consistantly
202 ! internal energy is incremented later in parent subroutine (mmain)
203 ! with total energy deformation which already includes plastic work
204 ! Edef = 0.5 * VOL * sum ( sig.eps_dot , i=1..6)
205 ! so internal energy and temperature remain consistant
206 ENDDO
207 END IF
208
209
210 RETURN
211 END
m49law
subroutine m49law(mat, pm, off, sig, epxe, theta, epd, cxx, df, d1, d2, d3, d4, d5, d6, rho0, dpdm, sigy, defp, dpla, espe, nel, jlag, jthe, fheat, vol)
Definition m49law.F:34
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21