OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
i21for3.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!|| i21for3 ../engine/source/interfaces/int21/i21for3.F
25!||--- called by ------------------------------------------------------
26!|| i21mainf ../engine/source/interfaces/int21/i21mainf.F
27!||--- calls -----------------------------------------------------
28!|| finter ../engine/source/tools/curve/finter.F
29!|| ibcoff ../engine/source/interfaces/interf/ibcoff.F
30!||====================================================================
31 SUBROUTINE i21for3(JLT ,NIN ,NOINT ,IBCC ,ICODT ,
32 2 FSAV ,GAP ,STIGLO ,VISC ,INACTI ,
33 3 MFROT ,IFQ ,IBAG ,IADM ,ICURV ,
34 4 STIF ,GAPV ,ITAB ,PENI ,ALPHA0 ,
35 5 IFPEN ,ICONTACT,RCONTACT,ACONTACT,PCONTACT,
36 6 NSVG ,X1 ,Y1 ,Z1 ,X2 ,
37 7 Y2 ,Z2 ,X3 ,Y3 ,Z3 ,
38 8 X4 ,Y4 ,Z4 ,XI ,YI ,
39 9 ZI ,VXI ,VYI ,VZI ,MSI ,
40 A VXM ,VYM ,VZM ,NX ,NY ,
41 B NZ ,PENE ,FXT ,FYT ,FZT ,
42 C FXN ,FYN ,FZN ,RCURVI ,ANGLMI ,
43 D PADM ,CAND_N_N, WEIGHT,IGAP ,GAP0 ,
44 E AREA0 ,PMAX ,IROT ,XG ,MXI ,
45 F MYI ,MZI ,STRI ,WXM ,WYM ,
46 G WZM ,XP ,YP ,ZP ,KT ,
47 H C ,ILEV ,FNI ,INTTH ,FHEAT ,
48 I EFRICT ,QFRIC ,IFRIC ,XFRIC ,TEMPI ,
49 J TEMPM ,NPC ,TF ,IX1 ,IX2 ,
50 K IX3 ,IX4 ,DT2T ,NELTST ,ITYPTST,
51 L KINET ,NISUB ,ISENSINT,FSAVPARIT,NFT ,
52 M IFTLIM ,PSKIDFLAG,PRATIO,PMAXSKID,INTEREFRIC,
53 N EFRIC_L ,FRICC ,FRIC_COEFS)
54C-----------------------------------------------
55C M o d u l e s
56C-----------------------------------------------
57
58C-----------------------------------------------
59C I m p l i c i t T y p e s
60C-----------------------------------------------
61#include "implicit_f.inc"
62#include "comlock.inc"
63C-----------------------------------------------
64C G l o b a l P a r a m e t e r s
65C-----------------------------------------------
66#include "mvsiz_p.inc"
67C-----------------------------------------------
68C C o m m o n B l o c k s
69C-----------------------------------------------
70#include "com01_c.inc"
71#include "com06_c.inc"
72#include "com08_c.inc"
73#include "scr05_c.inc"
74#include "scr11_c.inc"
75#include "sms_c.inc"
76#include "scr18_c.inc"
77#include "units_c.inc"
78#include "kincod_c.inc"
79C-----------------------------------------------
80C D u m m y A r g u m e n t s
81C-----------------------------------------------
82 INTEGER JLT, IBCC, INACTI, IBAG, NIN, NOINT, IADM,
83 . MFROT, IFQ, ICURV(3), IGAP, IROT, ILEV,INTTH,IFRIC,
84 . NELTST,ITYPTST,IFTLIM,
85 . ICODT(*), ITAB(*),IFPEN(*) ,ICONTACT(*), NPC(*),KINET(*)
86 INTEGER NSVG(MVSIZ),CAND_N_N(MVSIZ), WEIGHT(*),
87 . IX1(MVSIZ), IX2(MVSIZ), IX3(MVSIZ), IX4(MVSIZ), NISUB,
88 . ISENSINT(*),NFT,PSKIDFLAG
89 INTEGER , INTENT(IN) :: INTEREFRIC
90 my_real
91 . STIGLO, PENI(*), FSAV(*), TF(*),
92 . ALPHA0, GAP, VISC,DT2T,PMAXSKID ,
93 . FNI(MVSIZ),
94 . FXN(MVSIZ), FYN(MVSIZ), FZN(MVSIZ),
95 . FXT(MVSIZ), FYT(MVSIZ), FZT(MVSIZ)
96 my_real
97 . STIF(MVSIZ), GAPV(MVSIZ),
98 . VXI(MVSIZ),VYI(MVSIZ),VZI(MVSIZ),MSI(MVSIZ),
99 . X1(MVSIZ),Y1(MVSIZ),Z1(MVSIZ),
100 . X2(MVSIZ),Y2(MVSIZ),Z2(MVSIZ),
101 . X3(MVSIZ),Y3(MVSIZ),Z3(MVSIZ),
102 . X4(MVSIZ),Y4(MVSIZ),Z4(MVSIZ),
103 . XI(MVSIZ),YI(MVSIZ),ZI(MVSIZ),
104 . nx(mvsiz),ny(mvsiz),nz(mvsiz),pene(mvsiz),
105 . gap0(mvsiz), area0(mvsiz), pmax, ftmax,
106 . vxm, vym, vzm, wxm, wym, wzm,
107 . xp(mvsiz), yp(mvsiz), zp(mvsiz)
108 my_real
109 . rcurvi(*), rcontact(*), acontact(*),
110 . pcontact(*), padm, anglmi(*),
111 . xg(3), mxi(mvsiz), myi(mvsiz), mzi(mvsiz), stri(mvsiz),
112 . kt(mvsiz), c(mvsiz), fheat,efrict(mvsiz),qfric,
113 . tempi(mvsiz),tempm(mvsiz),xfric,
114 . fsavparit(nisub+1,11,*),pratio(mvsiz)
115 my_real , INTENT(INOUT) :: efric_l(mvsiz)
116 my_real , INTENT(IN) :: fric_coefs(mvsiz,10), fricc(mvsiz)
117C-----------------------------------------------
118C L o c a l V a r i a b l e s
119C-----------------------------------------------
120 INTEGER I, IG, JG , NI,
121 . IBCM,IBCS
122 my_real
123 . FXI(MVSIZ), FYI(MVSIZ), FZI(MVSIZ),
124 . XMU(MVSIZ),
125 . VX(MVSIZ), VY(MVSIZ), VZ(MVSIZ), VN(MVSIZ), DTMI(MVSIZ),
126 . AA,
127 . V2, ALPHA, BETA,
128 . FX, FY, FZ, FT, FN, FTN,
129 . ECONTT, ECONVT, FS2,ECONTDT,
130 . fsav1, fsav2, fsav3, fsav4, fsav5, fsav6, fsav8,
131 . fsav9, fsav10, fsav11, fsav12, fsav13, fsav14, fsav15,fsav25,
132 . vv,area,p,vv1,vv2,dmu,
133 . dt1inv, vis, pa, plin, fs, qfrict,dtmi0, mas2,beta1
134 my_real
135 . prec
136 my_real
137 .
138 . penx(mvsiz),stif0(mvsiz)
139 my_real
140 . impx,impy,impz,xx,yy,zz
141 my_real
142 . finter
143 EXTERNAL finter
144C
145C-----------------------------------------------
146 IF (iresp==1) THEN
147 prec = fiveem4
148 ELSE
149 prec = em10
150 ENDIF
151 IF(dt1>zero)THEN
152 dt1inv = one/dt1
153 ELSE
154 dt1inv =zero
155 ENDIF
156 IF(pskidflag >0) THEN
157 DO i=1,jlt
158 pratio(i) = zero
159 ENDDO
160 ENDIF
161 efric_l(1:jlt) = zero
162 efrict(1:jlt) = zero
163C---------------------
164C PENE INITIALE
165C---------------------
166 IF(ilev==0)THEN
167 IF(igap<=1)THEN
168 IF(inacti==5)THEN
169 DO i=1,jlt
170C REDUCTION OF THE INITIAL PENETRATION
171 peni(cand_n_n(i))=min(peni(cand_n_n(i)),
172 . ((one-fiveem2)*peni(cand_n_n(i))+fiveem2*pene(i)) )
173C SUBTRACTION OF THE INITIAL PENALTY FROM THE PENALTY AND THE GAP
174 pene(i)=max(zero,pene(i)-peni(cand_n_n(i)))
175 IF( pene(i)==zero ) stif(i) = zero
176 gapv(i)=gapv(i)-peni(cand_n_n(i))
177 ENDDO
178 ELSEIF(inacti==6)THEN
179 DO i=1,jlt
180C REDUCTION OF THE INITIAL PENALTY
181 peni(cand_n_n(i))=min(peni(cand_n_n(i)),
182 . ( (one-fiveem2)*peni(cand_n_n(i))
183 . +fiveem2*(pene(i)+fiveem2*(gapv(i)-pene(i)))) )
184C SUBTRACTION OF THE INITIAL PENALTY FROM THE PENALTY AND THE GAP
185 pene(i)=max(zero,pene(i)-peni(cand_n_n(i)))
186 IF( pene(i)==zero ) stif(i) = zero
187 gapv(i)=gapv(i)-peni(cand_n_n(i))
188 ENDDO
189 ELSE
190 DO i=1,jlt
191 IF( pene(i)==zero ) stif(i) = zero
192 ENDDO
193 ENDIF
194 ELSE
195 IF(inacti==5)THEN
196 DO i=1,jlt
197C REDUCTION OF THE INITIAL PENALTY
198 penx(i)=pene(i)
199 IF(penx(i) > max(gapv(i)-gap0(i),zero))
200 . penx(i) = penx(i)-max(gapv(i)-gap0(i),zero)
201 peni(cand_n_n(i))=min(peni(cand_n_n(i)),
202 . ((one-fiveem2)*peni(cand_n_n(i))+fiveem2*penx(i)) )
203C SUBTRACTION OF THE INITIAL PENALTY FROM THE PENALTY AND THE GAP
204 pene(i)=max(zero,pene(i)-peni(cand_n_n(i)))
205 IF( pene(i)==zero ) stif(i) = zero
206 gapv(i)=gapv(i)-peni(cand_n_n(i))
207 gap0(i)=gap0(i)-peni(cand_n_n(i))
208 ENDDO
209 ELSEIF(inacti==6)THEN
210 DO i=1,jlt
211C REDUCTION OF THE INITIAL PENALTY
212 penx(i)=pene(i)
213 IF(penx(i) > max(gapv(i)-gap0(i),zero))
214 . penx(i) = penx(i)-max(gapv(i)-gap0(i),zero)
215 peni(cand_n_n(i))=min(peni(cand_n_n(i)),
216 . ( (one-fiveem2)*peni(cand_n_n(i))
217 . +fiveem2*(penx(i)+fiveem2*(gap0(i)-penx(i)))) )
218C SUBTRACTION OF THE INITIAL PENALTY FROM THE PENALTY AND THE GAP
219 pene(i)=max(zero,pene(i)-peni(cand_n_n(i)))
220 IF( pene(i)==zero ) stif(i) = zero
221 gapv(i)=gapv(i)-peni(cand_n_n(i))
222 gap0(i)=gap0(i)-peni(cand_n_n(i))
223 ENDDO
224 ELSE
225 DO i=1,jlt
226 IF( pene(i)==zero ) stif(i) = zero
227 ENDDO
228 ENDIF
229 END IF
230 ELSE
231C
232C---- ILEV=1
233 IF(inacti==6)THEN
234 DO i=1,jlt
235C REDUCTION OF THE INITIAL PENALTY
236 peni(cand_n_n(i))=min(peni(cand_n_n(i)),
237 . ( (one-fiveem2)*peni(cand_n_n(i))
238 . +fiveem2*(pene(i)+fiveem2*abs(gapv(i)-pene(i)))) )
239C SUBTRACTION OF THE INITIAL PENALTY FROM THE PENALTY AND THE GAP
240 pene(i)=max(zero,pene(i)-peni(cand_n_n(i)))
241 IF( pene(i)==zero .AND.
242 . ( ifpen(cand_n_n(i))/=1 .OR.tt==zero ) ) stif(i) = zero
243 ENDDO
244 ELSE
245 DO i=1,jlt
246C REDUCTION OF THE INITIAL PENALTY
247 peni(cand_n_n(i))=min(peni(cand_n_n(i)),
248 . ((one-fiveem2)*peni(cand_n_n(i))+fiveem2*pene(i)) )
249C SUBTRACTION OF THE INITIAL PENALTY FROM THE PENALTY AND THE GAP
250 pene(i)=max(zero,pene(i)-peni(cand_n_n(i)))
251 IF( pene(i)==zero .AND.
252 . (ifpen(cand_n_n(i))/=1 .OR.
253 . (inacti==5.AND.tt==zero) ) ) stif(i) = zero
254 ENDDO
255
256 END IF
257 END IF
258C
259 DO i=1,jlt
260 stif0(i) = stif(i)
261 ENDDO
262C
263C-------------------------------------------
264C FNI + STIF
265C---------------------------------
266 econtt = zero
267 econvt = zero
268 qfrict = zero
269 econtdt = zero
270 IF(igap<=1)THEN
271 DO i=1,jlt
272 ig=nsvg(i)
273 IF(stiglo<=zero)THEN
274 econtt = econtt - half*stiglo*stif(i)*pene(i)**2
275 . * weight(ig)
276 stif(i) = -stiglo*stif(i)
277
278 ELSEIF(stif(i)/=zero)THEN
279 econtt = econtt + stiglo**pene(i)**2 * weight(ig)
280 stif(i) = stiglo
281 ENDIF
282 fni(i)= - stif(i) * pene(i)
283 END DO
284 ELSE
285 DO i=1,jlt
286 ig=nsvg(i)
287 IF(stiglo<=zero)THEN
288 stif(i) = -stiglo*stif(i)
289 ELSEIF(stif(i)/=zero)THEN
290 stif(i) = stiglo
291 ENDIF
292 pa = area0(i)*pmax
293 IF(stif(i)*max(gapv(i)-gap0(i),zero) <= pa)THEN
294 fni(i)= - stif(i)*pene(i)
295 econtt = econtt - half * pene(i) * fni(i)* weight(ig)
296 ELSE
297 fni(i)= - stif(i)*max(pene(i)-max(gapv(i)-gap0(i),zero),zero)
298 . - min(pa,stif(i)*pene(i))
299 plin = -fni(i)/max(em20,stif(i))
300 econtt = econtt + weight(ig)*
301 . ( max(pene(i)-plin,zero)*area0(i)*pmax - half *plin *fni(i) )
302 END IF
303 END DO
304 END IF
305C---------------------------------
306C DAMPING
307C---------------------------------
308 DO i=1,jlt
309 xx=xp(i)-xg(1)
310 yy=yp(i)-xg(2)
311 zz=zp(i)-xg(3)
312 vx(i)=vxi(i)-(vxm + wym*zz-wzm*yy)
313 vy(i)=vyi(i)-(vym + wzm*xx-wxm*zz)
314 vz(i)=vzi(i)-(vzm + wxm*yy-wym*xx)
315 vn(i) = nx(i)*vx(i) + ny(i)*vy(i) + nz(i)*vz(i)
316 ENDDO
317C-------------------------------------------------
318C
319C IF(KDTINT==0.AND.IDTMINS/=2)THEN
320 IF(idtmins/=2.AND.idtmins_int==0)THEN
321 DO i=1,jlt
322 vis = visc * sqrt(two * stif(i) * msi(i))
323 stif(i) = stif(i) + vis *dt1inv
324 fni(i) = fni(i) + vis * vn(i)
325 ig=nsvg(i)
326 econtdt = econtdt + vis * vn(i) * vn(i) * dt1 * weight(ig)
327 ENDDO
328 ELSE
329 DO i=1,jlt
330 c(i) = visc * sqrt(two * stif(i) * msi(i))
331 kt(i)= stif(i)
332 stif(i) = stif(i) + c(i) *dt1inv
333 fni(i) = fni(i) + c(i) * vn(i)
334 ig=nsvg(i)
335 econtdt = econtdt + c(i) * vn(i) * vn(i) * dt1 * weight(ig)
336 ENDDO
337 END IF
338C---------------------------------
339C CALCULATION OF THE NORMAL FORCE
340C---------------------------------
341 DO i=1,jlt
342 fxn(i)=fni(i)*nx(i)
343 fyn(i)=fni(i)*ny(i)
344 fzn(i)=fni(i)*nz(i)
345 END DO
346C---------------------------------
347C SAUVEGARDE DE L'IMPULSION NORMALE
348C---------------------------------
349 fsav1 = zero
350 fsav2 = zero
351 fsav3 = zero
352 fsav8 = zero
353 fsav9 = zero
354 fsav10= zero
355 fsav11= zero
356 DO i=1,jlt
357 ig=nsvg(i)
358 impx=fxn(i)*dt12*weight(ig)
359 impy=fyn(i)*dt12*weight(ig)
360 impz=fzn(i)*dt12*weight(ig)
361 fsav1 =fsav1 +impx
362 fsav2 =fsav2 +impy
363 fsav3 =fsav3 +impz
364 fsav8 =fsav8 +abs(impx)
365 fsav9 =fsav9 +abs(impy)
366 fsav10=fsav10+abs(impz)
367 fsav11=fsav11+fni(i)*dt12
368 ENDDO
369c
370#include "lockon.inc"
371 fsav(1)=fsav(1)+fsav1
372 fsav(2)=fsav(2)+fsav2
373 fsav(3)=fsav(3)+fsav3
374 fsav(8)=fsav(8)+fsav8
375 fsav(9)=fsav(9)+fsav9
376 fsav(10)=fsav(10)+fsav10
377 fsav(11)=fsav(11)+fsav11
378#include "lockoff.inc"
379c
380 IF(isensint(1)/=0) THEN
381 DO i=1,jlt
382 ig=nsvg(i)
383 fsavparit(1,1,i+nft) = fxn(i)*weight(ig)
384 fsavparit(1,2,i+nft) = fyn(i)*weight(ig)
385 fsavparit(1,3,i+nft) = fzn(i)*weight(ig)
386 ENDDO
387 ENDIF
388C---------------------------------
389C NEW FRICTION MODELS
390C---------------------------------
391
392 IF (mfrot==0) THEN
393C--- Coulomb friction
394 DO i=1,jlt
395 xmu(i) = fricc(i)
396 ENDDO
397 ELSEIF (mfrot==1) THEN
398C--- Viscous friction
399 DO i=1,jlt
400C attention : normal <> normal to the elt
401 aa = nx(i)*vx(i) + ny(i)*vy(i) + nz(i)*vz(i)
402 v2 = (vx(i) - nx(i)*aa)**2
403 . + (vy(i) - ny(i)*aa)**2
404 . + (vz(i) - nz(i)*aa)**2
405 vv = sqrt(max(em30,v2))
406C
407C should be the current nodal surface instead of the initial one.
408 area = area0(i)
409 p = -fni(i)/area
410 xmu(i) = fricc(i)+ (fric_coefs(i,1) + fric_coefs(i,4)*p ) * p
411 . +(fric_coefs(i,2) + fric_coefs(i,3)*p) * vv + fric_coefs(i,5)*v2
412 xmu(i) = max(xmu(i),em30)
413 ENDDO
414 ELSEIF(mfrot==2)THEN
415C--- Darmstad LAW
416 DO i=1,jlt
417 aa = nx(i)*vx(i) + ny(i)*vy(i) + nz(i)*vz(i)
418 v2 = (vx(i) - nx(i)*aa)**2
419 . + (vy(i) - ny(i)*aa)**2
420 . + (vz(i) - nz(i)*aa)**2
421 vv = sqrt(max(em30,v2))
422C
423C should be the current nodal surface instead of the initial one.
424 area = area0(i)
425 p = -fni(i)/area
426 xmu(i) = fricc(i)
427 . + fric_coefs(i,1)*exp(fric_coefs(i,2)*vv)*p*p
428 . + fric_coefs(i,3)*exp(fric_coefs(i,4)*vv)*p
429 . + fric_coefs(i,5)*exp(fric_coefs(i,6)*vv)
430 xmu(i) = max(xmu(i),em30)
431 ENDDO
432 ELSEIF (mfrot==3) THEN
433C--- Renard LAW
434 DO i=1,jlt
435 aa = nx(i)*vx(i) + ny(i)*vy(i) + nz(i)*vz(i)
436 v2 = (vx(i) - nx(i)*aa)**2
437 . + (vy(i) - ny(i)*aa)**2
438 . + (vz(i) - nz(i)*aa)**2
439 vv = sqrt(max(em30,v2))
440 IF(vv>=0.AND.vv<=fric_coefs(i,5)) THEN
441 dmu = fric_coefs(i,3)-fric_coefs(i,1)
442 vv1 = vv / fric_coefs(i,5)
443 xmu(i) = fric_coefs(i,1)+ dmu*vv1*(two-vv1)
444 ELSEIF(vv>fric_coefs(i,5).AND.vv<fric_coefs(i,6)) THEN
445 dmu = fric_coefs(i,4)-fric_coefs(i,3)
446 vv1 = (vv - fric_coefs(i,5))/(fric_coefs(i,6)-fric_coefs(i,5))
447 xmu(i) = fric_coefs(i,3)+ dmu * (three-two*vv1)*vv1**2
448 ELSE
449 dmu = fric_coefs(i,2)-fric_coefs(i,4)
450 vv2 = (vv - fric_coefs(i,6))**2
451 xmu(i) = fric_coefs(i,2) - dmu / (one + dmu*vv2)
452 ENDIF
453 xmu(i) = max(xmu(i),em30)
454 ENDDO
455 ELSEIF(mfrot==4)THEN
456C--- Exponential decay model
457 DO i=1,jlt
458 aa = nx(i)*vx(i) + ny(i)*vy(i) + nz(i)*vz(i)
459 v2 = (vx(i) - nx(i)*aa)**2
460 . + (vy(i) - ny(i)*aa)**2
461 . + (vz(i) - nz(i)*aa)**2
462 vv = sqrt(max(em30,v2))
463 xmu(i) = fric_coefs(i,1)
464 . + (fric_coefs(i,1)-fricc(i))*exp(-fric_coefs(i,2)*vv)
465 xmu(i) = max(xmu(i),em30)
466 ENDDO
467 ENDIF
468C------------------
469C TANGENT FORCE CALCULATION
470C------------------
471 fsav4 = zero
472 fsav5 = zero
473 fsav6 = zero
474 fsav12= zero
475 fsav13= zero
476 fsav14= zero
477 fsav15= zero
478 fsav25= zero
479C---------------------------------
480C INCREMENTAL (STIFFNESS) FORMULATION
481C---------------------------------
482 IF (ifq==13) THEN
483 alpha = max(one,alpha0*dt12)
484 ELSE
485 alpha = alpha0
486 ENDIF
487 IF(iftlim == 0 ) THEN ! if flag to limit tangential force is ON : FT<= YIELD/(S*sqrt(3))
488 ftmax = pmax
489 ELSE
490 ftmax = ep30
491 ENDIF
492 DO i=1,jlt
493 fx = stif0(i)*vx(i)*dt12
494 fy = stif0(i)*vy(i)*dt12
495 fz = stif0(i)*vz(i)*dt12
496
497 fx = fxt(i) + alpha*fx
498 fy = fyt(i) + alpha*fy
499 fz = fzt(i) + alpha*fz
500
501 ftn = fx*nx(i) + fy*ny(i) + fz*nz(i)
502
503 fx = fx - ftn*nx(i)
504 fy = fy - ftn*ny(i)
505 fz = fz - ftn*nz(i)
506
507 ft = fx*fx + fy*fy + fz*fz
508 ft = max(ft,em30)
509
510 fn = fni(i)*fni(i)
511
512 fn = sqrt(fn)
513 ft = sqrt(ft)
514
515 beta1 = xmu(i)*fn/ft
516 beta = min(one,beta1)
517
518C should be the current nodal surface instead of the initial one.
519 area = area0(i)
520 fs = ftmax/sqrt(three)*area
521 beta = min(beta,fs/ft)
522 IF(pskidflag >0) THEN
523 beta1 = min(beta1,fs/ft)
524 IF(beta1 <= one) THEN
525 fs2 = pmaxskid/sqrt(three)*area
526 pratio(i) = min(one,beta*ft/fs2)
527 ENDIF
528 ENDIF
529C
530 fxt(i) = fx * beta
531 fyt(i) = fy * beta
532 fzt(i) = fz * beta
533C------- total force
534 fxi(i)=fxn(i) + fxt(i)
535 fyi(i)=fyn(i) + fyt(i)
536 fzi(i)=fzn(i) + fzt(i)
537 ig=nsvg(i)
538C---------------------------------
539C CONTACT ENERGY CALCULATION
540C---------------------------------
541 efric_l(i) = dt1*(vx(i)*fxt(i)+vy(i)*fyt(i)+vz(i)*fzt(i))
542 . * weight(ig)
543 econvt = econvt + efric_l(i)
544
545 IF( intth > 0 .AND.beta/=zero) THEN
546 efrict(i) = (fx-fxt(i))*fxt(i) + (fy-fyt(i))*fyt(i) +
547 . (fz-fzt(i))*fzt(i) ! FRICTIONAL ENERGY
548 efrict(i) = efrict(i)/stif0(i)
549 qfrict = qfrict + efrict(i)
550 ENDIF
551 ENDDO
552C---------------------------------
553 DO i=1,jlt
554 ig=nsvg(i)
555 impx=fxt(i)*dt12*weight(ig)
556 impy=fyt(i)*dt12*weight(ig)
557 impz=fzt(i)*dt12*weight(ig)
558 fsav4 =fsav4 +impx
559 fsav5 =fsav5 +impy
560 fsav6 =fsav6 +impz
561 impx=fxi(i)*dt12
562 impy=fyi(i)*dt12
563 impz=fzi(i)*dt12
564 fsav12=fsav12+abs(impx)
565 fsav13=fsav13+abs(impy)
566 fsav14=fsav14+abs(impz)
567 fsav15=fsav15+sqrt(impx*impx+impy*impy+impz*impz)
568 ENDDO
569c
570#include "lockon.inc"
571 fsav(4) = fsav(4) + fsav4
572 fsav(5) = fsav(5) + fsav5
573 fsav(6) = fsav(6) + fsav6
574 fsav(12) = fsav(12) + fsav12
575 fsav(13) = fsav(13) + fsav13
576 fsav(14) = fsav(14) + fsav14
577 fsav(15) = fsav(15) + fsav15
578 fsav(25) = fsav(25) + fheat*qfrict
579 fsav(26) = fsav(26) + econtt
580 fsav(27) = fsav(27) + econvt- fheat*qfrict
581 fsav(28) = fsav(28) + econtdt
582#include "lockoff.inc"
583c
584 IF(isensint(1)/=0) THEN
585 DO i=1,jlt
586 ig=nsvg(i)
587 fsavparit(1,4,i+nft) = fxt(i)*weight(ig)
588 fsavparit(1,5,i+nft) = fyt(i)*weight(ig)
589 fsavparit(1,6,i+nft) = fzt(i)*weight(ig)
590 ENDDO
591 ENDIF
592C---------------------------------
593#include "lockon.inc"
594 econtv = econtv + econvt ! Frictional Energy
595 econt = econt + econtt ! Elastic Energy
596 econtd = econtd + econtdt ! Damping Energy
597 IF (intth/=0) THEN
598 qfric = qfric + fheat*qfrict ! FRICTIONAL HEAT ADDED TO INTERNAL ENERGY
599 econtv = econtv - fheat*qfrict ! FRICTIONAL HEAT REMOVED FROM CONTACT ENERGY
600 ENDIF
601#include "lockoff.inc"
602C---------------------------------
603 IF(irot/=0)THEN
604 DO i=1,jlt
605 xx=xp(i)-xg(1)
606 yy=yp(i)-xg(2)
607 zz=zp(i)-xg(3)
608 mxi(i) =yy*fzi(i)-zz*fyi(i)
609 myi(i) =zz*fxi(i)-xx*fzi(i)
610 mzi(i) =xx*fyi(i)-yy*fxi(i)
611 stri(i)= stif(i)*(xx*xx+yy*yy+zz*zz)
612 END DO
613 END IF
614C---------------------------------
615C
616c IF(IDTMIN(10)==1.OR.IDTMIN(10)==2.OR.
617c . IDTMIN(10)==5.OR.IDTMIN(10)==6)THEN
618 IF(idtmin(10)==1)THEN
619 dtmi0 = ep20
620 DO i=1,jlt
621 dtmi(i) = ep20
622 mas2 = two * msi(i)
623 jg = nsvg(i)
624 IF(mas2>zero.AND.stif(i)>zero.AND.
625 . irb(kinet(jg))==0.AND.irb2(kinet(jg))==0)THEN
626 dtmi(i) = min(dtmi(i),dtfac1(10)*sqrt(mas2/stif(i)))
627 ENDIF
628 dtmi0 = min(dtmi0,dtmi(i))
629 ENDDO
630
631 IF(dtmi0<=dtmin1(10))THEN
632 DO i=1,jlt
633 IF(dtmi(i)<=dtmin1(10))THEN
634 jg = nsvg(i)
635 ni = itab(jg)
636 IF(idtmin(10)==1)THEN
637#include "lockon.inc"
638 WRITE(iout,'(A,E12.4,A,I10,A,E12.4,A)')
639 . ' **WARNING MINIMUM TIME STEP ',dtmi(i),
640 . ' IN INTERFACE ',noint,'(DTMIN =',dtmin1(10),')'
641 WRITE(iout,'(A,I10)') ' SECONDARY NODE : ',ni
642c WRITE(IOUT,'(A,4I10)')' MAIN NODES : ',
643c . ITAB(IX1(I)),ITAB(IX2(I)),ITAB(IX3(I)),ITAB(IX4(I))
644#include "lockoff.inc"
645 tstop = tt
646 IF ( istamping == 1) THEN
647 WRITE(istdo,'(A)')'The run encountered a problem in an in
648 .terface Type 21.'
649 WRITE(istdo,'(A)')'Maximum penetration may be reached'
650 WRITE(istdo,'(A)')'You may need to check if contact normals
651 .of tools are oriented toward the blank,'
652 WRITE(iout, '(A)')'The run encountered a problem in an in
653 .terface Type 21.'
654 WRITE(iout, '(A)')'Maximum penetration may be reached'
655 WRITE(iout, '(A)')'You may need to check if contact normals
656 .of tools are oriented toward the blank,'
657 ENDIF
658 ENDIF
659 ENDIF
660 ENDDO
661 ENDIF
662 ENDIF
663C-----------------------------------------------------
664 IF(ibag/=0.OR.iadm/=0)THEN
665 DO i=1,jlt
666 IF(pene(i)/=zero)THEN
667 jg = nsvg(i)
668 icontact(jg)=1
669 ENDIF
670 ENDDO
671 ENDIF
672 IF(iadm/=0)THEN
673 DO i=1,jlt
674 jg = nsvg(i)
675#include "lockon.inc"
676 rcontact(jg)=min(rcontact(jg),rcurvi(i))
677#include "lockoff.inc"
678 END DO
679 END IF
680 IF(iadm>=2)THEN
681 DO i=1,jlt
682 jg = nsvg(i)
683#include "lockon.inc"
684 pcontact(jg)=max(pcontact(jg),pene(i)/(padm*gapv(i)))
685 acontact(jg)=min(acontact(jg),anglmi(i))
686#include "lockoff.inc"
687 END DO
688 END IF
689C
690 IF(ibcc==0) RETURN
691C
692 DO 400 i=1,jlt
693 IF(pene(i)==zero)GOTO 400
694 ibcm = ibcc / 8
695 ibcs = ibcc - 8 * ibcm
696 IF(ibcs>0) THEN
697 ig=nsvg(i)
698 CALL ibcoff(ibcs,icodt(ig))
699 ENDIF
700 400 CONTINUE
701C
702 RETURN
703 END
area
subroutine area(d1, x, x2, y, y2, eint, stif0)
Definition hm_read_prop32.F:567
i21for3
subroutine i21for3(jlt, nin, noint, ibcc, icodt, fsav, gap, stiglo, visc, inacti, mfrot, ifq, ibag, iadm, icurv, stif, gapv, itab, peni, alpha0, ifpen, icontact, rcontact, acontact, pcontact, nsvg, x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4, xi, yi, zi, vxi, vyi, vzi, msi, vxm, vym, vzm, nx, ny, nz, pene, fxt, fyt, fzt, fxn, fyn, fzn, rcurvi, anglmi, padm, cand_n_n, weight, igap, gap0, area0, pmax, irot, xg, mxi, myi, mzi, stri, wxm, wym, wzm, xp, yp, zp, kt, c, ilev, fni, intth, fheat, efrict, qfric, ifric, xfric, tempi, tempm, npc, tf, ix1, ix2, ix3, ix4, dt2t, neltst, ityptst, kinet, nisub, isensint, fsavparit, nft, iftlim, pskidflag, pratio, pmaxskid, interefric, efric_l, fricc, fric_coefs)
Definition i21for3.F:54
ibcoff
subroutine ibcoff(ibc, icodt)
Definition ibcoff.F:44
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
stri
character *8 function stri(n)
Definition stri.F:24