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,dydx,
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, J1, IG, J, JG , K0,NBINTER,K1S,K,IL,IE, NN, NI,
121 . NA1,NA2,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 . VNX, VNY, VNZ, AA, S2,
127 . V2, FF, ALPHA, BETA,
128 . FX, FY, FZ, FT, FN, FMAX, FTN,
129 . ECONTT, ECONVT, FS2,ECONTDT,
130 . fsav1, fsav2, fsav3, fsav4, fsav5, fsav6, fsav7, fsav8,
131 . fsav9, fsav10, fsav11, fsav12, fsav13, fsav14, fsav15,fsav25,
132 . vv,ax1,ax2,ay1,ay2,az1,az2,ax,ay,az,area,p,vv1,vv2,dmu,
133 . dt1inv, vis, pa, plin, fs, qfrict, dti,dti2,dtmi0, mas2,beta1
134 my_real
135 . prec
136 my_real
137 . st1(mvsiz),st2(mvsiz),st3(mvsiz),st4(mvsiz),stv(mvsiz),
138 . penx(mvsiz),stif0(mvsiz),frict(mvsiz)
139 my_real
140 . impx,impy,impz,xx,yy,zz,th
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 DE LA PENE INITIALE
171 peni(cand_n_n(i))=min(peni(cand_n_n(i)),
172 . ((one-fiveem2)*peni(cand_n_n(i))+fiveem2*pene(i)) )
173C SOUSTRACTION DE LA PENE INITIALE A LA PENE ET AU 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 DE LA PENE INITIALE
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 SOUSTRACTION DE LA PENE INITIALE A LA PENE ET AU 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 DE LA PENE INITIALE
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 SOUSTRACTION DE LA PENE INITIALE A LA PENE ET AU 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 DE LA PENE INITIALE
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 SOUSTRACTION DE LA PENE INITIALE A LA PENE ET AU 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 DE LA PENE INITIALE
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 SOUSTRACTION DE LA PENE INITIALE A LA PENE ET AU 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 DE LA PENE INITIALE
247 peni(cand_n_n(i))=min(peni(cand_n_n(i)),
248 . ((one-fiveem2)*peni(cand_n_n(i))+fiveem2*pene(i)) )
249C SOUSTRACTION DE LA PENE INITIALE A LA PENE ET AU 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 CALCUL DE LA FORCE NORMALE
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 : normale <> normale a l'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
417C attention : normale <> normale a l'elt
418 aa = nx(i)*vx(i) + ny(i)*vy(i) + nz(i)*vz(i)
419 v2 = (vx(i) - nx(i)*aa)**2
420 . + (vy(i) - ny(i)*aa)**2
421 . + (vz(i) - nz(i)*aa)**2
422 vv = sqrt(max(em30,v2))
423C
424C should be the current nodal surface instead of the initial one.
425 area = area0(i)
426 p = -fni(i)/area
427 xmu(i) = fricc(i)
428 . + fric_coefs(i,1)*exp(fric_coefs(i,2)*vv)*p*p
429 . + fric_coefs(i,3)*exp(fric_coefs(i,4)*vv)*p
430 . + fric_coefs(i,5)*exp(fric_coefs(i,6)*vv)
431 xmu(i) = max(xmu(i),em30)
432 ENDDO
433 ELSEIF (mfrot==3) THEN
434C--- Renard LAW
435 DO i=1,jlt
436C attention : normale <> normale a l'elt
437 aa = nx(i)*vx(i) + ny(i)*vy(i) + nz(i)*vz(i)
438 v2 = (vx(i) - nx(i)*aa)**2
439 . + (vy(i) - ny(i)*aa)**2
440 . + (vz(i) - nz(i)*aa)**2
441 vv = sqrt(max(em30,v2))
442 IF(vv>=0.AND.vv<=fric_coefs(i,5)) THEN
443 dmu = fric_coefs(i,3)-fric_coefs(i,1)
444 vv1 = vv / fric_coefs(i,5)
445 xmu(i) = fric_coefs(i,1)+ dmu*vv1*(two-vv1)
446 ELSEIF(vv>fric_coefs(i,5).AND.vv<fric_coefs(i,6)) THEN
447 dmu = fric_coefs(i,4)-fric_coefs(i,3)
448 vv1 = (vv - fric_coefs(i,5))/(fric_coefs(i,6)-fric_coefs(i,5))
449 xmu(i) = fric_coefs(i,3)+ dmu * (three-two*vv1)*vv1**2
450 ELSE
451 dmu = fric_coefs(i,2)-fric_coefs(i,4)
452 vv2 = (vv - fric_coefs(i,6))**2
453 xmu(i) = fric_coefs(i,2) - dmu / (one + dmu*vv2)
454 ENDIF
455 xmu(i) = max(xmu(i),em30)
456 ENDDO
457 ELSEIF(mfrot==4)THEN
458C--- Exponential decay model
459 DO i=1,jlt
460 aa = nx(i)*vx(i) + ny(i)*vy(i) + nz(i)*vz(i)
461 v2 = (vx(i) - nx(i)*aa)**2
462 . + (vy(i) - ny(i)*aa)**2
463 . + (vz(i) - nz(i)*aa)**2
464 vv = sqrt(max(em30,v2))
465 xmu(i) = fric_coefs(i,1)
466 . + (fric_coefs(i,1)-fricc(i))*exp(-fric_coefs(i,2)*vv)
467 xmu(i) = max(xmu(i),em30)
468 ENDDO
469 ENDIF
470C------------------
471C TANGENT FORCE CALCULATION
472C------------------
473 fsav4 = zero
474 fsav5 = zero
475 fsav6 = zero
476 fsav12= zero
477 fsav13= zero
478 fsav14= zero
479 fsav15= zero
480 fsav25= zero
481C---------------------------------
482C INCREMENTAL (STIFFNESS) FORMULATION
483C---------------------------------
484 IF (ifq==13) THEN
485 alpha = max(one,alpha0*dt12)
486 ELSE
487 alpha = alpha0
488 ENDIF
489 IF(iftlim == 0 ) THEN ! if flag to limit tangential force is ON : FT<= YIELD/(S*sqrt(3))
490 ftmax = pmax
491 ELSE
492 ftmax = ep30
493 ENDIF
494 DO i=1,jlt
495 fx = stif0(i)*vx(i)*dt12
496 fy = stif0(i)*vy(i)*dt12
497 fz = stif0(i)*vz(i)*dt12
498
499 fx = fxt(i) + alpha*fx
500 fy = fyt(i) + alpha*fy
501 fz = fzt(i) + alpha*fz
502
503 ftn = fx*nx(i) + fy*ny(i) + fz*nz(i)
504
505 fx = fx - ftn*nx(i)
506 fy = fy - ftn*ny(i)
507 fz = fz - ftn*nz(i)
508
509 ft = fx*fx + fy*fy + fz*fz
510 ft = max(ft,em30)
511
512 fn = fni(i)*fni(i)
513
514 fn = sqrt(fn)
515 ft = sqrt(ft)
516
517 beta1 = xmu(i)*fn/ft
518 beta = min(one,beta1)
519
520C should be the current nodal surface instead of the initial one.
521 area = area0(i)
522 fs = ftmax/sqrt(three)*area
523 beta = min(beta,fs/ft)
524 IF(pskidflag >0) THEN
525 beta1 = min(beta1,fs/ft)
526 IF(beta1 <= one) THEN
527 fs2 = pmaxskid/sqrt(three)*area
528 pratio(i) = min(one,beta*ft/fs2)
529 ENDIF
530 ENDIF
531C
532 fxt(i) = fx * beta
533 fyt(i) = fy * beta
534 fzt(i) = fz * beta
535C------- total force
536 fxi(i)=fxn(i) + fxt(i)
537 fyi(i)=fyn(i) + fyt(i)
538 fzi(i)=fzn(i) + fzt(i)
539 ig=nsvg(i)
540C---------------------------------
541C CONTACT ENERGY CALCULATION
542C---------------------------------
543 efric_l(i) = dt1*(vx(i)*fxt(i)+vy(i)*fyt(i)+vz(i)*fzt(i))
544 . * weight(ig)
545 econvt = econvt + efric_l(i)
546
547 IF( intth > 0 .AND.beta/=zero) THEN
548 efrict(i) = (fx-fxt(i))*fxt(i) + (fy-fyt(i))*fyt(i) +
549 . (fz-fzt(i))*fzt(i) ! FRICTIONAL ENERGY
550 efrict(i) = efrict(i)/stif0(i)
551 qfrict = qfrict + efrict(i)
552 ENDIF
553 ENDDO
554C---------------------------------
555 DO i=1,jlt
556 ig=nsvg(i)
557 impx=fxt(i)*dt12*weight(ig)
558 impy=fyt(i)*dt12*weight(ig)
559 impz=fzt(i)*dt12*weight(ig)
560 fsav4 =fsav4 +impx
561 fsav5 =fsav5 +impy
562 fsav6 =fsav6 +impz
563 impx=fxi(i)*dt12
564 impy=fyi(i)*dt12
565 impz=fzi(i)*dt12
566 fsav12=fsav12+abs(impx)
567 fsav13=fsav13+abs(impy)
568 fsav14=fsav14+abs(impz)
569 fsav15=fsav15+sqrt(impx*impx+impy*impy+impz*impz)
570 ENDDO
571c
572#include "lockon.inc"
573 fsav(4) = fsav(4) + fsav4
574 fsav(5) = fsav(5) + fsav5
575 fsav(6) = fsav(6) + fsav6
576 fsav(12) = fsav(12) + fsav12
577 fsav(13) = fsav(13) + fsav13
578 fsav(14) = fsav(14) + fsav14
579 fsav(15) = fsav(15) + fsav15
580 fsav(25) = fsav(25) + fheat*qfrict
581 fsav(26) = fsav(26) + econtt
582 fsav(27) = fsav(27) + econvt- fheat*qfrict
583 fsav(28) = fsav(28) + econtdt
584#include "lockoff.inc"
585c
586 IF(isensint(1)/=0) THEN
587 DO i=1,jlt
588 ig=nsvg(i)
589 fsavparit(1,4,i+nft) = fxt(i)*weight(ig)
590 fsavparit(1,5,i+nft) = fyt(i)*weight(ig)
591 fsavparit(1,6,i+nft) = fzt(i)*weight(ig)
592 ENDDO
593 ENDIF
594C---------------------------------
595#include "lockon.inc"
596 econtv = econtv + econvt ! Frictional Energy
597 econt = econt + econtt ! Elastic Energy
598 econtd = econtd + econtdt ! Damping Energy
599 IF (intth/=0) THEN
600 qfric = qfric + fheat*qfrict ! FRICTIONAL HEAT ADDED TO INTERNAL ENERGY
601 econtv = econtv - fheat*qfrict ! FRICTIONAL HEAT REMOVED FROM CONTACT ENERGY
602 ENDIF
603#include "lockoff.inc"
604C---------------------------------
605 IF(irot/=0)THEN
606 DO i=1,jlt
607 xx=xp(i)-xg(1)
608 yy=yp(i)-xg(2)
609 zz=zp(i)-xg(3)
610 mxi(i) =yy*fzi(i)-zz*fyi(i)
611 myi(i) =zz*fxi(i)-xx*fzi(i)
612 mzi(i) =xx*fyi(i)-yy*fxi(i)
613 stri(i)= stif(i)*(xx*xx+yy*yy+zz*zz)
614 END DO
615 END IF
616C---------------------------------
617C
618c IF(IDTMIN(10)==1.OR.IDTMIN(10)==2.OR.
619c . IDTMIN(10)==5.OR.IDTMIN(10)==6)THEN
620 IF(idtmin(10)==1)THEN
621 dtmi0 = ep20
622 DO i=1,jlt
623 dtmi(i) = ep20
624 mas2 = two * msi(i)
625 jg = nsvg(i)
626 IF(mas2>zero.AND.stif(i)>zero.AND.
627 . irb(kinet(jg))==0.AND.irb2(kinet(jg))==0)THEN
628 dtmi(i) = min(dtmi(i),dtfac1(10)*sqrt(mas2/stif(i)))
629 ENDIF
630 dtmi0 = min(dtmi0,dtmi(i))
631 ENDDO
632
633 IF(dtmi0<=dtmin1(10))THEN
634 DO i=1,jlt
635 IF(dtmi(i)<=dtmin1(10))THEN
636 jg = nsvg(i)
637 ni = itab(jg)
638 IF(idtmin(10)==1)THEN
639#include "lockon.inc"
640 WRITE(iout,'(A,E12.4,A,I10,A,E12.4,A)')
641 . ' **WARNING MINIMUM TIME STEP ',dtmi(i),
642 . ' IN INTERFACE ',noint,'(DTMIN =',dtmin1(10),')'
643 WRITE(iout,'(A,I10)') ' SECONDARY NODE : ',ni
644c WRITE(IOUT,'(A,4I10)')' MAIN NODES : ',
645c . ITAB(IX1(I)),ITAB(IX2(I)),ITAB(IX3(I)),ITAB(IX4(I))
646#include "lockoff.inc"
647 tstop = tt
648 IF ( istamping == 1) THEN
649 WRITE(istdo,'(A)')'The run encountered a problem in an in
650 .terface Type 21.'
651 WRITE(istdo,'(A)')'Maximum penetration may be reached'
652 WRITE(istdo,'(A)')'You may need to check if contact normals
653 .of tools are oriented toward the blank,'
654 WRITE(iout, '(A)')'The run encountered a problem in an in
655 .terface Type 21.'
656 WRITE(iout, '(A)')'Maximum penetration may be reached'
657 WRITE(iout, '(A)')'You may need to check if contact normals
658 .of tools are oriented toward the blank,'
659 ENDIF
660 ENDIF
661 ENDIF
662 ENDDO
663 ENDIF
664 ENDIF
665C-----------------------------------------------------
666 IF(ibag/=0.OR.iadm/=0)THEN
667 DO i=1,jlt
668 IF(pene(i)/=zero)THEN
669 jg = nsvg(i)
670 icontact(jg)=1
671 ENDIF
672 ENDDO
673 ENDIF
674 IF(iadm/=0)THEN
675 DO i=1,jlt
676 jg = nsvg(i)
677#include "lockon.inc"
678 rcontact(jg)=min(rcontact(jg),rcurvi(i))
679#include "lockoff.inc"
680 END DO
681 END IF
682 IF(iadm>=2)THEN
683 DO i=1,jlt
684 jg = nsvg(i)
685#include "lockon.inc"
686 pcontact(jg)=max(pcontact(jg),pene(i)/(padm*gapv(i)))
687 acontact(jg)=min(acontact(jg),anglmi(i))
688#include "lockoff.inc"
689 END DO
690 END IF
691C
692 IF(ibcc==0) RETURN
693C
694 DO 400 i=1,jlt
695 IF(pene(i)==zero)GOTO 400
696 ibcm = ibcc / 8
697 ibcs = ibcc - 8 * ibcm
698 IF(ibcs>0) THEN
699 ig=nsvg(i)
700 CALL ibcoff(ibcs,icodt(ig))
701 ENDIF
702 400 CONTINUE
703C
704 RETURN
705 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