OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
spdens.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!|| spdens ../engine/source/elements/sph/spdens.F
25!||--- called by ------------------------------------------------------
26!|| forintp ../engine/source/elements/forintp.F
27!||--- calls -----------------------------------------------------
28!|| weight1 ../engine/source/elements/sph/weight.F
29!||--- uses -----------------------------------------------------
30!|| sphbox ../engine/share/modules/sphbox.F
31!||====================================================================
32 SUBROUTINE spdens(
33 1 X ,V ,MS ,SPBUF ,ITAB ,
34 2 KXSP ,IXSP ,NOD2SP ,ISPSYM ,XSPSYM ,
35 3 VSPSYM ,IPARG ,WA ,WACOMP )
36C-----------------------------------------------
37C M o d u l e s
38C-----------------------------------------------
39 USE sphbox
40C-----------------------------------------------
41C I m p l i c i t T y p e s
42C-----------------------------------------------
43#include "implicit_f.inc"
44C-----------------------------------------------
45C C o m m o n B l o c k s
46C-----------------------------------------------
47#include "vect01_c.inc"
48#include "com08_c.inc"
49#include "sphcom.inc"
50#include "param_c.inc"
51C-----------------------------------------------
52C D u m m y A r g u m e n t s
53C-----------------------------------------------
54 INTEGER KXSP(NISP,*),IXSP(KVOISPH,*),NOD2SP(*),ITAB(*),
55 . ISPSYM(NSPCOND,*),IPARG(NPARG,*)
56 my_real
57 . x(3,*) ,v(3,*) ,ms(*) ,
58 . spbuf(nspbuf,*) ,xspsym(3,*) ,vspsym(3,*) ,
59 . wa(kwasph,*) ,wacomp(16,*)
60C-----------------------------------------------
61C L o c a l V a r i a b l e s
62C-----------------------------------------------
63 INTEGER I,N,INOD,JNOD,J,NVOIS,M,
64 . NVOISS,SM,JS,NC,NS,NN
65 my_real
66 . xi,yi,zi,di,rhoi,xj,yj,zj,dj,rhoj,dij,
67 . vxi,vyi,vzi,vxj,vyj,vzj,
68 . vj,vjx,vjy,vjz,
69 . drho,wght,wgrad(3),divv,
70 . muij,mumax,
71 . fact,wun,wvx,wvy,wvz,wrho,
72 . dxx,dxy,dxz,dyx,dyy,dyz,dzx,dzy,dzz,dt1d2,
73 . exx,exy,exz,eyx,eyy,eyz,ezx,ezy,ezz,
74 . rotvx,rotvy,rotvz,rotv,
75 . alphai,alphaxi,alphayi,alphazi,
76 . betaxi,betayi,betazi,
77 . betaxxi,betayxi,betazxi,
78 . betaxyi,betayyi,betazyi,
79 . betaxzi,betayzi,betazzi,
80 . betax,wgrdx,wgrdy,wgrdz
81C-------------------------------------------
82 DO i=lft,llt
83 n =nft+i
84 wa(1,n)=zero
85 wa(2,n)=zero
86 wa(3,n)=zero
87 wa(4,n)=zero
88 wa(5,n)=zero
89 wa(6,n)=zero
90 wa(7,n)=zero
91 wa(8,n)=zero
92 wa(9,n)=zero
93C characteristic length = diameter.
94 wa(11,n)=spbuf(1,n)
95 wa(12,n)=zero
96 ENDDO
97C-----------------------------------------------
98C Calcul des densites et tenseurs de deformations.
99C-----------------------------------------------
100 dt1d2=0.5*dt1
101 DO 10 i=lft,llt
102 n =nft+i
103 IF(kxsp(2,n)<=0)GOTO 10
104 inod =kxsp(3,n)
105 xi=x(1,inod)
106 yi=x(2,inod)
107 zi=x(3,inod)
108 di=spbuf(1,n)
109 vxi=v(1,inod)
110 vyi=v(2,inod)
111 vzi=v(3,inod)
112 rhoi =wa(10,n)
113C-----
114C for DT,SPCEL stability computation:
115 mumax=zero
116C
117 rotvx=zero
118 rotvy=zero
119 rotvz=zero
120C------
121 alphai=wacomp(1,n)
122C BETAXI=WACOMP(2,N)
123C BETAYI=WACOMP(3,N)
124C BETAZI=WACOMP(4,N)
125 alphaxi=wacomp( 5,n)
126 alphayi=wacomp( 6,n)
127 alphazi=wacomp( 7,n)
128 betaxxi=wacomp( 8,n)
129 betayxi=wacomp( 9,n)
130 betazxi=wacomp(10,n)
131 betaxyi=wacomp(11,n)
132 betayyi=wacomp(12,n)
133 betazyi=wacomp(13,n)
134 betaxzi=wacomp(14,n)
135 betayzi=wacomp(15,n)
136 betazzi=wacomp(16,n)
137C------
138 nvois=kxsp(4,n)
139 DO j=1,nvois
140 jnod=ixsp(j,n)
141 IF(jnod>0)THEN
142 m=nod2sp(jnod)
143 xj=x(1,jnod)
144 yj=x(2,jnod)
145 zj=x(3,jnod)
146 dj=spbuf(1,m)
147 dij=0.5*(di+dj)
148 CALL weight1(xi,yi,zi,xj,yj,zj,dij,wght,wgrad)
149 vxj=v(1,jnod)
150 vyj=v(2,jnod)
151 vzj=v(3,jnod)
152 rhoj=wa(10,m)
153 vj=spbuf(12,m)/max(em20,rhoj)
154 ELSE
155 nn = -jnod
156 xj=xsphr(3,nn)
157 yj=xsphr(4,nn)
158 zj=xsphr(5,nn)
159 dj=xsphr(2,nn)
160 dij=0.5*(di+dj)
161 CALL weight1(xi,yi,zi,xj,yj,zj,dij,wght,wgrad)
162 vxj=xsphr(9,nn)
163 vyj=xsphr(10,nn)
164 vzj=xsphr(11,nn)
165 rhoj=xsphr(7,nn)
166 vj=xsphr(8,nn)/max(em20,rhoj)
167 END IF
168 wgrdx=wgrad(1)*alphai+wght*alphaxi
169 . +wgrad(1)*betaxxi+wgrad(2)*betaxyi+wgrad(3)*betaxzi
170 wgrdy=wgrad(2)*alphai+wght*alphayi
171 . +wgrad(1)*betayxi+wgrad(2)*betayyi+wgrad(3)*betayzi
172 wgrdz=wgrad(3)*alphai+wght*alphazi
173 . +wgrad(1)*betazxi+wgrad(2)*betazyi+wgrad(3)*betazzi
174! Old order1 correction
175! BETAX=ONE +BETAXI*(XI-XJ)+BETAYI*(YI-YJ)+BETAZI*(ZI-ZJ)
176! WGRDX=WGRAD(1)*ALPHAI*BETAX
177! . +WGHT*(ALPHAXI*BETAX+ALPHAI*
178! . (BETAXXI*(XI-XJ)+BETAYXI*(YI-YJ)+BETAZXI*(ZI-ZJ)+BETAXI))
179! WGRDY=WGRAD(2)*ALPHAI*BETAX
180! . +WGHT*(ALPHAYI*BETAX+ALPHAI*
181! . (BETAXYI*(XI-XJ)+BETAYYI*(YI-YJ)+BETAZYI*(ZI-ZJ)+BETAYI))
182! WGRDZ=WGRAD(3)*ALPHAI*BETAX
183! . +WGHT*(ALPHAZI*BETAX+ALPHAI*
184! . (BETAXZI*(XI-XJ)+BETAYZI*(YI-YJ)+BETAZZI*(ZI-ZJ)+BETAZI))
185 wgrad(1)=wgrdx
186 wgrad(2)=wgrdy
187 wgrad(3)=wgrdz
188 vjx=vj*(vxi-vxj)
189 vjy=vj*(vyi-vyj)
190 vjz=vj*(vzi-vzj)
191C
192C stores DXX,DYY,DZZ,DXY,DYZ,DXZ,DYX,DZY,DZX
193C 2nd order correction :
194 dxx=-vjx*wgrad(1)
195 dyy=-vjy*wgrad(2)
196 dzz=-vjz*wgrad(3)
197 dxy=-vjx*wgrad(2)
198 dyz=-vjy*wgrad(3)
199 dxz=-vjx*wgrad(3)
200 dyx=-vjy*wgrad(1)
201 dzy=-vjz*wgrad(2)
202 dzx=-vjz*wgrad(1)
203C--------
204 exy = dxy
205 . -dt1d2*(dxx*dxy+dyx*dyy+dzx*dzy)
206 eyz = dyz
207 . -dt1d2*(dyy*dyz+dzy*dzz+dxy*dxz)
208 exz = dxz
209 . -dt1d2*(dzz*dzx+dxz*dxx+dyz*dyx)
210 eyx = dyx
211 . -dt1d2*(dxx*dxy+dyx*dyy+dzx*dzy)
212 ezy = dzy
213 . -dt1d2*(dyy*dyz+dzy*dzz+dxy*dxz)
214 ezx = dzx
215 . -dt1d2*(dzz*dzx+dxz*dxx+dyz*dyx)
216 exx = dxx
217 . -dt1d2*(dxx*dxx+dyx*dyx+dzx*dzx)
218 eyy = dyy
219 . -dt1d2*(dyy*dyy+dzy*dzy+dxy*dxy)
220 ezz = dzz
221 . -dt1d2*(dzz*dzz+dxz*dxz+dyz*dyz)
222 wa(1,n)=wa(1,n)+exx
223 wa(2,n)=wa(2,n)+eyy
224 wa(3,n)=wa(3,n)+ezz
225c if(n==315)print*,'vj,VYJ,=',n,vj,VYJ
226 wa(4,n)=wa(4,n)+exy
227 wa(5,n)=wa(5,n)+eyz
228 wa(6,n)=wa(6,n)+exz
229 wa(7,n)=wa(7,n)+eyx
230 wa(8,n)=wa(8,n)+ezy
231 wa(9,n)=wa(9,n)+ezx
232C--------
233C for stability computation into material routines:
234 muij=(vxi-vxj)*(xi-xj)
235 . +(vyi-vyj)*(yi-yj)
236 . +(vzi-vzj)*(zi-zj)
237 muij=dij*muij/
238 . ((xi-xj)*(xi-xj)
239 . +(yi-yj)*(yi-yj)
240C . +(ZI-ZJ)*(ZI-ZJ))
241 . +(zi-zj)*(zi-zj)+em02*dij*dij)
242 mumax=max(mumax,-muij)
243C--------
244C for artificial viscosity computation :
245 rotvx=rotvx+vjy*wgrad(3)-vjz*wgrad(2)
246 rotvy=rotvy+vjz*wgrad(1)-vjx*wgrad(3)
247 rotvz=rotvz+vjx*wgrad(2)-vjy*wgrad(1)
248 END DO
249C------
250C partie symetrique.
251 nvoiss=kxsp(6,n)
252 DO j=kxsp(5,n)+1,kxsp(5,n)+nvoiss
253 js=ixsp(j,n)
254 IF(js>0)THEN
255 sm=js/(nspcond+1)
256 nc=mod(js,nspcond+1)
257 js=ispsym(nc,sm)
258 xj =xspsym(1,js)
259 yj =xspsym(2,js)
260 zj =xspsym(3,js)
261 vxj=vspsym(1,js)
262 vyj=vspsym(2,js)
263 vzj=vspsym(3,js)
264 dj =spbuf(1,sm)
265 dij =half*(di+dj)
266 rhoj=wa(10,sm)
267 CALL weight1(xi,yi,zi,xj,yj,zj,dij,wght,wgrad)
268 vj=spbuf(12,sm)/max(em20,rhoj)
269 ELSE
270 sm=-js/(nspcond+1)
271 nc=mod(-js,nspcond+1)
272 js=ispsymr(nc,sm)
273 xj =xspsym(1,js)
274 yj =xspsym(2,js)
275 zj =xspsym(3,js)
276 vxj=vspsym(1,js)
277 vyj=vspsym(2,js)
278 vzj=vspsym(3,js)
279 dj =xsphr(2,sm)
280 dij =half*(di+dj)
281 rhoj=xsphr(7,sm)
282 CALL weight1(xi,yi,zi,xj,yj,zj,dij,wght,wgrad)
283 vj=xsphr(8,sm)/max(em20,rhoj)
284 END IF
285C
286 wgrdx=wgrad(1)*alphai+wght*alphaxi
287 . +wgrad(1)*betaxxi+wgrad(2)*betaxyi+wgrad(3)*betaxzi
288 wgrdy=wgrad(2)*alphai+wght*alphayi
289 . +wgrad(1)*betayxi+wgrad(2)*betayyi+wgrad(3)*betayzi
290 wgrdz=wgrad(3)*alphai+wght*alphazi
291 . +wgrad(1)*betazxi+wgrad(2)*betazyi+wgrad(3)*betazzi
292! Old order1 correction
293! BETAX=ONE + BETAXI*(XI-XJ)+BETAYI*(YI-YJ)+BETAZI*(ZI-ZJ)
294! WGRDX=WGRAD(1)*ALPHAI*BETAX
295! . +WGHT*(ALPHAXI*BETAX+ALPHAI*
296! . (BETAXXI*(XI-XJ)+BETAYXI*(YI-YJ)+BETAZXI*(ZI-ZJ)+BETAXI))
297! WGRDY=WGRAD(2)*ALPHAI*BETAX
298! . +WGHT*(ALPHAYI*BETAX+ALPHAI*
299! . (BETAXYI*(XI-XJ)+BETAYYI*(YI-YJ)+BETAZYI*(ZI-ZJ)+BETAYI))
300! WGRDZ=WGRAD(3)*ALPHAI*BETAX
301! . +WGHT*(ALPHAZI*BETAX+ALPHAI*
302! . (BETAXZI*(XI-XJ)+BETAYZI*(YI-YJ)+BETAZZI*(ZI-ZJ)+BETAZI))
303 wgrad(1)=wgrdx
304 wgrad(2)=wgrdy
305 wgrad(3)=wgrdz
306 vjx=vj*(vxi-vxj)
307 vjy=vj*(vyi-vyj)
308 vjz=vj*(vzi-vzj)
309C
310C stores DXX,DYY,DZZ,DXY,DYZ,DXZ,DYX,DZY,DZX
311C 2nd order correction :
312 dxx=-vjx*wgrad(1)
313 dyy=-vjy*wgrad(2)
314 dzz=-vjz*wgrad(3)
315 dxy=-vjx*wgrad(2)
316 dyz=-vjy*wgrad(3)
317 dxz=-vjx*wgrad(3)
318 dyx=-vjy*wgrad(1)
319 dzy=-vjz*wgrad(2)
320 dzx=-vjz*wgrad(1)
321C--------
322 exy = dxy
323 . -dt1d2*(dxx*dxy+dyx*dyy+dzx*dzy)
324 eyz = dyz
325 . -dt1d2*(dyy*dyz+dzy*dzz+dxy*dxz)
326 exz = dxz
327 . -dt1d2*(dzz*dzx+dxz*dxx+dyz*dyx)
328 eyx = dyx
329 . -dt1d2*(dxx*dxy+dyx*dyy+dzx*dzy)
330 ezy = dzy
331 . -dt1d2*(dyy*dyz+dzy*dzz+dxy*dxz)
332 ezx = dzx
333 . -dt1d2*(dzz*dzx+dxz*dxx+dyz*dyx)
334 exx = dxx
335 . -dt1d2*(dxx*dxx+dyx*dyx+dzx*dzx)
336 eyy = dyy
337 . -dt1d2*(dyy*dyy+dzy*dzy+dxy*dxy)
338 ezz = dzz
339 . -dt1d2*(dzz*dzz+dxz*dxz+dyz*dyz)
340 wa(1,n)=wa(1,n)+exx
341 wa(2,n)=wa(2,n)+eyy
342 wa(3,n)=wa(3,n)+ezz
343 wa(4,n)=wa(4,n)+exy
344 wa(5,n)=wa(5,n)+eyz
345 wa(6,n)=wa(6,n)+exz
346 wa(7,n)=wa(7,n)+eyx
347 wa(8,n)=wa(8,n)+ezy
348 wa(9,n)=wa(9,n)+ezx
349C-------
350C for stability computation into material routines:
351 muij=(vxi-vxj)*(xi-xj)
352 . +(vyi-vyj)*(yi-yj)
353 . +(vzi-vzj)*(zi-zj)
354 muij=dij*muij/
355 . ((xi-xj)*(xi-xj)
356 . +(yi-yj)*(yi-yj)
357C . +(ZI-ZJ)*(ZI-ZJ))
358 . +(zi-zj)*(zi-zj)+em02*dij*dij)
359 mumax=max(mumax,-muij)
360C-------
361C for artificial viscosity computation :
362 rotvx=rotvx+vjy*wgrad(3)-vjz*wgrad(2)
363 rotvy=rotvy+vjz*wgrad(1)-vjx*wgrad(3)
364 rotvz=rotvz+vjx*wgrad(2)-vjy*wgrad(1)
365 END DO
366C------
367C for stability computation into material routines:
368 wa(12,n)=mumax
369C characteristic length = diameter.
370 wa(11,n)=spbuf(1,n)
371C------
372C actual density:
373 divv = wa(1,n)+wa(2,n)+wa(3,n)
374 drho =-divv*rhoi
375c if(drho/=0)print*,'DRHO=',DRHO
376 spbuf(2,n) = max(em20,rhoi+drho*dt1)
377C------
378C for h adaptation + artificial viscosity computation :
379 wa(13,n)=divv
380 rotv =sqrt(rotvx*rotvx+rotvy*rotvy+rotvz*rotvz)
381 wa(14,n)=rotv
382C------
383 10 CONTINUE
384C-----------------------------------------------
385 RETURN
386 END
max
#define max(a, b)
Definition macros.h:21
sphbox::ispsymr
integer, dimension(:,:), allocatable ispsymr
Definition sphbox.F:93
spdens
subroutine spdens(x, v, ms, spbuf, itab, kxsp, ixsp, nod2sp, ispsym, xspsym, vspsym, iparg, wa, wacomp)
Definition spdens.F:36
weight1
subroutine weight1(xi, yi, zi, xj, yj, zj, h, w, wgrad)
Definition weight.F:79