OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
spforcp.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!|| spforcp ../engine/source/elements/sph/spforcp.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 spforcp(
33 1 PM ,GEO ,X ,V ,MS ,
34 2 SPBUF ,ITAB ,PLD ,BUFMAT ,BUFGEO ,
35 3 PARTSAV ,FSAV ,DT2T ,IPARG ,NPC ,
36 4 KXSP ,IXSP ,NOD2SP ,NELTST ,ITYPTST ,
37 5 IPART ,IPARTSP ,ISPCOND ,XFRAME ,ISPSYM ,
38 6 XSPSYM ,VSPSYM ,WA ,WASIGSM ,WACOMP ,
39 7 WSMCOMP ,WASPACT ,WAR ,STAB ,WFEXT)
40C-----------------------------------------------
41C M o d u l e s
42C-----------------------------------------------
43 USE sphbox
44C-----------------------------------------------
45C I m p l i c i t T y p e s
46C-----------------------------------------------
47#include "implicit_f.inc"
48#include "comlock.inc"
49C-----------------------------------------------
50C C o m m o n B l o c k s
51C-----------------------------------------------
52#include "vect01_c.inc"
53#include "com06_c.inc"
54#include "com08_c.inc"
55#include "sphcom.inc"
56#include "param_c.inc"
57#include "scr02_c.inc"
58#include "scr17_c.inc"
59#include "task_c.inc"
60C-----------------------------------------------
61C D u m m y A r g u m e n t s
62C-----------------------------------------------
63 INTEGER KXSP(NISP,*),IXSP(KVOISPH,*),NOD2SP(*),
64 . IPART(LIPART1,*) ,IPARTSP(*), NPC(*), IPARG(NPARG,*),
65 . NELTST,ITYPTST,ITAB(*) ,ISPCOND(NISPCOND,*),
66 . ISPSYM(NSPCOND,*), WASPACT(*)
67 my_real
68 . X(3,*) ,V(3,*) ,MS(*) ,
69 . PM(NPROPM,*), GEO(NPROPG,*),
70 . bufmat(*) ,bufgeo(*) ,pld(*) ,fsav(nthvki,*) ,
71 . spbuf(nspbuf,*) ,partsav(npsav,*) ,dt2t ,
72 . xframe(nxframe,*) ,xspsym(3,*) ,vspsym(3,*), wa(kwasph,*),
73 . wasigsm(6,*), wacomp(16,*), wsmcomp(6,*), war(10,*), stab(7,*)
74 DOUBLE PRECISION,INTENT(INOUT) :: WFEXT
75C-----------------------------------------------
76C L o c a l V a r i a b l e s
77C-----------------------------------------------
78 INTEGER N,INOD,JNOD,J,NVOIS,M,IPRT,IPROP,IMAT,I,
79 . NVOISS,IC,NC,IS,SM,JS,ISLIDE,SS,NSTAB,NN,KS,NR,
80 . K,JPERM(KVOISPH),IERROR
81 my_real
82 . xi,yi,zi,di,rhoi,xj,yj,zj,dj,rhoj,
83 . sxx,sxy,sxz,syy,syz,szz,
84 . txx,txy,txz,tyy,tyz,tzz,
85 . ax,ay,az,bx,by,bz,fx,fy,fz,
86 . wght,wgrad(3),
87 . dij,mm,
88 . vxi,vyi,vzi,vxj,vyj,vzj,muij,muij2,pij,ssp,
89 . fact,qa,qb,
90 . ox,oy,oz,ux,uy,uz,vx,vy,vz,wx,wy,wz,
91 . uxx,uxy,uxz,uyx,uyy,uyz,uzx,uzy,uzz,
92 . vxx,vxy,vxz,vyy,vyz,vzz,
93 . divvi,divvj,rotvi,rotvj,fi,fj,
94 . fvx,fvy,fvz,wvis,
95 . stif,pj,stij,sfac,dldt,l,cij,dzeta,wnorm,
96 . vi,vj,vij,
97 . drhoi,drhoj,ssp2i,ssp2j,stii,drhoidr,drhojdr,
98 . alphai,alphaxi,alphayi,alphazi,
99 . betaxi,betayi,betazi,
100 . betaxxi,betayxi,betazxi,
101 . betaxyi,betayyi,betazyi,
102 . betaxzi,betayzi,betazzi,
103 . alphaj,alphaxj,alphayj,alphazj,
104 . betaxj,betayj,betazj,
105 . betaxxj,betayxj,betazxj,
106 . betaxyj,betayyj,betazyj,
107 . betaxzj,betayzj,betazzj,
108 . betax,wgrdx,wgrdy,wgrdz,wgrd(3),
109 . voli,volj,dxij,nxij,
110 . cx,cy,cz,dx,dy,dz,tx,ty,tz,wfextt, ww, wi, wr, zstab
111C-----------------------------------------------
112 wfextt=zero
113 DO 10 i=lft,llt
114 n =nft+i
115 IF(kxsp(2,n)<=0.AND.isph2sol==0)GOTO 10
116 IF(kxsp(2,n)==0.AND.isph2sol/=0)GOTO 10
117 inod =kxsp(3,n)
118 xi=x(1,inod)
119 yi=x(2,inod)
120 zi=x(3,inod)
121 vxi=v(1,inod)
122 vyi=v(2,inod)
123 vzi=v(3,inod)
124 di =spbuf(1,n)
125 rhoi =spbuf(2,n)
126C
127 sxx=wa(1,n)
128 syy=wa(2,n)
129 szz=wa(3,n)
130 sxy=wa(4,n)
131 syz=wa(5,n)
132 sxz=wa(6,n)
133 iprt =ipartsp(n)
134 iprop=ipart(2,iprt)
135 qa = geo(14,iprop)
136 qb = geo(15,iprop)
137C
138C for artificial viscosity computation.
139 divvi=abs(wa(13,n))
140 rotvi=wa(14,n)
141 fi =divvi/max(em20,divvi+rotvi)
142C
143 alphai=wacomp(1,n)
144C BETAXI=WACOMP(2,N)
145C BETAYI=WACOMP(3,N)
146C BETAZI=WACOMP(4,N)
147 alphaxi=wacomp( 5,n)
148 alphayi=wacomp( 6,n)
149 alphazi=wacomp( 7,n)
150 betaxxi=wacomp( 8,n)
151 betayxi=wacomp( 9,n)
152 betazxi=wacomp(10,n)
153 betaxyi=wacomp(11,n)
154 betayyi=wacomp(12,n)
155 betazyi=wacomp(13,n)
156 betaxzi=wacomp(14,n)
157 betayzi=wacomp(15,n)
158 betazzi=wacomp(16,n)
159C----------------------------------
160C FORCES + ARTIFICIAL VISCOUS FORCES + STABILITY:
161C----------------------------------
162 nvois=kxsp(4,n)
163 DO j=1,nvois
164 jnod=ixsp(j,n)
165 IF(jnod>0)THEN
166 m=nod2sp(jnod)
167C
168C Solids to SPH, no interaction if both particles are inactive
169 IF(kxsp(2,n)<=0.AND.kxsp(2,m)<=0)cycle
170C
171 xj=x(1,jnod)
172 yj=x(2,jnod)
173 zj=x(3,jnod)
174 vxj=v(1,jnod)
175 vyj=v(2,jnod)
176 vzj=v(3,jnod)
177 dj =spbuf(1,m)
178 dij=half*(di+dj)
179 rhoj =spbuf(2,m)
180 txx=wa(1,m)
181 tyy=wa(2,m)
182 tzz=wa(3,m)
183 txy=wa(4,m)
184 tyz=wa(5,m)
185 txz=wa(6,m)
186 CALL weight1(xi,yi,zi,xj,yj,zj,dij,wght,wgrad)
187 wgrdx=wgrad(1)
188 wgrdy=wgrad(2)
189 wgrdz=wgrad(3)
190 wgrad(1)=wgrdx*alphai+wght*alphaxi
191 . +wgrdx*betaxxi+wgrdy*betaxyi+wgrdz*betaxzi
192 wgrad(2)=wgrdy*alphai+wght*alphayi
193 . +wgrdx*betayxi+wgrdy*betayyi+wgrdz*betayzi
194 wgrad(3)=wgrdz*alphai+wght*alphazi
195 . +wgrdx*betazxi+wgrdy*betazyi+wgrdz*betazzi
196! Old order1 correction
197! BETAX=1.+BETAXI*(XI-XJ)+BETAYI*(YI-YJ)+BETAZI*(ZI-ZJ)
198! WGRAD(1)=WGRDX*ALPHAI*BETAX
199! . +WGHT*(ALPHAXI*BETAX+ALPHAI*
200! . (BETAXXI*(XI-XJ)+BETAYXI*(YI-YJ)+BETAZXI*(ZI-ZJ)+BETAXI))
201! WGRAD(2)=WGRDY*ALPHAI*BETAX
202! . +WGHT*(ALPHAYI*BETAX+ALPHAI*
203! . (BETAXYI*(XI-XJ)+BETAYYI*(YI-YJ)+BETAZYI*(ZI-ZJ)+BETAYI))
204! WGRAD(3)=WGRDZ*ALPHAI*BETAX
205! . +WGHT*(ALPHAZI*BETAX+ALPHAI*
206! . BETAXZI*(XI-XJ)+BETAYZI*(YI-YJ)+BETAZZI*(ZI-ZJ)+BETAZI))
207C----------
208C noyau conjugue Grad[Wa(b)]
209 alphaj=wacomp(1,m)
210C BETAXJ=WACOMP(2,M)
211C BETAYJ=WACOMP(3,M)
212C BETAZJ=WACOMP(4,M)
213 alphaxj=wacomp( 5,m)
214 alphayj=wacomp( 6,m)
215 alphazj=wacomp( 7,m)
216 betaxxj=wacomp( 8,m)
217 betayxj=wacomp( 9,m)
218 betazxj=wacomp(10,m)
219 betaxyj=wacomp(11,m)
220 betayyj=wacomp(12,m)
221 betazyj=wacomp(13,m)
222 betaxzj=wacomp(14,m)
223 betayzj=wacomp(15,m)
224 betazzj=wacomp(16,m)
225C
226 wgrd(1)=-wgrdx*alphaj+wght*alphaxj
227 . -wgrdx*betaxxj-wgrdy*betaxyj-wgrdz*betaxzj
228 wgrd(2)=-wgrdy*alphaj+wght*alphayj
229 . -wgrdx*betayxj-wgrdy*betayyj-wgrdz*betayzj
230 wgrd(3)=-wgrdz*alphaj+wght*alphazj
231 . -wgrdx*betazxj-wgrdy*betazyj-wgrdz*betazzj
232C
233! Old order1 correction
234! BETAX=ONE +BETAXJ*(XJ-XI)+BETAYJ*(YJ-YI)+BETAZJ*(ZJ-ZI)
235! WGRD(1)=-WGRDX*ALPHAJ*BETAX
236! . +WGHT*(ALPHAXJ*BETAX+ALPHAJ*
237! . (BETAXXJ*(XJ-XI)+BETAYXJ*(YJ-YI)+BETAZXJ*(ZJ-ZI)+BETAXJ))
238! WGRD(2)=-WGRDY*ALPHAJ*BETAX
239! . +WGHT*(ALPHAYJ*BETAX+ALPHAJ*
240! . (BETAXYJ*(XJ-XI)+BETAYYJ*(YJ-YI)+BETAZYJ*(ZJ-ZI)+BETAYJ))
241! WGRD(3)=-WGRDZ*ALPHAJ*BETAX
242! . +WGHT*(ALPHAZJ*BETAX+ALPHAJ*
243! . BETAXZJ*(XJ-XI)+BETAYZJ*(YJ-YI)+BETAZZJ*(ZJ-ZI)+BETAZJ))
244C
245 ax=sxx*wgrad(1)+sxy*wgrad(2)+sxz*wgrad(3)
246 ay=sxy*wgrad(1)+syy*wgrad(2)+syz*wgrad(3)
247 az=sxz*wgrad(1)+syz*wgrad(2)+szz*wgrad(3)
248C BX=TXX*WGRAD(1)+TXY*WGRAD(2)+TXZ*WGRAD(3)
249C BY=TXY*WGRAD(1)+TYY*WGRAD(2)+TYZ*WGRAD(3)
250C BZ=TXZ*WGRAD(1)+TYZ*WGRAD(2)+TZZ*WGRAD(3)
251 bx=-(txx*wgrd(1)+txy*wgrd(2)+txz*wgrd(3))
252 by=-(txy*wgrd(1)+tyy*wgrd(2)+tyz*wgrd(3))
253 bz=-(txz*wgrd(1)+tyz*wgrd(2)+tzz*wgrd(3))
254 mm=spbuf(12,n)*spbuf(12,m)
255C--------
256 vi =spbuf(12,n)/max(em20,rhoi)
257 vj =spbuf(12,m)/max(em20,rhoj)
258 vij=vi*vj
259 fx=vij*(ax+bx)
260 fy=vij*(ay+by)
261 fz=vij*(az+bz)
262C--------
263 wi=zero
264 IF(stab(7,n)/=zero.AND.stab(7,m)/=zero)THEN
265 cx=stab(1,n)*wgrad(1)+stab(4,n)*wgrad(2)+stab(6,n)*wgrad(3)
266 cy=stab(4,n)*wgrad(1)+stab(2,n)*wgrad(2)+stab(5,n)*wgrad(3)
267 cz=stab(6,n)*wgrad(1)+stab(5,n)*wgrad(2)+stab(3,n)*wgrad(3)
268 dx=-(stab(1,m)*wgrd(1)+stab(4,m)*wgrd(2)+stab(6,m)*wgrd(3))
269 dy=-(stab(4,m)*wgrd(1)+stab(2,m)*wgrd(2)+stab(5,m)*wgrd(3))
270 dz=-(stab(6,m)*wgrd(1)+stab(5,m)*wgrd(2)+stab(3,m)*wgrd(3))
271C CALL WEIGHT0(XI,YI,ZI,XJ,YJ,ZJ,DIJ,WGHT)
272 ww=wght*dij*dij*dij
273 wr=half*(stab(7,n)+stab(7,m))
274 wi=ww*ww*ww*ww*wr
275 tx=vij*wi*(cx+dx)
276 ty=vij*wi*(cy+dy)
277 tz=vij*wi*(cz+dz)
278 fx=fx+tx
279 fy=fy+ty
280 fz=fz+tz
281 wfextt=wfextt+(tx*v(1,inod)+ty*v(2,inod)+tz*v(3,inod))*dt1
282 END IF
283C--------
284C artitifial viscosity.
285 dxij=(vxi-vxj)*(xi-xj)
286 . +(vyi-vyj)*(yi-yj)
287 . +(vzi-vzj)*(zi-zj)
288 nxij=(xi-xj)*(xi-xj)
289 . +(yi-yj)*(yi-yj)
290 . +(zi-zj)*(zi-zj)
291 muij=dij*dxij/(nxij+em02*dij*dij)
292 divvj=abs(wa(13,m))
293 rotvj=wa(14,m)
294 fj =divvj/max(em20,divvj+rotvj)
295 muij=min(muij,zero)
296 muij=muij*(fi+fj)*half
297 ssp=(wa(8,n)+wa(8,m))*half
298 muij2=muij*muij
299 pij =(qa*muij2-qb*ssp*muij)*two/max(em20,rhoi+rhoj)
300 fact=mm*pij
301C---------
302C FV(j,i)=-FV(i,j)
303 wgrdx=(wgrad(1)-wgrd(1))*half
304 wgrdy=(wgrad(2)-wgrd(2))*half
305 wgrdz=(wgrad(3)-wgrd(3))*half
306 fvx=-fact*wgrdx
307 fvy=-fact*wgrdy
308 fvz=-fact*wgrdz
309C
310 IF((nodadt/=0).OR.(i7kglo/=0))THEN
311C
312C nodal stability
313 dldt=abs(dxij)
314 l=sqrt(nxij)
315 dldt=max(em20,dldt/max(em20,l))
316C
317 volj=spbuf(12,m)/max(em20,rhoj)
318 drhoidr= volj
319 . * (wgrad(1)*wgrad(1)+wgrad(2)*wgrad(2)+wgrad(3)*wgrad(3))
320 ssp2i=wa(9,n)*wa(9,n)
321 stii = volj*spbuf(12,n)*ssp2i*drhoidr
322C
323 voli =spbuf(12,n)/max(em20,rhoi)
324 drhojdr= voli
325 . * (wgrd(1)*wgrd(1)+wgrd(2)*wgrd(2)+wgrd(3)*wgrd(3))
326 ssp2j=wa(9,m)*wa(9,m)
327 stij = voli*spbuf(12,m)*ssp2j*drhojdr
328C
329 stij=two*(stii+stij)
330C
331C K*=FV divise par DL (C=FV divise par V)
332 wnorm=sqrt(wgrdx*wgrdx+wgrdy*wgrdy+wgrdz*wgrdz)
333 cij =mm*abs(pij)*wnorm/dldt
334 dzeta=cij/max(em20,sqrt(two*stij*spbuf(12,n)))
335 sfac =sqrt(1.+dzeta*dzeta)-dzeta
336 stij =stij/max(em20,sfac*sfac)
337 wa(7,n)=wa(7,n)+stij*(one+wi)
338 ENDIF
339 ELSE ! cellule remote
340 nn = -jnod
341C
342C Solids to SPH, no interaction if both particles are inactive
343 IF(kxsp(2,n)<=0.AND.xsphr(13,nn)<=0)cycle
344 xj=xsphr(3,nn)
345 yj=xsphr(4,nn)
346 zj=xsphr(5,nn)
347 vxj=xsphr(9,nn)
348 vyj=xsphr(10,nn)
349 vzj=xsphr(11,nn)
350 dj =xsphr(2,nn)
351 dij=half*(di+dj)
352 rhoj =xsphr(7,nn)
353 txx=war(1,nn)
354 tyy=war(2,nn)
355 tzz=war(3,nn)
356 txy=war(4,nn)
357 tyz=war(5,nn)
358 txz=war(6,nn)
359 CALL weight1(xi,yi,zi,xj,yj,zj,dij,wght,wgrad)
360 wgrdx=wgrad(1)
361 wgrdy=wgrad(2)
362 wgrdz=wgrad(3)
363 wgrad(1)=wgrdx*alphai+wght*alphaxi
364 . +wgrdx*betaxxi+wgrdy*betaxyi+wgrdz*betaxzi
365 wgrad(2)=wgrdy*alphai+wght*alphayi
366 . +wgrdx*betayxi+wgrdy*betayyi+wgrdz*betayzi
367 wgrad(3)=wgrdz*alphai+wght*alphazi
368 . +wgrdx*betazxi+wgrdy*betazyi+wgrdz*betazzi
369! Old order1 correction
370! BETAX=1.+BETAXI*(XI-XJ)+BETAYI*(YI-YJ)+BETAZI*(ZI-ZJ)
371! WGRAD(1)=WGRDX*ALPHAI*BETAX
372! . +WGHT*(ALPHAXI*BETAX+ALPHAI*
373! . (betaxxi*(xi-xj)+betayxi*(yi-yj)+betazxi*(zi-zj)+betaxi))
374! WGRAD(2)=WGRDY*ALPHAI*BETAX
375! . +WGHT*(ALPHAYI*BETAX+ALPHAI*
376! . (BETAXYI*(XI-XJ)+BETAYYI*(YI-YJ)+BETAZYI*(ZI-ZJ)+BETAYI))
377! WGRAD(3)=WGRDZ*ALPHAI*BETAX
378! . +WGHT*(ALPHAZI*BETAX+ALPHAI*
379! . (BETAXZI*(XI-XJ)+BETAYZI*(YI-YJ)+BETAZZI*(ZI-ZJ)+BETAZI))
380C----------
381C noyau conjugue Grad[Wa(b)]
382 alphaj=wacompr(1,nn)
383C BETAXJ=WACOMPR(2,NN)
384C BETAYJ=WACOMPR(3,NN)
385C BETAZJ=WACOMPR(4,NN)
386 alphaxj=wacompr( 5,nn)
387 alphayj=wacompr( 6,nn)
388 alphazj=wacompr( 7,nn)
389 betaxxj=wacompr( 8,nn)
390 betayxj=wacompr( 9,nn)
391 betazxj=wacompr(10,nn)
392 betaxyj=wacompr(11,nn)
393 betayyj=wacompr(12,nn)
394 betazyj=wacompr(13,nn)
395 betaxzj=wacompr(14,nn)
396 betayzj=wacompr(15,nn)
397 betazzj=wacompr(16,nn)
398C
399 wgrd(1)=-wgrdx*alphaj+wght*alphaxj
400 . -wgrdx*betaxxj-wgrdy*betaxyj-wgrdz*betaxzj
401 wgrd(2)=-wgrdy*alphaj+wght*alphayj
402 . -wgrdx*betayxj-wgrdy*betayyj-wgrdz*betayzj
403 wgrd(3)=-wgrdz*alphaj+wght*alphazj
404 . -wgrdx*betazxj-wgrdy*betazyj-wgrdz*betazzj
405! Old order1 correction
406! BETAX=ONE +BETAXJ*(XJ-XI)+BETAYJ*(YJ-YI)+BETAZJ*(ZJ-ZI)
407! WGRD(1)=-WGRDX*ALPHAJ*BETAX
408! . +WGHT*(ALPHAXJ*BETAX+ALPHAJ*
409! . (BETAXXJ*(XJ-XI)+BETAYXJ*(YJ-YI)+BETAZXJ*(ZJ-ZI)+BETAXJ))
410! WGRD(2)=-WGRDY*ALPHAJ*BETAX
411! . +WGHT*(ALPHAYJ*BETAX+ALPHAJ*
412! . (BETAXYJ*(XJ-XI)+BETAYYJ*(YJ-YI)+BETAZYJ*(ZJ-ZI)+BETAYJ))
413! WGRD(3)=-WGRDZ*ALPHAJ*BETAX
414! . +WGHT*(ALPHAZJ*BETAX+ALPHAJ*
415! . (BETAXZJ*(XJ-XI)+BETAYZJ*(YJ-YI)+BETAZZJ*(ZJ-ZI)+BETAZJ))
416C
417 ax=sxx*wgrad(1)+sxy*wgrad(2)+sxz*wgrad(3)
418 ay=sxy*wgrad(1)+syy*wgrad(2)+syz*wgrad(3)
419 az=sxz*wgrad(1)+syz*wgrad(2)+szz*wgrad(3)
420C BX=TXX*WGRAD(1)+TXY*WGRAD(2)+TXZ*WGRAD(3)
421C BY=TXY*WGRAD(1)+TYY*WGRAD(2)+TYZ*WGRAD(3)
422C BZ=TXZ*WGRAD(1)+TYZ*WGRAD(2)+TZZ*WGRAD(3)
423 bx=-(txx*wgrd(1)+txy*wgrd(2)+txz*wgrd(3))
424 by=-(txy*wgrd(1)+tyy*wgrd(2)+tyz*wgrd(3))
425 bz=-(txz*wgrd(1)+tyz*wgrd(2)+tzz*wgrd(3))
426 mm=spbuf(12,n)*xsphr(8,nn)
427C--------
428 vi =spbuf(12,n)/max(em20,rhoi)
429 vj =xsphr(8,nn)/max(em20,rhoj)
430 vij=vi*vj
431 fx=vij*(ax+bx)
432 fy=vij*(ay+by)
433 fz=vij*(az+bz)
434C--------
435 wi=zero
436 IF(stab(7,n)/=zero.AND.stab(7,numsph+nn)/=zero)THEN
437 cx=stab(1,n)*wgrad(1)+stab(4,n)*wgrad(2)+stab(6,n)*wgrad(3)
438 cy=stab(4,n)*wgrad(1)+stab(2,n)*wgrad(2)+stab(5,n)*wgrad(3)
439 cz=stab(6,n)*wgrad(1)+stab(5,n)*wgrad(2)+stab(3,n)*wgrad(3)
440 nr=numsph+nn
441 dx=-(stab(1,nr)*wgrd(1)+stab(4,nr)*wgrd(2)+stab(6,nr)*wgrd(3))
442 dy=-(stab(4,nr)*wgrd(1)+stab(2,nr)*wgrd(2)+stab(5,nr)*wgrd(3))
443 dz=-(stab(6,nr)*wgrd(1)+stab(5,nr)*wgrd(2)+stab(3,nr)*wgrd(3))
444C CALL WEIGHT0(XI,YI,ZI,XJ,YJ,ZJ,DIJ,WGHT)
445 ww=wght*dij*dij*dij
446 wr=half*(stab(7,n)+stab(7,numsph+nn))
447 wi=ww*ww*ww*ww*wr
448 tx=vij*wi*(cx+dx)
449 ty=vij*wi*(cy+dy)
450 tz=vij*wi*(cz+dz)
451 fx=fx+tx
452 fy=fy+ty
453 fz=fz+tz
454 wfextt=wfextt+(tx*v(1,inod)+ty*v(2,inod)+tz*v(3,inod))*dt1
455 END IF
456C--------
457C artitifial viscosity.
458 dxij=(vxi-vxj)*(xi-xj)
459 . +(vyi-vyj)*(yi-yj)
460 . +(vzi-vzj)*(zi-zj)
461 nxij=(xi-xj)*(xi-xj)
462 . +(yi-yj)*(yi-yj)
463 . +(zi-zj)*(zi-zj)
464 muij=dij*dxij/(nxij+em02*dij*dij)
465C WA(13) <=> WAR(9) ; WA(14) <=> WAR(10)
466 divvj=abs(war(9,nn))
467 rotvj=war(10,nn)
468 fj =divvj/max(em20,divvj+rotvj)
469 muij=min(muij,zero)
470 muij=muij*(fi+fj)*half
471C WA(8) <=> WAR(7)
472 ssp=(wa(8,n)+war(7,nn))*half
473 muij2=muij*muij
474 pij =(qa*muij2-qb*ssp*muij)*two/max(em20,rhoi+rhoj)
475 fact=mm*pij
476C---------
477C FV(j,i)=-FV(i,j)
478 wgrdx=(wgrad(1)-wgrd(1))*half
479 wgrdy=(wgrad(2)-wgrd(2))*half
480 wgrdz=(wgrad(3)-wgrd(3))*half
481 fvx=-fact*wgrdx
482 fvy=-fact*wgrdy
483 fvz=-fact*wgrdz
484 IF((nodadt/=0).OR.(i7kglo/=0))THEN
485C
486C nodal stability
487 dldt=abs(dxij)
488 l=sqrt(nxij)
489 dldt=max(em20,dldt/max(em20,l))
490C
491 volj=xsphr(8,nn)/max(em20,rhoj)
492 drhoidr= volj
493 . * (wgrad(1)*wgrad(1)+wgrad(2)*wgrad(2)+wgrad(3)*wgrad(3))
494 ssp2i=wa(9,n)*wa(9,n)
495 stii = volj*spbuf(12,n)*ssp2i*drhoidr
496C
497 voli =spbuf(12,n)/max(em20,rhoi)
498 drhojdr= voli
499 . * (wgrd(1)*wgrd(1)+wgrd(2)*wgrd(2)+wgrd(3)*wgrd(3))
500C WA(9) <=> WAR(8)
501 ssp2j=war(8,nn)*war(8,nn)
502 stij = voli*xsphr(8,nn)*ssp2j*drhojdr
503C
504 stij=two*(stii+stij)
505C
506C K*=FV divise par DL (C=FV divise par V)
507 wnorm=sqrt(wgrdx*wgrdx+wgrdy*wgrdy+wgrdz*wgrdz)
508 cij =mm*abs(pij)*wnorm/dldt
509 dzeta=cij/max(em20,sqrt(two*stij*spbuf(12,n)))
510 sfac =sqrt(1.+dzeta*dzeta)-dzeta
511 stij =stij/max(em20,sfac*sfac)
512 wa(7,n)=wa(7,n)+stij*(one+wi)
513 ENDIF
514 END IF
515C
516 fx=fx+fvx
517 fy=fy+fvy
518 fz=fz+fvz
519 wa(10,n)=wa(10,n)+fx
520 wa(11,n)=wa(11,n)+fy
521 wa(12,n)=wa(12,n)+fz
522C
523C pour travail des forces de visc. artificielle (par cellule).
524 wvis = min(zero,
525 . half*(fvx*(vxi-vxj)+fvy*(vyi-vyj)+fvz*(vzi-vzj)))
526 spbuf(11,n)=spbuf(11,n)+wvis
527 ENDDO
528C----------------------------------
529 nvoiss=kxsp(6,n)
530 DO 200 j=kxsp(5,n)+1,kxsp(5,n)+nvoiss
531 js=ixsp(j,n)
532 IF(js>0)THEN
533 sm=js/(nspcond+1)
534C
535C Solids to SPH, no interaction if both particles are inactive
536 IF(kxsp(2,n)<=0.AND.kxsp(2,sm)<=0)cycle
537C
538 nc=mod(js,nspcond+1)
539 js=ispsym(nc,sm)
540 xj =xspsym(1,js)
541 yj =xspsym(2,js)
542 zj =xspsym(3,js)
543 vxj=vspsym(1,js)
544 vyj=vspsym(2,js)
545 vzj=vspsym(3,js)
546 dj =spbuf(1,sm)
547 dij =half*(di+dj)
548 rhoj=spbuf(2,sm)
549 jnod=kxsp(3,sm)
550C
551 txx=wasigsm(1,js)
552 tyy=wasigsm(2,js)
553 tzz=wasigsm(3,js)
554 txy=wasigsm(4,js)
555 tyz=wasigsm(5,js)
556 txz=wasigsm(6,js)
557 CALL weight1(xi,yi,zi,xj,yj,zj,dij,wght,wgrad)
558 wgrdx=wgrad(1)
559 wgrdy=wgrad(2)
560 wgrdz=wgrad(3)
561 wgrad(1)=wgrdx*alphai+wght*alphaxi
562 . +wgrdx*betaxxi+wgrdy*betaxyi+wgrdz*betaxzi
563 wgrad(2)=wgrdy*alphai+wght*alphayi
564 . +wgrdx*betayxi+wgrdy*betayyi+wgrdz*betayzi
565 wgrad(3)=wgrdz*alphai+wght*alphazi
566 . +wgrdx*betazxi+wgrdy*betazyi+wgrdz*betazzi
567! Old order1 correction
568! BETAX=ONE +BETAXI*(XI-XJ)+BETAYI*(YI-YJ)+BETAZI*(ZI-ZJ)
569! WGRAD(1)=WGRDX*ALPHAI*BETAX
570! . +WGHT*(ALPHAXI*BETAX+ALPHAI*
571! . (BETAXXI*(XI-XJ)+BETAYXI*(YI-YJ)+BETAZXI*(ZI-ZJ)+BETAXI))
572! WGRAD(2)=WGRDY*ALPHAI*BETAX
573! . +WGHT*(ALPHAYI*BETAX+ALPHAI*
574! . (BETAXYI*(XI-XJ)+BETAYYI*(YI-YJ)+BETAZYI*(ZI-ZJ)+BETAYI))
575! WGRAD(3)=WGRDZ*ALPHAI*BETAX
576! . +WGHT*(ALPHAZI*BETAX+ALPHAI*
577! . (BETAXZI*(XI-XJ)+BETAYZI*(YI-YJ)+BETAZZI*(ZI-ZJ)+BETAZI))
578C----------
579C noyau conjugue.
580 alphaj=wacomp(1,sm)
581C BETAXJ=WACOMP(2,SM)
582C BETAYJ=WACOMP(3,SM)
583C BETAZJ=WACOMP(4,SM)
584 betaxj=wsmcomp(1,js)
585 betayj=wsmcomp(2,js)
586 betazj=wsmcomp(3,js)
587C ALPHAXJ=WACOMP( 5,SM)
588C ALPHAYJ=WACOMP( 6,SM)
589C ALPHAZJ=WACOMP( 7,SM)
590 alphaxj=wsmcomp( 4,js)
591 alphayj=wsmcomp( 5,js)
592 alphazj=wsmcomp( 6,js)
593 betaxxj=wacomp( 8,sm)
594 betayxj=wacomp( 9,sm)
595 betazxj=wacomp(10,sm)
596 betaxyj=wacomp(11,sm)
597 betayyj=wacomp(12,sm)
598 betazyj=wacomp(13,sm)
599 betaxzj=wacomp(14,sm)
600 betayzj=wacomp(15,sm)
601 betazzj=wacomp(16,sm)
602 wgrd(1)=-wgrdx*alphaj+wght*alphaxj
603 . -wgrdx*betaxxj-wgrdy*betaxyj-wgrdz*betaxzj
604 wgrd(2)=-wgrdy*alphaj+wght*alphayj
605 . -wgrdx*betayxj-wgrdy*betayyj-wgrdz*betayzj
606 wgrd(3)=-wgrdz*alphaj+wght*alphazj
607 . -wgrdx*betazxj-wgrdy*betazyj-wgrdz*betazzj
608! Old order1 correction
609! BETAX=ONE +BETAXJ*(XJ-XI)+BETAYJ*(YJ-YI)+BETAZJ*(ZJ-ZI)
610! WGRD(1)=-WGRDX*ALPHAJ*BETAX
611! . +WGHT*(ALPHAXJ*BETAX+ALPHAJ*
612! . (BETAXXJ*(XJ-XI)+BETAYXJ*(YJ-YI)+BETAZXJ*(ZJ-ZI)+BETAXJ))
613! WGRD(2)=-WGRDY*ALPHAJ*BETAX
614! . +WGHT*(ALPHAYJ*BETAX+ALPHAJ*
615! . (BETAXYJ*(XJ-XI)+BETAYYJ*(YJ-YI)+BETAZYJ*(ZJ-ZI)+BETAYJ))
616! WGRD(3)=-WGRDZ*ALPHAJ*BETAX
617! . +WGHT*(ALPHAZJ*BETAX+ALPHAJ*
618! . (betaxzj*(xj-xi)+betayzj*(yj-yi)+betazzj*(zj-zi)+betazj))
619 ax=sxx*wgrad(1)+sxy*wgrad(2)+sxz*wgrad(3)
620 ay=sxy*wgrad(1)+syy*wgrad(2)+syz*wgrad(3)
621 az=sxz*wgrad(1)+syz*wgrad(2)+szz*wgrad(3)
622 bx=-(txx*wgrd(1)+txy*wgrd(2)+txz*wgrd(3))
623 by=-(txy*wgrd(1)+tyy*wgrd(2)+tyz*wgrd(3))
624 bz=-(txz*wgrd(1)+tyz*wgrd(2)+tzz*wgrd(3))
625 mm=spbuf(12,n)*spbuf(12,sm)
626C--------
627 vi =spbuf(12,n)/max(em20,rhoi)
628 vj =spbuf(12,sm)/max(em20,rhoj)
629 vij=vi*vj
630 fx=vij*(ax+bx)
631 fy=vij*(ay+by)
632 fz=vij*(az+bz)
633 wi=zero
634 IF(stab(7,n)/=zero.AND.stab(7,sm)/=zero)THEN
635 cx=stab(1,n)*wgrad(1)+stab(4,n)*wgrad(2)+stab(6,n)*wgrad(3)
636 cy=stab(4,n)*wgrad(1)+stab(2,n)*wgrad(2)+stab(5,n)*wgrad(3)
637 cz=stab(6,n)*wgrad(1)+stab(5,n)*wgrad(2)+stab(3,n)*wgrad(3)
638 ks=numsph+nsphr+js
639 dx=-(stab(1,ks)*wgrd(1)+stab(4,ks)*wgrd(2)+stab(6,ks)*wgrd(3))
640 dy=-(stab(4,ks)*wgrd(1)+stab(2,ks)*wgrd(2)+stab(5,ks)*wgrd(3))
641 dz=-(stab(6,ks)*wgrd(1)+stab(5,ks)*wgrd(2)+stab(3,ks)*wgrd(3))
642C CALL WEIGHT0(XI,YI,ZI,XJ,YJ,ZJ,DIJ,WGHT)
643 ww=wght*dij*dij*dij
644 wr=half*(stab(7,n)+stab(7,sm))
645 wi=ww*ww*ww*ww*wr
646 tx=vij*wi*(cx+dx)
647 ty=vij*wi*(cy+dy)
648 tz=vij*wi*(cz+dz)
649 fx=fx+tx
650 fy=fy+ty
651 fz=fz+tz
652 wfextt=wfextt+(tx*v(1,inod)+ty*v(2,inod)+tz*v(3,inod))*dt1
653 END IF
654C--------
655C artitifial viscosity.
656 dxij=(vxi-vxj)*(xi-xj)
657 . +(vyi-vyj)*(yi-yj)
658 . +(vzi-vzj)*(zi-zj)
659 nxij=(xi-xj)*(xi-xj)
660 . +(yi-yj)*(yi-yj)
661 . +(zi-zj)*(zi-zj)
662 muij=dij*dxij/(nxij+em02*dij*dij)
663 divvj=abs(wa(13,sm))
664 rotvj=wa(14,sm)
665 fj=divvj/max(em20,divvj+rotvj)
666 muij=min(muij,zero)
667 muij=muij*(fi+fj)*half
668 ssp=(wa(8,n)+wa(8,sm))*half
669 muij2=muij*muij
670 pij =(qa*muij2-qb*ssp*muij)*two/max(em20,rhoi+rhoj)
671 fact=mm*pij
672C---------
673 wgrdx=(wgrad(1)-wgrd(1))*half
674 wgrdy=(wgrad(2)-wgrd(2))*half
675 wgrdz=(wgrad(3)-wgrd(3))*half
676 fvx=-fact*wgrdx
677 fvy=-fact*wgrdy
678 fvz=-fact*wgrdz
679 IF((nodadt/=0).OR.(i7kglo/=0))THEN
680C nodal stability
681 dldt=abs(dxij)
682 l =sqrt(nxij)
683 dldt=max(em20,dldt/max(em20,l))
684C
685 volj=spbuf(12,sm)/max(em20,rhoj)
686 drhoidr= volj
687 . * (wgrad(1)*wgrad(1)+wgrad(2)*wgrad(2)+wgrad(3)*wgrad(3))
688 ssp2i=wa(9,n)*wa(9,n)
689 stii = volj*spbuf(12,n)*ssp2i*drhoidr
690C
691 voli =spbuf(12,n)/max(em20,rhoi)
692 drhojdr= voli
693 . * (wgrd(1)*wgrd(1)+wgrd(2)*wgrd(2)+wgrd(3)*wgrd(3))
694 ssp2j=wa(9,sm)*wa(9,sm)
695 stij = voli*spbuf(12,sm)*ssp2j*drhojdr
696C
697 stij=two*(stii+stij)
698C
699C K*=FV divise par DL (C=FV divise par V)
700 wnorm=sqrt(wgrdx*wgrdx+wgrdy*wgrdy+wgrdz*wgrdz)
701 cij =mm*abs(pij)*wnorm/dldt
702 dzeta=cij/max(em20,sqrt(two*stij*spbuf(12,n)))
703 sfac =sqrt(one +dzeta*dzeta)-dzeta
704 stij =stij/max(em20,sfac*sfac)
705 wa(7,n)=wa(7,n)+stij*(one+wi)
706 ENDIF
707 fx=fx+fvx
708 fy=fy+fvy
709 fz=fz+fvz
710 wa(10,n)=wa(10,n)+fx
711 wa(11,n)=wa(11,n)+fy
712 wa(12,n)=wa(12,n)+fz
713C
714C pour travail des forces de visc. artificielle (par cellule).
715 wvis = min(zero,
716 . half*(fvx*(vxi-vxj)+fvy*(vyi-vyj)+fvz*(vzi-vzj)))
717 spbuf(11,n)=spbuf(11,n)+wvis
718 ELSE ! particule symetrique de particule remote
719 sm=-js/(nspcond+1)
720C
721C Solids to SPH, no interaction if both particles are inactive
722 IF(kxsp(2,n)<=0.AND.xsphr(13,sm)<=0)cycle
723 nc=mod(-js,nspcond+1)
724 js=ispsymr(nc,sm)
725 xj =xspsym(1,js)
726 yj =xspsym(2,js)
727 zj =xspsym(3,js)
728 vxj=vspsym(1,js)
729 vyj=vspsym(2,js)
730 vzj=vspsym(3,js)
731 dj =xsphr(2,sm)
732 dij =half*(di+dj)
733 rhoj=xsphr(7,sm)
734C
735 txx=wasigsm(1,js)
736 tyy=wasigsm(2,js)
737 tzz=wasigsm(3,js)
738 txy=wasigsm(4,js)
739 tyz=wasigsm(5,js)
740 txz=wasigsm(6,js)
741 CALL weight1(xi,yi,zi,xj,yj,zj,dij,wght,wgrad)
742 wgrdx=wgrad(1)
743 wgrdy=wgrad(2)
744 wgrdz=wgrad(3)
745 wgrad(1)=wgrdx*alphai+wght*alphaxi
746 . +wgrdx*betaxxi+wgrdy*betaxyi+wgrdz*betaxzi
747 wgrad(2)=wgrdy*alphai+wght*alphayi
748 . +wgrdx*betayxi+wgrdy*betayyi+wgrdz*betayzi
749 wgrad(3)=wgrdz*alphai+wght*alphazi
750 . +wgrdx*betazxi+wgrdy*betazyi+wgrdz*betazzi
751! Old order1 correction
752! BETAX=ONE +BETAXI*(XI-XJ)+BETAYI*(YI-YJ)+BETAZI*(ZI-ZJ)
753! WGRAD(1)=WGRDX*ALPHAI*BETAX
754! . +WGHT*(ALPHAXI*BETAX+ALPHAI*
755! . (BETAXXI*(XI-XJ)+BETAYXI*(YI-YJ)+BETAZXI*(ZI-ZJ)+BETAXI))
756! WGRAD(2)=WGRDY*ALPHAI*BETAX
757! . +WGHT*(ALPHAYI*BETAX+ALPHAI*
758! . (BETAXYI*(XI-XJ)+BETAYYI*(YI-YJ)+BETAZYI*(ZI-ZJ)+BETAYI))
759! WGRAD(3)=WGRDZ*ALPHAI*BETAX
760! . +WGHT*(ALPHAZI*BETAX+ALPHAI*
761! . (BETAXZI*(XI-XJ)+BETAYZI*(YI-YJ)+BETAZZI*(ZI-ZJ)+BETAZI))
762C----------
763C noyau conjugue.
764 alphaj=wacompr(1,sm)
765C BETAXJ=WACOMPR(2,SM)
766C BETAYJ=WACOMPR(3,SM)
767C BETAZJ=WACOMPR(4,SM)
768 betaxj=wsmcomp(1,js)
769 betayj=wsmcomp(2,js)
770 betazj=wsmcomp(3,js)
771C ALPHAXJ=WACOMPR( 5,SM)
772C ALPHAYJ=WACOMPR( 6,SM)
773C ALPHAZJ=WACOMPR( 7,SM)
774 alphaxj=wsmcomp( 4,js)
775 alphayj=wsmcomp( 5,js)
776 alphazj=wsmcomp( 6,js)
777 betaxxj=wacompr( 8,sm)
778 betayxj=wacompr( 9,sm)
779 betazxj=wacompr(10,sm)
780 betaxyj=wacompr(11,sm)
781 betayyj=wacompr(12,sm)
782 betazyj=wacompr(13,sm)
783 betaxzj=wacompr(14,sm)
784 betayzj=wacompr(15,sm)
785 betazzj=wacompr(16,sm)
786C
787 wgrd(1)=-wgrdx*alphaj+wght*alphaxj
788 . -wgrdx*betaxxj-wgrdy*betaxyj-wgrdz*betaxzj
789 wgrd(2)=-wgrdy*alphaj+wght*alphayj
790 . -wgrdx*betayxj-wgrdy*betayyj-wgrdz*betayzj
791 wgrd(3)=-wgrdz*alphaj+wght*alphazj
792 . -wgrdx*betazxj-wgrdy*betazyj-wgrdz*betazzj
793! Old order1 correction
794! BETAX=ONE +BETAXJ*(XJ-XI)+BETAYJ*(YJ-YI)+BETAZJ*(ZJ-ZI)
795! WGRD(1)=-WGRDX*ALPHAJ*BETAX
796! . +WGHT*(ALPHAXJ*BETAX+ALPHAJ*
797! . (BETAXXJ*(XJ-XI)+BETAYXJ*(YJ-YI)+BETAZXJ*(ZJ-ZI)+BETAXJ))
798! WGRD(2)=-WGRDY*ALPHAJ*BETAX
799! . +WGHT*(ALPHAYJ*BETAX+ALPHAJ*
800! . (BETAXYJ*(XJ-XI)+BETAYYJ*(YJ-YI)+BETAZYJ*(ZJ-ZI)+BETAYJ))
801! WGRD(3)=-WGRDZ*ALPHAJ*BETAX
802! . +WGHT*(ALPHAZJ*BETAX+ALPHAJ*
803! . (BETAXZJ*(XJ-XI)+BETAYZJ*(YJ-YI)+BETAZZJ*(ZJ-ZI)+BETAZJ))
804 ax=sxx*wgrad(1)+sxy*wgrad(2)+sxz*wgrad(3)
805 ay=sxy*wgrad(1)+syy*wgrad(2)+syz*wgrad(3)
806 az=sxz*wgrad(1)+syz*wgrad(2)+szz*wgrad(3)
807 bx=-(txx*wgrd(1)+txy*wgrd(2)+txz*wgrd(3))
808 by=-(txy*wgrd(1)+tyy*wgrd(2)+tyz*wgrd(3))
809 bz=-(txz*wgrd(1)+tyz*wgrd(2)+tzz*wgrd(3))
810 mm=spbuf(12,n)*xsphr(8,sm)
811C--------
812 vi =spbuf(12,n)/max(em20,rhoi)
813 vj =xsphr(8,sm)/max(em20,rhoj)
814 vij=vi*vj
815 fx=vij*(ax+bx)
816 fy=vij*(ay+by)
817 fz=vij*(az+bz)
818 wi=zero
819 IF(stab(7,n)/=zero.AND.stab(7,numsph+sm)/=zero)THEN
820 cx=stab(1,n)*wgrad(1)+stab(4,n)*wgrad(2)+stab(6,n)*wgrad(3)
821 cy=stab(4,n)*wgrad(1)+stab(2,n)*wgrad(2)+stab(5,n)*wgrad(3)
822 cz=stab(6,n)*wgrad(1)+stab(5,n)*wgrad(2)+stab(3,n)*wgrad(3)
823 ks=numsph+nsphr+js
824 dx=-(stab(1,ks)*wgrd(1)+stab(4,ks)*wgrd(2)+stab(6,ks)*wgrd(3))
825 dy=-(stab(4,ks)*wgrd(1)+stab(2,ks)*wgrd(2)+stab(5,ks)*wgrd(3))
826 dz=-(stab(6,ks)*wgrd(1)+stab(5,ks)*wgrd(2)+stab(3,ks)*wgrd(3))
827C CALL WEIGHT0(XI,YI,ZI,XJ,YJ,ZJ,DIJ,WGHT)
828 ww=wght*dij*dij*dij
829 wr=half*(stab(7,n)+stab(7,numsph+sm))
830 wi=ww*ww*ww*ww*wr
831 tx=vij*wi*(cx+dx)
832 ty=vij*wi*(cy+dy)
833 tz=vij*wi*(cz+dz)
834 fx=fx+tx
835 fy=fy+ty
836 fz=fz+tz
837 wfextt=wfextt+(tx*v(1,inod)+ty*v(2,inod)+tz*v(3,inod))*dt1
838 END IF
839C--------
840C artitifial viscosity.
841 dxij=(vxi-vxj)*(xi-xj)
842 . +(vyi-vyj)*(yi-yj)
843 . +(vzi-vzj)*(zi-zj)
844 nxij=(xi-xj)*(xi-xj)
845 . +(yi-yj)*(yi-yj)
846 . +(zi-zj)*(zi-zj)
847 muij=dij*dxij/(nxij+em02*dij*dij)
848C WA(13) <=> WAR(9) ; WA(14) <=> WAR(10)
849 divvj=abs(war(9,sm))
850 rotvj=war(10,sm)
851 fj=divvj/max(em20,divvj+rotvj)
852 muij=min(muij,zero)
853 muij=muij*(fi+fj)*half
854C WA(8) <=> WAR(7)
855 ssp=(wa(8,n)+war(7,sm))*half
856 muij2=muij*muij
857 pij =(qa*muij2-qb*ssp*muij)*two/max(em20,rhoi+rhoj)
858 fact=mm*pij
859C---------
860 wgrdx=(wgrad(1)-wgrd(1))*half
861 wgrdy=(wgrad(2)-wgrd(2))*half
862 wgrdz=(wgrad(3)-wgrd(3))*half
863 fvx=-fact*wgrdx
864 fvy=-fact*wgrdy
865 fvz=-fact*wgrdz
866 IF((nodadt/=0).OR.(i7kglo/=0))THEN
867C nodal stability
868 dldt=abs(dxij)
869 l =sqrt(nxij)
870 dldt=max(em20,dldt/max(em20,l))
871C
872 volj=xsphr(8,sm)/max(em20,rhoj)
873 drhoidr= volj
874 . * (wgrad(1)*wgrad(1)+wgrad(2)*wgrad(2)+wgrad(3)*wgrad(3))
875 ssp2i=wa(9,n)*wa(9,n)
876 stii = volj*spbuf(12,n)*ssp2i*drhoidr
877C
878 voli =spbuf(12,n)/max(em20,rhoi)
879 drhojdr= voli
880 . * (wgrd(1)*wgrd(1)+wgrd(2)*wgrd(2)+wgrd(3)*wgrd(3))
881C WA(9) <=> WAR(8)
882 ssp2j=war(8,sm)*war(8,sm)
883 stij = voli*xsphr(8,sm)*ssp2j*drhojdr
884C
885 stij=two*(stii+stij)
886C
887C K*=FV divise par DL (C=FV divise par V)
888 wnorm=sqrt(wgrdx*wgrdx+wgrdy*wgrdy+wgrdz*wgrdz)
889 cij =mm*abs(pij)*wnorm/dldt
890 dzeta=cij/max(em20,sqrt(two*stij*spbuf(12,n)))
891 sfac =sqrt(one +dzeta*dzeta)-dzeta
892 stij =stij/max(em20,sfac*sfac)
893 wa(7,n)=wa(7,n)+stij*(one+wi)
894 ENDIF
895 fx=fx+fvx
896 fy=fy+fvy
897 fz=fz+fvz
898 wa(10,n)=wa(10,n)+fx
899 wa(11,n)=wa(11,n)+fy
900 wa(12,n)=wa(12,n)+fz
901C
902C pour travail des forces de visc. artificielle (par cellule).
903 wvis = min(zero,
904 . half*(fvx*(vxi-vxj)+fvy*(vyi-vyj)+fvz*(vzi-vzj)))
905 spbuf(11,n)=spbuf(11,n)+wvis
906 END IF
907 200 CONTINUE
908C-------
909 10 CONTINUE
910#include "lockon.inc"
911 wfext=wfext+wfextt
912#include "lockoff.inc"
913C----------------------------------
914 RETURN
915 END
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
sphbox::ispsymr
integer, dimension(:,:), allocatable ispsymr
Definition sphbox.F:93
sphbox::nsphr
integer nsphr
Definition sphbox.F:83
spforcp
subroutine spforcp(pm, geo, x, v, ms, spbuf, itab, pld, bufmat, bufgeo, partsav, fsav, dt2t, iparg, npc, kxsp, ixsp, nod2sp, neltst, ityptst, ipart, ipartsp, ispcond, xframe, ispsym, xspsym, vspsym, wa, wasigsm, wacomp, wsmcomp, waspact, war, stab, wfext)
Definition spforcp.F:40
weight1
subroutine weight1(xi, yi, zi, xj, yj, zj, h, w, wgrad)
Definition weight.F:79