OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps71pi.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!|| sigeps71pi ../engine/source/materials/mat/mat071/sigeps71pi.F
25!||--- called by ------------------------------------------------------
26!|| mulaw_ib ../engine/source/elements/beam/mulaw_ib.F
27!||====================================================================
28 SUBROUTINE sigeps71pi(
29 1 NEL ,NUPARAM ,UPARAM ,IPM ,IMAT ,
30 2 OFF ,DEPSXX ,DEPSXY ,DEPSXZ ,
31 3 SIGOXX ,SIGOXY ,SIGOXZ ,EPSXX ,EPSXY ,
32 4 EPSXZ ,SIGNXX ,SIGNXY ,SIGNXZ ,ETSE ,
33 5 NUVAR ,UVAR ,JTHE ,TEMP ,FM ,
34 7 TREPS)
35c Law for SMA (Shape memory alloy - NiTinol)
36c based on Auricchio 1997
37C-----------------------------------------------
38C I m p l i c i t T y p e s
39C-----------------------------------------------
40#include "implicit_f.inc"
41C---------+---------+---+---+--------------------------------------------
42C C o m m o n B l o c k s
43C-----------------------------------------------
44#include "param_c.inc"
45#include "com04_c.inc"
46C-----------------------------------------------
47C I N P U T A r g u m e n t s
48C-----------------------------------------------
49 INTEGER ,INTENT(IN) :: IMAT
50 INTEGER ,INTENT(IN) :: NEL,NUPARAM,NUVAR,JTHE
51 INTEGER ,DIMENSION(NPROPMI,NUMMAT) ,INTENT(IN) :: IPM
52 my_real ,DIMENSION(*) ,INTENT(IN) :: UPARAM
53 my_real ,DIMENSION(NEL) ,INTENT(IN) :: EPSXX,EPSXY,EPSXZ,
54 . DEPSXX,DEPSXY,DEPSXZ,SIGOXX,SIGOXY,SIGOXZ
55 my_real ,DIMENSION(NEL) ,INTENT(INOUT) :: temp
56C-----------------------------------------------
57C O U T P U T A r g u m e n t s
58C-----------------------------------------------
59 my_real ,DIMENSION(NEL) ,INTENT(OUT) :: signxx,signxy,signxz,etse
60 my_real ,DIMENSION(NEL) ,INTENT(INOUT) :: off
61 my_real ,DIMENSION(NEL,NUVAR) ,INTENT(INOUT):: uvar
62 my_real fm(nel),treps(nel)
63C-----------------------------------------------
64C L o c a l V a r i a b l e s
65C-----------------------------------------------
66 INTEGER :: I,IADBUF,KK
67 INTEGER ,DIMENSION(NEL) ::EFLAG
68 my_real :: fass,fsas,fasf,fsaf,rsas,rfas,fs0,rssa,rfsa,dfm, fss,scale_sig,
69 . sqdt,aa,bb,cc,fct,fctp,beta,dgt,gs,det,shfact,dfmsa,dfmas ,
70 . delta, sol1,sol2
71 my_real, DIMENSION(NEL) ::e,k,emart,g,g2, alpha,epsl,gm,km,
72 . yld_ass,yld_asf,yld_sas,tini,
73 . yld_saf,cas,csa,tsas,tfas, tssa,tfsa,fs,et,gt
74C-----------------------------------------------
75c UVAR 2= FS0
76c UVAR 3= FS0
77C=======================================================================
78
79 shfact = five_over_6
80 sqdt = sqrt(two/three)
81 iadbuf = ipm(7,imat)-1
82 DO i=1,nel
83 e(i) = uparam(iadbuf+1)
84 g(i) = uparam(iadbuf+3)
85 k(i) = uparam(iadbuf+4)
86 yld_ass(i) = uparam(iadbuf+6)
87 yld_asf(i) = uparam(iadbuf+7)
88 yld_sas(i) = uparam(iadbuf+8)
89 yld_saf(i) = uparam(iadbuf+9)
90 alpha(i) = uparam(iadbuf+10)
91 epsl(i) = uparam(iadbuf+11)
92 emart(i) = uparam(iadbuf+14)
93 eflag(i) = int(uparam(iadbuf+15))
94 gm(i) = uparam(iadbuf+16)
95 km(i) = uparam(iadbuf+17)
96 g2(i) = two*g(i)
97 cas(i) = uparam(iadbuf+18)
98 csa(i) = uparam(iadbuf+19)
99 tsas(i) = uparam(iadbuf+20)
100 tfas(i) = uparam(iadbuf+21)
101 tssa(i) = uparam(iadbuf+22)
102 tfsa(i) = uparam(iadbuf+23)
103 tini(i) = uparam(iadbuf+25)
104 ENDDO
105 ! Compute temperature - adiabatic conditions
106 IF (jthe == 0 ) THEN
107 DO i=1,nel
108 temp(i) = tini(i) ! + EINT(I)/ RHO0(I)/CP/MAX(EM15,VOLUME(I))
109 ENDDO
110 ENDIF
111C=======================================================================
112 DO i = 1,nel
113 !compute limits for start and end of transformation
114 rsas = yld_ass(i) -cas(i)*tsas(i) !start from aust to martensite
115 rfas = yld_asf(i) -cas(i)*tfas(i) !end from aust to martensite
116 rssa = yld_sas(i) -csa(i)*tssa(i) !start from martensite to austenite
117 rfsa = yld_saf(i) -csa(i)*tfsa(i) !end from martensite to austenite
118 IF (epsxx(i)< zero)THEN
119 scale_sig = (one + alpha(i))/(one - alpha(i))
120 scale_sig = (sqrt(two_third) + alpha(i))/(sqrt(two_third) - alpha(i))
121 rsas = yld_ass(i) * scale_sig -cas(i)*tsas(i) !start from aust to martensite compression
122 rfas = yld_asf(i) * scale_sig -cas(i)*tfas(i) !end from aust to martensite compression
123 rssa = yld_sas(i) * scale_sig -csa(i)*tssa(i) !start from martensite to austenite compression
124 rfsa = yld_saf(i) * scale_sig -csa(i)*tfsa(i) !end from martensite to austenite compression
125 epsl(i) = epsl(i) / scale_sig
126 ENDIF
127
128 IF (eflag(i) > zero)THEN ! young modulus dependent on martensite fraction
129 gt(i) = g(i) + fm(i) * (gm(i) - g(i)) !G_n
130 et(i) = e(i) + fm(i) * (emart(i) - e(i)) !E_n
131 gs = shfact*gt(i)
132 !trial stress (sigma_tr_n+1)
133 signxx(i) = et(i)*(epsxx(i) - epsl(i)*fm(i)* sign(one,epsxx(i)) )
134 signxy(i) = gs * epsxy(i)
135 signxz(i) = gs * epsxz(i)
136 etse(i) = one
137
138 beta = et(i)*epsl(i)* sign(one,epsxx(i))-(emart(i)- e(i))*epsxx(i)
139 . + (emart(i) - e(i)) *epsl(i)* sign(one,epsxx(i)) *fm(i)
140 dfmsa = zero
141 dfmas = zero
142 !---------------
143 !Check Austenite -----> martensite
144 !---------------
145 fs(i) = abs(signxx(i) ) -cas(i)*temp(i)! trial
146 fass = fs(i) - rsas
147 fasf = fs(i) - rfas
148 fs0 = uvar(i,2)
149 IF((fs(i) - fs0) > zero .AND. fass > zero.AND. fm(i) <= one )THEN
150 aa = (one - fm(i))*(emart(i) - e(i))*epsl(i)
151 bb = -(fs0 - rfas - (one - fm(i))*beta* sign(one,epsxx(i)))
152 cc = -(one - fm(i)) * (fs(i) - fs0)
153 dfmas = min(one , - (one - fm(i))*(fs(i) - fs0 ) / (fs0-rfas - (one - fm(i))*e(i)*epsl(i) ) )
154 DO kk = 1,3
155 fct = dfmas*dfmas *aa+ dfmas* bb + cc
156 fctp = two*dfmas *aa+ bb
157 dfmas = dfmas - fct / fctp
158 ENDDO
159 dfmas = min(one,dfmas )
160 ENDIF
161 !---------------
162 !Check marteniste -----> austenite
163 !---------------
164 fs(i) = abs(signxx(i) ) -csa(i)*temp(i)! trial
165 fsas = fs(i) - rssa
166 fsaf = fs(i) - rfsa
167 fs0 = uvar(i,3)
168 IF((fs(i) - fs0) < zero .AND. fsas < zero .AND. fm(i) >zero)THEN
169 aa = fm(i)*(emart(i) - e(i))*epsl(i)
170 bb = fs0-rfsa + fm(i)*beta* sign(one,epsxx(i))
171 cc = -fm(i) * (fs(i)-fs0 )
172 dfmsa = max(-one , fm(i) * ( fs(i) - fs0 ) / (fs0 - rfsa + fm(i)* e(i)*epsl(i) ))
173 DO kk = 1,3
174 fct = dfmsa*dfmsa *aa+ dfmsa* bb +cc
175 fctp = two*dfmsa *aa+ bb
176 dfmsa = dfmsa - fct / fctp
177 ENDDO
178 dfmsa = max(-one , dfmsa )
179 ENDIF
180 !--------------------------
181 !new martensite fraction increment
182 !--------------------------
183 dfm = dfmas + dfmsa
184
185 IF(dfm < zero .AND. fm(i) == zero) dfm = zero
186 ! --------------------------
187 ! UPDATE
188 ! --------------------------
189 dgt = dfm * (gm(i) - g(i))
190 det = dfm * (emart(i) - e(i))
191 signxx(i) = signxx(i) + det * ( epsxx(i) - epsl(i) * (fm(i) + dfm)* sign(one,epsxx(i)) )
192 . - et(i) * epsl(i) * dfm* sign(one,epsxx(i))
193
194
195 fs(i) = abs(signxx(i) )
196 fm(i) = fm(i) + dfm
197 fm(i) = max(zero,fm(i))
198 fm(i) = min(one ,fm(i))
199 treps(i) = fm(i) * sign(one,epsxx(i)) * epsl(i) !transformation strain for output
200 uvar(i,2) = fs(i) - cas(i)*temp(i)
201 uvar(i,3) = fs(i) - csa(i)*temp(i)
202
203C=======================================================================
204 ELSE ! constant young modulus
205 gs = shfact*g(i)
206 !trial stress(sigma_tr_n+1)
207 signxx(i) = e(i)*(epsxx(i) - epsl(i)*fm(i)* sign(one,epsxx(i))) ! deviator of stress
208 signxy(i) = gs * epsxy(i)
209 signxz(i) = gs * epsxz(i)
210
211 etse(i) = one
212 fs(i) = abs(signxx(i) ) - cas(i)*temp(i)! trial
213 !---------------
214 !Check Austenite -----> martensite
215 !---------------
216 fass = fs(i) - rsas
217 fasf = fs(i) - rfas
218 fs0 = uvar(i,2)
219 dfmsa = zero
220 dfmas = zero
221 IF((fs(i) - fs0) > zero .AND. fass > zero.AND. fm(i) < one )THEN
222 ! transformation to martensite
223 ! compute martensite phase fraction increment
224 dfmas = min(one , - (fs(i) - fs0 )*(one - fm(i)) / (fs0-rfas - (one - fm(i))*e(i)*epsl(i) ) )
225 ENDIF
226 !---------------
227 !Check marteniste -----> austenite
228 !---------------
229 fs(i) = abs(signxx(i) ) - csa(i)*temp(i)! trial
230 fsas = fs(i) - rssa
231 fsaf = fs(i) - rfsa
232 fs0 = uvar(i,3)
233 IF((fs(i) - fs0) < zero .AND. fsas < zero .AND. fm(i) >zero )THEN
234 ! transformation to austenite
235 ! compute martensite phase fraction increment
236 dfmsa = max(-one , fm(i) * ( fs(i) - fs0 ) / (fs0 - rfsa + fm(i) * e(i)*epsl(i) ))
237 ENDIF
238 !--------------------------
239 !new martensite fraction increment
240 !--------------------------
241 dfm = dfmas + dfmsa
242 IF(dfm < zero .AND. fm(i) == zero) dfm = zero
243 !--------------------------
244 !UPDATE
245 !--------------------------
246 signxx(i) = signxx(i) - e(i) * dfm * epsl(i) * sign(one,epsxx(i))
247 fs(i) = abs(signxx(i) )
248c
249 fm(i) = fm(i) + dfm
250 fm(i) = max(zero,fm(i))
251 fm(i) = min(one ,fm(i))
252
253 uvar(i,2) = fs(i) - cas(i)*temp(i)
254 uvar(i,3) = fs(i) - csa(i)*temp(i)
255 treps(i) = fm(i) * epsl(i) * sign(one,epsxx(i)) !transformation strain for output
256 ENDIF ! EFLAG
257 ENDDO
258 RETURN
259 END
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
sigeps71pi
subroutine sigeps71pi(nel, nuparam, uparam, ipm, imat, off, depsxx, depsxy, depsxz, sigoxx, sigoxy, sigoxz, epsxx, epsxy, epsxz, signxx, signxy, signxz, etse, nuvar, uvar, jthe, temp, fm, treps)
Definition sigeps71pi.F:35