Dyskusja wyników badań
7.2 Zastosowanie metody KRSR podczas analizy procesów zależnych od czasu
Procedura przeprowadzania dyskretyzacji problemu przy użyciu metody KRSR jest bardzo podobna do tej znanej dla metody różnic skończonych. Polega ona na intuicyjnie naturalnym zastępowaniu pochodnych funkcji w RRCz przez wyrażenia kwadraturowe (§3.3). Przekształcenie takie prowadzi do uzyskania układu wyrażeń algebraicznych, które polega na wykonaniu skończonej liczby mnożeń i sumowań. Operacje takie można bardzo łatwo wykonać na maszynie cyfrowej. Zastosowanie omawianej metody jest dużo prostsze niż na przykład metody elementów skończonych. Podobnie jak w przypadku metody RS, także w tym wypadku istnieją bardziej wyszukane sposoby aproksymacji warunków brzegowych, są one analogiczne do tych stosowanych w przypadku metody KR [1, 65]. W rozprawie tej jednak została przedstawiona najprostsza metoda wprowadzania warunków brzegowych. Mimo tego wynikiem obliczeń jest rozwiązanie numeryczne wysokiej dokładności.
Jako problem posiadający rozwiązanie dokładne wybrano problem stygnięcia nieskończonej płyty staliwnej. Porównanie wyników symulacji wykonanej przy użyciu metody KRSR pozwoliło zweryfikować dokładność rozwiązania uzyskanego dla problemu dynamicznego przy użyciu omawianej metody. Równania różniczkowe cząstkowe opisujące przewodzenie ciepła podczas stygnięcia nieskończonej płyty staliwnej zostały przekształcone do postaci dyskretnej przy użyciu metody KRSR (5.4) oraz RS (5.5). Układy równań opisujące przybliżony przebieg zmian temperatury w układzie zostały rozwiązane przy użyciu komputera. W obu przypadkach czas dyskretyzowany był jawnym schematem różnicowym Eulera, w przypadku obu metod zastosowano ten sam krok czasowy. Metoda RS jest powszechnie stosowana podczas numerycznego rozwiązywania równań różniczkowych cząstkowych. Stanowi podstawę wielu programów komercyjnych, w tym znanego systemu do numerycznego obliczania procesów odlewniczych: MAGMASoft®. Jak jednak widać (Rys. 5.1) rozwiązanie prezentowanego modelu znacząco odchyla się od wartości dokładnych. Związane jest to z własnością metody RS, która polega na tym, że stosunek długości kroku czasowego do przestrzennego ma znaczny wpływ na dokładność schematu. Zmiana ta jest większa w przypadku gęstszej siatki. Inaczej jest w przypadku rozwiązania uzyskanego metodą KRSR. Nie jest możliwe ogólne stwierdzenie dotyczące pogarszania się jakości
Strona 104 rozwiązania wraz ze wzrostem liczby węzłów siatki. Takie zachowanie jest obserwowane po 250 s w przypadku metody o parametrach (12; 16), jednak dla metody o parametrach (5; 9) obserwuje się zachowanie odwrotne - dokładność metody jest wyższa dla gęstszych siatek. Potwierdza to omówiony wcześniej wpływ doboru parametrów metody KRSR do gęstości wybranej siatki, dla której przeprowadzane są obliczenia. Wahania dokładności aproksymacji rozwiązania dla różnych siatek są nieznaczne. Wzrost błędu w trakcie obliczeń jest jednak ciągle obecny (Rys. 5.2). Podczas rozwiązywania problemów dynamicznych niedokładności, które popełniane były podczas aproksymacji pochodnej, w każdym węźle siatki, kumulują się zwielokrotnione dodatkowo poprzez błędy związane z zastosowaniem schematu dyskretyzacji czasu (§3.3.1). Jakkolwiek, wraz ze wzrostem liczby węzłów siatki wartość błędu popełnianego podczas dyskretyzacji metodą KRSR w centrum dziedziny wzrasta (Rys. 4.6) niedokładność w węzłach brzegowych maleje. Obserwowany wzrost niedokładności w centrum dziedziny, choć dotyczy licznego zbioru węzłów, jest nieznaczny w stosunku do poprawy jakości rozwiązania w węzłach brzegowych i bezpośrednio sąsiadujących z brzegiem. Ta właściwość schematów opartych na metodzie KRSR o zdefiniowanym, jak w tej rozprawie rozkładzie rzędów może chronić rozwiązanie przybliżone przed tak szybkim wzrostem wartości błędu jak jest to obserwowane w metodzie RS. Stosunkowo duża wartość niedokładności w pojedynczych węzłach znajdujących się na brzegu, może mieć dużo większy wpływ na niedokładność obliczeń rozwiązania procesów dynamicznych, niż na średnią wartość błędu w całej dziedzinie będącą wskaźnikiem jakości aproksymacji podczas analizy przybliżenia drugiej pochodnej. Wielkość o jaką rozwiązanie przybliżone różni się od dokładnego w punkcie brzegowym kumuluje się i zaburza wartości w węzłach sąsiadujących. Opisana propagacja błędu spowodowana jest tym, że w rozwiązaniach numerycznych w każdym kolejnym kroku czasowym dokonuje się przybliżone całkowanie równania o wartościach początkowych (czasem również brzegowych) wyznaczonych w sposób przybliżony w poprzednim kroku
(
krokach w przypadku schematów różnicowych wielokrokowych [77])
obliczeń numerycznych. Funkcja, którą opisują RRCz oraz jej przebieg w znacznym stopniu zależą od warunków początkowo brzegowych [3, 17, 76, 77]. Jeżeli wartości początkowe są już znacznie oddalone od wartości, które dałoby rozwiązanie dokładne, problem, który jest rozwiązywany będzie różnił się od rozwiązania modelu matematycznego. Dlatego też, rozwiązanie numeryczne, nawet jeżeli uzyskane metodą o wysokiej dokładności, będzie znacznie różniło się od rozwiązania dokładnego analizowanego problemu, gdyż dla błędnych warunkach początkowych będzie już rozwiązaniem całkiem innego problemu (Rys. 7.1). Taki wpływ może mieć niedokładność kumulująca się na brzegu. Modyfikuje ona analizowany problem do innego problemu, którego rozwiązanie oddala się od poszukiwanego, nawet pomimo aproksymacji wysokiej dokładności w kolejnych węzłach bliższych centrum.Strona 105
Rys. 7.1. Zmiana przebiegu rozwiązania równania różniczkowego zwyczajnego y x
dx dy
exp w wyniku nieznacznych zaburzeń warunku początkowego. Niebieskie strzałki, to wektory u styczne do krzywej będącej rozwiązaniem szczególnym tego równania różniczkowego wpunkcie (x0, y0) u=x1,expy0x0, gdzie x jest długością składowej wektora na osi argumentów. Krzywe zaznaczone na wykresie są wybranymi rozwiązaniami szczególnymi dla nieznacznie różniących się warunków początkowych. Wydruk za programem MAXIMA [112].
W świetle wcześniejszych spostrzeżeń, możliwe jest stwierdzenie, iż rozbieżności pomiędzy dokładną wartością rozwiązywanego modelu matematycznego, a rozwiązaniem numerycznym mogą być zależne od chwilowego charakteru rozwiązania. Analiza taka została przeprowadzona na podstawie danych zebranych podczas eksperymentu numerycznego i zgromadzonych na rysunku Rys. 5.3. Zestawienie szybkości stygnięcia, aktualnej wartości temperatury oraz szybkości narastania niedokładności potwierdza to stwierdzenie. Przebieg krzywej opisującej zmianę niedokładności metody KRSR związany jest z aktualnym przebiegiem krzywej obrazującej wartości gradientu temperatury w pobliżu analizowanego punktu. Przyspieszenie narastania niedokładności metody w ciągu pierwszych 50 s od rozpoczęcia obliczeń może być spowodowane błędami zaokrągleń i operacji arytmetycznych powstałych w trakcie realizacji formuł opisujących warunek symetrii (3.26), a także przesuwaniem się pierwszych efektów spadku temperatury od brzegu na którym zadany jest warunek trzeciego rodzaju (3.31). Ze względu na znaczną różnicę pomiędzy temperaturą na brzegu analizowanego układu, a temperaturą otoczenia, błąd popełniany w pierwszych krokach obliczeń na brzegu stycznym z otoczeniem będzie większy. Spowodowane jest to faktem, że w momencie rozpoczęcia obliczeń założony został jednorodny rozkład temperatury. Spadek temperatury w jednym, dwóch tylko węzłach brzegowych, powoduje
Strona 106 znaczne wygięcia wielomianów, opisujących w trakcie aproksymacji przebieg temperatury. W takim przypadku zmiany wypukłości funkcji interpolacyjnej wpływają na znak parzystych pochodnych wielomianu i powodują pojawienie się dodatkowych odchyleń, w punktach, w których temperatura jeszcze nie powinna się zmieniać zgodnie z modelem matematycznym. Obserwowany wzrost gradientu po 125 s wywołuje spadek przyspieszenia wartości błędu obliczeń numerycznych. Może być to związane z tym, że na brzegu stycznym z otoczeniem temperatura już spadła, a jej rozkład jest gładki w całej dziedzinie. Kolejnym powodem jest stabilizacja zachowania temperatury w układzie. Zaczyna ona spadać w każdym węźle, a szybkość z jaką wzrasta gradient jest funkcją gładką. Przypadkowe zaburzenia wartości temperatury spowodowane błędami aproksymacji są znikome w stosunku do zmian wywołanych stygnięciem. Omówiony powyżej proces stabilizacji przebiegu ewolucji problemu fizycznego i jego wpływu na charakter aproksymacji widoczny jest podczas dalszych obliczeń. Po 175 s szybkość zmiany niedokładności metody KRSR spada poniżej 4,5 10-8 1/s, a po 300 s poniżej 4 10-8 1/s. Wiąże się to z utrzymywaniem się podczas obliczeń tych samych czynników, które opisane były powyżej. Stabilizacja przyrostu błędu na bardzo niskim poziomie świadczy o tym, że metoda KRSR może znaleźć zastosowanie podczas symulacji procesów przewodzenia ciepła. Niskie wartości błędów w kolejnych krokach czasowych wpływają na wysoką dokładność obliczeń metodą kwadratur różniczkowych sterowanego rzędu w stosunku do metody różnic skończonych.
Metoda KRSR jest metodą dyskretyzacji pochodnych przestrzennych występujących w równaniach różniczkowych cząstkowych. Występujące w procesach związanych z odlewnictwem RRCz często są równaniami nieliniowymi. Oprócz dyskretyzacji pochodnych względem zmiennej opisującej zmianę czasu oraz pochodnych opisujących zmianę pola poszukiwanej wartości względem zmiennej przestrzennej konieczna jest aproksymacja członu nieliniowego (1.2). Czasem na drodze odpowiednich operacji matematycznych możliwe jest przekształcenie równania FK zawierającego funkcję źródła do równania liniowego. Tak jest w przypadku zastosowania metody zastępczej pojemności cieplnej strefy dwufazowej (§1.1.5). Zastosowanie jej w przypadku modelowania krzepnięcia odlewu sprowadziłoby rozwiązywany problem, do zagadnienia, które można zapisać w postaci dyskretnej z użyciem metody KRSR w sposób analogiczny do problemu prezentowanego w paragrafie §5.1. Podejście takie podczas krzepnięcia stopów w przedziale temperatury uniemożliwia kontrolę ilości ciepła wydzielonego podczas tego procesu. Wielu badaczy stosuje jednak metodę zastępczej pojemności cieplnej podczas symulacji procesów krzepnięcia osiągając dużą zgodność z wynikami pomiarów [3, 7, 20, 113]. Metoda ta jest jednak nieefektywna podczas modelowania stopów krzepnących w relatywnie dużych przedziałach temperatury. Z tego względu w niniejszej pracy została zastosowana odmienna metoda modelowania krzepnięcia, bazująca na metodach spotykanych w pracach naukowych związanych z omawianą tematyką (§6.2.2). Prezentowany w niniejszej rozprawie model procesów krzepnięcia wymaga znajomości funkcji udziału objętości zakrzepłej
fS w zależności od aktualnej temperatury w dyskretnym elemencie, odwzorowującym podobszar odlewu. Funkcja ta została wyznaczona doświadczalnie na podstawie krzywych stygnięcia odlewu płyty wykonanej ze stopu AZ91 (Tabela 5).
Strona 107 Analiza danych doświadczalnych w celu określenia pewnych parametrów, których znalezienie w inny sposób przysparza problemów jest często stosowana w praktyce inżynierskiej. Krzywa udziału objętości zakrzepłej zawiera w swoim przebiegu informację o kinetyce krzepnięcia (Rys. 6.6, Rys. 6.7). W pierwszym etapie zaznacza się wzrost szybkości krzepnięcia w wyniku zwiększania się udziału objętości fazy pierwotnej -Mg, w drugim krzywa niesie informację o zmianie kinetyki krzepnięcia, gdy pojawia się przemiana eutektyczna. Wyznaczona w trakcie eksperymentu ilość ciepła przemiany jest bliska wartości określanej przez innych badaczy [83, 103-105]. Stop AZ91 jest określony przez pewne graniczne udziały pierwiastków stopowych [114, 115], konkretne stopy mogą zatem różnić się składem, co za tym idzie parametrami termofizycznymi. Teoretyczna analiza parametrów termofizycznych przeprowadzona przez Kaufmanna i wsp. [106] przy użyciu programu Thermo-Calc wykazuje znaczne wahania wartości parametrów termofizycznych. Dlatego przeprowadzenie eksperymentów dla konkretnego stopu powinno prowadzić do podniesienia dokładności obliczeń numerycznych.
Zaproponowany model numeryczny (§6.2.3) pozwolił przygotować program komputerowy realizujący symulację procesu krzepnięcia nieskończonej płyty wykonanej ze stopu AZ91 (Tabela 5) w formie piaskowej (Dodatek A). W celu weryfikacji obliczeń wykonany został odlew płyty nieskończonej, której krzepnięcie odpowiada procesowi opisanemu przez model matematyczny (§6.1). Wyniki symulacji wykazują dużą zgodność z rezultatami eksperymentu (Rys. 6.8). Kształty krzywych stygnięcia są bardzo zbliżone. Świadczy to dobrym odwzorowaniu proporcji ciepła wydzielanego podczas krzepnięcia i przepływającego przez media. Zarówno przez płytę, jak do formy i dalej, przez otaczającą ją formę piaskową do otoczenia. Analogie przebiegu krzywych obserwowalne są również podczas porównań szybkości stygnięcia (Rys. 6.9), szybkości z jaką wydziela się ciepło przemiany (Rys. 6.11) oraz zmiany objętości zakrzepłej (Rys. 6.12) względem temperatury. Na podstawie wspomnianych wykresów można stwierdzić dużą zgodność wyznaczonego schematem, opartym na metodzie KRSR przepływu ciepła ze stanem rzeczywistym. Wymienione powyżej wykresy sporządzone zostały względem argumentu będącego temperaturą, ze względu na różnice między czasem krzepnięcia wyznaczonym numerycznie i zarejestrowanym podczas eksperymentu. Ta różnica wynosi około 100 s (Tabela 6). Czas krzepnięcia wyznaczony numerycznie jest dłuższy. Różnica taka może być związana ze sformułowaniem modelu numerycznego. Podczas symulacji zaniedbany został początkowy rozkład temperatury w odlewie oraz w formie po jej wypełnieniu. Dodatkowo parametry termofizyczne w dużej mierze zostały wzięte z bazy danych programu MAGMASoft® [83]. Zachowanie funkcji je opisujących może różnić się od rzeczywistego przebiegu dla stopu o składzie takim jak użyty podczas eksperymentu (Tabela 5). Założenie stałej temperatury w całej objętości odlewu i formy powoduje pojawienie się sytuacji, w której na brzegu następują znaczne zmiany temperatury (w formie bardzo szybkie nagrzewanie, w odlewnie chłodzenie). Tak znaczny spadek temperatury w punkcie brzegowym prowadzi do wypaczenia charakteru funkcji opisującej rozkład temperatury w kierunku centrum odlewu. Taka sytuacja ma miejsce na początku analizowanego procesu. Może ona być powodem przyspieszenia wzrostu niedokładności podczas pierwszych sekund obliczeń o czym było już
Strona 108 mówione podczas opisu i interpretacji wyników symulacji stygnięcia nieskończonej płyty staliwnej (Rys. 5.3). Dodatkowo, analogicznie jak w przypadku warunku brzegowego III rodzaju, podczas kilku pierwszych kroków czasowych te znaczne różnice między wartościami temperatury w kilku węzłach sąsiadujących z brzegiem pomiędzy odlewem, a formą mogą zaburzać płynną zmianę temperatury w węźle brzegowym, wyznaczoną z warunku brzegowego IV rodzaju (3.35).
Ważnym aspektem podczas modelowania krzepnięcia jest zachowanie ilości ciepła przemiany, która wydzieli się podczas symulacji. Obliczenia takie zostały przeprowadzone dla krzywej stygnięcia wyznaczonej numerycznie z zastosowaniem metody KRSR (Rys. 6.10). Wyznaczona wartość utajonego ciepła przemiany jest bardzo bliska wartości zadanej w modelu numerycznym (Tabela 7). Różnice między tymi dwiema wartościami są na tyle małe, że można je przypisać metodzie numerycznego całkowania pola powierzchni pod krzywą. Określona wzorem (6.24) operacja całkowania pola powierzchni pomiędzy teoretyczną krzywą zerową, a rzeczywistą szybkością stygnięcia pomnożonego przez wartość ciepła właściwego daje wartość przybliżoną. Dokładność aproksymacji w przypadku danych eksperymentalnych zależy od długości kroku próbkowania oraz dokładności z jaką rejestrowana jest temperatura. Na krzywej stygnięcia uzyskanej podczas eksperymentu często obserwowane są różnego rodzaju zaburzenia, które mogą zaburzać dokładność z jaką określona jest ilość ciepła wydzielonego podczas przemiany. W przypadku obliczeń numerycznych dokładność z jaką wyznaczonej jest L1 zależy od długości kroku całkowania. W omawianym przypadku kroki całkowania były różne w przypadku analizy danych eksperymentalnych oraz wyników obliczeń numerycznych, co ze względu na charakter przybliżonej metody całkowania musi wpływać bezpośrednio na wyznaczone wartości [18, 25]. Porównanie obliczonych ilości ciepła przemiany potwierdza zatem, że zaproponowany model uwzględnienia przemian mających miejsce podczas krzepnięcia stopu AZ91 prowadzi do uzyskania schematu dającego rozwiązanie wysokiej dokładności, a także pozwala na precyzyjną kontrolę ilości ciepła, jaka wydziela się podczas krzepnięcia.
Oparty na aproksymacji pochodnych przestrzennych występujących w RRCz metodą kwadratur różniczkowych sterowanego rzędu schemat numeryczny umożliwił wyznaczenie przebiegu zmian temperatury bardzo zbliżonego do zarejestrowanego w trakcie eksperymentu. W przypadku tym model matematyczny charakteryzuje się obecnością dużo bardziej skomplikowanych równań różniczkowych cząstkowych. Nieliniowość występująca w równaniu FK może skutecznie uniemożliwić przeprowadzenie obliczeń numerycznych. Przeprowadzony eksperyment numeryczny wykazał, że metoda KRSR może zostać zastosowana podczas symulacji procesów, do opisu których stosuje się nieliniowe RRCz. Ponadto rozważane dwie poddziedziny, z których jedna reprezentuje obszar odlewu, druga obszar formy utrudniają przeprowadzenie obliczeń numerycznych. Szczególnie w przypadku dużych różnic pomiędzy temperaturami początkowymi obu poddziedzin. W tym wypadku zastosowanie metody KRSR nie przysparza trudności, a przybliżone rozwiązanie wiernie oddaje przebieg analizowanego procesu. Wyniki omówionych powyżej analiz pozwalają
Strona 109 postawić metodę KRSR pośród metod numerycznych, które należy stosować do problemów występujących w odlewnictwie.