OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
damping.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!|| damping51 ../engine/source/assembly/damping.F
25!||--- called by ------------------------------------------------------
26!|| resol ../engine/source/engine/resol.F
27!||--- calls -----------------------------------------------------
28!|| damping_vref ../engine/source/assembly/damping_vref.F
29!|| damping_vref_compute_dampa ../engine/source/assembly/damping_vref_compute_dampa.F90
30!|| get_u_func ../engine/source/user_interface/ufunc.F
31!||--- uses -----------------------------------------------------
32!|| damping_vref_compute_dampa_mod ../engine/source/assembly/damping_vref_compute_dampa.f90
33!|| groupdef_mod ../common_source/modules/groupdef_mod.f
34!||====================================================================
35 SUBROUTINE damping51(
36 1 DIM ,V ,
37 2 VR ,A ,AR ,MS ,IN ,
38 3 DAMPR ,DAMP ,IGRNOD ,WEIGHT ,TAGSLV_RBY,
39 4 SKEW ,ICONTACT, I_DAMP_RDOF_TAB,NDAMP_VREL,
40 5 ID_DAMP_VREL,FR_DAMP_VREL,IPARIT,ISPMD,WFEXT)
41C-----------------------------------------------
42C M o d u l e s
43C-----------------------------------------------
44 USE groupdef_mod
45 use damping_vref_compute_dampa_mod
46C-----------------------------------------------
47C I m p l i c i t T y p e s
48C-----------------------------------------------
49#include "implicit_f.inc"
50#include "comlock.inc"
51C-----------------------------------------------
52C C o m m o n B l o c k s
53C-----------------------------------------------
54#include "com01_c.inc"
55#include "com04_c.inc"
56#include "com06_c.inc"
57#include "com08_c.inc"
58#include "sms_c.inc"
59#include "param_c.inc"
60C-----------------------------------------------
61C D u m m y A r g u m e n t s
62C-----------------------------------------------
63 INTEGER DIM,ITASK,WEIGHT(*),ICONTACT(*),I_DAMP_RDOF_TAB(*),TAGSLV_RBY(*)
64 INTEGER ,INTENT(IN) :: NDAMP_VREL,ID_DAMP_VREL(NDAMP_VREL)
65 INTEGER ,INTENT(IN) :: FR_DAMP_VREL(NSPMD+2,NDAMP_VREL),IPARIT,ISPMD
66 my_real V(3,*), VR(3,*), A(3,*), AR(3,*) ,MS(*), IN(*),DAMPR(NRDAMP,*), DAMP(DIM,*), SKEW(LSKEW,*)
67 my_real,DIMENSION(:,:), ALLOCATABLE :: VSKW,ASKW,DAMPSKW
68 DOUBLE PRECISION,INTENT(INOUT):: WFEXT
69C-----------------------------------------------
70 TYPE (GROUP_) , DIMENSION(NGRNOD) :: IGRNOD
71C-----------------------------------------------
72C L o c a l V a r i a b l e s
73C-----------------------------------------------
74 LOGICAL :: NEED_COMPUTE
75 INTEGER I,J,N,ND,IGR,ISK,FL_VREL,ID_RBY,FL_FREQ_RANGE,ITYPE
76 my_real, DIMENSION(3) :: DAMP_A,DAMP_B
77 my_real, DIMENSION(3) :: DA_,DW_LOCAL,OMEGA_,BETASDT_
78 my_real :: AA,DW,DA,BETASDT,OMEGA,FACTB,DAMPT,D_TSTART,D_TSTOP,
79 . v_ref(3,ndamp),a_ref(3,ndamp),v_refskw(3),a_refskw(3),
80 . dam(3),darm(3),dfm,get_u_func
81 EXTERNAL get_u_func
82C-----------------------------------------------
83C C = a M + b K
84C======================================================================|
85 ALLOCATE(vskw(3,numnod))
86 ALLOCATE(askw(3,numnod))
87 ALLOCATE(dampskw(3,numnod))
88
89 IF(idtmins==2.OR.idtmins_int/=0)GOTO 1000
90C
91C-----Computation of average velocity for /DAMP/VREF -------------------------------
92C
93 IF (ndamp_vrel > 0) THEN
94 CALL damping_vref(igrnod,v,a,ms,dampr,
95 . weight,v_ref,a_ref,ndamp_vrel,id_damp_vrel,
96 . fr_damp_vrel,ngrnod,numnod,ndamp,nrdamp,
97 . nspmd,iparit,ispmd)
98 ENDIF
99C-----------------------------------------------
100 dw = zero
101 DO nd=1,ndamp
102 igr = nint(dampr(2,nd))
103 isk = nint(dampr(15,nd))
104 factb = one
105! FACTB = DAMPR(16,ND) ! used for alpha of func
106 dampt = min(dt1,dt2)*factb
107 d_tstart = dampr(17,nd)
108 d_tstop = dampr(18,nd)
109 itype = nint(dampr(21,nd))
110 fl_vrel = 0
111 fl_freq_range = 0
112 SELECT CASE (itype)
113 CASE(2)
114 fl_vrel = 1
115 CASE(3)
116 fl_freq_range = 1
117 END SELECT
118 id_rby = nint(dampr(25,nd))
119C
120 IF (tt>=d_tstart .AND. tt<=d_tstop) THEN
121C
122 IF(isk<=1)THEN
123C-----------------------------------------------------------------------------------
124C---------- Damping in global coordinates system -----------------------------------
125C-----------------------------------------------------------------------------------
126C
127 IF ((fl_vrel == 0).AND.(fl_freq_range == 0)) THEN
128C-------------Damping using global velocities
129 IF (dampr(19,nd)==0) THEN
130C-------------------------------------------------
131 ! --------------
132 damp_a(1) = dampr(3,nd)
133 damp_b(1) = dampr(4,nd)
134 betasdt_(1)= -min(damp_b(1),dampt)*dt1/max(dt1*dt1,em30)
135 omega_(1) = one/ (one + half * damp_a(1) * dt1)
136 ! --------------
137 damp_a(2) = dampr(5,nd)
138 damp_b(2) = dampr(6,nd)
139 betasdt_(2)= -min(damp_b(2),dampt)*dt1/max(dt1*dt1,em30)
140 omega_(2) = one/ (one + half * damp_a(2) * dt1)
141 ! --------------
142 damp_a(3) = dampr(7,nd)
143 damp_b(3) = dampr(8,nd)
144 betasdt_(3)= -min(damp_b(3),dampt)*dt1/max(dt1*dt1,em30)
145 omega_(3) = one/ (one + half * damp_a(3) * dt1)
146 ! --------------
147
148 DO n=1,igrnod(igr)%NENTITY
149 i=igrnod(igr)%ENTITY(n)
150 IF(tagslv_rby(i)==0) THEN
151 ! --------------
152 da_(1) = a(1,i) - damp_a(1)*v(1,i) - betasdt_(1) *(a(1,i) - damp(1,i))
153 da_(1) = da_(1) * omega_(1) - a(1,i)
154 damp(1,i) = a(1,i)
155 a(1,i) = a(1,i) + da_(1)
156 ! --------------
157 da_(2) = a(2,i) - damp_a(2)*v(2,i) - betasdt_(2) *(a(2,i) - damp(2,i))
158 da_(2) = da_(2) * omega_(2) - a(2,i)
159 damp(2,i) = a(2,i)
160 a(2,i) = a(2,i) + da_(2)
161 ! --------------
162 da_(3) = a(3,i) - damp_a(3)*v(3,i) - betasdt_(3) *(a(3,i) - damp(3,i))
163 da_(3) = da_(3) * omega_(3) - a(3,i)
164 damp(3,i) = a(3,i)
165 a(3,i) = a(3,i) + da_(3)
166 ! --------------
167 dw_local(1) = da_(1)*(v(1,i)+half*a(1,i)*dt1)
168 dw_local(2) = da_(2)*(v(2,i)+half*a(2,i)*dt1)
169 dw_local(3) = da_(3)*(v(3,i)+half*a(3,i)*dt1)
170 dw =dw+ms(i)*(dw_local(1)+dw_local(2)+dw_local(3))*dt12*weight(i)
171 ENDIF
172 ENDDO
173 ENDIF
174C
175 IF(iroddl/=0)THEN
176 ! --------------
177 damp_a(1) = dampr(9,nd)
178 damp_b(1) = dampr(10,nd)
179 betasdt_(1)= -min(damp_b(1),dampt)*dt1/max(dt1*dt1,em30)
180 omega_(1) = one/ (one + half * damp_a(1) * dt1)
181 ! --------------
182 damp_a(2) = dampr(11,nd)
183 damp_b(2) = dampr(12,nd)
184 betasdt_(2)= -min(damp_b(2),dampt)*dt1/max(dt1*dt1,em30)
185 omega_(2) = one/ (one + half * damp_a(2) * dt1)
186 ! --------------
187 damp_a(3) = dampr(13,nd)
188 damp_b(3) = dampr(14,nd)
189 betasdt_(3)= -min(damp_b(3),dampt)*dt1/max(dt1*dt1,em30)
190 omega_(3) = one/ (one + half * damp_a(3) * dt1)
191 ! --------------
192
193 DO n=1,igrnod(igr)%NENTITY
194 i=igrnod(igr)%ENTITY(n)
195 IF(tagslv_rby(i)==0) THEN
196C----- Damping seulement pour les noeuds en contact avec temporisation -----
197 need_compute = .true.
198 IF (dampr(19,nd)/=0) THEN
199 IF (icontact(i)/=0) i_damp_rdof_tab(i) = dampr(19,nd)
200 IF (i_damp_rdof_tab(i)==0) need_compute = .false.
201 ENDIF
202
203 IF(need_compute) THEN
204 ! --------------
205 da_(1) = ar(1,i) - damp_a(1)*vr(1,i) - betasdt_(1) *(ar(1,i)-damp(4,i))
206 da_(1) = da_(1) * omega_(1) - ar(1,i)
207 damp(4,i) = ar(1,i)
208 ar(1,i) = ar(1,i) + da_(1)
209 ! --------------
210 da_(2) = ar(2,i) - damp_a(2)*vr(2,i) - betasdt_(2) *(ar(2,i)-damp(5,i))
211 da_(2) = da_(2) * omega_(2) - ar(2,i)
212 damp(5,i) = ar(2,i)
213 ar(2,i) = ar(2,i) + da_(2)
214 ! --------------
215 da_(3) = ar(3,i) - damp_a(3)*vr(3,i) - betasdt_(3) *(ar(3,i)-damp(6,i))
216 da_(3) = da_(3) * omega_(3) - ar(3,i)
217 damp(6,i) = ar(3,i)
218 ar(3,i) = ar(3,i) + da_(3)
219 ! --------------
220 dw_local(1) = da_(1)*(vr(1,i)+half*ar(1,i)*dt1)
221 dw_local(2) = da_(2)*(vr(2,i)+half*ar(2,i)*dt1)
222 dw_local(3) = da_(3)*(vr(3,i)+half*ar(3,i)*dt1)
223 dw = dw+in(i)*(dw_local(1)+dw_local(2)+dw_local(3)) *dt12*weight(i)
224 ! --------------
225 IF (dampr(19,nd)/=0) i_damp_rdof_tab(i)=i_damp_rdof_tab(i)-1
226 ENDIF
227 ENDIF
228 ENDDO
229 END IF
230C
231 ELSEIF ((fl_vrel == 1).AND.(id_rby == 0)) THEN
232C-------------Damping using average relative velocity - trans dof only ----
233C
234C---------------alpha = Cdamp_M / dt--------------
235 call damping_vref_compute_dampa(nd,ndamp,nrdamp,dampr,dt1,tt,damp_a)
236!
237 omega_(1) = one/ (one + half * damp_a(1) * dt1)
238 omega_(2) = one/ (one + half * damp_a(2) * dt1)
239 omega_(3) = one/ (one + half * damp_a(3) * dt1)
240C-------------------------------------------------
241
242 DO n=1,igrnod(igr)%NENTITY
243 i=igrnod(igr)%ENTITY(n)
244 IF(tagslv_rby(i)==0) THEN
245 da_(1) = a(1,i)-damp_a(1)*(v(1,i)-v_ref(1,nd))
246 . +half*damp_a(1)*dt1*a_ref(1,nd)
247 da_(1) = da_(1) * omega_(1) - a(1,i)
248 damp(1,i) = a(1,i)
249 a(1,i) = a(1,i) + da_(1)
250 dw =dw+ms(i)*da_(1)*(v(1,i)+half*a(1,i)*dt1)*dt12*weight(i)
251 ENDIF
252 ENDDO
253C-------------------------------------------------
254
255 DO n=1,igrnod(igr)%NENTITY
256 i=igrnod(igr)%ENTITY(n)
257 IF(tagslv_rby(i)==0) THEN
258 da_(2) = a(2,i)-damp_a(2)*(v(2,i)-v_ref(2,nd))
259 . +half*damp_a(2)*dt1*a_ref(2,nd)
260 da_(2) = da_(2) * omega_(2) - a(2,i)
261 damp(2,i) = a(2,i)
262 a(2,i) = a(2,i) + da_(2)
263 dw =dw+ms(i)*da_(2)*(v(2,i)+half*a(2,i)*dt1)*dt12*weight(i)
264 ENDIF
265 ENDDO
266C-------------------------------------------------
267 DO n=1,igrnod(igr)%NENTITY
268 i=igrnod(igr)%ENTITY(n)
269 IF(tagslv_rby(i)==0) THEN
270 da_(3) = a(3,i)-damp_a(3)*(v(3,i)-v_ref(3,nd))
271 . +half*damp_a(3)*dt1*a_ref(3,nd)
272 da_(3) = da_(3) * omega_(3) - a(3,i)
273 damp(3,i) = a(3,i)
274 a(3,i) = a(3,i) + da_(3)
275 dw =dw+ms(i)*da_(3)*(v(3,i)+half*a(3,i)*dt1)*dt12*weight(i)
276 ENDIF
277 ENDDO
278C-------------------------------------------------
279 ENDIF
280C
281 ELSE
282C-----------------------------------------------------------------------------------
283C---------- Damping in local coordinates system ------------------------------------
284C-----------------------------------------------------------------------------------
285C
286 IF ((fl_vrel == 0).AND.(fl_freq_range == 0)) THEN
287C-------------Damping using global velocities
288 IF (dampr(19,nd)==0) THEN
289C-------------------------------------------------
290#include "vectorize.inc"
291 DO n=1,igrnod(igr)%NENTITY
292 i=igrnod(igr)%ENTITY(n)
293 vskw(1,i)= skew(1,isk)*v(1,i)+skew(2,isk)*v(2,i)+skew(3,isk)*v(3,i)
294 vskw(2,i)= skew(4,isk)*v(1,i)+skew(5,isk)*v(2,i)+skew(6,isk)*v(3,i)
295 vskw(3,i)= skew(7,isk)*v(1,i)+skew(8,isk)*v(2,i)+skew(9,isk)*v(3,i)
296 askw(1,i)= skew(1,isk)*a(1,i)+skew(2,isk)*a(2,i)+skew(3,isk)*a(3,i)
297 askw(2,i)= skew(4,isk)*a(1,i)+skew(5,isk)*a(2,i)+skew(6,isk)*a(3,i)
298 askw(3,i)= skew(7,isk)*a(1,i)+skew(8,isk)*a(2,i)+skew(9,isk)*a(3,i)
299 dampskw(1,i)= skew(1,isk)*damp(1,i)+skew(2,isk)*damp(2,i)+skew(3,isk)*damp(3,i)
300 dampskw(2,i)= skew(4,isk)*damp(1,i)+skew(5,isk)*damp(2,i)+skew(6,isk)*damp(3,i)
301 dampskw(3,i)= skew(7,isk)*damp(1,i)+skew(8,isk)*damp(2,i)+skew(9,isk)*damp(3,i)
302 END DO
303 damp_a(1) = dampr(3,nd)
304 damp_b(1) = dampr(4,nd)
305 betasdt_(1)= -min(damp_b(1),dampt)*dt1/max(dt1*dt1,em30)
306 omega_(1) = one/ (one + half * damp_a(1) * dt1)
307
308 damp_a(2) = dampr(5,nd)
309 damp_b(2) = dampr(6,nd)
310 betasdt_(2)= -min(damp_b(2),dampt)*dt1/max(dt1*dt1,em30)
311 omega_(2) = one/ (one + half * damp_a(2) * dt1)
312
313 damp_a(3) = dampr(7,nd)
314 damp_b(3) = dampr(8,nd)
315 betasdt_(3)= -min(damp_b(3),dampt)*dt1/max(dt1*dt1,em30)
316 omega_(3) = one/ (one + half * damp_a(3) * dt1)
317
318 DO n=1,igrnod(igr)%NENTITY
319 i=igrnod(igr)%ENTITY(n)
320 IF(tagslv_rby(i)==0) THEN
321 da_(1) = askw(1,i) - damp_a(1)*vskw(1,i) - betasdt_(1) *(askw(1,i) - dampskw(1,i))
322 da_(1) = da_(1) * omega_(1) - askw(1,i)
323 dampskw(1,i) = askw(1,i)
324 askw(1,i) = askw(1,i) + da_(1)
325 dw_local(1)= da_(1)*(vskw(1,i)+half*askw(1,i)*dt1)
326
327 da_(2) = askw(2,i) - damp_a(2)*vskw(2,i)- betasdt_(2) *(askw(2,i) - dampskw(2,i))
328 da_(2) = da_(2) * omega_(2) - askw(2,i)
329 dampskw(2,i) = askw(2,i)
330 askw(2,i) = askw(2,i) + da_(2)
331 dw_local(2)= da_(2)*(vskw(2,i)+half*askw(2,i)*dt1)
332
333 da_(3) = askw(3,i) - damp_a(3)*vskw(3,i) - betasdt_(3) *(askw(3,i) - dampskw(3,i))
334 da_(3) = da_(3) * omega_(3) - askw(3,i)
335 dampskw(3,i) = askw(3,i)
336 askw(3,i) = askw(3,i) + da_(3)
337 dw_local(3)= da_(3)*(vskw(3,i)+half*askw(3,i)*dt1)
338
339 dw =dw +ms(i)*(dw_local(1)+dw_local(2)+dw_local(3))*dt12*weight(i)
340 ENDIF
341 ENDDO
342
343#include "vectorize.inc"
344 DO n=1,igrnod(igr)%NENTITY
345 i=igrnod(igr)%ENTITY(n)
346 a(1,i)= skew(1,isk)*askw(1,i)
347 . +skew(4,isk)*askw(2,i)
348 . +skew(7,isk)*askw(3,i)
349 a(2,i)= skew(2,isk)*askw(1,i)
350 . +skew(5,isk)*askw(2,i)
351 . +skew(8,isk)*askw(3,i)
352 a(3,i)= skew(3,isk)*askw(1,i)
353 . +skew(6,isk)*askw(2,i)
354 . +skew(9,isk)*askw(3,i)
355 damp(1,i)= skew(1,isk)*dampskw(1,i)
356 . +skew(4,isk)*dampskw(2,i)
357 . +skew(7,isk)*dampskw(3,i)
358 damp(2,i)= skew(2,isk)*dampskw(1,i)
359 . +skew(5,isk)*dampskw(2,i)
360 . +skew(8,isk)*dampskw(3,i)
361 damp(3,i)= skew(3,isk)*dampskw(1,i)
362 . +skew(6,isk)*dampskw(2,i)
363 . +skew(9,isk)*dampskw(3,i)
364 END DO
365 ENDIF
366
367 IF(iroddl/=0)THEN
368#include "vectorize.inc"
369 DO n=1,igrnod(igr)%NENTITY
370 i=igrnod(igr)%ENTITY(n)
371 vskw(1,i)= skew(1,isk)*vr(1,i)
372 . +skew(2,isk)*vr(2,i)
373 . +skew(3,isk)*vr(3,i)
374 vskw(2,i)= skew(4,isk)*vr(1,i)
375 . +skew(5,isk)*vr(2,i)
376 . +skew(6,isk)*vr(3,i)
377 vskw(3,i)= skew(7,isk)*vr(1,i)
378 . +skew(8,isk)*vr(2,i)
379 . +skew(9,isk)*vr(3,i)
380 askw(1,i)= skew(1,isk)*ar(1,i)
381 . +skew(2,isk)*ar(2,i)
382 . +skew(3,isk)*ar(3,i)
383 askw(2,i)= skew(4,isk)*ar(1,i)
384 . +skew(5,isk)*ar(2,i)
385 . +skew(6,isk)*ar(3,i)
386 askw(3,i)= skew(7,isk)*ar(1,i)
387 . +skew(8,isk)*ar(2,i)
388 . +skew(9,isk)*ar(3,i)
389 dampskw(1,i)= skew(1,isk)*damp(4,i)
390 . +skew(2,isk)*damp(5,i)
391 . +skew(3,isk)*damp(6,i)
392 dampskw(2,i)= skew(4,isk)*damp(4,i)
393 . +skew(5,isk)*damp(5,i)
394 . +skew(6,isk)*damp(6,i)
395 dampskw(3,i)= skew(7,isk)*damp(4,i)
396 . +skew(8,isk)*damp(5,i)
397 . +skew(9,isk)*damp(6,i)
398 END DO
399
400 damp_a(1) = dampr(9,nd)
401 damp_b(1) = dampr(10,nd)
402 betasdt_(1)= -min(damp_b(1),dampt)*dt1/max(dt1*dt1,em30)
403 omega_(1) = one/ (one + half * damp_a(1) * dt1)
404
405 damp_a(2) = dampr(11,nd)
406 damp_b(2) = dampr(12,nd)
407 betasdt_(2)= -min(damp_b(2),dampt)*dt1/max(dt1*dt1,em30)
408 omega_(2) = one/ (one + half * damp_a(2) * dt1)
409
410 damp_a(3) = dampr(13,nd)
411 damp_b(3) = dampr(14,nd)
412 betasdt_(3)= -min(damp_b(3),dampt)*dt1/max(dt1*dt1,em30)
413 omega_(3) = one/ (one + half * damp_a(3) * dt1)
414
415 DO n=1,igrnod(igr)%NENTITY
416 i=igrnod(igr)%ENTITY(n)
417 IF(tagslv_rby(i)==0) THEN
418C----- Damping seulement pour les noeuds en contact avec temporisation -----
419 need_compute = .true.
420 IF (dampr(19,nd)/=0) THEN
421 IF (icontact(i)/=0) i_damp_rdof_tab(i) = dampr(19,nd)
422 IF (i_damp_rdof_tab(i)==0) need_compute=.false.
423 ENDIF
424C--------------------------------------------------------
425 IF(need_compute) THEN
426 da_(1) = askw(1,i) - damp_a(1)*vskw(1,i)
427 . - betasdt_(1) *(askw(1,i) - dampskw(1,i))
428 da_(1) = da_(1) * omega_(1) - askw(1,i)
429 dampskw(1,i) = askw(1,i)
430 askw(1,i) = askw(1,i) + da_(1)
431 dw_local(1)=da_(1)*(vskw(1,i)+half*askw(1,i)*dt1)
432
433 da_(2) = askw(2,i) - damp_a(2)*vskw(2,i)
434 . - betasdt_(2) *(askw(2,i) - dampskw(2,i))
435 da_(2) = da_(2) * omega_(2) - askw(2,i)
436 dampskw(2,i) = askw(2,i)
437 askw(2,i) = askw(2,i) + da_(2)
438 dw_local(2) =da_(2)*(vskw(2,i)+half*askw(2,i)*dt1)
439
440 da_(3) = askw(3,i) - damp_a(3)*vskw(3,i)
441 . - betasdt_(3) *(askw(3,i) - dampskw(3,i))
442 da_(3) = da_(3) * omega_(3) - askw(3,i)
443 dampskw(3,i) = askw(3,i)
444 askw(3,i) = askw(3,i) + da_(3)
445 dw_local(3)= da_(3)*(vskw(3,i)+half*askw(3,i)*dt1)
446
447 dw =dw +in(i)*(dw_local(1)+dw_local(2)+dw_local(3))*dt12*weight(i)
448 IF (dampr(19,nd)/=0) i_damp_rdof_tab(i)=i_damp_rdof_tab(i)-1
449 ENDIF
450 ENDIF
451 ENDDO
452 ENDIF
453#include "vectorize.inc"
454 DO n=1,igrnod(igr)%NENTITY
455 i=igrnod(igr)%ENTITY(n)
456 ar(1,i)= skew(1,isk)*askw(1,i)
457 . +skew(4,isk)*askw(2,i)
458 . +skew(7,isk)*askw(3,i)
459 ar(2,i)= skew(2,isk)*askw(1,i)
460 . +skew(5,isk)*askw(2,i)
461 . +skew(8,isk)*askw(3,i)
462 ar(3,i)= skew(3,isk)*askw(1,i)
463 . +skew(6,isk)*askw(2,i)
464 . +skew(9,isk)*askw(3,i)
465 damp(4,i)= skew(1,isk)*dampskw(1,i)
466 . +skew(4,isk)*dampskw(2,i)
467 . +skew(7,isk)*dampskw(3,i)
468 damp(5,i)= skew(2,isk)*dampskw(1,i)
469 . +skew(5,isk)*dampskw(2,i)
470 . +skew(8,isk)*dampskw(3,i)
471 damp(6,i)= skew(3,isk)*dampskw(1,i)
472 . +skew(6,isk)*dampskw(2,i)
473 . +skew(9,isk)*dampskw(3,i)
474 END DO
475C
476 ELSEIF ((fl_vrel == 1).AND.(id_rby == 0)) THEN
477C-------------Damping using average relative velocity - trans dof only ----
478C
479#include "vectorize.inc"
480 DO n=1,igrnod(igr)%NENTITY
481 i=igrnod(igr)%ENTITY(n)
482 vskw(1,i)= skew(1,isk)*v(1,i)+skew(2,isk)*v(2,i)+skew(3,isk)*v(3,i)
483 vskw(2,i)= skew(4,isk)*v(1,i)+skew(5,isk)*v(2,i)+skew(6,isk)*v(3,i)
484 vskw(3,i)= skew(7,isk)*v(1,i)+skew(8,isk)*v(2,i)+skew(9,isk)*v(3,i)
485 askw(1,i)= skew(1,isk)*a(1,i)+skew(2,isk)*a(2,i)+skew(3,isk)*a(3,i)
486 askw(2,i)= skew(4,isk)*a(1,i)+skew(5,isk)*a(2,i)+skew(6,isk)*a(3,i)
487 askw(3,i)= skew(7,isk)*a(1,i)+skew(8,isk)*a(2,i)+skew(9,isk)*a(3,i)
488 END DO
489C
490C---------------Computation of ref velocity in local skew
491C
492 v_refskw(1)= skew(1,isk)*v_ref(1,nd)+skew(2,isk)*v_ref(2,nd)
493 . +skew(3,isk)*v_ref(3,nd)
494 v_refskw(2)= skew(4,isk)*v_ref(1,nd)+skew(5,isk)*v_ref(2,nd)
495 . +skew(6,isk)*v_ref(3,nd)
496 v_refskw(3)= skew(7,isk)*v_ref(1,nd)+skew(8,isk)*v_ref(2,nd)
497 . +skew(9,isk)*v_ref(3,nd)
498 a_refskw(1)= skew(1,isk)*a_ref(1,nd)+skew(2,isk)*a_ref(2,nd)
499 . +skew(3,isk)*a_ref(3,nd)
500 a_refskw(2)= skew(4,isk)*a_ref(1,nd)+skew(5,isk)*a_ref(2,nd)
501 . +skew(6,isk)*a_ref(3,nd)
502 a_refskw(3)= skew(7,isk)*a_ref(1,nd)+skew(8,isk)*a_ref(2,nd)
503 . +skew(9,isk)*a_ref(3,nd)
504C
505C---------------alpha = Cdamp_M / dt--------------
506 call damping_vref_compute_dampa(nd,ndamp,nrdamp,dampr,dt1,tt,damp_a)
507!
508 omega_(1) = one/ (one + half * damp_a(1) * dt1)
509 omega_(2) = one/ (one + half * damp_a(2) * dt1)
510 omega_(3) = one/ (one + half * damp_a(3) * dt1)
511C-------------------------------------------------
512 DO n=1,igrnod(igr)%NENTITY
513 i=igrnod(igr)%ENTITY(n)
514 IF(tagslv_rby(i)==0) THEN
515 da_(1) = askw(1,i)-damp_a(1)*(vskw(1,i)-v_refskw(1))
516 . +half*damp_a(1)*dt1*a_refskw(1)
517 da_(1) = da_(1) * omega_(1) - askw(1,i)
518 askw(1,i) = askw(1,i) + da_(1)
519 dw_local(1)=da_(1)*(vskw(1,i)+half*askw(1,i)*dt1)
520C
521 da_(2) = askw(2,i)-damp_a(2)*(vskw(2,i)-v_refskw(2))
522 . +half*damp_a(2)*dt1*a_refskw(2)
523 da_(2) = da_(2) * omega_(2) - askw(2,i)
524 askw(2,i) = askw(2,i) + da_(2)
525 dw_local(2)=da_(2)*(vskw(2,i)+half*askw(2,i)*dt1)
526C
527 da_(3) = askw(3,i)-damp_a(3)*(vskw(3,i)-v_refskw(3))
528 . +half*damp_a(3)*dt1*a_refskw(3)
529 da_(3) = da_(3) * omega_(3) - askw(3,i)
530 askw(3,i) = askw(3,i) + da_(3)
531 dw_local(3)=da_(3)*(vskw(3,i)+half*askw(3,i)*dt1)
532C
533 dw =dw +ms(i)*(dw_local(1)+dw_local(2)+dw_local(3))*dt12*weight(i)
534 ENDIF
535 ENDDO
536C-------------------------------------------------
537#include "vectorize.inc"
538 DO n=1,igrnod(igr)%NENTITY
539 i=igrnod(igr)%ENTITY(n)
540 a(1,i)= skew(1,isk)*askw(1,i)+skew(4,isk)*askw(2,i)
541 . +skew(7,isk)*askw(3,i)
542 a(2,i)= skew(2,isk)*askw(1,i)+skew(5,isk)*askw(2,i)
543 . +skew(8,isk)*askw(3,i)
544 a(3,i)= skew(3,isk)*askw(1,i)+skew(6,isk)*askw(2,i)
545 . +skew(9,isk)*askw(3,i)
546 END DO
547C
548 ENDIF
549C
550 END IF ! isk<=1
551C
552 ENDIF ! tt>tstart...
553 ENDDO ! do nd=1,ndamp
554C
555#include "lockon.inc"
556 wfext = wfext + dw
557#include "lockoff.inc"
558 RETURN
559C-----------------------------------------------
560C AMS
561C-----------------------------------------------
562 1000 CONTINUE
563 dw = zero
564 DO nd=1,ndamp
565 igr = nint(dampr(2,nd))
566 isk = nint(dampr(15,nd))
567 factb = dampr(16,nd)
568 dampt = min(dt1,dt2)*factb
569 d_tstart = dampr(17,nd)
570 d_tstop = dampr(18,nd)
571 fl_freq_range = nint(dampr(31,nd))
572 IF ((tt>=d_tstart .AND. tt<=d_tstop).AND.(fl_freq_range==0)) THEN
573 IF(isk<=1)THEN
574C-------------------------------------------------
575 IF(iroddl/=0)THEN
576 dampa = dampr(9,nd)
577 dampb = dampr(10,nd)
578 betasdt= -min(dampb,dampt)*dt1/max(dt1*dt1,em30)
579 omega = one/ (one + half * dampa * dt1)
580
581 DO n=1,igrnod(igr)%NENTITY
582 i=igrnod(igr)%ENTITY(n)
583 IF(tagslv_rby(i)/=0) cycle
584C----- Damping seulement pour les noeuds en contact avec temporisation -----
585 IF (dampr(19,nd)/=0) THEN
586 IF (icontact(i)/=0) i_damp_rdof_tab(i) = dampr(19,nd)
587 IF (i_damp_rdof_tab(i)==0) GOTO 351
588 ENDIF
589C--------------------------------------------------------
590 da = ar(1,i) - dampa*vr(1,i)
591 . - betasdt *(ar(1,i)-damp(4,i))
592 da = da * omega - ar(1,i)
593 damp(4,i) = ar(1,i)
594 ar(1,i) = ar(1,i) + da
595 dw = dw+in(i)*da*(vr(1,i)+half*ar(1,i)*dt1)*dt12*weight(i)
596351 CONTINUE
597 ENDDO
598 dampa = dampr(11,nd)
599 dampb = dampr(12,nd)
600 betasdt= -min(dampb,dampt)*dt1/max(dt1*dt1,em30)
601 omega = one/ (one + half * dampa * dt1)
602
603 DO n=1,igrnod(igr)%NENTITY
604 i=igrnod(igr)%ENTITY(n)
605 IF(tagslv_rby(i)/=0) cycle
606C----- Damping seulement pour les noeuds en contact avec temporisation -----
607 IF (dampr(19,nd)/=0) THEN
608 IF (i_damp_rdof_tab(i)==0) GOTO 451
609 ENDIF
610C--------------------------------------------------------
611 da = ar(2,i) - dampa*vr(2,i)
612 . - betasdt *(ar(2,i)-damp(5,i))
613 da = da * omega - ar(2,i)
614 damp(5,i) = ar(2,i)
615 ar(2,i) = ar(2,i) + da
616 dw = dw+in(i)*da*(vr(2,i)+half*ar(2,i)*dt1)*dt12*weight(i)
617451 CONTINUE
618 ENDDO
619 dampa = dampr(13,nd)
620 dampb = dampr(14,nd)
621 betasdt= -min(dampb,dampt)*dt1/max(dt1*dt1,em30)
622 omega = one/ (one + half * dampa * dt1)
623
624 DO n=1,igrnod(igr)%NENTITY
625 i=igrnod(igr)%ENTITY(n)
626 IF(tagslv_rby(i)/=0) cycle
627C----- Damping seulement pour les noeuds en contact avec temporisation -----
628 IF (dampr(19,nd)/=0) THEN
629 IF (i_damp_rdof_tab(i)==0) GOTO 551
630 i_damp_rdof_tab(i)=i_damp_rdof_tab(i)-1
631 ENDIF
632C--------------------------------------------------------
633 da = ar(3,i) - dampa*vr(3,i)
634 . - betasdt *(ar(3,i)-damp(6,i))
635 da = da * omega - ar(3,i)
636 damp(6,i) = ar(3,i)
637 ar(3,i) = ar(3,i) + da
638 dw = dw+in(i)*da*(vr(3,i)+half*ar(3,i)*dt1)*dt12*weight(i)
639551 CONTINUE
640 ENDDO
641 END IF
642 ELSE
643 IF(iroddl/=0)THEN
644#include "vectorize.inc"
645 DO n=1,igrnod(igr)%NENTITY
646 i=igrnod(igr)%ENTITY(n)
647 vskw(1,i)= skew(1,isk)*vr(1,i)
648 . +skew(2,isk)*vr(2,i)
649 . +skew(3,isk)*vr(3,i)
650 vskw(2,i)= skew(4,isk)*vr(1,i)
651 . +skew(5,isk)*vr(2,i)
652 . +skew(6,isk)*vr(3,i)
653 vskw(3,i)= skew(7,isk)*vr(1,i)
654 . +skew(8,isk)*vr(2,i)
655 . +skew(9,isk)*vr(3,i)
656 askw(1,i)= skew(1,isk)*ar(1,i)
657 . +skew(2,isk)*ar(2,i)
658 . +skew(3,isk)*ar(3,i)
659 askw(2,i)= skew(4,isk)*ar(1,i)
660 . +skew(5,isk)*ar(2,i)
661 . +skew(6,isk)*ar(3,i)
662 askw(3,i)= skew(7,isk)*ar(1,i)
663 . +skew(8,isk)*ar(2,i)
664 . +skew(9,isk)*ar(3,i)
665 dampskw(1,i)= skew(1,isk)*damp(4,i)
666 . +skew(2,isk)*damp(5,i)
667 . +skew(3,isk)*damp(6,i)
668 dampskw(2,i)= skew(4,isk)*damp(4,i)
669 . +skew(5,isk)*damp(5,i)
670 . +skew(6,isk)*damp(6,i)
671 dampskw(3,i)= skew(7,isk)*damp(4,i)
672 . +skew(8,isk)*damp(5,i)
673 . +skew(9,isk)*damp(6,i)
674 END DO
675 dampa = dampr(9,nd)
676 dampb = dampr(10,nd)
677 betasdt= -min(dampb,dampt)*dt1/max(dt1*dt1,em30)
678 omega = one/ (one + half * dampa * dt1)
679
680 DO n=1,igrnod(igr)%NENTITY
681 i=igrnod(igr)%ENTITY(n)
682 IF(tagslv_rby(i)/=0) cycle
683C----- Damping seulement pour les noeuds en contact avec temporisation -----
684 IF (dampr(19,nd)/=0) THEN
685 IF (icontact(i)/=0) i_damp_rdof_tab(i) = dampr(19,nd)
686 IF (i_damp_rdof_tab(i)==0) GOTO 751
687 ENDIF
688C--------------------------------------------------------
689 da = askw(1,i) - dampa*vskw(1,i)
690 . - betasdt *(askw(1,i) - dampskw(1,i))
691 da = da * omega - askw(1,i)
692 dampskw(1,i) = askw(1,i)
693 askw(1,i) = askw(1,i) + da
694 dw =dw
695 . +in(i)*da*(vskw(1,i)+half*askw(1,i)*dt1)*dt12*weight(i)
696751 CONTINUE
697 ENDDO
698 dampa = dampr(11,nd)
699 dampb = dampr(12,nd)
700 betasdt= -min(dampb,dampt)*dt1/max(dt1*dt1,em30)
701 omega = one/ (one + half * dampa * dt1)
702
703 DO n=1,igrnod(igr)%NENTITY
704 i=igrnod(igr)%ENTITY(n)
705 IF(tagslv_rby(i)/=0) cycle
706C----- Damping seulement pour les noeuds en contact avec temporisation -----
707 IF (dampr(19,nd)/=0) THEN
708 IF (i_damp_rdof_tab(i)==0) GOTO 851
709 ENDIF
710C--------------------------------------------------------
711 da = askw(2,i) - dampa*vskw(2,i)
712 . - betasdt *(askw(2,i) - dampskw(2,i))
713 da = da * omega - askw(2,i)
714 dampskw(2,i) = askw(2,i)
715 askw(2,i) = askw(2,i) + da
716 dw =dw
717 . +in(i)*da*(vskw(2,i)+half*askw(2,i)*dt1)*dt12*weight(i)
718851 CONTINUE
719 ENDDO
720 dampa = dampr(13,nd)
721 dampb = dampr(14,nd)
722 betasdt= -min(dampb,dampt)*dt1/max(dt1*dt1,em30)
723 omega = one/ (one + half * dampa * dt1)
724
725 DO n=1,igrnod(igr)%NENTITY
726 i=igrnod(igr)%ENTITY(n)
727 IF(tagslv_rby(i)/=0) cycle
728C----- Damping seulement pour les noeuds en contact avec temporisation -----
729 IF (dampr(19,nd)/=0) THEN
730 IF (i_damp_rdof_tab(i)==0) GOTO 951
731 i_damp_rdof_tab(i)=i_damp_rdof_tab(i)-1
732 ENDIF
733C--------------------------------------------------------
734 da = askw(3,i) - dampa*vskw(3,i)
735 . - betasdt *(askw(3,i) - dampskw(3,i))
736 da = da * omega - askw(3,i)
737 dampskw(3,i) = askw(3,i)
738 askw(3,i) = askw(3,i) + da
739 dw =dw
740 . +in(i)*da*(vskw(3,i)+half*askw(3,i)*dt1)*dt12*weight(i)
741951 CONTINUE
742 ENDDO
743#include "vectorize.inc"
744 DO n=1,igrnod(igr)%NENTITY
745 i=igrnod(igr)%ENTITY(n)
746 IF(tagslv_rby(i)/=0) cycle
747 ar(1,i)= skew(1,isk)*askw(1,i)
748 . +skew(4,isk)*askw(2,i)
749 . +skew(7,isk)*askw(3,i)
750 ar(2,i)= skew(2,isk)*askw(1,i)
751 . +skew(5,isk)*askw(2,i)
752 . +skew(8,isk)*askw(3,i)
753 ar(3,i)= skew(3,isk)*askw(1,i)
754 . +skew(6,isk)*askw(2,i)
755 . +skew(9,isk)*askw(3,i)
756 damp(4,i)= skew(1,isk)*dampskw(1,i)
757 . +skew(4,isk)*dampskw(2,i)
758 . +skew(7,isk)*dampskw(3,i)
759 damp(5,i)= skew(2,isk)*dampskw(1,i)
760 . +skew(5,isk)*dampskw(2,i)
761 . +skew(8,isk)*dampskw(3,i)
762 damp(6,i)= skew(3,isk)*dampskw(1,i)
763 . +skew(6,isk)*dampskw(2,i)
764 . +skew(9,isk)*dampskw(3,i)
765 END DO
766 END IF
767 END IF
768 ENDIF
769 ENDDO
770C
771#include "lockon.inc"
772 wfext = wfext + dw
773#include "lockoff.inc"
774 DEALLOCATE(vskw)
775 DEALLOCATE(askw)
776 DEALLOCATE(dampskw)
777
778C
779 RETURN
780 END
781C
782!||====================================================================
783!|| damping44 ../engine/source/assembly/damping.F
784!||--- called by ------------------------------------------------------
785!|| resol ../engine/source/engine/resol.F
786!||--- uses -----------------------------------------------------
787!|| groupdef_mod ../common_source/modules/groupdef_mod.F
788!||====================================================================
789 SUBROUTINE damping44(
790 1 DIM ,V ,
791 2 VR ,A ,AR ,MS ,IN ,
792 3 DAMPR ,DAMP ,IGRNOD ,WEIGHT ,TAGSLV_RBY,
793 4 WFEXT)
794C-----------------------------------------------
795C M o d u l e s
796C-----------------------------------------------
797 USE groupdef_mod
798C-----------------------------------------------
799C I m p l i c i t T y p e s
800C-----------------------------------------------
801#include "implicit_f.inc"
802#include "comlock.inc"
803C-----------------------------------------------
804C C o m m o n B l o c k s
805C-----------------------------------------------
806#include "com01_c.inc"
807#include "com04_c.inc"
808#include "com06_c.inc"
809#include "com08_c.inc"
810C-----------------------------------------------
811C D u m m y A r g u m e n t s
812C-----------------------------------------------
813 INTEGER DIM,ITASK,WEIGHT(*),TAGSLV_RBY(*)
814 my_real V(3,*), VR(3,*), A(3,*), AR(3,*) ,MS(*), IN(*),DAMPR(4,*), DAMP(DIM,*)
815 DOUBLE PRECISION, INTENT(INOUT) :: WFEXT
816C-----------------------------------------------
817 TYPE (GROUP_) , DIMENSION(NGRNOD) :: IGRNOD
818C-----------------------------------------------
819C L o c a l V a r i a b l e s
820C-----------------------------------------------
821 INTEGER I,J,N,ND,IGR
822 my_real betasdt,aa,da,dw,omega
823C-----------------------------------------------
824C C = a M + b K
825C======================================================================|
826 dw = zero
827 DO nd=1,ndamp
828 igr = nint(dampr(2,nd))
829 dampa = dampr(3,nd)
830 dampb = dampr(4,nd)
831 betasdt= -min(dampb,dt1,dt2)*dt1/max(dt1*dt1,em30)
832 omega = one/ (one + half * dampa * dt1)
833 DO j=1,3
834 IF(iroddl==0)THEN
835 DO n=1,igrnod(igr)%NENTITY
836 i=igrnod(igr)%ENTITY(n)
837 IF(tagslv_rby(i)/=0) cycle
838 da = a(j,i) - dampa*v(j,i) - betasdt *(a(j,i) - damp(j,i))
839 da = da * omega - a(j,i)
840 damp(j,i) = a(j,i)
841 a(j,i) = a(j,i) + da
842 dw =dw+ms(i)*da*(v(j,i)+half*a(j,i)*dt1)*dt12*weight(i)
843 ENDDO
844 ELSE
845 DO n=1,igrnod(igr)%NENTITY
846 i=igrnod(igr)%ENTITY(n)
847 IF(tagslv_rby(i)/=0) cycle
848 da = a(j,i) - dampa*v(j,i) - betasdt *(a(j,i) - damp(j,i))
849 da = da * omega - a(j,i)
850 damp(j,i) = a(j,i)
851 a(j,i) = a(j,i) + da
852 dw =dw+ms(i)*da*(v(j,i)+half*a(j,i)*dt1)*dt12*weight(i)
853 da = ar(j,i) - dampa*vr(j,i)
854 . - betasdt *(ar(j,i)-damp(j+3,i))
855 da = da * omega - ar(j,i)
856 damp(j+3,i) = ar(j,i)
857 ar(j,i) = ar(j,i) + da
858 dw = dw+in(i)*da*(vr(j,i)+half*ar(j,i)*dt1)*dt12*weight(i)
859 ENDDO
860 ENDIF
861 ENDDO
862 ENDDO
863C
864#include "lockon.inc"
865 wfext = wfext + dw
866#include "lockoff.inc"
867C
868 RETURN
869 END
870C--------------------------
871C DAMPING alpha M + beta K
872C--------------------------
873!||====================================================================
874!|| damping ../engine/source/assembly/damping.F
875!||--- called by ------------------------------------------------------
876!|| resol ../engine/source/engine/resol.F
877!||--- uses -----------------------------------------------------
878!|| groupdef_mod ../common_source/modules/groupdef_mod.F
879!||====================================================================
880 SUBROUTINE damping(NODFT,NODLT,V ,VR,A ,AR ,DAMP,MS,IN,
881 . IGRNOD,DIM,ITASK,WEIGHT,TAGSLV_RBY, WFEXT)
882C-----------------------------------------------
883C M o d u l e s
884C-----------------------------------------------
885 USE groupdef_mod
886C-----------------------------------------------
887C I m p l i c i t T y p e s
888C-----------------------------------------------
889#include "implicit_f.inc"
890#include "comlock.inc"
891C-----------------------------------------------
892C C o m m o n B l o c k s
893C-----------------------------------------------
894#include "com01_c.inc"
895#include "com04_c.inc"
896#include "com06_c.inc"
897#include "com08_c.inc"
898C-----------------------------------------------
899C D u m m y A r g u m e n t s
900C-----------------------------------------------
901 integer
902 . nodft,nodlt,dim,itask,
903 . weight(*),tagslv_rby(*)
904 my_real v(3,*) ,vr(3,*),a(3,*) ,ar(3,*) ,damp(dim,*),ms(*),in(*)
905 DOUBLE PRECISION,INTENT(INOUT) :: WFEXT
906C-----------------------------------------------
907 TYPE (GROUP_) , DIMENSION(NGRNOD) :: IGRNOD
908C-----------------------------------------------
909C L o c a l V a r i a b l e s
910C-----------------------------------------------
911 INTEGER I,J,N
912 my_real BETASDT,AA,DA,DW,OMEGA
913C-----------------------------------------------
914C
915C C = a M + b K
916C
917C
918C Fint(1) = Fint0(0) + K v(1/2) dt1 + C v(1)
919C Fint0(1) = Fint0(0) + K v(1/2) dt1
920C Fint(1) = Fint0(1) + C v(1)
921C Fint(1) = Fint0(1) + [a M + b K] v(1)
922C
923C a0(1) = a0(0) - K/M v(1/2) dt1
924C a(1) = a0(1) - [a + b K/M] v(1)
925C
926C a(1) = a0(1) - a v(1) - b/dt1 K/M v(1) dt1
927C a(1) = a0(1) - a (v(1/2) + a(1) dt1/2) - b K/M (v(1/2) + a(1) dt1/2)
928C a(1) = a0(1) - a v(1/2) - a a(1) dt1/2 - b/dt1 [a0(1) - a0(0)]
929C - b K/M a(1) dt1/2
930C
931C a(1) ~= a0(1) - a v(1/2) - a a(1) dt1/2 - b/dt1 [a0(1) - a0(0)]
932C a(1)[1+a dt1/2] ~= a0(1) - a v(1/2) - b/dt1 [a0(1) - a0(0)]
933C
934C a(1) = [a0(1) - a v(1/2) - b/dt1 [a0(1) - a0(0)]]/[1+a dt1/2]
935C da = a(1)-a0(1)
936C da = [a0(1) - a v(1/2) - b/dt1 [a0(1) - a0(0)]]/[1+a dt1/2] - a0(1)
937C-----------------------------------------------
938C
939C a(1) ~= [1-a dt1/2- b/dt1]a0(1) - a v(1/2) + b/dt1 a0(0)
940C
941C-----------------------------------------------
942
943 dw = zero
944C
945 IF(dampb>=0.0)THEN
946 betasdt= -min(dampb,dt1,dt2)*dt1/max(dt1*dt1,em30)
947 ELSE
948 betasdt= dampb
949 ENDIF
950 omega = one/ (one + half * dampa * dt1)
951C
952 IF(idampg==0)THEN
953 DO j=1,3
954 IF(iroddl==0)THEN
955 DO i=nodft,nodlt
956 IF(tagslv_rby(i)/=0) cycle
957 da = a(j,i) - dampa*v(j,i) - betasdt *(a(j,i) - damp(j,i))
958 da = da * omega - a(j,i)
959 damp(j,i) = a(j,i)
960 a(j,i) = a(j,i) + da
961 dw=dw+weight(i)*(ms(i)*da*(v(j,i)+half*a(j,i)*dt1)*dt12)
962 ENDDO
963 ELSE
964 DO i=nodft,nodlt
965 IF(tagslv_rby(i)/=0) cycle
966 da = a(j,i) - dampa*v(j,i) - betasdt *(a(j,i) - damp(j,i))
967 da = da * omega - a(j,i)
968 damp(j,i) = a(j,i)
969 a(j,i) = a(j,i) + da
970 dw = dw+ms(i)*da*(v(j,i)+half*a(j,i)*dt1)*dt12*weight(i)
971 da = ar(j,i) - dampa*vr(j,i) - betasdt *(ar(j,i)-damp(j+3,i))
972 da = da * omega - ar(j,i)
973 damp(j+3,i) = ar(j,i)
974 ar(j,i) = ar(j,i) + da
975 dw =dw+weight(i)*in(i)*da*(vr(j,i)+half*ar(j,i)*dt1)*dt12
976 ENDDO
977 ENDIF
978 ENDDO
979 ELSEIF(itask==0)THEN
980 DO j=1,3
981 IF(iroddl==0)THEN
982 DO n=1,igrnod(idampg)%NENTITY
983 i=igrnod(idampg)%ENTITY(n)
984 IF(tagslv_rby(i)/=0) cycle
985 da = a(j,i) - dampa*v(j,i) - betasdt *(a(j,i) - damp(j,i))
986 da = da * omega - a(j,i)
987 damp(j,i) = a(j,i)
988 a(j,i) = a(j,i) + da
989 dw =dw+weight(i)*ms(i)*da*(v(j,i)+half*a(j,i)*dt1)*dt12
990 ENDDO
991 ELSE
992 DO n=1,igrnod(idampg)%NENTITY
993 i=igrnod(idampg)%ENTITY(n)
994 IF(tagslv_rby(i)/=0) cycle
995 da = a(j,i) - dampa*v(j,i) - betasdt *(a(j,i) - damp(j,i))
996 da = da * omega - a(j,i)
997 damp(j,i) = a(j,i)
998 a(j,i) = a(j,i) + da
999 dw = dw+ms(i)*da*(v(j,i)+half*a(j,i)*dt1)*dt12*weight(i)
1000 da = ar(j,i) - dampa*vr(j,i) - betasdt *(ar(j,i)-damp(j+3,i))
1001 da = da * omega - ar(j,i)
1002 damp(j+3,i) = ar(j,i)
1003 ar(j,i) = ar(j,i) + da
1004 dw =dw+weight(i)*in(i)*da*(vr(j,i)+half*ar(j,i)*dt1)*dt12
1005 ENDDO
1006 ENDIF
1007 ENDDO
1008 ENDIF
1009C
1010#include "lockon.inc"
1011c EDAMP = EDAMP + DW
1012 wfext = wfext + dw
1013#include "lockoff.inc"
1014C
1015 RETURN
1016 END
my_real
#define my_real
Definition cppsort.cpp:32
damping
subroutine damping(nodft, nodlt, v, vr, a, ar, damp, ms, in, igrnod, dim, itask, weight, tagslv_rby, wfext)
Definition damping.F:882
damping44
subroutine damping44(dim, v, vr, a, ar, ms, in, dampr, damp, igrnod, weight, tagslv_rby, wfext)
Definition damping.F:794
damping51
subroutine damping51(dim, v, vr, a, ar, ms, in, dampr, damp, igrnod, weight, tagslv_rby, skew, icontact, i_damp_rdof_tab, ndamp_vrel, id_damp_vrel, fr_damp_vrel, iparit, ispmd, wfext)
Definition damping.F:41
damping_vref
subroutine damping_vref(igrnod, v, a, ms, dampr, weight, v_ref, a_ref, ndamp_vrel, id_damp_vrel, fr_damp_vrel, ngrnod, numnod, ndamp, nrdamp, nspmd, iparit, ispmd)
Definition damping_vref.F:38
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
groupdef_mod
Definition groupdef_mod.F:662