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

Go to the source code of this file.

Functions/Subroutines

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)

Function/Subroutine Documentation

◆ m49law()

subroutine m49law ( integer, dimension(mvsiz) mat,
pm,
off,
sig,
epxe,
theta,
epd,
cxx,
df,
d1,
d2,
d3,
d4,
d5,
d6,
rho0,
dpdm,
sigy,
defp,
dpla,
espe,
integer nel,
integer, intent(in) jlag,
integer, intent(in) jthe,
intent(inout) fheat,
intent(in) vol )
Parameters
[in]jlagflag for lagrangian framework
[in]jtheflag for heat equation

Definition at line 28 of file m49law.F.

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
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