35
36
37
38 USE elbufdef_mod
39
40
41
42#include "implicit_f.inc"
43
44
45
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"
54
55
56
57 INTEGER ITAB(*),WEIGHT(*),IXR(NIXR,*),
58 . IPART(LIPART1,*),IPARTR(*),IGEO(NPROPGI,*),NPBY(NNPBY,*),
59 . IPARG(NPARG,*)
60
61 my_real stifn(*), stifr(*),ms(*) ,in(*),x(3,*),
62 . dmas,diner,geo(npropg,*)
63
64 TYPE(ELBUF_STRUCT_), TARGET, DIMENSION(NGROUP) :: ELBUF_TAB
65
66
67
68 INTEGER I,J,K,M1,M2,IG,IGTYP,N1,N2,KAD,ITYP,NG,JFT,JLT,NEL,
69 . ,FLAG,NFT,NUVAR,JNTYP,IRB1,IRB2,FLAG_S,FLAG_PR,NV
70 my_real mass,iner1,iner2,km1,krm1,km2,krm2,xx,kx1,kx2
71 . xl,kxmax,krmax,kxmin1,kxmin2,kx,kr,dtcg,scf,get_u_geo,
72 . xx1,xx2,dta,mass1,mass2
73
74 TYPE(G_BUFEL_),POINTER :: GBUF
75
76 EXTERNAL get_u_geo
77
78
79 flag_pr = 0
80
81 dtc = dt2
82 dta = dtmin1(11)
83
84 IF (tt==0) THEN
85
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"
96 ELSE
97 dtc = dta
98 ENDIF
99 ENDIF
100
102
103 DO ng=1,ngroup
104 ityp = iparg(5,ng)
105 nel = iparg(2,ng)
106 nft = iparg(3,ng)
107 jft = 1
109 gbuf => elbuf_tab(ng)%GBUF
110 IF (ityp /= 6) GOTO 250
111
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
124
125 IF (flag_pr==0) THEN
126 WRITE(iout,100)
127 WRITE(iout,200)
128 flag_pr = 1
129 ENDIF
130
131 n1 = ixr(2,j)
132 n2 = ixr(3,j)
133
134 irb1 = nint(gbuf%VAR(nv + 37))
135 IF (irb1 > 0) THEN
136 m1 = npby(1,irb1)
137 ELSE
138 m1 = n1
139 ENDIF
140
141 irb2 = nint(gbuf%VAR(nv + 38))
142 IF (irb2 > 0) THEN
143 m2 = npby(1,irb2)
144 ELSE
145 m2 = n2
146 ENDIF
147
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)
155
156
157 mass1 = ms(m1)
158 iner1 = in(m1)
159 km1 = stifn(m1)
160 krm1 = stifr(m1)
161
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
171
172
173 mass2 = ms(m2)
174 iner2 = in(m2)
175 km2 = stifn(m2)
176 krm2 = stifr(m2)
177
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
186
187 kxmax =
min(kx1,kx2,kxmax)
189
190
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
194
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
202
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
210
211 IF (flag_s==1) THEN
212 kr =
max(kr,1.3*kx*(xx+xl))
213 ENDIF
214
215 kx = scf*kx
216 kr = scf*kr
217
218 WRITE(iout,'(4X,I10,5X,I2,8X,1PE11.4,8X,1PE11.4)')
219 . ixr(nixr,j),jntyp,kx,kr
220
221 gbuf%VAR(nv + 16) = kx
222 gbuf%VAR(nv + 17) = kr
223
224
225 IF (jntyp==1) THEN
226 gbuf%VAR(nv + 18) = kx
227 gbuf%VAR(nv + 19) = kx
228 gbuf%VAR(nv + 20) = kx
229
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
236
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
242
243 ELSEIF (jntyp==4) THEN
244 gbuf%VAR(nv + 18) = kx
245 gbuf%VAR(nv + 31) = kr
246 gbuf%VAR(nv + 32) = kr
247
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
253
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
260
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
266
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
275
276 ENDIF
277 ENDIF
278
279 ENDDO
280250 CONTINUE
281 ENDDO
282
283 WRITE(iout,'( )')
284
285
286 ENDIF
287
288 RETURN
289100 FORMAT(/,
290 & 1X,'automatic stiffness computation
for joints
'/)
291200 FORMAT(1X,'joint
id',11X,'type
',6X,'knn
',16X,'knr
')
292
293300 FORMAT(1X,'warning trans. stiff. imposed by stiff. on rbody')
294400 FORMAT(1X,'warning rot. stiff. imposed by stiff. on rbody')
295
2961300 FORMAT(1X,'warning trans. stiff. imposed by stiff. on connected structures')
2971400 FORMAT(1X,'warning')
for(i8=*sizetab-1;i8 >=0;i8--)