OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
joint_block_stiffness.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!|| joint_block_stiffness ../engine/source/elements/joint/joint_block_stiffness.F
25!||--- called by ------------------------------------------------------
26!|| resol ../engine/source/engine/resol.F
27!||--- calls -----------------------------------------------------
28!|| arret ../engine/source/system/arret.F
29!||--- uses -----------------------------------------------------
30!|| elbufdef_mod ../common_source/modules/mat_elem/elbufdef_mod.F90
31!||====================================================================
32 SUBROUTINE joint_block_stiffness(ITAB,MS,IN,STIFN,STIFR,
33 1 WEIGHT ,IXR,IPART,X,
34 2 IPARTR,IGEO,GEO,NPBY,IPARG,ELBUF_TAB,DMAS,DINER)
35C-----------------------------------------------
36C M o d u l e s
37C-----------------------------------------------
38 USE elbufdef_mod
39C----6---------------------------------------------------------------7---------8
40C I m p l i c i t T y p e s
41C-----------------------------------------------
42#include "implicit_f.inc"
43C-----------------------------------------------
44C G l o b a l P a r a m e t e r s
45C-----------------------------------------------
46#include "com08_c.inc"
47#include "com01_c.inc"
48#include "param_c.inc"
49#include "scr02_c.inc"
50#include "scr17_c.inc"
51#include "scr18_c.inc"
52#include "cong2_c.inc"
53#include "units_c.inc"
54C-----------------------------------------------------------------
55C D u m m y A r g u m e n t s
56C-----------------------------------------------
57 INTEGER ITAB(*),WEIGHT(*),IXR(NIXR,*),
58 . IPART(LIPART1,*),IPARTR(*),IGEO(NPROPGI,*),NPBY(NNPBY,*),
59 . iparg(nparg,*)
60C REAL
61 my_real stifn(*), stifr(*),ms(*) ,in(*),x(3,*),
62 . dmas,diner,geo(npropg,*)
63C
64 TYPE(elbuf_struct_), TARGET, DIMENSION(NGROUP) :: ELBUF_TAB
65C-----------------------------------------------
66C L o c a l V a r i a b l e s
67C-----------------------------------------------
68 INTEGER I,J,K,M1,M2,IG,IGTYP,N1,N2,KAD,ITYP,NG,JFT,JLT,NEL,
69 . NB8,FLAG,NFT,NUVAR,JNTYP,IRB1,IRB2,FLAG_S,FLAG_PR,NV
70 my_real mass,iner1,iner2,km1,krm1,km2,krm2,xx,kx1,kx2,kr1,dtc,alpha,
71 . xl,kxmax,krmax,kxmin1,kxmin2,kx,kr,dtcg,scf,get_u_geo,
72 . xx1,xx2,dta,mass1,mass2
73C
74 TYPE(g_bufel_),POINTER :: GBUF
75C-----------------------------------------------
76 EXTERNAL get_u_geo
77C----------------------------------------------------------
78
79 flag_pr = 0
80
81 dtc = dt2
82 dta = dtmin1(11)
83
84 IF (tt==0) THEN
85C----------------------------------------------------------
86 IF (nodadt>0) THEN
87 dtc = max(dtc,dtmin1(11)/dtfac1(11))
88 IF (dtmx>em20) dtc = min(dtc,dtmx)
89 ENDIF
90
91 IF (dtc>90000) THEN
92 IF (dta==0) THEN
93 print *,"ERROR NO TARGET TIME STEP DT=",dtc
94 print *,"STIFFNESS CAN NOT BE COMPUTED"
95 CALL arret(2)
96 ELSE
97 dtc = dta
98 ENDIF
99 ENDIF
100
101 alpha = 4.1
102
103 DO ng=1,ngroup
104 ityp = iparg(5,ng)
105 nel = iparg(2,ng)
106 nft = iparg(3,ng)
107 jft = 1
108 jlt = min(nvsiz,nel)
109 gbuf => elbuf_tab(ng)%GBUF
110 IF (ityp /= 6) GOTO 250
111C--------> Boucle sur les elements ressort -------
112 DO i=jft,jlt
113 j = i + nft
114 ig = ipart(2,ipartr(j))
115 igtyp = igeo(11,ig)
116 nuvar = nint(geo(25,ig))
117 nv = nuvar*(i-1) + 1
118 IF (igtyp==45) THEN
119 scf = get_u_geo(11,ig)
120 jntyp = nint(get_u_geo(1,ig))
121 flag = nint(get_u_geo(10,ig))
122 flag_s = 0
123 IF (flag==0) THEN
124C
125 IF (flag_pr==0) THEN
126 WRITE(iout,100)
127 WRITE(iout,200)
128 flag_pr = 1
129 ENDIF
130C
131 n1 = ixr(2,j)
132 n2 = ixr(3,j)
133C
134 irb1 = nint(gbuf%VAR(nv + 37))
135 IF (irb1 > 0) THEN
136 m1 = npby(1,irb1)
137 ELSE
138 m1 = n1
139 ENDIF
140C
141 irb2 = nint(gbuf%VAR(nv + 38))
142 IF (irb2 > 0) THEN
143 m2 = npby(1,irb2)
144 ELSE
145 m2 = n2
146 ENDIF
147C
148 xl = ((x(1,n1)-x(1,n2))**2)+((x(2,n1)-x(2,n2))**2)
149 . +((x(3,n1)-x(3,n2))**2)
150 xx1 = ((x(1,n1)-x(1,m1))**2)+((x(2,n1)-x(2,m1))**2)
151 . +((x(3,n1)-x(3,m1))**2)
152 xx2 = ((x(1,n2)-x(1,m2))**2)+((x(2,n2)-x(2,m2))**2)
153 . +((x(3,n2)-x(3,m2))**2)
154 xx = max(xx1,xx2)
155C
156C--------> Calcul raideur cote noeud 1 -------
157 mass1 = ms(m1)
158 iner1 = in(m1)
159 km1 = stifn(m1)
160 krm1 = stifr(m1)
161C
162 kx1 = (2*mass1/(alpha*dtc*dtc)) - km1
163 IF (iner1 > zero) THEN
164 kx2 = 0.8*(iner1/(alpha*dtc*dtc)- krm1)/(max(em20,(xx+xl)))
165 kr = iner1/(alpha*dtc*dtc)- krm1
166 ELSE
167 kx2 = ep30
168 kr = zero
169 ENDIF
170 kxmax = min(kx1,kx2)
171C
172C--------> Calcul raideur cote noeud 2 -------
173 mass2 = ms(m2)
174 iner2 = in(m2)
175 km2 = stifn(m2)
176 krm2 = stifr(m2)
177C
178 kx1 = (2*mass2/(alpha*dtc*dtc)) - km2
179 IF (iner2 > zero) THEN
180 kx2 = 0.8*(iner2/(alpha*dtc*dtc)- krm2)/(max(em20,(xx+xl)))
181 kr1 = iner2/(alpha*dtc*dtc)- krm2
182 ELSE
183 kx2 = ep30
184 kr1 = zero
185 ENDIF
186C
187 kxmax = min(kx1,kx2,kxmax)
188 kr = min(kr,kr1)
189
190C--------> Calcul raideur finale et affectation-------
191 kx = max(kxmax,2*km1,2*km2)
192 IF ((kx - kxmax)>1e-8) flag_s = 1
193 IF ((iner1 == zero).OR.(iner2 == zero)) THEN
194C-- Raideur en rotation mise zero si kjoint attach un noeud de solid --
195 kr = zero
196 krmax = kr
197 flag_s = 0
198 ELSE
199 krmax = kr
200 kr = max(kr,2*krm1,2*krm2)
201 ENDIF
202C
203 IF ((irb1 > 0) .AND.(irb2 > 0)) THEN
204 IF ((kx - kxmax)>1e-8) WRITE(iout,300)
205 IF ((kr - krmax)>1e-8) WRITE(iout,400)
206 ELSE
207 IF ((kx - kxmax)>1e-8) WRITE(iout,1300)
208 IF ((kr - krmax)>1e-8) WRITE(iout,1400)
209 ENDIF
210C
211 IF (flag_s==1) THEN
212 kr = max(kr,1.3*kx*(xx+xl))
213 ENDIF
214C
215 kx = scf*kx
216 kr = scf*kr
217C
218 WRITE(iout,'(4X,I10,5X,I2,8X,1PE11.4,8X,1PE11.4)')
219 . ixr(nixr,j),jntyp,kx,kr
220C
221 gbuf%VAR(nv + 16) = kx
222 gbuf%VAR(nv + 17) = kr
223
224C--------> Spherical Joint-------
225 IF (jntyp==1) THEN
226 gbuf%VAR(nv + 18) = kx
227 gbuf%VAR(nv + 19) = kx
228 gbuf%VAR(nv + 20) = kx
229C--------> Revolutional Joint-------
230 ELSEIF (jntyp==2) THEN
231 gbuf%VAR(nv + 18) = kx
232 gbuf%VAR(nv + 19) = kx
233 gbuf%VAR(nv + 20) = kx
234 gbuf%VAR(nv + 31) = kr
235 gbuf%VAR(nv + 32) = kr
236C--------> Cylindrical Joint-------
237 ELSEIF (jntyp==3) THEN
238 gbuf%VAR(nv + 19) = kx
239 gbuf%VAR(nv + 20) = kx
240 gbuf%VAR(nv + 31) = kr
241 gbuf%VAR(nv + 32) = kr
242C--------> Planar Joint-------
243 ELSEIF (jntyp==4) THEN
244 gbuf%VAR(nv + 18) = kx
245 gbuf%VAR(nv + 31) = kr
246 gbuf%VAR(nv + 32) = kr
247C--------> Universal Joint-------
248 ELSEIF (jntyp==5) THEN
249 gbuf%VAR(nv + 18) = kx
250 gbuf%VAR(nv + 19) = kx
251 gbuf%VAR(nv + 20) = kx
252 gbuf%VAR(nv + 30) = kr
253C--------> Translational Joint-------
254 ELSEIF (jntyp==6) THEN
255 gbuf%VAR(nv + 19) = kx
256 gbuf%VAR(nv + 20) = kx
257 gbuf%VAR(nv + 30) = kr
258 gbuf%VAR(nv + 31) = kr
259 gbuf%VAR(nv + 32) = kr
260C--------> Oldham Joint-------
261 ELSEIF (jntyp==7) THEN
262 gbuf%VAR(nv + 18) = kx
263 gbuf%VAR(nv + 30) = kr
264 gbuf%VAR(nv + 31) = kr
265 gbuf%VAR(nv + 32) = kr
266C--------> Rigid Joint-------
267 ELSEIF (jntyp==8) THEN
268 gbuf%VAR(nv + 18) = kx
269 gbuf%VAR(nv + 19) = kx
270 gbuf%VAR(nv + 20) = kx
271 gbuf%VAR(nv + 30) = kr
272 gbuf%VAR(nv + 31) = kr
273 gbuf%VAR(nv + 32) = kr
274 ENDIF ! IF (JNTYP)
275C
276 ENDIF ! IF (FLAG==0)
277 ENDIF ! IF (IGTYP==45)
278C
279 ENDDO
280250 CONTINUE
281 ENDDO
282
283 WRITE(iout,'( )')
284C----------------------------------------------------------
285
286 ENDIF
287C
288 RETURN
289100 FORMAT(/,
290 & 1x,'AUTOMATIC STIFFNESS COMPUTATION FOR JOINTS'/)
291200 FORMAT(1x,'JOINT ID',11x,'TYPE',6x,'KNN',16x,'KNR')
292C
293300 FORMAT(1x,'WARNING TRANS. STIFF. IMPOSED BY STIFF. ON RBODY')
294400 FORMAT(1x,'WARNING ROT. STIFF. IMPOSED BY STIFF. ON RBODY')
295C
2961300 FORMAT(1x,'WARNING TRANS. STIFF. IMPOSED BY STIFF. ON CONNECTED STRUCTURES')
2971400 FORMAT(1x,'WARNING ROT. STIFF. IMPOSED BY STIFF. ON CONNECTED STRUCTURES')
298 END
my_real
#define my_real
Definition cppsort.cpp:32
alpha
#define alpha
Definition eval.h:35
joint_block_stiffness
subroutine joint_block_stiffness(itab, ms, in, stifn, stifr, weight, ixr, ipart, x, ipartr, igeo, geo, npby, iparg, elbuf_tab, dmas, diner)
Definition joint_block_stiffness.F:35
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
arret
subroutine arret(nn)
Definition arret.F:87