43
44
45
46 USE sensor_mod
47
48
49
50#include "implicit_f.inc"
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93#include "com04_c.inc"
94#include "tabsiz_c.inc"
95
96
97
98 INTEGER ,INTENT(IN) :: NEL
99 INTEGER ,INTENT(IN) :: NFUNC
100 INTEGER ,INTENT(IN) :: NUPARAM
101 INTEGER ,INTENT(IN) :: NIPARAM
102 INTEGER ,INTENT(IN) :: NUVAR
103 INTEGER ,INTENT(IN) :: NSENSOR
104 my_real ,
INTENT(IN) :: time,timestep
105 INTEGER ,DIMENSION(NFUNC) ,INTENT(IN) :: IFUNC
106 INTEGER ,DIMENSION(SNPC) ,INTENT(IN) :: NPF
107 my_real ,
DIMENSION(STF) ,
INTENT(IN) :: tf
108 INTEGER ,DIMENSION(NIPARAM) ,INTENT(IN) :: IPARAM
109 my_real ,
DIMENSION(NUPARAM) ,
INTENT(IN) :: uparam
110 my_real ,
DIMENSION(NEL) ,
INTENT(IN) ::
area,rho0,thkly,offgg,
111 . depsxx,depsyy,depsxy,depsyz,depszx,epsxx,epsyy,epsxy,epsyz,epszx,
112 . sigoxx,sigoyy,sigoxy,sigoyz,sigozx,shf,aldt
113 TYPE (SENSOR_STR_) ,DIMENSION(NSENSOR) ,INTENT(IN) :: SENSOR_TAB
114
115
116
117 my_real ,
DIMENSION(NEL) ,
INTENT(OUT) :: tan_phi,etse,
118 . signxx,signyy,signxy,signyz,signzx,
119 . sigvxx,sigvyy,sigvxy,soundsp,viscmax
120
121
122
123 my_real ,
DIMENSION(NEL,NUVAR) ,
INTENT(INOUT) :: uvar
124
125
126
127 INTEGER :: I,ITER,NITER,ISENS,IFUNCC,IFUNCT,IFUNCS,IFUNF1,IFUNF2
128 INTEGER ,DIMENSION(NEL) :: IADF1,IADF2,IAD1,IAD2,IAD3,ILEN1,ILEN2,ILEN3,
129 . ILENF1,ILENF2,IPOS,IPOS1,IPOS2
130 my_real :: kmax,gmax,gfrot,gsh,stiff,logc,logt,kflex,kflex1,kflex2,
131 . mass,dyc,dyt,dc0,dt0,h0,hc0,ht0,visce,viscg,lc0,lt0,hdc,hdt,
132 . tfrot,sigg,ec2,et2,sinn,tan2,damp,v1,v2,dtinv,etc,ett,trace,
133 . zerostress,dsig,lmin,tstart,phi1,phi2,j11,j12,j21,j22,det
135 my_real ,
DIMENSION(NEL) :: sigc,sigt,ec,et,lc,lt,fn,angle,kg,gxy,
136 . dtang,tfold,cvisc,cvist,flexc,flext,flexf,hc,ht,dc,dt,dcc,dtt,
137 . fc,ft,fpc,fpt,dflxc,dflxt
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160 isens = iparam(1)
161
162 lc0 = one
163 lt0 = one
164 dc0 = uparam(1)
165 dt0 = uparam(2)
166 hc0 = uparam(3)
167 ht0 = uparam(4)
168 kflex = uparam(5)
169 kflex1 = uparam(6)
170 kflex2 = uparam(7)
171 zerostress = uparam(8)
172 kmax = uparam(10)
173 gmax = uparam(11)
174 yfac(1) = uparam(12)
175 yfac(2) = uparam(13)
176 yfac(3) = uparam(14)
177
178 ifuncc = ifunc(1)
179 ifunct = ifunc(2)
180 ifuncs = ifunc(3)
181 ifunf1 = ifunc(4)
182 ifunf2 = ifunc(5)
183
184 gfrot = gmax
185 gsh = shf(1)
186 visce = em02
187 viscg = em02
188
189 IF (zerostress > zero .and. isens > 0) THEN
190 tstart = sensor_tab(isens)%TSTART
191 ELSE
192 tstart = zero
193 ENDIF
194 dtinv = timestep/
max(timestep**2,em20)
195
196 h0 = hc0 + ht0
197 niter = 10
198
199#include "vectorize.inc"
200 DO i = 1,nel
201 tfold(i) = uvar(i,9)
202 mass = rho0(i)*
area(i)*thkly(i)*fourth
203 cvisc(i) = sqrt(mass*kflex1)*dtinv*third
204 cvist(i) = sqrt(mass*kflex2)*dtinv*third
205 ENDDO
206
207
208#include "vectorize.inc"
209 DO I = 1,NEL
210 ETC = UVAR(I,4) + DEPSXX(I)
211 ETT = UVAR(I,5) + DEPSYY(I)
212 UVAR(I,4) = ETC
213 UVAR(I,5) = ETT
214 EC(I) = EXP(ETC) - ONE ! eng strain dir 1
215 ET(I) = EXP(ETT) - ONE ! eng strain dir 2
216 LC(I) = LC0 * (ONE + EC(I)) ! element length dir 1
217 LT(I) = LT0 * (ONE + ET(I)) ! element length dir 2
218 ENDDO
219
220 DO I = 1,NEL
221 BETA(I) = ONE ! increment scale factor for Newton iterations
222 PH01(I) = EP10
223 PH02(I) = EP10
224 ENDDO
225
226#include "vectorize.inc"
227 DO I = 1,NEL
228 IPOS(I) = 1
229 IAD1(I) = NPF(IFUNCC) / 2 + 1
230 IAD2(I) = NPF(IFUNCT) / 2 + 1
231 IAD3(I) = NPF(IFUNCS) / 2 + 1
232 ILEN1(I) = NPF(IFUNCC+1) / 2 - IAD1(I) - IPOS(I)
233 ILEN2(I) = NPF(IFUNCT+1) / 2 - IAD2(I) - IPOS(I)
234 ILEN3(I) = NPF(IFUNCS+1) / 2 - IAD3(I) - IPOS(I)
235 END DO
236 IF (IFUNF1 > 0) THEN
237 DO I = 1,NEL
238 IADF1(I) = NPF(IFUNF1) / 2 + 1
239 ILENF1(I) = NPF(IFUNF1+1) / 2 - IADF1(I) - IPOS(I)
240 END DO
241 END IF
242 IF (IFUNF2 > 0) THEN
243 DO I = 1,NEL
244 IADF2(I) = NPF(IFUNF2) / 2 + 1
245 ILENF2(I) = NPF(IFUNF2+1) / 2 - IADF2(I) - IPOS(I)
246 END DO
247 END IF
248
249
250
251 EPSFC(1:NEL) = UVAR(1:NEL,7) ! eng strain flex dir 1
252 EPSFT(1:NEL) = UVAR(1:NEL,8) ! eng strain flex dir 2
253
254 DO ITER = 1,NITER
255
256 EPSF (1:NEL) =(HC0 * EPSFC(1:NEL) + HT0 * EPSFT(1:NEL)) / H0 ! eng strain coupling flex spring
257 FLEXF(1:NEL) = KFLEX * EPSF(1:NEL)
258 DO I = 1,NEL
259 XC(I) = EPSFC(I)
260 XT(I) = EPSFT(I)
261 HC(I) = HC0 * (EPSFC(I) + ONE) ! length of flex spring dir 1
262 HT(I) = HT0 * (EPSFT(I) + ONE) ! length of flex spring dir 2
263 DC(I) = SQRT(LC(I)**2 + HC(I)**2) ! fiber length in dir 1
264 DT(I) = SQRT(LT(I)**2 + HT(I)**2) ! fiber length in dir 2
265 DCC(I) = DC(I) - DC0
266 DTT(I) = DT(I) - DT0
267 END DO
268 IPOS1(1:NEL) = 1
269 IPOS2(1:NEL) = 1
270
271 CALL VINTER2(TF,IAD1,IPOS1,ILEN1,NEL,DCC,FPC,FC) ! fiber force dir 1
272 CALL VINTER2(TF,IAD2,IPOS2,ILEN2,NEL,DTT,FPT,FT) ! fiber force dir 2
273
274 ! flex force dir 1
275 IF (IFUNF1 > 0) THEN
276 IPOS(1:NEL) = 1
277 CALL VINTER2(TF,IADF1,IPOS,ILENF1,NEL,EPSFC,DFLXC,FLEXC)
278 FLEXC(1:NEL) = FLEXC(1:NEL) * KFLEX1
279 DFLXC(1:NEL) = DFLXC(1:NEL) * KFLEX1
280 ELSE
281 DO I = 1,NEL
282 FLEXC(I) = KFLEX1 * LOG(EPSFC(I) + ONE)
283 DFLXC(I) = KFLEX1 / (EPSFC(I) + ONE)
284 END DO
285 END IF
286
287 IF (IFUNF2 > 0) THEN
288 IPOS(1:NEL) = 1
289 CALL VINTER2(TF,IADF2,IPOS,ILENF2,NEL,EPSFT,DFLXT,FLEXT)
290 FLEXT(1:NEL) = FLEXT(1:NEL) * KFLEX2
291 DFLXT(1:NEL) = DFLXT(1:NEL) * KFLEX2
292 ELSE
293 DO I = 1,NEL
294 FLEXT(I) = KFLEX2 * LOG(EPSFT(I) + ONE)
295 DFLXT(I) = KFLEX2 / (EPSFT(I) + ONE)
296 END DO
297 END IF
298
299#include "vectorize.inc"
300 DO I = 1,NEL
301 DYC = EPSFC(I) - UVAR(I,7)
302 DYT = EPSFT(I) - UVAR(I,8)
303 HDC = HC(I) / DC(I)
304 HDT = HT(I) / DT(I)
305 PHI1 = FC(I) * HDC + FLEXC(I) + FLEXF(I) + CVISC(I)*DYC
306 PHI2 = FT(I) * HDT + FLEXT(I) + FLEXF(I) + CVIST(I)*DYT
307 J12 = KFLEX * HT0 / H0
308 J21 = KFLEX * HC0 / H0
309 J11 = J12 + FPC(I)*HDC*HC0 + DFLXC(I) + CVISC(I)
310 J22 = J21 + FPT(I)*HDT*HT0 + DFLXT(I) + CVIST(I)
311 DET = J11 * J22 - J12 * J21
312
313 EPSFC(I) = EPSFC(I) - BETA(I) * (J22 * PHI1 - J12 * PHI2) / DET
314 EPSFT(I) = EPSFT(I) + BETA(I) * (J12 * PHI1 - J11 * PHI2) / DET
315
316 EPSFC(I) = MAX(EPSFC(I), EM04 - ONE)
317 EPSFT(I) = MAX(EPSFT(I), EM04 - ONE)
318
319.and. IF (ABS(PHI1) > PH01(I) ABS(PHI2) > PH02(I)) THEN
320 EPSFC(I) = XC(I)
321 EPSFT(I) = XT(I)
322 BETA(I) = BETA(I) * HALF
323 BETA(I) = MAX(BETA(I), EM02)
324 END IF
325 PH01(I) = ABS(PHI1)
326 PH02(I) = ABS(PHI2)
327 END DO
328
329 END DO ! ITER
330
331#include "vectorize.inc"
332 DO I = 1,NEL
333 SIGC(I) = FC(I) * LC(I) / DC(I)
334 SIGT(I) = FT(I) * LT(I) / DT(I)
335 UVAR(I,7) = EPSFC(I)
336 UVAR(I,8) = EPSFT(I)
337 UVAR(I,1) = SIGC(I)
338 UVAR(I,2) = SIGT(I)
339 UVAR(I,15) = DC(I)
340 UVAR(I,16) = DT(I)
341 FN(I) = FLEXF(I) ! normal force for friction in shear
342 ENDDO !i,nel
343
344
345
346
347
348
349
350
351#include "vectorize.inc"
352 DO I = 1,NEL
353 TRACE = EXP(EPSXX(I) + EPSYY(I))! =exp(tr)
354 EC2 = MAX(TRACE / (EC(I) + ONE), EM6)
355 ET2 = MAX(TRACE / (ET(I) + ONE), EM6)
356
357 SIGNXX(I) = SIGC(I) / EC2
358 SIGNYY(I) = SIGT(I) / ET2
359 ENDDO !i,nel
360
361
362
363#include "vectorize.inc"
364 DO I = 1,NEL
365 TAN_PHI(I)= DEPSXY(I)
366 ANGLE(I) = ATAN(TAN_PHI(I))*HUNDRED80/PI
367 IPOS(I) = 1
368 ENDDO
369
370 CALL VINTER2(TF,IAD3,IPOS,ILEN3,NEL,ANGLE,GXY,SIGNXY) ! shear stress =f(angle)
371
372#include "vectorize.inc"
373 DO I = 1,NEL
374 TAN2 = TAN_PHI(I)**2
375 KG(I) = GXY(I) * TAN2 * YFAC(3)
376
377 DAMP = SQRT(RHO0(I)*AREA(I)*THKLY(I)*HALF)
378 V1 = VISCE*DAMP*SQRT(KMAX)
379 V2 = VISCE*DAMP*SQRT(KMAX)
380 SIGVXX(I) = DTINV*(DEPSXX(I))*V1
381 SIGVYY(I) = DTINV*(DEPSYY(I))*V2
382
383 IF (FN(I) > ZERO) THEN
384 TFROT = TWO_THIRD*VISCG*FN(I)*(HC0+HT0)/(LC(I)+LT(I))
385 DTANG(I) = DEPSXY(I) - TAN_PHI(I)
386 SIGG = TFOLD(I) + GFROT*DTANG(I)
387 IF (ABS(SIGG) > TFROT) THEN
388 SIGVXY(I) = SIGN(TFROT,SIGG)
389 ELSE
390 SIGVXY(I) = SIGG
391 ENDIF
392 ENDIF
393
394 SINN = TAN_PHI(I) / SQRT(ONE + TAN2)
395 STIFF = KMAX*(ONE+SINN) + GMAX
396 LMIN = MIN(DC(I)/DC0,DT(I)/DT0)*UVAR(I,14)
397 SOUNDSP(I) = SQRT(STIFF/(RHO0(I)))*ALDT(I)/LMIN
398 VISCMAX(I) = MAX(V1,V2)
399 ETSE(I) = ONE
400 ENDDO
401
402#include "vectorize.inc"
403 DO I = 1,NEL
404 UVAR(I,3) = SIGNXY(I) + SIGVXY(I)
405 TAN_PHI(I)= DEPSXY(I)
406 UVAR(I,6) = DEPSXY(I)
407 UVAR(I,9) = SIGVXY(I)
408
409 SIGNYZ(I) = SIGOYZ(I) + GSH * DEPSYZ(I)
410 SIGNZX(I) = SIGOZX(I) + GSH * DEPSZX(I)
411 ENDDO
412
413
414
415 IF (ZEROSTRESS > ZERO)THEN
416 IF (TIME <= TSTART) THEN
417#include "vectorize.inc"
418 DO I = 1,NEL
419 UVAR(I,11) = SIGNXX(I)
420 UVAR(I,12) = SIGNYY(I)
421 UVAR(I,13) = SIGNXY(I)
422 SIGNXX(I) = ZERO
423 SIGNYY(I) = ZERO
424 SIGNXY(I) = ZERO
425 ENDDO
426 ELSE
427#include "vectorize.inc"
428 DO I = 1,NEL
429 DSIG = SIGNXX(I) - SIGOXX(I) - UVAR(I,11)
430.AND. IF((UVAR(I,11) > ZERO)(DSIG < ZERO))THEN
431 UVAR(I,11) = MAX(ZERO,UVAR(I,11)+ZEROSTRESS*DSIG)
432.AND. ELSEIF((UVAR(I,11) < ZERO)(DSIG > ZERO))THEN
433 UVAR(I,11) = MIN(ZERO,UVAR(I,11)+ZEROSTRESS*DSIG)
434 ENDIF
435 DSIG = SIGNYY(I) - SIGOYY(I) - UVAR(I,12)
436.AND. IF((UVAR(I,12) > ZERO)(DSIG < ZERO))THEN
437 UVAR(I,12) = MAX(ZERO,UVAR(I,12)+ZEROSTRESS*DSIG)
438.AND. ELSEIF((UVAR(I,12) < ZERO)(DSIG > ZERO))THEN
439 UVAR(I,12) = MIN(ZERO,UVAR(I,12)+ZEROSTRESS*DSIG)
440 ENDIF
441 DSIG = SIGNXY(I) - SIGOXY(I) - UVAR(I,13)
442.AND. IF((UVAR(I,13) > ZERO)(DSIG < ZERO))THEN
443 UVAR(I,13) = MAX(ZERO,UVAR(I,13)+ZEROSTRESS*DSIG)
444.AND. ELSEIF((UVAR(I,13) < ZERO)(DSIG > ZERO))THEN
445 UVAR(I,13) = MIN(ZERO,UVAR(I,13)+ZEROSTRESS*DSIG)
446 ENDIF
447 SIGNXX(I) = SIGNXX(I) - UVAR(I,11)
448 SIGNYY(I) = SIGNYY(I) - UVAR(I,12)
449 SIGNXY(I) = SIGNXY(I) - UVAR(I,13)
450 ENDDO
451 ENDIF
452 ENDIF
453
454 DO I = 1,NEL
455 SOUNDSP(I) = SQRT(KMAX*TWO/(RHO0(I)))
456 ENDDO
457
458 RETURN
subroutine area(d1, x, x2, y, y2, eint, stif0)