SPEED
Functions/Subroutines
SYSCOMMON_SUBROUTINES.f90 File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine central_difference (ndof, m, m_inv, c, dt, u, u1, u0, p, flag_minv)
 Computes the displacement of the oscillator through central difference scheme.
 
subroutine matinv (a, b, n)
 
subroutine ksteel02 (props, s, e, de, et, statev, spd, yield, ideath)
 

Function/Subroutine Documentation

◆ central_difference()

subroutine central_difference ( integer, intent(in)  ndof,
real*8, dimension(ndof, ndof), intent(in)  m,
real*8, dimension(ndof, ndof), intent(inout)  m_inv,
real*8, dimension(ndof, ndof), intent(in)  c,
real*8, intent(in)  dt,
real*8, dimension(ndof), intent(inout)  u,
real*8, dimension(ndof), intent(in)  u1,
real*8, dimension(ndof), intent(in)  u0,
real*8, dimension(ndof), intent(in)  p,
integer, intent(inout)  flag_minv 
)

Computes the displacement of the oscillator through central difference scheme.

Author
Aline Herlin, Srihari
Date
November, 2020
Version
1.0
Parameters
[in]mmass of the system
[in]cdamping of system
[in]dTtime step for system
[in]U1,U0displacement at time instants n and n-1
[in]Pequivalent force [-Mu''-Fint]
[out]Udisplacement at time instant n+1

Definition at line 31 of file SYSCOMMON_SUBROUTINES.f90.

32
33 implicit none
34
35 integer, intent(in) :: NDOF
36 integer, intent(inout) :: flag_Minv
37 integer :: idof, j
38
39 real*8, dimension(NDOF, NDOF), intent(in) :: m, c
40 real*8, dimension(NDOF), intent(in) :: u1, u0, p
41 real*8, intent(in) :: dt
42
43 real*8, dimension(NDOF, NDOF), intent(inout) :: m_inv
44 real*8, dimension(NDOF), intent(inout) :: u
45 real*8, dimension(NDOF, NDOF) :: m1, m2, m3, m4
46 real*8, dimension(NDOF) :: t
47 real*8 :: x, y
48
49 ! Check if for some cases we need to change damping matrix (C); then we can pre define these variables (m1, m2, m3, m4); may be we dont need M,M1_inv,C matrices later
50 u=0;
51
52 m1=m/dt/dt + c/2./dt !!! coeff for u(n+1) -> 1/M_inv
53 m2=-2.*m/dt/dt !!! coeff for u(n)
54 m3=m/dt/dt-c/2./dt !!! coeff for u(n-1)
55
56 if (flag_minv.eq.0) then
57 if (ndof.gt.1) then
58 call matinv(m1, m4, ndof)
59 elseif (ndof.eq.1) then
60 m4 = 1/m1
61 endif
62 m_inv = m4
63 flag_minv = 1
64 else
65 m4 = m_inv
66 endif
67
68 do idof= 1,ndof
69 x=0; y=0;
70 do j=1, ndof;
71 x = x + m2(idof,j)*u1(j);
72 y = y + m3(idof,j)*u0(j);
73 enddo
74 t(idof) = p(idof) -x -y
75 enddo
76
77 !For SDOF
78 !T=P-m2*U1-m3*U0
79 !U=T/m1
80
81 do idof=1,ndof
82 do j=1,ndof
83 u(idof) = u(idof) + m_inv(idof,j)*t(j)
84 enddo
85 enddo
86
87 return
matinv
subroutine matinv(a, b, n)
Definition SYSCOMMON_SUBROUTINES.f90:92

References matinv().

Referenced by sdof_shear_model().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ ksteel02()

subroutine ksteel02 ( real*8, dimension(10)  props,
real*8  s,
real*8  e,
real*8  de,
real*8  et,
real*8, dimension(11)  statev,
real*8  spd,
integer  yield,
integer  ideath 
)

Definition at line 149 of file SYSCOMMON_SUBROUTINES.f90.

150!
151 implicit none
152 real*8 e0,sy0,eta,mu,gama,esoft,alpha,beta,a_k,omega
153 real*8 emax,emin,ert,srt,erc,src,ehc,eh1,dt,dc,eu
154 real*8 de,s,e,s0,et,e_unload,sign,sy,evs,eve,epeak,smax,max
155 real*8 sres,eres,x,e_slip,s_slip,e_close,s_close,srel,et1
156 real*8 smin,spd,strain_end
157
158 !real*8 mu
159 integer kon, yield, IDeath !ndof
160 !real*8 M(ndof,ndof)!
161 real*8 props(10), statev(11)
162
163 e0 = props(1)
164 sy0 = props(2)
165 eta = props(3)
166 mu = props(4)
167 gama= props(5)
168 esoft=props(6)
169 alpha= props(7)
170 beta = props(8)
171 a_k= props(9)
172 omega= props(10)
173 emax = statev(1) !maximum strain
174 emin = statev(2) !minimum strain
175 ert = statev(3) !strain at load reversal toward tension
176 srt = statev(4)
177 erc = statev(5) !strain at load reversal toward compression
178 src = statev(6)
179 kon = nint(statev(7)) !
180 ehc = statev(8) !effective cummulative hysteresis energy
181 eh1 = statev(9) !hysteresis energy in a half cycle
182 dt = statev(10) !damage index for tension
183 dc = statev(11) !damage index for compression
184
185 eu = mu * sy0/e0 !characteristic ultimate strain
186
187
188 if(omega<=0.) omega=0.5
189 if(a_k<0.) a_k=0.
190 if(eta<=0.) eta=1.d-6;
191 if(esoft>=0.) esoft=-1.d-6
192
193 if (kon.eq.0) then
194 emax = sy0/e0
195 emin = -beta*sy0/e0
196 if (de.ge.0.0) then
197 kon = 1
198 else
199 kon = 2
200 end if
201 else if ((kon.eq.1).and.(de.lt.0.0)) then !Load reversal
202 kon = 2
203 if (s.gt.0.0) then
204 erc = e
205 src = s
206 end if
207 ehc = ehc + eh1 * (erc / eu ) ** 2.0
208 eh1 = 0.0 !a new half cycle is to begin
209 if (e.gt.emax) emax = e
210 else if ((kon.eq.2).and.(de.gt.0.0)) then !Load reversal
211 kon = 1
212 if (s.lt.0.0) then
213 ert = e
214 srt = s
215 endif
216 ehc = ehc + eh1 * (ert / eu ) ** 2.0
217 eh1 = 0.0 !a new half cycle is to begin
218 if (e.lt.emin) emin = e
219 end if
220!
221 s0=s
222 s = s + e0 * de
223 et = e0
224
225 if(a_k>0.) then
226 if(s0>0.) e_unload=e0*(abs(emax/(sy0/e0)))**(-a_k)
227 if(s0<0.) e_unload=e0*(abs(emin/(sy0/e0)))**(-a_k)
228 if(e_unload<0.1*e0) e_unload=0.1*e0
229 else
230 e_unload=e0
231 end if
232
233 if(s0*de<0.) then
234 s = s0 + e_unload * de
235 et = e_unload
236 if(s*s0<0.) then
237 de=de-s0/e_unload
238 s0=1.d-6*sy0*sign(1.d0,s)
239 et=e0
240 end if
241 end if
242
243
244 if ( de .ge. 0.0 .and. s0>=0.) then
245 sy = (1.0 - dt) * sy0
246 !loading envelope
247 ! Hardening
248 if(e+de>sy/e0) then
249 evs = max( sy + ( e + de - sy/e0) * eta * e0, 0.)
250 eve = eta * e0
251 if (s .ge. evs) then
252 s = evs
253 et = eve
254 yield=1
255 end if
256 end if
257 ! Softening
258 epeak=sy/e0+(alpha-1.)*sy/e0/eta
259 if(e+0.5*de>epeak) then
260 evs=max(sy*alpha+esoft*e0*(e+de-epeak),0.0*sy)
261 if(sy*alpha+esoft*e0*(e+de-epeak)<=0.) ideath=1 ! Complete damage
262 eve=esoft*e0
263 if (s .ge. evs) then
264 s = evs
265 et = eve
266 yield=1
267 end if
268 end if
269
270 !reloading envelope
271 smax = max(sy, sy + (emax - sy/e0) * eta * e0)
272 if(emax>epeak) then
273 smax=max(sy*alpha+esoft*e0*(emax-epeak),0.0*sy)
274 end if
275 sres = 0.02 * smax
276 eres = ert - (srt - sres) / e_unload
277
278 x=emax-smax/e0
279 e_slip=gama*emax+(1.-gama)*x
280 s_slip=smax*gama
281 e_close=e_slip*omega
282 s_close=(e_close-eres)/(e_slip-eres) * (s_slip-sres) + sres
283
284 if (eres .le. emax - smax / e0) then
285 if(e+0.5*de<e_close) then
286 srel = (e+de-eres)/(e_slip-eres) * (s_slip-sres) + sres
287 et1=(s_slip-sres)/(e_slip-eres)
288 else
289 srel=(e+de-e_close)/(emax-e_close)*(smax-s_close)+s_close
290 et1 = (smax - s_close) / (emax - e_close)
291 end if
292 if (s .gt. srel) then
293 s = max( srel, 0.)
294 et = et1
295 end if
296 end if
297
298 elseif ( de .lt. 0.0 .and. s0<0. ) then
299 sy = (1.0 - dc) * sy0 *beta
300 !loading envelope
301 ! Hardening
302 if(e+de<-sy/e0) then
303 evs = min(-sy + ( e + de + sy/e0) * eta * e0,-0.0*sy)
304 eve = eta * e0
305 if (s .le. evs) then
306 s = evs
307 et = eve
308 yield=1
309 end if
310 end if
311 ! Softening
312 epeak=-sy/e0-(alpha-1.)*sy/e0/eta
313 if(e+0.5*de<epeak) then
314 evs=min(-sy*alpha+esoft*e0*(e+de-epeak),-0.*sy)
315 if(-sy*alpha+esoft*e0*(e+de-epeak)>=0.) ideath=1
316 eve=esoft*e0
317 if (s .le. evs) then
318 s = evs
319 et = eve
320 yield=1
321 end if
322 end if
323
324 !reloading envelope
325 smin = min(-sy, -sy + (emin + sy/e0) * eta * e0)
326 if(emin<epeak) then
327 smin=min(-sy*alpha+esoft*e0*(emin-epeak),0.)
328 end if
329 sres = 0.02 * smin
330 eres = erc - (src - sres) / e_unload
331
332 x=emin-smin/e0
333 e_slip=gama*emin+(1.-gama)*x
334 s_slip=smin*gama
335 e_close=e_slip*omega
336 s_close=(e_close-eres)/(e_slip-eres) * (s_slip-sres) + sres
337
338 if (eres .ge. emin - smin / e0) then
339 if(e+0.5*de>e_close) then
340 srel = (e+de-eres)/(e_slip-eres) * (s_slip-sres) + sres
341 et1=(s_slip-sres)/(e_slip-eres)
342 else
343 srel=(e+de-e_close)/(emin-e_close)*(smin-s_close)+s_close
344 et1 = (smin - s_close) / (emin - e_close)
345 end if
346 if (s .lt. srel) then
347 s = min(srel, 0.)
348 et = et1
349 end if
350 end if
351 end if
352
353
354 if (et.ne.e0 .and. et.ne. e_unload) then
355 spd = spd + s * de
356 eh1 = eh1 + s * de
357 if ( s .ge. 0.0 ) then
358 dc = min(ehc /(3.0 * beta* sy0 * eu), 0.7)
359 else
360 dt = min(ehc /(3.0 * sy0 * eu), 0.7)
361 end if
362 end if
363
364 x=max(sy0, beta*sy0)
365 epeak=x/e0+(alpha-1.)*x/e0/eta
366 strain_end=epeak+abs(x/(esoft*e0))
367 x=max(abs(emax),abs(emin),abs(e+de))
368
369 if (ideath==1) then
370 s=0;
371 et=1e-6*e0
372 end if
373 !
374 statev(1) = emax
375 statev(2) = emin
376 statev(3) = ert
377 statev(4) = srt
378 statev(5) = erc
379 statev(6) = src
380 statev(7) = kon
381 statev(8) = ehc
382 statev(9) = eh1
383 statev(10) = dt
384 statev(11) = dc
385 return

Referenced by sdof_shear_model().

Here is the caller graph for this function:

◆ matinv()

subroutine matinv ( real(8), dimension(n,n), intent(in)  a,
real(8), dimension(n,n)  b,
integer(4), intent(in)  n 
)

Definition at line 91 of file SYSCOMMON_SUBROUTINES.f90.

92
93 implicit none
94
95 integer(4), intent(in) :: N
96 integer(4) :: IS(N),JS(N)
97 integer(4) :: I, J, K
98
99 real(8), intent (in) :: A(N,N)
100 real(8) :: B(N,N)
101 real(8) :: D, T
102
103 b=a
104 do k=1,n
105 d=0.0d0
106
107 do i=k,n
108 do j=k,n
109 if(abs(b(i,j))>d) then
110 d=abs(b(i,j)); is(k)=i; js(k)=j
111 end if
112 enddo
113 enddo
114
115 do j=1,n
116 t=b(k,j); b(k,j)=b(is(k),j); b(is(k),j)=t
117 enddo
118
119 do i=1,n
120 t=b(i,k); b(i,k)=b(i,js(k)); b(i,js(k))=t
121 enddo
122
123 b(k,k)=1.d0/b(k,k);
124 do j=1,n; if(j.NE.k) b(k,j)=b(k,j)*b(k,k); enddo
125
126 do i=1,n
127 if(i.NE.k) then
128 do j=1,n; if(j.NE.k) b(i,j)=b(i,j)-b(i,k)*b(k,j); enddo
129 endif
130 enddo
131
132 do i=1,n; if(i.NE.k) b(i,k)=-b(i,k)*b(k,k); enddo
133 enddo
134
135 do k=n,1,-1
136 do j=1,n
137 t=b(k,j); b(k,j)=b(js(k),j); b(js(k),j)=t
138 enddo
139 do i=1,n
140 t=b(i,k); b(i,k)=b(i,is(k)); b(i,is(k))=t
141 enddo
142 enddo
143
144 return

Referenced by central_difference(), and make_sdof_system().

Here is the caller graph for this function: