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, scale_sig,
69 . sqdt,aa,bb,cc,fct,fctp,beta,dgt,gs,det,shfact,dfmsa,dfmas
70 my_real, DIMENSION(NEL) ::e,k,emart,g,g2, alpha,epsl,gm,km,
71 . yld_ass,yld_asf,yld_sas,tini,
72 . yld_saf,cas,csa,tsas,tfas, tssa,tfsa,fs,et,gt
73C-----------------------------------------------
74c UVAR 2= FS0
75c UVAR 3= FS0
76C=======================================================================
77
78 shfact = five_over_6
79 sqdt = sqrt(two/three)
80 iadbuf = ipm(7,imat)-1
81 DO i=1,nel
82 e(i) = uparam(iadbuf+1)
83 g(i) = uparam(iadbuf+3)
84 k(i) = uparam(iadbuf+4)
85 yld_ass(i) = uparam(iadbuf+6)
86 yld_asf(i) = uparam(iadbuf+7)
87 yld_sas(i) = uparam(iadbuf+8)
88 yld_saf(i) = uparam(iadbuf+9)
89 alpha(i) = uparam(iadbuf+10)
90 epsl(i) = uparam(iadbuf+11)
91 emart(i) = uparam(iadbuf+14)
92 eflag(i) = int(uparam(iadbuf+15))
93 gm(i) = uparam(iadbuf+16)
94 km(i) = uparam(iadbuf+17)
95 g2(i) = two*g(i)
96 cas(i) = uparam(iadbuf+18)
97 csa(i) = uparam(iadbuf+19)
98 tsas(i) = uparam(iadbuf+20)
99 tfas(i) = uparam(iadbuf+21)
100 tssa(i) = uparam(iadbuf+22)
101 tfsa(i) = uparam(iadbuf+23)
102 tini(i) = uparam(iadbuf+25)
103 ENDDO
104 ! Compute temperature - adiabatic conditions
105 IF (jthe == 0 ) THEN
106 DO i=1,nel
107 temp(i) = tini(i) ! + EINT(I)/ RHO0(I)/CP/MAX(EM15,VOLUME(I))
108 ENDDO
109 ENDIF
110C=======================================================================
111 DO i = 1,nel
112 !compute limits for start and end of transformation
113 rsas = yld_ass(i) -cas(i)*tsas(i) !start from aust to martensite
114 rfas = yld_asf(i) -cas(i)*tfas(i) !end from aust to martensite
115 rssa = yld_sas(i) -csa(i)*tssa(i) !start from martensite to austenite
116 rfsa = yld_saf(i) -csa(i)*tfsa(i) !end from martensite to austenite
117 IF (epsxx(i)< zero)THEN
118 scale_sig = (one + alpha(i))/(one - alpha(i))
119 scale_sig = (sqrt(two_third) + alpha(i))/(sqrt(two_third) - alpha(i))
120 rsas = yld_ass(i) * scale_sig -cas(i)*tsas(i) !start from aust to martensite compression
121 rfas = yld_asf(i) * scale_sig -cas(i)*tfas(i) !end from aust to martensite compression
122 rssa = yld_sas(i) * scale_sig -csa(i)*tssa(i) !start from martensite to austenite compression
123 rfsa = yld_saf(i) * scale_sig -csa(i)*tfsa(i) !end from martensite to austenite compression
124 epsl(i) = epsl(i) / scale_sig
125 ENDIF
126
127 IF (eflag(i) > zero)THEN ! young modulus dependent on martensite fraction
128 gt(i) = g(i) + fm(i) * (gm(i) - g(i)) !G_n
129 et(i) = e(i) + fm(i) * (emart(i) - e(i)) !E_n
130 gs = shfact*gt(i)
131 !trial stress (sigma_tr_n+1)
132 signxx(i) = et(i)*(epsxx(i) - epsl(i)*fm(i)* sign(one,epsxx(i)) )
133 signxy(i) = gs * epsxy(i)
134 signxz(i) = gs * epsxz(i)
135 etse(i) = one
136
137 beta = et(i)*epsl(i)* sign(one,epsxx(i))-(emart(i)- e(i))*epsxx(i)
138 . + (emart(i) - e(i)) *epsl(i)* sign(one,epsxx(i)) *fm(i)
139 dfmsa = zero
140 dfmas = zero
141 !---------------
142 !Check Austenite -----> martensite
143 !---------------
144 fs(i) = abs(signxx(i) ) -cas(i)*temp(i)! trial
145 fass = fs(i) - rsas
146 fasf = fs(i) - rfas
147 fs0 = uvar(i,2)
148 IF((fs(i) - fs0) > zero .AND. fass > zero.AND. fm(i) <= one )THEN
149 aa = (one - fm(i))*(emart(i) - e(i))*epsl(i)
150 bb = -(fs0 - rfas - (one - fm(i))*beta* sign(one,epsxx(i)))
151 cc = -(one - fm(i)) * (fs(i) - fs0)
152 dfmas = min(one , - (one - fm(i))*(fs(i) - fs0 ) / (fs0-rfas - (one - fm(i))*e(i)*epsl(i) ) )
153 DO kk = 1,3
154 fct = dfmas*dfmas *aa+ dfmas* bb + cc
155 fctp = two*dfmas *aa+ bb
156 dfmas = dfmas - fct / fctp
157 ENDDO
158 dfmas = min(one,dfmas )
159 ENDIF
160 !---------------
161 !Check marteniste -----> austenite
162 !---------------
163 fs(i) = abs(signxx(i) ) -csa(i)*temp(i)! trial
164 fsas = fs(i) - rssa
165 fsaf = fs(i) - rfsa
166 fs0 = uvar(i,3)
167 IF((fs(i) - fs0) < zero .AND. fsas < zero .AND. fm(i) >zero)THEN
168 aa = fm(i)*(emart(i) - e(i))*epsl(i)
169 bb = fs0-rfsa + fm(i)*beta* sign(one,epsxx(i))
170 cc = -fm(i) * (fs(i)-fs0 )
171 dfmsa = max(-one , fm(i) * ( fs(i) - fs0 ) / (fs0 - rfsa + fm(i)* e(i)*epsl(i) ))
172 DO kk = 1,3
173 fct = dfmsa*dfmsa *aa+ dfmsa* bb +cc
174 fctp = two*dfmsa *aa+ bb
175 dfmsa = dfmsa - fct / fctp
176 ENDDO
177 dfmsa = max(-one , dfmsa )
178 ENDIF
179 !--------------------------
180 !new martensite fraction increment
181 !--------------------------
182 dfm = dfmas + dfmsa
183
184 IF(dfm < zero .AND. fm(i) == zero) dfm = zero
185 ! --------------------------
186 ! UPDATE
187 ! --------------------------
188 dgt = dfm * (gm(i) - g(i))
189 det = dfm * (emart(i) - e(i))
190 signxx(i) = signxx(i) + det * ( epsxx(i) - epsl(i) * (fm(i) + dfm)* sign(one,epsxx(i)) )
191 . - et(i) * epsl(i) * dfm* sign(one,epsxx(i))
192
193
194 fs(i) = abs(signxx(i) )
195 fm(i) = fm(i) + dfm
196 fm(i) = max(zero,fm(i))
197 fm(i) = min(one ,fm(i))
198 treps(i) = fm(i) * sign(one,epsxx(i)) * epsl(i) !transformation strain for output
199 uvar(i,2) = fs(i) - cas(i)*temp(i)
200 uvar(i,3) = fs(i) - csa(i)*temp(i)
201
202C=======================================================================
203 ELSE ! constant young modulus
204 gs = shfact*g(i)
205 !trial stress(sigma_tr_n+1)
206 signxx(i) = e(i)*(epsxx(i) - epsl(i)*fm(i)* sign(one,epsxx(i))) ! deviator of stress
207 signxy(i) = gs * epsxy(i)
208 signxz(i) = gs * epsxz(i)
209
210 etse(i) = one
211 fs(i) = abs(signxx(i) ) - cas(i)*temp(i)! trial
212 !---------------
213 !Check Austenite -----> martensite
214 !---------------
215 fass = fs(i) - rsas
216 fasf = fs(i) - rfas
217 fs0 = uvar(i,2)
218 dfmsa = zero
219 dfmas = zero
220 IF((fs(i) - fs0) > zero .AND. fass > zero.AND. fm(i) < one )THEN
221 ! transformation to martensite
222 ! compute martensite phase fraction increment
223 dfmas = min(one , - (fs(i) - fs0 )*(one - fm(i)) / (fs0-rfas - (one - fm(i))*e(i)*epsl(i) ) )
224 ENDIF
225 !---------------
226 !Check marteniste -----> austenite
227 !---------------
228 fs(i) = abs(signxx(i) ) - csa(i)*temp(i)! trial
229 fsas = fs(i) - rssa
230 fsaf = fs(i) - rfsa
231 fs0 = uvar(i,3)
232 IF((fs(i) - fs0) < zero .AND. fsas < zero .AND. fm(i) >zero )THEN
233 ! transformation to austenite
234 ! compute martensite phase fraction increment
235 dfmsa = max(-one , fm(i) * ( fs(i) - fs0 ) / (fs0 - rfsa + fm(i) * e(i)*epsl(i) ))
236 ENDIF
237 !--------------------------
238 !new martensite fraction increment
239 !--------------------------
240 dfm = dfmas + dfmsa
241 IF(dfm < zero .AND. fm(i) == zero) dfm = zero
242 !--------------------------
243 !UPDATE
244 !--------------------------
245 signxx(i) = signxx(i) - e(i) * dfm * epsl(i) * sign(one,epsxx(i))
246 fs(i) = abs(signxx(i) )
247c
248 fm(i) = fm(i) + dfm
249 fm(i) = max(zero,fm(i))
250 fm(i) = min(one ,fm(i))
251
252 uvar(i,2) = fs(i) - cas(i)*temp(i)
253 uvar(i,3) = fs(i) - csa(i)*temp(i)
254 treps(i) = fm(i) * epsl(i) * sign(one,epsxx(i)) !transformation strain for output
255 ENDIF ! EFLAG
256 ENDDO
257 RETURN
258 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