Přednáška cvičení 1.10.2002

Co Vás chci naučit:

Analýzu systémů

1

Univerzální programy

Vlastní program řešení teplotního pole. Metoda sítí. Program v Power Fortranu (terminal server). Porovnání C, Fortran

Diskrétní systémy, bilance (PRO II, ASPEN, …) Charakteristika: Aparát je chápán jako BLACK BOX: - integrální bilance (uzlů – aparátů)- algebraické rov.(stacionární), obyč.dif.rov. (nestacionární)- řešení Newtonovou metodou nebo metodami Runge Kutta

Kontinuální systémy – proudění, vícefázové systémy, teplo, pružnost, ,..,Charakteristika: Analýza procesů, které se odehrávají v aparátech WHITE BOX- diferenciální bilance - parciální diferenciální rovnice (CFD transportní děje), variační principy (pružnost)- řešení metodu sítí, kontrolních objemů, konečných prvků, hraničních prvků,

spektrální metody

Komerční programypružnost pevnost, teplotní poleCOSMOS/M, ANSYS, ABAQUS, ADINA, NASTRAN…CFD (Com.Fluid Dynam.)FLUENT, CFX

Freeware SSIIM (3D water and sediment flow) LCA codes for astrophysical simulation and visualization PHOENICS/CHAM shareware (multi-phase flow, N-S, combustion) NASA LeRC (LPDF2D) Collection of CFD codes for ship design Software for 3-D Fluid Interfaces (TVD/AC) GASDYN (Euler, PVM) ASME software collection Hanley Innovations PHASES Engineering Solutions DIFFPACK

vlastní programyFEMINA

Speciální programy a Postprocessing

Komerční programyExcel, GNUplot

Vlastní programy: havarie MitasTUPLEX foulingTemelín chlazeníSpirálový filtrVakuové chlazeníKalcinační pecKolimované detektoryMikrovlnný ohřevTok tixotropní kapaliny

Problém odhadu, kterou kategorii použít – zda zvolit komerční program nebo psát svůj vlastní

Příklad: Zpracování potravin za vysokého tlaku, počáteční řešení programem COSMOS – problémy při zadávání konstitutivní rovnice vody za vysokého tlaku

www. icemcfd. com /cfd/CFD_codes.html

ACRi ARSoftware (TEP: a combustion analysis tool for windows) COSMIC NASA Fluent Inc. (FLUENT v5, FIDAP, POLYFLOW, GAMBIT, TGrid, Icepak, Airpak, MixSim) Flowtech Int. AB (SHIPFLOW: analysis of flow around ships) Fluid Dynamics International, Inc. (FIDAP) AEA Technology (CFX: 3D fluid flow/heat transfer code) ICEM CFD (ICEM CFD, Icepak) KIVA (reactive flows) CFD Research Corporation (ACE: reactive flows) Computational Dynamics Ltd. (STAR-CD) Analytical Methods, Inc. (VSAERO, USAERO, OMNI3D) AeroSoft, Inc. (GASP and GUST) Ithaca Combustion Enterprises (PDF2DS) Flow Science, Inc. (FLOW3D) ALGOR, Inc. (ALGOR) Engineering Mechanics Research Corp. (NISA) Reaction Engineering International (BANFF/GLACIER) Combustion Dynamics Ltd. (SuperSTATE) AVL List Gmbh. (FIRE) IBM Corp. catalogue (30 positions) Sun Microsystems catalogue (70 positions) Cray Research catalogue (100 positions) Silicon Graphics, Inc. catalogue (75 positions) Pointwise, Inc. (Gridgen - structured grids) Simulog (N3S Finite Element code, MUSCL) Directory of CFD codes on IBM supercomputer environment ANSYS, Inc. (FLOTRAN) Flomercis Inc. (FLOTHERM) Computational Mechanics Corporation Computational Mechanics Company, Inc. (COMCO) KASIMIR (shock tube simulation program) Livermore Software Technology Corporation (LS-DYNA3D) Advanced Combustion Eng. Research Center (PCGC, FBED) NUMECA International s.a. (FINE, FINE/Turbo, FINE/Aero, IGG, IGG/Autogrid) Computational Engineering International., Inc. (EnSight, ...) Blocon Software Agency (HEAT2, HEAT3) Adaptive Research Corp. (CFD2000) Unicom Technology Systems (VORSTAB-PC) Incinerator Consultants Incorporated (ICI) PHOENICS/CHAM (multi-phase flow, N-S, combustion) Innovative Aerodynamic Technologies (LAMDA) XYZ Scientific Applications, Inc. (TrueGrid) South Bay Simulations, Inc. (SPLASH) PHASES Engineering Solutions Amtec Engineering, Inc. (INCA, Tecplot) Engineering Sciences, Inc. (UNIC) Catalpa Research, Inc. (TIGER) Swansea NS codes (LAM2D, TURB) Engineering Systems International S.A. (PAM-FLOW, PAM-FLUID) Daat Research Corp. (COOLIT) Flomerics Inc. (FLOVENT) Innovative Research, Inc. Centric Engineering Systems, Inc. (SPECTRUM) Blue Ridge Numerics, Inc. WinPipeD Exa Corporation (PowerFLOW) Polyflow s.a. Flow Pro

Computational Aerodynamics Systems Co. Tahoe Design Software ADINA-F YFLOW PSW Advanced Visual Systems Flo++ KSNIS Flowcode Concert SMARTFIRE VISCOUS Polydynamics Cullimore and Ring Technologies, Inc. (SINDA/FLUINT, SINAPS) Linflow (ANKER - ZEMER ENGINEERING) PFDReaction Airfoil Analysis Institute of Computational Continuum Mechanics GmbH CFD++ RADIOSS-CFD VECTIS MAYA Simulation Compass Arena Flow

Příklad informací o KIWAKIVA (modeling of reactive flows) for more information search OSTI WAIS gate with keyword KIVAArticle: 576 of sci.physics.computational.fluid-dynamics From Robin B. Lake / rbl@hal.cwru.eduDate: 10 Jul 1994 14:52:30 GMT

In article <2vo5bn$5jl@news.u.washington.edu> hinkey@aa.washington.edu (John Hinkey) writes:>>Does anyone know of any CFD codes capable of >modeling a two-phase detonation, i.e. a >gaseous oxidizer with a solid or liquid fuel>suspended in it?>There is a code "KIVA" that BP uses for such explosive environments asinternal combustion engine cylinders. In that it can handle chemicaltransformations, it may be suitable.

In that BP has reorganized and eliminated its US Research, I can't giveyou the name of any competent staff at their Sunbury-on-Thames ResearchCentre in the UK. You might try to contact Dr. Wendell Mills, EngineeringComputing Inc., who occupies space at BP's former Warrensville ResearchCenter in Cleveland. Call the central BP number in Cleveland. Can't helpmuch, as they've changed the phone numbers recently (and probably thelocks, too!).

Rob LakeEnvironmental Modeling Inc.rbl@hal.cwru.edu

See also KIVA page provided by Fluid Dynamics Group (T-3) at LANL.

2

Vedení tepla – relativně snadno řešitelný problém.

Fourierova Kirchhoffova rovnice (u,v složky rychlosti)

q=0 kartézský, q=1 cylindrický souřadný systém

Speciální případ kartézský souřadný systém, bez proudění a konstantní vlastnosti

Stacionární případ a bez vnitřního zdroje tepla – Laplaceova rovnice

Diskretizace náhrada derivací diferencemi

Diferenční náhrada Laplaceovy rovnice je tedy

a speciální případ pro hx=hy=h .

Okrajové podmínky 1.druhu (Dirichlet)

3

ii-1 i+1

h

Přednáška 7.10.02

2.druhu (Neaumann)

3.druhu (Newton)

Konec přednášky 2.10.2002

Nahrazení první derivace symetrickou diferencí. Zavedení fiktivního bodu mimo hranici oblasti

Eliminace fiktivní hodnoty z Laplaceovy rovnice

Speciální případ – osa symetrie A=0 – nulový tepelný tok ( )

a speciální případ pro hx=hy=h

Přibližně ( )

4

ii-1 i+1

h

Příklad programu

Elektrolyzér

Rozložení elektrického potenciálu Laplaceova rovnice

Elektrody U=constStěna (isolant) nulová normálová derivace potenciálu.

Hustota výkonu

Postup zpracování v režimu Terminal Server (Citrix)1. stáhnout FEMINA.ZIP2. rozbalit do adresáře femina3. překopírovat MINIPF.LIB F.BAT FC.BAT FPSVARS.BAT

Program překládaný příkazem F THERM.FORfl32 /MW /4Yb %1 /link minipf.lib

Základy fortranu

Příkladinteger preal a(10)write(*,*)’ Nazdar Rudolfe, zadej pocet bodu’read(*,*)pdo i=1,p

write(*,*)’ zadej ‘,i,’ –tou hodnotu vektoru’ read(*,*)a(i) enddo end

Jména proměnných, procedur, funkcí musí začínat písmenem (nebo $)Typy proměnných Fortran REAL INTEGER COMPLEX LOGICAL CHARACTERdeklarace REAL I,VEC(10),A(10,20) (indexy od 1) implicitní typy první písmeno I,J,K,L,M,N integer

5

1.fáze – rozložení elektrického potenciálu

2.fáze – teplotní pole

U=1U=0

0yU

0yU

COMMON /jméno bloku/ seznam proměnných

Příklad: FEMINA PARAMETER (MAXAUX=70)

COMMON /FEM/IAUX(MAXAUX),NPT,NCR,NSF,…

Příklad přenosu hodnot mezi moduly

common /fem/…,npt,……… npt=2 npt=2 call proc call procp(npt) write(*,*)npt write(*,*)npt end end

subroutine proc subroutine procp(np) common /fem/…,np,… np=np+1 np=np+1 end end

Přiřazovací příkaz (každý příkaz na nové řádce), C příkazy jsou oddělovány ;proměnná op=výraz je ekvivalentní s výrazem proměnná=proměnná op (výraz) s=s+1 (C s++, nebo s+=1)res=res+abs(a-b) (C res+=abs(a-b))a(i,j)=sin(x) (C a[i][j]=sin(x))operátory +-*/ ** ( ) (C nemá umocňování, operátor % modulo)

Inkrementace / dekrementace v C++proměnná - hodnota proměnné je nejprve zvětšena o jedničku a teprve potom je tato nová hodnota vrácena jako hodnota aritmetického výrazu proměnná++ - nejprve je vrácena původní hodnota aritmetického výrazu a teprve potom je hodnota proměnné zvětšena o jedničku.

Podmínky a cyklydo while(podmínka) while(podmínka)… {…}end do ;

do i=1,n for(i=1;i<=n,i++)… {…}end do ;

if(podmínka)then if(podmínka)… {…};else else… {…};endif

Podmínkarelační operátory .LT. .GT. .LE. .GE. .EQ. .NE.logické operátory .AND. .OR. .NOT.(v jazyce C && || == <>)

Ternární logické výrazy v C

6

Vazba skutečné a formální parametry

Sdílení operační paměti

Chceme-li vyhodnotit výraz v závislosti na nějaké logické podmínce, použijeme tzv. podmíněný výraz, který v C se zapisuje pomocí tzv. ternárního operátoru ? : Chceme-li např. vypočítat hodnotu sin(x) / x pro x různé od nuly a hodnotu 1 pro x = 0 , píšeme x ? sin(x)/x : 1Chceme-li realizovat výpočet maxima ze dvou zadaných hodnot, můžeme psát (a > b) ? a : b

Příkazy vstupu a výstupuread(unit,format)a,b,… scanfwrite(unit,format)a,b,… printfZnak '\n' je interpretován jako přechod na novou řádku. Kombinace %d, %f jsou konverze určené pro výstup hodnot typu int a double. (obecně %[příznaky][šířka][.přesnost][modifikátor]konverze)

Operace se souboryFILE *soubor; ...soubor=fopen("DATA.TXT","r"); ...fclose(soubor); ...

7

ELEKROLYZER ANSI C#include <stdio.h>#include <math.h> #define MAX 101main(){int m,me,i,j,iter; double u[MAX][MAX],e,w,n,s,res,ua; printf("\n\nDimension M:"); scanf("%d",&m); me=m/2; printf("\nElectrode at %i M=%i",me,m); for(i=1;i<=m;i++) for(j=1;j<=m;j++) { if(i=me&&j<=me) u[i][j]=1.; else u[i][j]=0.; }; iter=0; res=1.; while(iter<50000&&res>1e-6) {res=0.;iter++; for(i=2;i<=m;i++) {for(j=1;j<=m;j++) if(i!=me||j>me) {e=u[i!=m?i+1:i-1][j]; w=u[i!=1?i-1:i+1][j]; n=u[i][j!=m?j+1:j-1]; s=u[i][j!=1?j-1:j+1]; ua=(e+w+n+s)/4; res+=abs(u[i][j]-ua); u[i][j]=ua; }; }; res/=(m*m); } printf("\n\nRes =%f iter=%i \n",res,iter);} 560 znakůANSI FORTRAN parameter (MAX=101) dimension u(max,max) real n read(*,*)m me=m/2

do i=1,m do j=1,m if(i.ne.me.or.j.gt.me)then u(i,j)=0 else u(i,j)=1 endif enddo enddo res=1 iter=0do while(res.gt.1e-6.and.iter.le.100000) res=0 do i=2,m do j=1,m if(i.eq.m)then e=u(i-1,j) else e=u(i+1,j) endif if(j.eq.1)then s=u(i,j+1) else s=u(i,j-1) endif if(j.eq.m)then n=u(i,j-1) else n=u(i,j+1) endif w=u(i-1,j) if(i.ne.me.or.j.gt.me)then ua=(w+e+n+s)/4 res=res+abs(ua-u(i,j)) u(i,j)=ua endif enddo enddo res=res/(m*m) iter=iter+1 enddo write(*,*)'Iter=',iter,' res=',res end

492 znaků

8

9

SGI-ORIGIN (sobota-malé zatížení)

0

10

20

30

40

50

60

70

80

90

50 60 70 80 90 100N

t [s]

C-SGIFortran-SGIFortran-GIN

Řešení programem FEMINAC* Opening elec date:07.10.02 at 10:29PT 1,0,0;PT 2,1,0;PT 3,1,1;PT 4,0,1;PT 5,.5,0;PT 6,.5,.5;SF4PT 1,1,2,3,4;CR2PT 5,5,6;MSF 1,50,50,1,1,4;NFCR 5,1,-1,1,1,1;NFCR 4,1,-1,0,0,0;EGROUP 1,0,0,0;THERMA 0,1,1;

Počet Gaussovych bodů 1x1 síť 50x50 čtyřúhelníkové prvky. Počet Gaussových uzlů 2x2

Q=elementy, 3x3 Gaussovy body. Triangles 1 Gaussův uzel

10

1

6

5

4 3

2

CR 5CR 4

Q 1 Q 2

Q 3 T 1

11

Trojúhelníková sit pouze 10 x 10tj. 200 elementu, s jedním Gaussovým uzlem.

Čtyřúhelníková sit pouze 10 x 10tj. 200 elementu, s jedním Gaussovým uzlem.

Čtyřúhelníková sit pouze 10 x 10tj. 200 elementu, s 2 x 2 Gaussovými uzly.

Q 1 Q 2

T 1

THERM0.FORc stacionarni teplotni pole. Okrajove podminky Dirichlet. Bez vystupu. use msflib dimension vektor(10000) write(*,*)' Nazdar Rudolfe, zadej H,NX,NY' read(*,*)h,nx,ny write(*,*)' Nx=',nx,' Ny=',ny,' H=',h call therm(vektor,nx,ny,h) end

subroutine therm(a,nx,ny,h) dimension a(nx,ny) a=0. hx=h/(nx-1) hy=h/(ny-1) do i=1,nxc Okrajova podminka Dirichletova pro Y=H x=(i-1)*hx a(i,ny)=100*sin(x) enddo do iter=1,100c Fixni pocet iteraci do i=2,nx-1 do j=2,ny-1 a(i,j)= / ((a(i-1,j)+a(i+1,j))/hx**2+(a(i,j-1)+a(i,j+1))/hy**2)/ / (2./hx**2+2./hy**2) enddo enddo write(*,*)' Iteration=',iter enddo end

Některé užitečné procedury z knihovny MINIPF.LIB

SUBROUTINE GMFSTW(INDS,ITYPE,IWINS,IFONT,IHEIGH,LABELX,LABELY, / CXMI,CXMA,CYMI,CYMA,XMI,XMA,YMI,YMA) C NASTAVENI SOUR. SYSTEMU : C INDS = 1,2,3,4,5,6 C ITYPE= 0 slouzi k vykreslovani vrstevnic (vymezuje se prostor pro zobrazeni skaly vrstevnic), nezaokrouhluji se rozsahy os. C = 1 aritmeticke osy se zaokrouhlovanim zadaneho rozsahu. C = 2 log-log, lisi se jen jinym pomerem velikosti vnitrniho a vnejsiho okna. C IWIN = 1,2,...40 child window C IFONT= 1,2,3,4,5 font (Courier, Arial, Times, Sans Serif, Small) C IHEIGH vyska fontu (v pixelech) C LABELX,LABELY *16 popisy os

SUBROUTINE GMFAXE(IND,ITYPC,ICOLI4)C TRANSFORMACE A OSY ITYP=0,4 LOG - LOG C ITYP=1,5 LOG - ARITM C ITYP=2,6 ARITM - LOG C ITYP=3,7 ARITM -ARITM C ITYPC>=10 VYBARVENI OKRAJE BARVOU ITYPC/10.

SUBROUTINE GMFRCTR(IND,LABEL,NLAB,VAL,M,N, / B1MI,B1MA,B2MI,B2MA,VMIN,VMAX,IVMIN,IVMAX,ICOL) C CONTOURS RECTANGULARS NET. Spojite spektrum barev, ale navic se barvou ICOL kresli i cary - vrstevnice. C IND - INDEX OKNA, LABEL (TEXT NLAB ZNAKU) C VAL(M,N) - FUNKCNI HODNOTY V OBDELNIKOVE MATICI C B1MI,... - OBDELNIK (X,Y), OHRANICENY BARVOU ICOL C VMIN,VMAX- ROZSAH VAL, IVMIN,IVMAX ODPOVIDAJICI KODY BAREV (0-2850) C ICOL - BARVA VRSTEVNIC A RAMECKU SUBROUTINE GMFPLD(XMI,DX,N,Y,ICOLR4)C Průběh EKVIDISTANTNI FUNKCE C XMI POCATECNI HODNOTA ARGUMENTU, DX PRIRUSTEK, N POCET BODU C Y(N) VEKTOR HODNOT, ICOLR4-BARVA CARY SUBROUTINE GMFPLF(N,X,Y,ICOL,IFIL)C FILL CURVE (N, X(N),Y(N), ICOL) BY COLOUR IFIL C CUBIC SPLINES - STRAIGHT LINE CLOSE THE CURVE

SUBROUTINE GMFPL(N,X,Y,ICOLR4)

12

T=100 sin x

T=0

C POLYLINE X(N),Y(N) barvou ICOLR4. THERM1.FORc Nadrelaxacni faktor a testovani residua. Grafika vrstevnice. use msflib dimension vektor(10000) write(*,*)' Nazdar Rudolfe, zadej H,NX,NY,omega'c omega je nadrelaxacni faktor read(*,*)h,nx,ny,omega write(*,*)' Nx=',nx,' Ny=',ny,' H=',h call therm(vektor,nx,ny,h,omega) end

subroutine therm(a,nx,ny,h,omega) dimension a(nx,ny)c label pouzijeme pro hlavicku obrazku character*11 label a=0. hx=h/(nx-1) hy=h/(ny-1) do i=1,nx x=(i-1)*hx a(i,ny)=100*sin(x) enddo do iter=1,1000

res=0 do i=2,nx-1 do j=2,ny-1 aux= / ((a(i-1,j)+a(i+1,j))/hx**2+(a(i,j-1)+a(i,j+1))/hy**2)/ /

(2./hx**2+2./hy**2) res=res+abs(a(i,j)-aux)

a(i,j)=(1-omega)*a(i,j)+omega*aux enddo enddo res=res/(nx*ny) write(*,*)' iter=',iter,' resid=',res if(res.le.1e-4)exit enddoc grafika call wopen(10,'teploty',30,90,10,20,480,320) call gmfstw(1,0,10,3,14,'x','y',0.,1.,0.,1.,0.,h,0.,h) call gmfaxe(1,7,15) write(label,'('' Niter='',i4)')iter call gmfrctr(1,label,11,a,nx,ny,0.,h,0.,h,0.,100.,0,800,0) end

THERM2.FORc doplneni okrajove podminky 3.druhu a kresleni grafu T(y). Vysledky do souboru. use msflib dimension vektor(10000) write(*,*)' Nazdar Rudolfe, zadej H,NX,NY,omega,alfa,te' read(*,*)h,nx,ny,omega,alfa,te write(*,*)' Nx=',nx,' Ny=',ny,' H=',h call therm(vektor,nx,ny,h,omega,alfa,te) end

subroutine therm(a,nx,ny,h,omega,alfa,te) dimension a(nx,ny),gy(100) character*36 label a=0. hx=h/(nx-1) hy=h/(ny-1) do i=1,nx x=(i-1)*hx a(i,ny)=100*sin(x) enddo do iter=1,1000

res=0 do j=2,ny-1 do i=2,nx if(i.eq.nx)then

13

c okrajova podminka tretiho druhu aux=(alfa*te+a(i-1,j)/hx+hx*(a(i,j-1)+a(i,j+1))/hy**2)/ / (alfa+1/hx+hx/hy**2) elsec vnitrni body oblasti aux= / ((a(i-1,j)+a(i+1,j))/hx**2+(a(i,j-1)+a(i,j+1))/hy**2)/

/ (2./hx**2+2./hy**2) endif

res=res+abs(a(i,j)-aux)c superrelaxacni metoda a(i,j)=(1-omega)*a(i,j)+omega*aux

enddo enddo

res=res/(nx*ny) write(*,*)' iter=',iter,' resid=',res if(res.le.1e-4)exit

enddocall wopen(10,'teploty',30,90,0,20,480,320)call gmfstw(1,0,10,3,14,'x','y',0.,1.,0.,1.,0.,h,0.,h)call gmfaxe(1,7,15)write(label,'('' Niter='',i4,'' alfa='',e10.3,'' Te='',f5.0)')

/iter,alfa,tecall gmfrctr(1,label,36,a,nx,ny,0.,h,0.,h,0.,100.,0,800,0)call wopen(20,'teploty',30,90,1,21,480,320)

tmax=0 do i=1,ny tmax=amax1(tmax,a(nx,i)) gy(i)=a(nx,i) enddo

call gmfstw(2,1,20,3,14,'y','T',0.,1.,0.,1.,0.,h,0.,tmax)call gmfaxe(2,3,15)call gmfpld(0.,hy,ny,gy,12)

c ulozeni vysledku do souboru open(1,file='therm2.out') do i=1,nx write(1,'(i3,100e10.3)')i,(a(i,j),j=1,ny) enddo end

Úkol pro studenty:Nalézt optimální nadrelaxační faktor. Doplnit procedurou pro hledání minima.

14

Nestacionární problém vedení tepla

Fourierova Kirchhoffova rovnice

Diskretizace

Odhad velikosti časového krokuUvažujme případ s nulovým zdrojovým členem a s nulovými teplotami v čase t s výjimkou jediné Tij=1.

Pro by nová teplota byla nulová (a teploty na všech dalších hladinách také nulové).

Pro větší hodnoty t by se objevily oscilace.

Program pro nestacionární vedení tepla – explicitní methoda

15

tento člen by měl být vyčíslen na časové hladině t+t plně implicitní metoda. Stabilní, ale časově náročné.

tento člen je na časové hladině t. Explicitní metoda. Podmíněně stabilní, musí být malý časový krok.

H=1 Nx=Ny=21 a=1e-6 dt=600 ze stabilitního kriteria dt=625c Nestacionarni vedeni tepla use msflib dimension vold(10000),vnew(10000) write(*,*)' Nazdar Rudolfe, zadej H,NX,NY,Ntime,dt' read(*,*)h,nx,ny,nt,dt write(*,*)' Nx=',nx,' Ny=',ny,' H=',h call therm(vold,vnew,nx,ny,h,nt,dt) end

subroutine therm(aold,anew,nx,ny,h,nt,dt)c assuming dimension aold(nx,ny),anew(nx,ny) character*36 labelc teplotni vodivost konstantni a=1e-6 a=1e-6c nulova pocatecni podminka aold=0 hx=h/(nx-1) hy=h/(ny-1) do i=1,nx x=(i-1)*hx aold(i,ny)=100*sin(x) anew(i,ny)=aold(i,ny) enddo call wopen(10,'teploty',30,90,0,20,480,320) call gmfstw(1,0,10,3,14,'x','y',0.,1.,0.,1.,0.,h,0.,h) call gmfaxe(1,7,15) do it=1,nt do j=2,ny-1 do i=2,nx-1c vnitrni body oblasti- explicitni metoda aux=aold(i,j)*(1/dt-2*a*(1/hx**2+1/hy**2))+ / a*((aold(i-1,j)+aold(i+1,j))/hx**2+ / (aold(i,j-1)+aold(i,j+1))/hy**2) anew(i,j)=aux*dt if(abs(anew(i,j)).gt.200.)stop

enddoenddo

if(mod(it,nt/10).eq.0)then write(label,'('' Nt='',i4,'' dt='',e10.3)')it,dt call gmfrctr(1,label,22,anew,nx,ny,0.,h,0.,h,0.,100.,0,800,0) endif do j=2,ny-1 do i=2,nx-1 aold(i,j)=anew(i,j)

enddoenddo

enddo endÚkol pro studenty: Naprogramovat plně implicitní metodu

16

ADI Implicitní metoda střídavých směrů

x-implicitní

y-implicitní

Diskretizacex-implicitní

y-implicitní

Kdybychom hledali jen stacionární řešení byla by optimální volba časového kroku

,

Pro hrubý odhad relativní velikosti časových půlkroků lze využít informaci o rychlosti šíření teplotní poruchy odpovídající konduktivnímu přenosu tepla

čemuž odpovídá poměr časových půlkroků

17

c Nestacionarni vedeni tepla ADI use msflib dimension vold(10000),vnew(10000) write(*,*)' Nazdar Rudolfe, zadej H,NX,NY,Ntime,dt' read(*,*)h,nx,ny,nt,dt write(*,*)' Nx=',nx,' Ny=',ny,' H=',h call therm(vold,vnew,nx,ny,h,nt,dt) end

subroutine therm(aold,anew,nx,ny,h,nt,dt)c assuming dimension aold(nx,ny),anew(nx,ny) dimension aa(100),bb(100),cc(100),dd(100),uu(100) character*36 labelc teplotni vodivost konstantni a=1e-6 a=1e-6c nulova pocatecni podminka anew=0c polovicni casovy krok x,y dtx=dt/2 dty=dt/2 hx=h/(nx-1) hy=h/(ny-1) do i=1,nx x=(i-1)*hx anew(i,ny)=100*sin(x) enddo aold=anew call wopen(10,'teploty',30,90,0,20,480,320) call gmfstw(1,0,10,3,14,'x','y',0.,1.,0.,1.,0.,h,0.,h) call gmfaxe(1,7,15) do it=1,nt atxx=a*dtx/hx**2 atxy=a*dtx/hy**2 atyy=a*dty/hy**2 atyx=a*dty/hx**2c x-implicitni do j=2,ny-1 aa(1)=0 bb(1)=1 cc(1)=0 aa(nx)=0

bb(nx)=1 cc(nx)=0 dd(1)=anew(1,j) dd(nx)=anew(nx,j) do i=2,nx-1 aa(i)=-atxx bb(i)=1+2*atxx cc(i)=-atxx

dd(i)=anew(i,j)+atxy*(anew(i,j+1)-2*anew(i,j)+anew(i,j-1)) enddo call tridag(aa,bb,cc,dd,uu,nx) do i=2,nx-1 aold(i,j)=uu(i) enddo enddoc y-implicitni do i=2,nx-1 aa(1)=0 bb(1)=1 cc(1)=0 aa(ny)=0 bb(ny)=1 cc(ny)=0 dd(1)=anew(i,1) dd(ny)=anew(i,ny) do j=2,ny-1 aa(j)=-atyy bb(j)=1+2*atyy cc(j)=-atyy

dd(j)=aold(i,j)+atyx*(aold(i+1,j)-2*aold(i,j)+aold(i-1,j)) enddo call tridag(aa,bb,cc,dd,uu,ny) do j=2,ny-1 anew(i,j)=uu(j) enddo

enddo if(mod(it,nt/10).eq.0)then write(label,'('' Nt='',i4,'' dt='',e10.3)')it,dt call gmfrctr(1,label,22,anew,nx,ny,0.,h,0.,h,0.,100.,0,800,0) endif enddo end

18

Řešení programem FEMINA

C* Opening test date:16.10.02 at 10:05PT 1,0,0;PT 2,1,0;PT 3,1,1;PT 4,0,1;SF4PT 1,1,2,3,4;EGROUP 1,0,3,0;MPROP 1,.600E+00,0,4200,0,998,0,.400E-01,0,.210E+12,0,.280E+00,0,.100E-02,0,0,0,0,0,0,0,0,0RCONST 1,1,0,0,0,0,0,0,0,0,0;MSF 1,20,20,1,1,4;FUNDEF -2,100*SIN(XX);NFCR 1,1,-1,0,0,0;NFCR 2,1,-1,0,0,0;NFCR 4,1,-1,0,0,0;NFCR 3,1,-2,1,1,1;INITIA 1,0;THERMA 0,1,1;EGROUP 1,1,3,0;INITIA 1,0;THERMA 0,100,1000;

Stacionární teplotní pole 20 x 20 Q3

Nestacionární teplotní pole 20 x 20 Q3 dt=1000 Nt=100

19

Okrajová podmínka 3.druhu je asi nejkomplikovanější

C* Opening test date:16.10.02 at 13:28PT 1,0,0;PT 2,1,0;PT 3,1,1;PT 4,0,1;SF4PT 1,1,2,3,4;EGROUP 1,0,3,0;MPROP 1,.600E+00,0,4200,0,998,0,.400E-01,0,.210E+12,0,.280E+00,0,.100E-02,0,0,0,0,0,0,0,0,0

RCONST 1,1,0,0,0,20,0,0,0,0,0;MSF 1,20,20,1,1,4;FUNDEF -2,100*SIN(XX);NFCR 1,1,-1,0,0,0;NFCR 4,1,-1,0,0,0;NFCR 3,1,-2,1,1,1;

NFCR 2,1,21,200,200,200;THERMA 0,1,0;

GCR 2,temp vykreslení průběhu teplot podél křivky číslo 2

20

Konvektivní přenos tepla

Povšimněte si analogie s rovnicí

můžeme použít implicitní nebo explicitní metodu řešení (souřadnice z hraje roli evoluční proměnné – času – vývoj teplotního profilu).

Stabilitní analýza

, speciálně pro 1D problém

Metoda kontrolních objemůUvažujme jen jednorozměrný problémDivergentní formulace

Pro nestlačitelnou kapalinu je tato rovnice totožná s konzervativní formulací

Pro konstantní u,a

Centrální formule

Stabilita

TW=0, TE=1 TP bude kladné jen pro Peh<2.

Protiproudá formule (upwind) 1.řád přesnosti

21

z

x

y

EWW W P

ew

pro u>0

, a TP bude vždy kladné.

Protiproudá formule (upwind) 2.řád přesnostiLineární extrapolace

Protiproudá formule QUICKkvadratický polynom

stabilita

Řád přesnostiTaylorův rozvoj řešení kolem bodu P

Dosazení do centrální formule

splňují všechny

kvadratické polynomy

Dosazení do upwind prvního řádu

22

Falešná difuze

Dosazení do upwind 2.řádu

Nabízí se možnost kombinovat predikci centrální formule a upwindu – kompenzace třetích derivací rozvoje.

Dosazení do formule QUICK

Doporučení

QUICK pro rovnice hybnostihybridní schemata pro k,,c.

23