37
38
39
40#include "implicit_f.inc"
41
42
43
44#include "param_c.inc"
45#include "scr17_c.inc"
46#include "com04_c.inc"
47#include "com08_c.inc"
48#include "units_c.inc"
49#include "comlock.inc"
50
51
52
53 INTEGER ,INTENT(IN) :: IMAT,NGL(NEL),IPT
54 INTEGER ,INTENT(IN) :: NEL,NUPARAM,NUVAR,
55 . NFUNC,IFUNC(NFUNC),NPF(*),NVARTMP
56 INTEGER ,DIMENSION(NPROPMI,NUMMAT) ,INTENT(IN) :: IPM
57 my_real ,
DIMENSION(NEL) ,
INTENT(IN) :: epst,
58 . depsxx,depsxy,depsxz,sigoxx,sigoxy,sigoxz
59 my_real ,
DIMENSION(*) ,
INTENT(IN) :: uparam
61 . tf(*)
62
63
64
65 my_real ,
DIMENSION(NEL) ,
INTENT(OUT) :: signxx,signxy,signxz,etse
66 my_real ,
DIMENSION(NEL) ,
INTENT(INOUT) :: pla,off,epsp
67 INTEGER :: VARTMP(NEL,NVARTMP)
68 my_real ,
DIMENSION(NEL,NUVAR) ,
INTENT(INOUT):: uvar
69 my_real ,
DIMENSION(NEL) ,
INTENT(INOUT) :: sigy
70
71
72
73 INTEGER :: I,IADBUF
75 INTEGER, DIMENSION(NEL) :: ICC,ISRATE,VFLAG,IAD,IPOS,ILEN
76 my_real,
DIMENSION(NEL) :: e,nu,g,g3,yld,yldmax,epmax,epdr,
77 . epsr1,epsr2,ca,cb,cn,cp,asrate,yscale,dfdpla,dpla
78
79 shfact = five_over_6
80 epif = zero
81
82 iadbuf = ipm(7,imat)-1
83 DO i=1,nel
84 e(i) = uparam(iadbuf+1)
85 nu(i) = uparam(iadbuf+2)
86 ca(i) = uparam(iadbuf+3)
87 yldmax(i)= uparam(iadbuf+4)
88 epmax(i) = uparam(iadbuf+5)
89 epsr1(i) = uparam(iadbuf+6)
90 epsr2(i) = uparam(iadbuf+7)
91 cb(i) = uparam(iadbuf+8)
92 cn(i) = uparam(iadbuf+9)
93 icc(i) = nint(uparam(iadbuf+10))
94 epdr(i) = uparam(iadbuf+11)
95 epif =
max(epif,epdr(i))
96 cp(i) = uparam(iadbuf+12)
97 g(i) = uparam(iadbuf+16)
98 g3(i) = uparam(iadbuf+18)
99 israte(i)= nint(uparam(iadbuf+13))
100 asrate(i)= uparam(iadbuf+14)
101 vflag(i) = nint(uparam(iadbuf+23))
102 yscale(i)= uparam(iadbuf+24)
103 IF (vflag(i) == 1) THEN
104 epsp(i) = uvar(i,1)
105 ENDIF
106 dpla(i) = zero
107 ENDDO
108
109
110
111 DO i = 1,nel
112 gs = shfact*g(i)
113 signxx(i) = sigoxx(i) + e(i)*depsxx(i)
114 signxy(i) = sigoxy(i) + gs*depsxy(i)
115 signxz(i) = sigoxz(i) + gs*depsxz(i)
116 etse(i) = one
117 ENDDO
118
119
120
121 IF (nfunc > 0) THEN
122 ipos(1:nel) = vartmp(1:nel,1)
123 iad(1:nel) = npf(ifunc(1)) / 2 + 1
124 ilen(1:nel) = npf(ifunc(1)+1) / 2 - iad(1:nel) - ipos(1:nel)
125 CALL vinter(tf,iad,ipos,ilen,nel,pla,dfdpla,yld)
126 vartmp(1:nel,1) = ipos(1:nel)
127 ENDIF
128
129
130
131 DO i = 1,nel
132 IF (nfunc > 0) THEN
133 yld(i) = yscale(i)*yld(i)
134 ELSE
135 yld(i) = ca(i)
136 ENDIF
137 ENDDO
138
139
140
141 DO i = 1,nel
142 IF (pla(i) > zero) THEN
143 IF (nfunc > 0) THEN
144 yld(i) = yscale(i)*yld(i)
145 ELSE
146 yld(i) = ca(i) + cb(i)*exp(cn(i)*log(pla(i)))
147 ENDIF
148 ENDIF
149 ENDDO
150
151
152
153 IF (epif > zero) THEN
154 DO i = 1,nel
155 IF (epdr(i) > zero) THEN
156 frate = one + (epsp(i)*epdr(i))**cp(i)
157 IF (icc(i)== 1) yldmax(i) = yldmax(i) * frate
158 IF ((nfunc > 0) .AND. (ca(i) /= zero)) THEN
159 IF (pla(i)>zero) THEN
160 yld(i) = yld(i) + (ca(i) + cb(i)*exp(cn(i)*log(pla(i))))*(frate-one)
161 ELSE
162 yld(i) = yld(i) + ca(i)*(frate-one)
163 ENDIF
164 ELSE
165 yld(i) = yld(i) * frate
166 ENDIF
167 ENDIF
168 ENDDO
169 ENDIF
170
171
172
173 DO i = 1,nel
174 yld(i) =
min(yld(i),yldmax(i))
175 sigy(i)= yld(i)
176 svm = signxx(i)**2 + three*(signxy(i)**2 + signxz(i)**2)
177 IF (svm > yld(i)**2) THEN
178 svm = sqrt(svm)
179 r =
min( one, yld(i) / svm)
180 signxx(i) = signxx(i)*r
181 signxy(i) = signxy(i)*r
182 signxz(i) = signxz(i)*r
183 dpla(i) = off(i)*svm*(one - r) / e(i)
184 pla(i) = pla(i) + off(i)*svm*(one - r) / e(i)
185 ENDIF
186 ENDDO
187
188
189
190 DO i=1,nel
191 IF (off(i) < em01) off(i) = zero
192 IF (off(i) < one) off(i) = off(i)*four_over_5
193 ENDDO
194
195
196
197 DO i = 1,nel
198 IF (off(i) == one) THEN
199 dmg = one
200 IF (epst(i) > epsr1(i)) THEN
201 dmg = (epsr2(i) - epst(i)) / (epsr2(i) - epsr1(i))
203 signxx(i) = signxx(i)*dmg
204 signxy(i) = signxy(i)*dmg
205 signxz(i) = signxz(i)*dmg
206 ENDIF
207
208 IF (dmg == zero .or. pla(i) >= epmax(i)) THEN
209 off(i) = four_over_5
210
211 WRITE(iout, 1000) ngl(i),ipt
212 WRITE(istdo,1000) ngl(i),ipt
213
214 ENDIF
215 IF (vflag(i) == 1) THEN
217 epsdot = dpla(i)/
max(em20,dt1)
218 epsp(i) =
alpha*epsdot + (one -
alpha)*epsp(i)
219 uvar(i,1) = epsp(i)
220 ENDIF
221
222 ENDIF
223 ENDDO
224
225
226 1000 FORMAT(5x,' FAILURE BEAM ELEMENT NUMBER',i3,', INTEGRATION POINT NUMBER ',i3)
227
228 RETURN
subroutine vinter(tf, iad, ipos, ilen, nel, x, dydx, y)