1! 2! CalculiX - A 3-dimensional finite element program 3! Copyright (C) 1998-2021 Guido Dhondt 4! 5! This program is free software; you can redistribute it and/or 6! modify it under the terms of the GNU General Public License as 7! published by the Free Software Foundation(version 2); 8! 9! 10! This program is distributed in the hope that it will be useful, 11! but WITHOUT ANY WARRANTY; without even the implied warranty of 12! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13! GNU General Public License for more details. 14! 15! You should have received a copy of the GNU General Public License 16! along with this program; if not, write to the Free Software 17! Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. 18! 19 subroutine frdfluid(co,nk,konf,ipkonf,lakonf,nef,vold, 20 & kode,time,ielmat,matname,filab,inum,ntrans,inotr,trab,mi, 21 & istep,stn,qfn,nactdohinv,xmach,xkappa,physcon,xturb, 22 & coel,vel,cofa,vfa,nface) 23! 24! stores the results in frd format 25! 26 implicit none 27! 28 character*1 c 29 character*3 m1,m2,m3,m4,m5 30 character*5 p0,p1,p2,p3,p4,p5,p6,p8,p10,p11,p12 31 character*8 lakonf(*),date,newclock,fmat 32 character*10 clock 33 character*20 newdate 34 character*80 matname(*) 35 character*87 filab(*) 36 character*132 text 37! 38 integer konf(*),nk,nef,kode,i,j,ipkonf(*),indexe,inum(*),mi(*), 39 & one,ielmat(mi(3),*),null,inotr(2,*),ntrans,nout,istep, 40 & nactdohinv(*),nface 41! 42 real*8 co(3,*),vold(0:mi(2),*),time,pi,oner,trab(7,*),xturb(2,*), 43 & a(3,3),stn(6,*),qfn(3,*),xmach(*),xkappa(*),physcon(*), 44 & coel(3,*),vel(nef,0:7),vfa(0:7,*),cofa(3,*) 45! 46 save nout 47! 48c do i=1,nef 49c write(*,*) 'frdfluid vel',i,coel(2,i),vel(i,1)/coel(2,i) 50c enddo 51c do i=1,nface 52c write(*,*) 'frdfluid vfa',i,cofa(2,i),vfa(1,i)/cofa(2,i) 53c enddo 54 kode=kode+1 55 pi=4.d0*datan(1.d0) 56! 57 c='C' 58! 59 m1=' -1' 60 m2=' -2' 61 m3=' -3' 62 m4=' -4' 63 m5=' -5' 64! 65 p0=' 0' 66 p1=' 1' 67 p2=' 2' 68 p3=' 3' 69 p4=' 4' 70 p5=' 5' 71 p6=' 6' 72 p8=' 8' 73 p10=' 10' 74 p11=' 11' 75 p12=' 12' 76! 77 fmat(1:8)='(e12.5) ' 78! 79 null=0 80 one=1 81 oner=1.d0 82 83! first time something is written in the frd-file: store 84! computational metadata, the nodal coordinates and the 85! topology */ 86! 87 if(kode.eq.1) then 88! 89 write(13,'(a5,a1)') p1,c 90 call date_and_time(date,clock) 91 newdate(1:20)=' ' 92 newdate(1:2)=date(7:8) 93 newdate(3:3)='.' 94 if(date(5:6).eq.'01') then 95 newdate(4:11)='january.' 96 newdate(12:15)=date(1:4) 97 elseif(date(5:6).eq.'02') then 98 newdate(4:12)='february.' 99 newdate(13:16)=date(1:4) 100 elseif(date(5:6).eq.'03') then 101 newdate(4:9)='march.' 102 newdate(10:13)=date(1:4) 103 elseif(date(5:6).eq.'04') then 104 newdate(4:9)='april.' 105 newdate(10:13)=date(1:4) 106 elseif(date(5:6).eq.'05') then 107 newdate(4:7)='may.' 108 newdate(8:11)=date(1:4) 109 elseif(date(5:6).eq.'06') then 110 newdate(4:8)='june.' 111 newdate(9:12)=date(1:4) 112 elseif(date(5:6).eq.'07') then 113 newdate(4:8)='july.' 114 newdate(9:12)=date(1:4) 115 elseif(date(5:6).eq.'08') then 116 newdate(4:10)='august.' 117 newdate(11:14)=date(1:4) 118 elseif(date(5:6).eq.'09') then 119 newdate(4:13)='september.' 120 newdate(14:17)=date(1:4) 121 elseif(date(5:6).eq.'10') then 122 newdate(4:11)='october.' 123 newdate(12:15)=date(1:4) 124 elseif(date(5:6).eq.'11') then 125 newdate(4:12)='november.' 126 newdate(13:16)=date(1:4) 127 elseif(date(5:6).eq.'12') then 128 newdate(4:12)='december.' 129 newdate(13:16)=date(1:4) 130 endif 131 newclock(1:2)=clock(1:2) 132 newclock(3:3)=':' 133 newclock(4:5)=clock(3:4) 134 newclock(6:6)=':' 135 newclock(7:8)=clock(5:6) 136 write(13,'(a5,''UUSER'')') p1 137 write(13,'(a5,''UDATE'',14x,a20)') p1,newdate 138 write(13,'(a5,''UTIME'',14x,a8)') p1,newclock 139 write(13,'(a5,''UHOST'')') p1 140 write(13,'(a5,''UPGM CalculiX'')') p1 141 write(13,'(a5,''UDIR'')') p1 142 write(13,'(a5,''UDBN'')') p1 143! 144! storing the coordinates of the nodes 145! 146 write(13,'(a5,a1,67x,i1)') p2,c,one 147! 148 nout=0 149 do i=1,nk 150 if(inum(i).le.0) cycle 151 nout=nout+1 152 write(13,100) m1,i,(co(j,i),j=1,3) 153 enddo 154! 155 write(13,'(a3)') m3 156! 157! storing the element topology 158! 159 write(13,'(a5,a1,67x,i1)') p3,c,one 160! 161 do i=1,nef 162! 163 if(ipkonf(i).lt.0) cycle 164 indexe=ipkonf(i) 165 if(lakonf(i)(4:4).eq.'2') then 166 if((lakonf(i)(7:7).eq.' ').or. 167 & (lakonf(i)(7:7).eq.'H')) then 168 write(13,'(a3,i10,3a5)') m1,nactdohinv(i),p4,p0, 169 & matname(ielmat(1,i))(1:5) 170 write(13,'(a3,10i10)') m2,(konf(indexe+j),j=1,10) 171 write(13,'(a3,10i10)') m2,(konf(indexe+j),j=11,12), 172 & (konf(indexe+j),j=17,19),konf(indexe+20), 173 & (konf(indexe+j),j=13,16) 174 elseif(lakonf(i)(7:7).eq.'B') then 175 write(13,'(a3,i10,3a5)')m1,nactdohinv(i),p12,p0, 176 & matname(ielmat(1,i))(1:5) 177 write(13,'(a3,3i10)') m2,konf(indexe+21),konf(indexe+23), 178 & konf(indexe+22) 179 else 180 write(13,'(a3,i10,3a5)')m1,nactdohinv(i),p10,p0, 181 & matname(ielmat(1,i))(1:5) 182 write(13,'(a3,8i10)') m2,(konf(indexe+20+j),j=1,8) 183 endif 184 elseif(lakonf(i)(4:4).eq.'8') then 185 write(13,'(a3,i10,3a5)') m1,nactdohinv(i),p1,p0, 186 & matname(ielmat(1,i))(1:5) 187 write(13,'(a3,8i10)') m2,(konf(indexe+j),j=1,8) 188 elseif(lakonf(i)(4:5).eq.'10') then 189 write(13,'(a3,i10,3a5)') m1,nactdohinv(i),p6,p0, 190 & matname(ielmat(1,i))(1:5) 191 write(13,'(a3,10i10)') m2,(konf(indexe+j),j=1,10) 192 elseif(lakonf(i)(4:4).eq.'4') then 193 write(13,'(a3,i10,3a5)') m1,nactdohinv(i),p3,p0, 194 & matname(ielmat(1,i))(1:5) 195 write(13,'(a3,4i10)') m2,(konf(indexe+j),j=1,4) 196 elseif(lakonf(i)(4:5).eq.'15') then 197 if((lakonf(i)(7:7).eq.' ')) then 198 write(13,'(a3,i10,3a5)') m1,nactdohinv(i),p5,p0, 199 & matname(ielmat(1,i))(1:5) 200 write(13,'(a3,10i10)') m2,(konf(indexe+j),j=1,9), 201 & konf(indexe+13) 202 write(13,'(a3,5i10)') m2,(konf(indexe+j),j=14,15), 203 & (konf(indexe+j),j=10,12) 204 else 205 write(13,'(a3,i10,3a5)') m1,nactdohinv(i),p8,p0, 206 & matname(ielmat(1,i))(1:5) 207 write(13,'(a3,6i10)') m2,(konf(indexe+15+j),j=1,6) 208 endif 209 elseif(lakonf(i)(4:4).eq.'6') then 210 write(13,'(a3,i10,3a5)') m1,nactdohinv(i),p2,p0, 211 & matname(ielmat(1,i))(1:5) 212 write(13,'(a3,6i10)') m2,(konf(indexe+j),j=1,6) 213 elseif(lakonf(i)(1:1).eq.'D') then 214 if((konf(indexe+1).eq.0).or.(konf(indexe+3).eq.0)) cycle 215 write(13,'(a3,i10,3a5)')m1,nactdohinv(i),p12,p0, 216 & matname(ielmat(1,i))(1:5) 217 write(13,'(a3,3i10)') m2,konf(indexe+1),konf(indexe+3), 218 & konf(indexe+2) 219 elseif(lakonf(i)(1:1).eq.'E') then 220 write(13,'(a3,i10,3a5)')m1,nactdohinv(i),p11,p0, 221 & matname(ielmat(1,i))(1:5) 222 write(13,'(a3,2i10)') m2,(konf(indexe+j),j=1,2) 223 endif 224! 225 enddo 226! 227 write(13,'(a3)') m3 228! 229 endif 230! 231 if(filab(34)(1:4).eq.'VF ') then 232 text=' 1PSTEP' 233 write(text(25:36),'(i12)') kode 234 write(13,'(a132)') text 235! 236 text= 237 & ' 100CL .00000E+00 3 1' 238 text(75:75)='1' 239 write(text(25:36),'(i12)') nout 240 write(text(8:12),'(i5)') 100+kode 241 write(text(13:24),fmat) time 242 write(text(59:63),'(i5)') kode 243 write(13,'(a132)') text 244 text=' -4 V3DF 4 1' 245 write(13,'(a132)') text 246 text=' -5 V1 1 2 1 0' 247 write(13,'(a132)') text 248 text=' -5 V2 1 2 2 0' 249 write(13,'(a132)') text 250 text=' -5 V3 1 2 3 0' 251 write(13,'(a132)') text 252 text=' -5 ALL 1 2 0 0 1ALL' 253 write(13,'(a132)') text 254! 255 if((ntrans.eq.0).or.(filab(21)(6:6).eq.'G')) then 256 do i=1,nk 257 if(inum(i).le.0) cycle 258 write(13,100) m1,i,(vold(j,i),j=1,3) 259 enddo 260 else 261 do i=1,nk 262 if(inum(i).le.0) cycle 263 if(inotr(1,i).eq.0) then 264 write(13,100) m1,i,(vold(j,i),j=1,3) 265 else 266 call transformatrix(trab(1,inotr(1,i)),co(1,i),a) 267 write(13,100) m1,i, 268 & vold(1,i)*a(1,1)+vold(2,i)*a(2,1)+vold(3,i)*a(3,1), 269 & vold(1,i)*a(1,2)+vold(2,i)*a(2,2)+vold(3,i)*a(3,2), 270 & vold(1,i)*a(1,3)+vold(2,i)*a(2,3)+vold(3,i)*a(3,3) 271 endif 272 enddo 273 endif 274! 275 write(13,'(a3)') m3 276 endif 277! 278 if(filab(35)(1:4).eq.'PSF ') then 279 text=' 1PSTEP' 280 write(text(25:36),'(i12)') kode 281 write(13,'(a132)') text 282! 283 text= 284 & ' 100CL .00000E+00 3 1' 285 text(75:75)='1' 286 write(text(25:36),'(i12)') nout 287 write(text(8:12),'(i5)') 100+kode 288 write(text(13:24),fmat) time 289 write(text(59:63),'(i5)') kode 290 write(13,'(a132)') text 291 text=' -4 PS3DF 1 1' 292 write(13,'(a132)') text 293 text=' -5 PS 1 1 0 0' 294 write(13,'(a132)') text 295! 296 do i=1,nk 297 if(inum(i).le.0) cycle 298 write(13,100) m1,i,vold(4,i) 299 enddo 300! 301 write(13,'(a3)') m3 302 endif 303! 304 if(filab(36)(1:4).eq.'TSF ') then 305 text=' 1PSTEP' 306 write(text(25:36),'(i12)') kode 307 write(13,'(a132)') text 308! 309 text= 310 & ' 100CL .00000E+00 3 1' 311 text(75:75)='1' 312 write(text(25:36),'(i12)') nout 313 write(text(8:12),'(i5)') 100+kode 314 write(text(13:24),fmat) time 315 write(text(59:63),'(i5)') kode 316 write(13,'(a132)') text 317 text=' -4 TS3DF 1 1' 318 write(13,'(a132)') text 319 text=' -5 TS 1 1 0 0' 320 write(13,'(a132)') text 321! 322 do i=1,nk 323 if(inum(i).le.0) cycle 324 write(13,100) m1,i,vold(0,i) 325 enddo 326! 327 write(13,'(a3)') m3 328 endif 329! 330! xkappa(i) contains kappa (cp/cv) 331! xmach(i) contains the Mach number 332! 333 if(filab(23)(1:4).eq.'MACH') then 334 text=' 1PSTEP' 335 write(text(25:36),'(i12)') kode 336 write(13,'(a132)') text 337! 338 text= 339 & ' 100CL .00000E+00 3 1' 340 text(75:75)='1' 341 write(text(25:36),'(i12)') nout 342 write(text(8:12),'(i5)') 100+kode 343 write(text(13:24),fmat) time 344 write(text(59:63),'(i5)') kode 345 write(13,'(a132)') text 346 text=' -4 M3DF 1 1' 347 write(13,'(a132)') text 348 text=' -5 MACH 1 1 0 0' 349 write(13,'(a132)') text 350! 351 do i=1,nk 352 if(inum(i).le.0) cycle 353 write(13,100) m1,i,xmach(i) 354 enddo 355! 356 write(13,'(a3)') m3 357 endif 358! 359 if(filab(38)(1:4).eq.'TTF ') then 360 text=' 1PSTEP' 361 write(text(25:36),'(i12)') kode 362 write(13,'(a132)') text 363! 364 text= 365 & ' 100CL .00000E+00 3 1' 366 text(75:75)='1' 367 write(text(25:36),'(i12)') nout 368 write(text(8:12),'(i5)') 100+kode 369 write(text(13:24),fmat) time 370 write(text(59:63),'(i5)') kode 371 write(13,'(a132)') text 372 text=' -4 TT3DF 1 1' 373 write(13,'(a132)') text 374 text=' -5 TT 1 1 0 0' 375 write(13,'(a132)') text 376! 377 do i=1,nk 378 if(inum(i).le.0) cycle 379 write(13,100) m1,i, 380 & vold(0,i)*(1.d0+(xkappa(i)-1.d0)/2*xmach(i)**2) 381 enddo 382! 383 write(13,'(a3)') m3 384 endif 385! 386 if(filab(37)(1:4).eq.'PTF ') then 387 text=' 1PSTEP' 388 write(text(25:36),'(i12)') kode 389 write(13,'(a132)') text 390! 391 text= 392 & ' 100CL .00000E+00 3 1' 393 text(75:75)='1' 394 write(text(25:36),'(i12)') nout 395 write(text(8:12),'(i5)') 100+kode 396 write(text(13:24),fmat) time 397 write(text(59:63),'(i5)') kode 398 write(13,'(a132)') text 399 text=' -4 PT3DF 1 1' 400 write(13,'(a132)') text 401 text=' -5 PT 1 1 0 0' 402 write(13,'(a132)') text 403! 404 do i=1,nk 405 if(inum(i).le.0) cycle 406 write(13,100) m1,i,vold(4,i)* 407 & (1.d0+(xkappa(i)-1.d0)/2*xmach(i)**2) 408 & **(xkappa(i)/(xkappa(i)-1.d0)) 409 enddo 410! 411 write(13,'(a3)') m3 412 endif 413! 414! storing the total stresses in the nodes 415! 416 if(filab(39)(1:4).eq.'SF ') then 417 text=' 1PSTEP' 418 write(text(25:36),'(i12)') kode 419 write(text(49:60),'(i12)') istep 420 write(13,'(a132)') text 421! 422 text= 423 & ' 100CL .00000E+00 3 1' 424 text(75:75)='1' 425 write(text(25:36),'(i12)') nout 426 write(text(8:12),'(i5)') 100+kode 427 write(text(13:24),fmat) time 428 write(text(59:63),'(i5)') kode 429 write(13,'(a132)') text 430 text=' -4 STRESS 6 1' 431 write(13,'(a132)') text 432 text=' -5 SXX 1 4 1 1' 433 write(13,'(a132)') text 434 text=' -5 SYY 1 4 2 2' 435 write(13,'(a132)') text 436 text=' -5 SZZ 1 4 3 3' 437 write(13,'(a132)') text 438 text=' -5 SXY 1 4 1 2' 439 write(13,'(a132)') text 440 text=' -5 SYZ 1 4 2 3' 441 write(13,'(a132)') text 442 text=' -5 SZX 1 4 3 1' 443 write(13,'(a132)') text 444 do i=1,nk 445 if(inum(i).le.0) cycle 446 write(13,100) m1,i,(stn(j,i)-vold(4,i),j=1,3), 447 & stn(4,i),stn(6,i),stn(5,i) 448 enddo 449 write(13,'(a3)') m3 450 endif 451! 452! storing the viscous stresses in the nodes 453! 454 if(filab(41)(1:4).eq.'SVF ') then 455 text=' 1PSTEP' 456 write(text(25:36),'(i12)') kode 457 write(13,'(a132)') text 458! 459 text= 460 & ' 100CL .00000E+00 3 1' 461 text(75:75)='1' 462 write(text(25:36),'(i12)') nout 463 write(text(8:12),'(i5)') 100+kode 464 write(text(13:24),fmat) time 465 write(text(59:63),'(i5)') kode 466 write(13,'(a132)') text 467 text=' -4 VSTRES 6 1' 468 write(13,'(a132)') text 469 text=' -5 SXX 1 4 1 1' 470 write(13,'(a132)') text 471 text=' -5 SYY 1 4 2 2' 472 write(13,'(a132)') text 473 text=' -5 SZZ 1 4 3 3' 474 write(13,'(a132)') text 475 text=' -5 SXY 1 4 1 2' 476 write(13,'(a132)') text 477 text=' -5 SYZ 1 4 2 3' 478 write(13,'(a132)') text 479 text=' -5 SZX 1 4 3 1' 480 write(13,'(a132)') text 481 do i=1,nk 482 if(inum(i).le.0) cycle 483 write(13,100) m1,i,(stn(j,i),j=1,4), 484 & stn(6,i),stn(5,i) 485 enddo 486 write(13,'(a3)') m3 487 endif 488! 489! storing the heat flux in the nodes 490! 491 if(filab(40)(1:4).eq.'HFLF') then 492 text=' 1PSTEP' 493 write(text(25:36),'(i12)') kode 494 write(13,'(a132)') text 495! 496 text= 497 & ' 100CL .00000E+00 3 1' 498 text(75:75)='1' 499 write(text(25:36),'(i12)') nout 500 write(text(8:12),'(i5)') 100+kode 501 write(text(13:24),fmat) time 502 write(text(59:63),'(i5)') kode 503 write(13,'(a132)') text 504 text=' -4 FLUX 4 1' 505 write(13,'(a132)') text 506 text=' -5 F1 1 2 1 0' 507 write(13,'(a132)') text 508 text=' -5 F2 1 2 2 0' 509 write(13,'(a132)') text 510 text=' -5 F3 1 2 3 0' 511 write(13,'(a132)') text 512 text=' -5 ALL 1 2 0 0 1ALL' 513 write(13,'(a132)') text 514 do i=1,nk 515 if(inum(i).le.0) cycle 516 write(13,100) m1,i,(qfn(j,i),j=1,3) 517 enddo 518 write(13,'(a3)') m3 519 endif 520! 521 if(filab(24)(1:4).eq.'CP ') then 522 text=' 1PSTEP' 523 write(text(25:36),'(i12)') kode 524 write(13,'(a132)') text 525! 526 text= 527 & ' 100CL .00000E+00 3 1' 528 text(75:75)='1' 529 write(text(25:36),'(i12)') nout 530 write(text(8:12),'(i5)') 100+kode 531 write(text(13:24),fmat) time 532 write(text(59:63),'(i5)') kode 533 write(13,'(a132)') text 534 text=' -4 CP3DF 1 1' 535 write(13,'(a132)') text 536 text=' -5 CP 1 1 0 0' 537 write(13,'(a132)') text 538! 539 do i=1,nk 540 if(inum(i).le.0) cycle 541 write(13,100) m1,i,(vold(4,i)-physcon(6))*2.d0/ 542 & (physcon(7)*physcon(5)**2) 543 enddo 544! 545 write(13,'(a3)') m3 546 endif 547! 548 if(filab(25)(1:4).eq.'TURB') then 549 text=' 1PSTEP' 550 write(text(25:36),'(i12)') kode 551 write(13,'(a132)') text 552! 553 text= 554 & ' 100CL .00000E+00 3 1' 555 text(75:75)='1' 556 write(text(25:36),'(i12)') nout 557 write(text(8:12),'(i5)') 100+kode 558 write(text(13:24),fmat) time 559 write(text(59:63),'(i5)') kode 560 write(13,'(a132)') text 561 text=' -4 TURB3DF 4 1' 562 write(13,'(a132)') text 563 text=' -5 K 1 1 0 0' 564 write(13,'(a132)') text 565 text=' -5 OM 1 1 0 0' 566 write(13,'(a132)') text 567 text=' -5 NU_T 1 1 0 0' 568 write(13,'(a132)') text 569 text=' -5 EPS 1 1 0 0' 570 write(13,'(a132)') text 571! 572 do i=1,nk 573 if(inum(i).le.0) cycle 574 write(13,100) m1,i,xturb(1,i),xturb(2,i), 575 & xturb(1,i)/xturb(2,i),xturb(1,i)*xturb(2,i) 576 enddo 577! 578 write(13,'(a3)') m3 579 endif 580! 581 100 format(a3,i10,1p,6e12.5) 582! 583 return 584 end 585 586 587