1 Zakład Mechaniki i Fizyki PłynówTomasz Michałek Metoda oceny wiarygodności symulacji numerycznych przepływów lepkich i termicznych Instytut Podstawowych Problemów Techniki Polska Akademia Nauk Zakład Mechaniki i Fizyki Płynów

2 Cel pracy: Opracowanie metody oceny wiarygodności symulacji numerycznych dla przepływów lepkich i termicznych. Zastosowano analizę wrażliwości rozwiązań numerycznych do określenia istotności parametrów fizycznych przy projektowaniu wzorców eksperymentalnych. Wzorce eksperymentalne są niezbędne w procesie potwierdzenia słuszności przyjętych założeń modelowych w symulacjach numerycznych (proces walidacji). Projektowanie wzorców eksperymentalnych: 1. Wyznaczenie kluczowych (istotnych) parametrów dla konfiguracji eksperymentalnych przy określaniu poprawności fizycznej modelu (walidacji). 2. Określenie dokładności pomiarów eksperymentalnych stanowiących poziom wiarygodności w procesie walidacji.

3 Równania Naviera - Stokes’aPrzedmiot badań: Równania Naviera - Stokes’a sprzężone z równaniem transportu ciepła Konwekcja naturalna Konwekcja wymuszona Warunki brzegowe: Warunki początkowe:

4 Wiarygodność symulacji numerycznychWeryfikacja Walidacja Weryfikacja programu (kodu) Weryfikacja obliczeń Konfiguracje laboratoryjne (testowe) Rzeczywiste konfiguracje Weryfikacja aparatu numerycznego Walidacja modelu wzorce numeryczne pomiary eksperymentalne wzorce eksperymentalne obliczenia numeryczne Analiza wrażliwości wzorce eksperymentalne

5 Wiarygodność symulacji numerycznychWeryfikacja Walidacja Weryfikacja programu (kodu) Weryfikacja obliczeń Konfiguracje laboratoryjne (testowe) Rzeczywiste konfiguracje metoda sztucznych rozwiązań (ang. method of manufactured solution) [Roache] rozwiązania analityczne wzorce numeryczne [Ghia, de Vahl Davis, Le Quere,…] Ekstrapolacja Richardsona (1904) Uogólniona ekstrapolacja Richardsona [Stern at all.] Współczynnik zbieżności na siatce (GCI) [Roache] Problemy jednostkowe (unit problems) Wzorce eksperymentalne (benchmark cases) Uproszczone konf. (simplified/partial flow path) [Sindir et al.] Analiza wrażliwości

6 ANALIZA WRAŻLIWOŚCI ze względu na parametry wejścioweRozw. numeryczne PARAMETRY WEJŚCIOWE: Warunki brzegowe Warunki początkowe Własności materiałowe MODEL WYNIK: Określenie wrażliwości funkcjonałów na zmiany parametrów wejściowych Zestaw kluczowych (istotnych) parametrów 2. Określenie niezbędnej dokładności pomiarów

7 Symulacje numeryczne: Pomiary eksperymentalne: Weryfikacja rozwiązań:Metodologia badań: Symulacje numeryczne: Pomiary eksperymentalne: Metoda różnic skończonych (SOLVSTR) Metoda bezsiatkowa (SOLVMEF) Metoda objętości skończonych (FLUENT) Pole prędkości (cyfrowa anemometria obrazowa) Pole temperatury (cyfrowa termometria obrazowa) Pomiary punktowe temperatury (termopary) Pomiary właściwości termo-fizycznych Weryfikacja rozwiązań: Ekstrapolacja Richardsona GCI Określanie rzędu zbieżności rozwiązań

8 Weryfikacja obliczeń numerycznychRe = 100, Przepływ w zagłębieniu (ang. moving lid driven cavity) Test na niezależność rozw. od siatki obl. Siatka min V2 GCI p  33x33 -0,222240 65x65 -0,241399 8 % 2.5% 129x129 -0,246071 2 % 2,03 0.6 % 257x257 -0,247051 0.3 % 2,25 0.13 % 513x513 -0,247027 0.01% 5,35 0.003% Ra = 105, Konwekcja naturalna (poziomy gradient temperatury) Th Test na niezależność rozw. od siatki obl. Siatka Max V2 GCI p  25x25 50x50 9% 2.9% 100x100 0.6% 3.94 0.2% Tc

9 Wzorzec zdefiniowany przez Ghia et all. [1980]Re = 100 Siatka min V2 GCI p  33x33 -0,222240 65x65 -0,241399 8 % 2.5% 129x129 -0,246071 2 % 2,03 0.6 % 257x257 -0,247051 0.3 % 2,25 0.13 % 513x513 -0,247027 0.01% 5,35 0.003% Re = 5000 Siatka Min V2 GCI p  33x33 -0,320896 65x65 -0,613754 47% 13% 129x129 -0,477550 28% 1,10 10% 257x257 -0,551242 0,88 4%

10 Wzorzec zdefiniowany przez G. de Vahl Davisa [1983]Ra = 103 Ra = 104 Ra = 105

11 Rozwiązanie wzorcowe dla wysokich liczb Rayleigha Le Quere [1991]

12 WZORZEC NUMERYCZNY do weryfikacji obliczeń z nieliniowym członem objętościowym2D lepki, nieściśliwy przepływ wody Równania Naviera – Stokesa sprzężone z równaniem transportu ciepła poprzez nieliniowy człon sił masowych opisujący zależność gęstości wody od temperatury Różnica temperatur ΔT = 10ºC Programy: Th = 10C Tc = 0C SOLVSTR (FDM) SOLVMEF (MEF) FLUENT (FVM) FIDAP (FEM) FRECON(FDM) Ra = 1.5 · Pr = 13.31

13 PROCEDURA WERYFIKACJIV Rozwiązanie wzorcowe Michalek T., Kowalewski T.A., Sarler B. ”Natural Convection for Anomalous Density Variation of Water: Numerical Benchmark” Progress in Computational Fluid Dynamics, 5 (3-5),pp ,2005 FRECON3V (FRE) FLUENT 6.1. (FLU) FIDAP (FID) SOLVSTR (STR) Oszacowanie dokładności rozwiązań poprzez porównywanie profili (nie wartości punktowych)

14 Metody eksperymentalne:Pomiar pól prędkości : Cyfrowa anemometria obrazowa (PIV) Pomiar pól temperatury: Cyfrowa termometria obrazowa (PIT) Wizualizacja pól prędkości Punktowe pomiary temperatury F(t0) Korelacja F(t0+t)

15 Układ pomiarowy: TH > TC light sheet

16 Naczynie pomiarowe: + pomiary kontrolnePrzekrój centralny TE1 TE2 ALUMINIUM WALL PLEXIGLASS WALL T7 T10 T14 T15 Th TL TP Tc

17 WYNIKI POMIARÓW EKSPERYMENTALNYCH zastosowanie dwóch zawiesin ciekłokrystalicznychPIT - Pomiary temperatury Ra = 1.5*106 Pr = 11.78 PIV – Pomiary prędkości Th = 10 C Tc = 0 C

18 WZORZEC EKSPERYMETALNY Ra = 1. 5x106 Pr = 11WZORZEC EKSPERYMETALNY Ra = 1.5x106 Pr = Profile temperatury i prędkości 2D Pole temperatury 2D Pole prędkości TH = 10 C Tc = 0 C T wzdłuż Y = 0.5L V wzdłuż Y = 0.5L V wzdłuż X = 0.9L

19 Oszacowanie niepewności pomiarów UD:PIV pomiar pól prędkości N = długość serii UD : Mix C Zakres temp. [C] Hue Kolor UD[C] 0.0 3.0 0.11 0.18 Czerwony 1.0 3.5 0.25 Żółty 0.5 3.9 0.48 Zielony 8.0 0.66 Niebieski BM100 5.5 6.4 0.12 0.28 6.5 0.35 7.5 0.55 9.5 0.70 1.5 PIT Pomiar pól temperatury dwie zawiesiny ciekłokrystaliczne

20 PRZYKŁADY WALIDACJI SYM. NUM.:PROCEDURA WALIDACJI Błąd porównania Metryka walidacji PRZYKŁADY WALIDACJI SYM. NUM.: Symulacja numeryczna przepływu konwekcyjnego wody (poziomy gradient temperatury), Ra = 1.5x106 Pr = 11.78 Symulacja numeryczna przepływu konwekcyjnego wody (poziomy gradient temperatury), Ra = 3x107 Pr = 9.53 Symulacja numeryczna przepływu konwekcyjnego glikolu polietylenowego (pionowy gradient temperatury) Ra~105 Pr~103

21 ANALIZA WRAŻLIWOŚCI ze względu na parametry wejścioweRozw. numeryczne PARAMETRY WEJŚCIOWE: Warunki brzegowe TH, TC, Text, Q1, Q2, Q3 Warunki początkowe Tinit. ,vinit Własności materiałowe ,,,,cp MODEL WYNIK: Określenie wrażliwości funkcjonałów na zmiany parametrów wejściowych Zestaw kluczowych (istotnych) parametrów dla rozważanej konfiguracji 2. Określenie niezbędnej dokładności pomiarów w celu przeprowadzenia walidacji obl. num.

22 ZASTOSOWANIE ANALIZY WRAŻLIWOŚCIParametry wejściowe: TH, TC, Text, 1, 2, 3 T0 , v0 , , cp, ,  (wody) , cp, , (plexi) , cp, , (aluminium) MODEL Zakres zmian parametrów: dla temp. TH , TC, 1 C , Text 2 C własności materiałowe (na post. tablic) dla wsp. przejmowania ciepła na podstawie pomiarów Funkcjonały:

23 PRZYKŁAD 1 : Wyniki analizy wrażliwościIstotne parametry: temperatura Tc wsp. przejmowania ciepła  wsp. lepkości 

24 PRZYKŁAD 1: Ra = 1.5*106 Pr ~ 10 Symulacja A (T),(T),cp(T)wł. mat. zależne od temp. (T),(T),cp(T) Symulacja B stałe wartości wł. mat. ,,cp = const. Symulacja C poziome ścianki adiabatyczne pionowe izotermiczne, stałe wł. mat. Pole temperatury Pole prędkości

25 PRZYKŁAD3: WALIDACJA ILOŚCIOWE PORÓWNANIETemp. wzdłuż Y = 0.5L V wzdłuż Y = 0.5L

26 Pomiary eksperymentalne: Ra ~ 3.0x107Pr 1 3*107 9.53 2 1.5 *108 7.01 3 1.8*108 4 4.4*108 5.41 Th = 18.0 C Tc = 4.0 C PIV Th = 23.2 C Tc = 9.0 C

27 Dodatkowe pomiary eksperymentalneRa = 3.107 Ra =

28 PRZYKŁAD 2 : Wyniki analizy wrażliwościIstotne parametry: temperatura Tc lepkość  Zmiana współczynnika przejmowania ciepła dla pomiarów temp. w górnej i dolnej ścianie naczynia

29 Pomiary eksperymentalne: Ra = 1.5x108Pr 1 3*107 9.53 2 1.5 *108 7.01 3 1.8*108 4 4.4*108 5.41 Th = 27.3 C Tc = 6.8 C PIV PIT with two TLCs Th = 27.2 C Tc = 6.8 C

30 Pomiary eksperymentalne: Ra = 1.8x108Pr 1 3*107 9.53 2 1.5 *108 7.01 3 1.8*108 4 4.4*108 5.41 Th = 36.4 C Tc = 10.2 C PIV PIT with two TLCs Tc = 10.2 C Th = 36.4 C

31 Pomiary eksperymentalne: Ra = 4.4x108Pr 1 3*107 9.53 2 1.5 *108 7.01 3 1.8*108 4 4.4*108 5.41 Th = 45.8 C Tc = 14.2 C PIV PIT with two TLCs Tc = 14.0 C Th = 45.8 C

32 ZASTOSOWANIE ANALIZY WRAŻLIWOŚCIKluczowe parametry dla badanej konfiguracji eksperymentalnej Zmiana TC o  1 C powoduje zmianę wartości funkcjonału || . ||3 o 46 %, || . ||2 o 19% , || . ||1 o 24 % Zmiana TH o  1 C powoduje zmianę wartości funkcjonału || . ||1 o 36 %, || . ||2 o 12 % , || . ||3 o 11 % Zmiana 1 o  10 Wm-2K-1 powoduje zmianę wartości funkcjonału || . ||1 o 22 % , || . ||6 o 10 % , || . ||7 o 12 % Zmiana lepkości  o  10 % powoduje zmianę wartości funkcjonału || . ||2 o 7 %, || . ||3 o 9 % Zmiana 2 i 3 o  100 Wm-2K-1 powoduje zmianę wartości funkcjonałów || . ||1 i || . ||3 o 4 % i 5 %. … Oszacowanie dokładności pomiarów na potrzeby walidacji Parametr Th 27.8% 11.4 % 11.2 % 8.7 % - 5.9 % 6.2 % 4.5 % Tc 21.2% 17.6 % 12.0 % 9.4 % Text 11.8% 3.4 % 3.0 % 3.1 % 3.2 % 7.7 % 1 11.6% 3.9 % 6.0 % 6.9 % 24.2 % 2 3.20% 3 5.60%  7.2 %  4.60% 10.3 % cp 11.8 % 10.7 %  5.40% 13.6 %

33 Symulacja numeryczna (SOLVSTR)PRZYKŁAD 2 : Ra ~ 3 x 107 Eksperyment Eksperyment Symulacja numeryczna (SOLVSTR) Funkcjonał D UD S USN E UV T7 18,22 0,48 17,99 0,07 0,23 0,49 T10 17,76 0,63 17,17 0,59 Umin -0,66 0,24 -0,65 0,01 Umax 0,69 0,65 0,04 Vmin -2,60 -2,40 0,09 0,20 0,26 Vmax 2,42 2,40 0,02 VP1 -2,48 0,58 -1,99 VP2 -1,85 0,42 -1,71 0,14 UP3 -0,24 -0,22 VP3 -0,75 0,21 -1,05 0,30 UP4 -0,58 -0,39 0,19 UP5 -0,60 0,16 -0,42 0,18 FD (SOLVSTR) Warunek nie jest spełniony

34 Sym. Numeryczna (FLUENT)PRZYKŁAD 2 : Ra ~ 3 x 107 Eksperyment Eksperyment Sym. Numeryczna (FLUENT) Funkcjonał D UD S USN E UV T12 18,67 0,38 18,92 0,02 0,25 T16 4,05 3,83 0,22 T7 18,22 0,48 18,39 0,17 T10 17,76 0,63 17,64 0,12 Umin -0,66 0,24 -0,73 0,01 0,07 Umax 0,69 0,68 Vmin -2,6 -2,22 0,05 Vmax 2,42 2,22 0,20 VP1 -2,48 0,58 -1,99 0,49 VP2 -1,85 0,42 -1,77 0,08 UP3 -0,24 0,09 -0,29 VP3 -0,75 0,21 -1,29 0,54 UP4 -0,58 0,14 -0,4 0,18 UP5 -0,6 0,16 -0,42 FV (Fluent) Warunek nie jest spełniony

35 PRZYKŁAD 3 : Analiza wrażliwościTc Parametry wejściowe: TC, Text,  T0 , v0 , , cp, ,  (PEG) T4 T3 Text Zakres zmian parametrów: dla temp. TC 1 C , Text 1 C własności materiałowe (na post. tablic) dla wsp. przejmowania ciepła na podstawie pomiarów T2 Funkcjonały: Kowalewski T.A., Cybulski A., Michałek T., Kowalczyk M. „Laboratoryjne wzorce do walidacji programów odlewniczych” Prace IPPT, 2005,

36 PRZYKŁAD 3 : Wyniki analizy wrażliwości[Wm-2K] Pole pionowej składowej prędkości [C] Eksperyment Sym. Num.

37 PRZYKŁAD 3 : Procedura walidacjiSym. Num. (  = const.) Eksperyment Symulacja Numeryczna Zmienna D UD S USN E UV T1 [K] 307,2 0,5 307,6 0,01 0,40 0,50 T2 [K] 306,8 306,5 0,02 0,30 T3 [K] 306,0 0,80 T4 [K] 305,5 2,0 302,9 0,20 2,60 2,01 ||F||5 = Wmin [mm/s] -0,965 0,1 -1,754 0,79 0,10 ||F||6 = Wmax [mm/s] 0,140 0,317 0,18 Sym. Num. (  = (T) ) Eksperyment Symulacja Numeryczna Zmienna D UD S USN E UV T1 [K] 307,2 0,5 307,3 0,01 0,10 0,50 T2 [K] 306,8 306,5 0,02 0,30 T3 [K] 306,3 T4 [K] 305,5 2,0 303,5 0,20 2,00 2,01 ||F||5 = Wmin [mm/s] -0,965 0,1 -0,963 0,00 ||F||6 = Wmax [mm/s] 0,140 0,125

38 PODSUMOWANIE I WNIOSKIZaproponowano metodę oceny wiarygodności symulacji numerycznych opartą na badaniu wrażliwości rozwiązań numerycznych Zdefiniowano wzorzec numeryczny dla przepływów konwekcyjnych z nieliniowym członem wypornościowym (Ra = 1.5x106 , Pr = 11.78) Wykorzystano zaproponowaną metodę do zaprojektowania referencyjnego eksperymentu Zdefiniowano wzorzec eksperymentalny w oparciu o przeprowadzone pomiary doświadczalne Przeprowadzono walidację symulacji numerycznych w oparciu o wyniki eksperymentalne z wykorzystaniem zaproponowanej metody