OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
rivet1.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!|| rivet1 ../engine/source/elements/rivet/rivet1.F
25!||--- called by ------------------------------------------------------
26!|| resol ../engine/source/engine/resol.F
27!||====================================================================
28 SUBROUTINE rivet1(MS ,IN ,A ,AR ,X ,
29 . IXRT ,RIVET ,GEO ,V ,VR ,
30 . ITASK )
31C-----------------------------------------------
32C I m p l i c i t T y p e s
33C-----------------------------------------------
34#include "implicit_f.inc"
35#include "comlock.inc"
36C-----------------------------------------------
37C G l o b a l P a r a m e t e r s
38C-----------------------------------------------
39#include "mvsiz_p.inc"
40C-----------------------------------------------
41C C o m m o n B l o c k s
42C-----------------------------------------------
43#include "com04_c.inc"
44#include "com08_c.inc"
45#include "units_c.inc"
46#include "param_c.inc"
47#include "scr11_c.inc"
48#include "task_c.inc"
49C-----------------------------------------------
50C D u m m y A r g u m e n t s
51C-----------------------------------------------
52 INTEGER IXRT(4,*), ITASK
53C REAL
54 my_real
55 . ms(*), in(*), a(3,*), ar(3,*), x(3,*), rivet(nrivf,*),
56 . geo(npropg,*), v(3,*), vr(3,*)
57C-----------------------------------------------
58C L o c a l V a r i a b l e s
59C-----------------------------------------------
60 INTEGER I,IG,N1,N2,IROT,IMOD,RUPT,LFT,LLT,PROC,
61 . JS, NN, IS, NIND1, NIND2, J, L,
62 . ind1(mvsiz), ind2(mvsiz)
63
64C REAL
65 my_real
66 . da1(3),da2(3),vt1(3),vt2(3),
67 . dmx2, dt12m1, fn2, ft2,
68 . xm, amx, amy, amz, an, an2, anx, any, anz,
69 . atx, aty, atz, at2, alp,mass,masm1,iner,inm1,i1,i2,
70 . xcdg,ycdg,zcdg,xx1,yy1,zz1,xx2,yy2,zz2, ww,wt,
71 . vx1,vy1,vz1,vx2,vy2,vz2,vmx1,vmy1,vmz1,vmx2,vmy2,vmz2,
72 . vx,vy,vz,vrx1,vrx2,vry1,vry2,vrz1,vrz2,vxx,vyy,vzz,
73 . ax,ay,az, off, s,
74 . xn(mvsiz), yn(mvsiz), zn(mvsiz), dx2(mvsiz), enrot_l
75C-----------------------------------------------
76 rupt = 0
77 enrot_l = zero
78 proc = ispmd+1
79 lft = 1 + itask*nrivet / nthread
80 llt = (itask+1)*nrivet / nthread
81 js = lft-1
82 DO l=lft,llt,nvsiz
83 nn = min(nvsiz,llt-js)
84 nind1 = 0
85 nind2 = 0
86 DO is = 1 , nn
87 i = js+is
88C test si rivet a traiter sur le proc
89 IF(ixrt(2,i)>0)THEN
90C----------------------------
91C RUPTURE
92C----------------------------
93 off = rivet(1,i)
94 IF (off/=zero) THEN
95 ig=ixrt(1,i)
96 n1=ixrt(2,i)
97 n2=ixrt(3,i)
98 dmx2= geo(3,ig)
99 imod=nint(geo(5,ig))
100 xn(is)=x(1,n2)-x(1,n1)
101 yn(is)=x(2,n2)-x(2,n1)
102 zn(is)=x(3,n2)-x(3,n1)
103 dx2(is)=xn(is)**2+yn(is)**2+zn(is)**2
104 IF(dx2(is)>dmx2)THEN
105 off=off-em01
106 IF(off<=zero)THEN
107 rivet(1,i)=-one
108 rupt = 1
109 ELSE
110 rivet(1,i) = off
111 ENDIF
112 ENDIF
113 IF (off>zero) THEN
114 IF(imod==1) THEN
115 nind1 = nind1 + 1
116 ind1(nind1) = is
117 ELSE
118 nind2 = nind2 + 1
119 ind2(nind2) = is
120 ENDIF
121 ENDIF
122 ENDIF
123 ENDIF
124 ENDDO
125C
126C-----------------------------------------------
127C RIGID BODY formulation
128C-----------------------------------------------
129#include "vectorize.inc"
130 DO j = 1, nind1
131 is = ind1(j)
132 i = js+is
133 ig=ixrt(1,i)
134 n1=ixrt(2,i)
135 n2=ixrt(3,i)
136 off = rivet(1,i)
137 fn2 = geo(1,ig)
138 ft2 = geo(2,ig)
139 dmx2= geo(3,ig)
140 irot=nint(geo(4,ig))
141 dt12m1=one/dt12
142C----------------------------
143C TRANSLATION DU CDG
144C----------------------------
145 mass = ms(n1)+ms(n2)
146 masm1= one / mass
147 xcdg=(x(1,n1)*ms(n1)+x(1,n2)*ms(n2))*masm1
148 ycdg=(x(2,n1)*ms(n1)+x(2,n2)*ms(n2))*masm1
149 zcdg=(x(3,n1)*ms(n1)+x(3,n2)*ms(n2))*masm1
150C
151 vx1= v(1,n1)+a(1,n1)*dt12
152 vy1= v(2,n1)+a(2,n1)*dt12
153 vz1= v(3,n1)+a(3,n1)*dt12
154 vx2= v(1,n2)+a(1,n2)*dt12
155 vy2= v(2,n2)+a(2,n2)*dt12
156 vz2= v(3,n2)+a(3,n2)*dt12
157C
158 vmx1= vx1*ms(n1)
159 vmy1= vy1*ms(n1)
160 vmz1= vz1*ms(n1)
161 vmx2= vx2*ms(n2)
162 vmy2= vy2*ms(n2)
163 vmz2= vz2*ms(n2)
164C
165 vx = (vmx1+vmx2)*masm1
166 vy = (vmy1+vmy2)*masm1
167 vz = (vmz1+vmz2)*masm1
168C
169 IF(irot==0) THEN
170C----------------------------
171C CALCUL FORCES
172C----------------------------
173 ax = (-a(1,n1)+(vx-v(1,n1))*dt12m1)*ms(n1)
174 ay = (-a(2,n1)+(vy-v(2,n1))*dt12m1)*ms(n1)
175 az = (-a(3,n1)+(vz-v(3,n1))*dt12m1)*ms(n1)
176C
177 IF(dx2(is)>em15)THEN
178C POINTS NON CONFONDUS CRITERE FN ET FT
179 s = one/sqrt(dx2(is))
180 xn(is) =xn(is)*s
181 yn(is) =yn(is)*s
182 zn(is) =zn(is)*s
183 an =ax*xn(is)+ay*yn(is)+az*zn(is)
184 an2=an**2
185 anx=an*xn(is)
186 any=an*yn(is)
187 anz=an*zn(is)
188 atx=ax-anx
189 aty=ay-any
190 atz=az-anz
191 at2=(atx**2+aty**2+atz**2)
192 ELSE
193C POINTS CONFONDUS CRITERE UNIQUEMENT SUR FN
194 an2=(ax**2+ay**2+az**2)
195 at2=0.
196 ENDIF
197C
198 alp=sqrt((an2/fn2)+(at2/ft2))
199 alp=off / max(alp,one)
200 ax=alp*ax
201 ay=alp*ay
202 az=alp*az
203C----------------------------
204C CALCUL ACCELERATIONS
205C----------------------------
206 a(1,n1)=a(1,n1)+ax/ms(n1)
207 a(2,n1)=a(2,n1)+ay/ms(n1)
208 a(3,n1)=a(3,n1)+az/ms(n1)
209 a(1,n2)=a(1,n2)-ax/ms(n2)
210 a(2,n2)=a(2,n2)-ay/ms(n2)
211 a(3,n2)=a(3,n2)-az/ms(n2)
212 ELSE
213C----------------------------
214C ROTATION DU CDG
215C----------------------------
216 xx1=x(1,n1)-xcdg
217 yy1=x(2,n1)-ycdg
218 zz1=x(3,n1)-zcdg
219 xx2=x(1,n2)-xcdg
220 yy2=x(2,n2)-ycdg
221 zz2=x(3,n2)-zcdg
222C
223 i1 = (xx1*xx1+yy1*yy1+zz1*zz1)*ms(n1) + in(n1)
224 i2 = (xx2*xx2+yy2*yy2+zz2*zz2)*ms(n2) + in(n2)
225 iner = i1 + i2
226 inm1 = one/iner
227C
228 vrx1= vr(1,n1)+ar(1,n1)*dt12
229 vry1= vr(2,n1)+ar(2,n1)*dt12
230 vrz1= vr(3,n1)+ar(3,n1)*dt12
231 vrx2= vr(1,n2)+ar(1,n2)*dt12
232 vry2= vr(2,n2)+ar(2,n2)*dt12
233 vrz2= vr(3,n2)+ar(3,n2)*dt12
234C
235 vxx = (vrx1*in(n1)+yy1*vmz1-zz1*vmy1
236 . + vrx2*in(n2)+yy2*vmz2-zz2*vmy2)*inm1
237 vyy = (vry1*in(n1)+zz1*vmx1-xx1*vmz1
238 . + vry2*in(n2)+zz2*vmx2-xx2*vmz2)*inm1
239 vzz = (vrz1*in(n1)+xx1*vmy1-yy1*vmx1
240 . + vrz2*in(n2)+xx2*vmy2-yy2*vmx2)*inm1
241C----------------------------
242C CALCUL FORCES
243C----------------------------
244 vt1(1) = zz1*vyy - yy1*vzz
245 vt1(2) = xx1*vzz - zz1*vxx
246 vt1(3) = yy1*vxx - xx1*vyy
247 vt2(1) = zz2*vyy - yy2*vzz
248 vt2(2) = xx2*vzz - zz2*vxx
249 vt2(3) = yy2*vxx - xx2*vyy
250C
251 ax = (-a(1,n2)+(vx+vt2(1)-v(1,n2))*dt12m1)*ms(n2)
252 ay = (-a(2,n2)+(vy+vt2(2)-v(2,n2))*dt12m1)*ms(n2)
253 az = (-a(3,n2)+(vz+vt2(3)-v(3,n2))*dt12m1)*ms(n2)
254C
255 ax = (-a(1,n1)+(vx+vt1(1)-v(1,n1))*dt12m1)*ms(n1)
256 ay = (-a(2,n1)+(vy+vt1(2)-v(2,n1))*dt12m1)*ms(n1)
257 az = (-a(3,n1)+(vz+vt1(3)-v(3,n1))*dt12m1)*ms(n1)
258C
259 IF(dx2(is)>em15)THEN
260C POINTS NON CONFONDUS CRITERE FN ET FT
261 s = one/sqrt(dx2(is))
262 xn(is) =xn(is)*s
263 yn(is) =yn(is)*s
264 zn(is) =zn(is)*s
265 an =ax*xn(is)+ay*yn(is)+az*zn(is)
266 an2=an**2
267 anx=an*xn(is)
268 any=an*yn(is)
269 anz=an*zn(is)
270 atx=ax-anx
271 aty=ay-any
272 atz=az-anz
273 at2=(atx**2+aty**2+atz**2)
274 ELSE
275C POINTS CONFONDUS CRITERE UNIQUEMENT SUR FN
276 an2=(ax**2+ay**2+az**2)
277 at2=zero
278 an = sqrt(an2)
279 ENDIF
280C
281 alp=sqrt((an2/fn2)+(at2/ft2))
282 alp=off / max(alp,one)
283 ax=alp*ax
284 ay=alp*ay
285 az=alp*az
286C----------------------------
287C CALCUL ACCELERATIONS
288C----------------------------
289 da1(1) = half*dt2*dt12m1*(vyy*vt1(3) - vzz*vt1(2))
290 da1(2) = half*dt2*dt12m1*(vzz*vt1(1) - vxx*vt1(3))
291 da1(3) = half*dt2*dt12m1*(vxx*vt1(1) - vyy*vt1(2))
292 da2(1) = half*dt2*dt12m1*(vyy*vt2(3) - vzz*vt2(2))
293 da2(2) = half*dt2*dt12m1*(vzz*vt2(1) - vxx*vt2(3))
294 da2(3) = half*dt2*dt12m1*(vxx*vt2(1) - vyy*vt2(2))
295C
296 a(1,n1)=a(1,n1)+ax/ms(n1)+da1(1)
297 a(2,n1)=a(2,n1)+ay/ms(n1)+da1(2)
298 a(3,n1)=a(3,n1)+az/ms(n1)+da1(3)
299 a(1,n2)=a(1,n2)-ax/ms(n2)+da2(1)
300 a(2,n2)=a(2,n2)-ay/ms(n2)+da2(2)
301 a(3,n2)=a(3,n2)-az/ms(n2)+da2(3)
302 rivet(2, i) = an
303 rivet(3, i) = sqrt(at2)
304 rivet(4, i) = ax
305 rivet(5, i) = ay
306 rivet(6, i) = az
307 rivet(7, i) = a(1,n1)*ms(n1)
308 rivet(8 ,i) = a(2,n1)*ms(n1)
309 rivet(9 ,i) = a(3,n1)*ms(n1)
310C
311 amx=-ar(1,n1)+(vxx-vr(1,n1))*dt12m1
312 amy=-ar(2,n1)+(vyy-vr(2,n1))*dt12m1
313 amz=-ar(3,n1)+(vzz-vr(3,n1))*dt12m1
314 ar(1,n1)=ar(1,n1)+amx*alp
315 ar(2,n1)=ar(2,n1)+amy*alp
316 ar(3,n1)=ar(3,n1)+amz*alp
317 rivet(10,i) = ar(1,n1)*i1
318 rivet(11,i) = ar(2,n1)*i1
319 rivet(12,i) = ar(3,n1)*i1
320C
321 amx=-ar(1,n2)+(vxx-vr(1,n2))*dt12m1
322 amy=-ar(2,n2)+(vyy-vr(2,n2))*dt12m1
323 amz=-ar(3,n2)+(vzz-vr(3,n2))*dt12m1
324 ar(1,n2)=ar(1,n2)+amx*alp
325 ar(2,n2)=ar(2,n2)+amy*alp
326 ar(3,n2)=ar(3,n2)+amz*alp
327C
328 rivet(10,i) = ar(1,n2)*i2
329 rivet(11,i) = ar(2,n2)*i2
330 rivet(12,i) = ar(3,n2)*i2
331 rivet(13,i) = zero
332C correction energie cinetique de rotation
333 ww = vxx**2+vyy**2+vzz**2
334 wt = (vyy*zn(is)-vzz*yn(is))**2
335 . + (vzz*xn(is)-vxx*zn(is))**2
336 . + (vxx*yn(is)-vyy*xn(is))**2
337 enrot_l = enrot_l + half*iner*(ww-wt)
338 ENDIF
339 ENDDO
340C-----------------------------------------------
341C OLD RIGID LINK formulation
342C-----------------------------------------------
343#include "vectorize.inc"
344 DO j = 1, nind2
345 is = ind2(j)
346 i = js+is
347 ig=ixrt(1,i)
348 n1=ixrt(2,i)
349 n2=ixrt(3,i)
350 off = rivet(1,i)
351 fn2 = geo(1,ig)
352 ft2 = geo(2,ig)
353 irot=nint(geo(4,ig))
354 xm=ms(n1)*ms(n2)/(ms(n1)+ms(n2))
355 amx=(a(1,n2)-a(1,n1)+(v(1,n2)-v(1,n1))/dt12)*xm
356 amy=(a(2,n2)-a(2,n1)+(v(2,n2)-v(2,n1))/dt12)*xm
357 amz=(a(3,n2)-a(3,n1)+(v(3,n2)-v(3,n1))/dt12)*xm
358 IF(dx2(is)>em15)THEN
359C POINTS NON CONFONDUS CRITERE FN ET FT
360 s = one/sqrt(dx2(is))
361 xn(is) =xn(is)*s
362 yn(is) =yn(is)*s
363 zn(is) =zn(is)*s
364 an =amx*xn(is)+amy*yn(is)+amz*zn(is)
365 an2=an**2
366 anx=an*xn(is)
367 any=an*yn(is)
368 anz=an*zn(is)
369 atx=amx-anx
370 aty=amy-any
371 atz=amz-anz
372 at2=(atx**2+aty**2+atz**2)
373 ELSE
374C POINTS CONFONDUS CRITERE UNIQUEMENT SUR FN
375 an2=(amx**2+amy**2+amz**2)
376 at2=zero
377 ENDIF
378 alp=sqrt((an2/fn2)+(at2/ft2))
379 alp=off / max(alp,one)
380 amx=alp*amx
381 amy=alp*amy
382 amz=alp*amz
383 a(1,n1)=a(1,n1)+amx/ms(n1)
384 a(2,n1)=a(2,n1)+amy/ms(n1)
385 a(3,n1)=a(3,n1)+amz/ms(n1)
386 a(1,n2)=a(1,n2)-amx/ms(n2)
387 a(2,n2)=a(2,n2)-amy/ms(n2)
388 a(3,n2)=a(3,n2)-amz/ms(n2)
389 IF(irot==1)THEN
390 inm1=one/(in(n1)+in(n2))
391 ar(1,n1)=(ar(1,n1)*in(n1)+ar(1,n2)*in(n2))*inm1
392 ar(2,n1)=(ar(2,n1)*in(n1)+ar(2,n2)*in(n2))*inm1
393 ar(3,n1)=(ar(3,n1)*in(n1)+ar(3,n2)*in(n2))*inm1
394 ar(1,n2)=ar(1,n1)
395 ar(2,n2)=ar(2,n1)
396 ar(3,n2)=ar(3,n1)
397 ENDIF
398C-----------------------------------------------
399 ENDDO
400 js = js + nn
401 ENDDO
402C
403#include "lockon.inc"
404 enrot = enrot + enrot_l
405#include "lockoff.inc"
406 IF (rupt==1) THEN
407 DO i=lft,llt
408 IF(rivet(1,i)==-one)THEN
409 rivet(1,i)=zero
410#include "lockon.inc"
411 WRITE(istdo,*)' FAILURE OF RIVET',ixrt(4,i)
412 WRITE(iout,*) ' FAILURE OF RIVET',ixrt(4,i)
413#include "lockoff.inc"
414 ENDIF
415 ENDDO
416 ENDIF
417C
418 RETURN
419 END
my_real
#define my_real
Definition cppsort.cpp:32
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
rivet1
subroutine rivet1(ms, in, a, ar, x, ixrt, rivet, geo, v, vr, itask)
Definition rivet1.F:31