Artykuł metodologiczny

Symulacja przepływu kondensacyjnego i wymiany ciepła w spiralnych wymiennikach ciepła dla nieazeotropowych mieszanin węglowodorów

51 wyświetleń

⸱

DOI:

10.3791/71595

⸱

8 września 2026

W tym artykule

Podsumowanie

W niniejszym artykule przedstawiono protokół numerycznego symulowania wymiany ciepła podczas kondensacji oraz charakterystyki przepływu nieazeotropowych mieszanin węglowodorów w wymiennikach ciepła o spiralnym zwoju. Metoda ta pozwala na ocenę warunków operacyjnych i roboczych w celu przewidzenia współczynników wymiany ciepła oraz spadków ciśnienia.

Streszczenie

Jako kluczowy element procesu skraplania gazu ziemnego, spiralne wymienniki ciepła odgrywają istotną rolę w produkcji LNG. Celem niniejszej pracy było kompleksowe zrozumienie charakterystyki przepływu kondensacyjnego i wymiany ciepła nieazeotropowych mieszanin węglowodorów wewnątrz rurek spiralnych. Przeprowadzone w badaniu wysokoprecyzyjne dochodzenie numeryczne oparto na ściśle określonym schemacie postępowania, który obejmował tworzenie geometrii za pomocą profesjonalnego narzędzia do modelowania, generowanie siatki w dedykowanym oprogramowaniu do siatkowania, obliczenia numeryczne w komercyjnym solverze (z uwzględnieniem monitorowania zbieżności w czasie rzeczywistym) oraz ilościowe przetwarzanie wyników. Zintegrowane podejście to zapewnia wysoką wierność otrzymanego modelu numerycznego. Maksymalne odchylenia od klasycznych danych eksperymentalnych (dane eksperymentalne Neeraasa) wyniosły poniżej 15% dla współczynnika wymiany ciepła i poniżej 10% dla gradientu spadku ciśnienia tarcia. Wyniki symulacji wykazują, że zmiana okresów i amplitud falowania powoduje podobne trendy oscylacyjne w procesie wymiany ciepła, wykazując zarówno efekty usprawnienia, jak i pogorszenia. W szczególności okres falowania zmienia wydajność wymiany ciepła o ±20%, podczas gdy amplituda falowania wpływa na nią w zakresie ±10%.

Wprowadzenie

Gaz ziemny, jako stosunkowo czyste paliwo kopalne, emituje podczas spalania znacznie mniej dwutlenku węgla i innych zanieczyszczeń niż węgiel i ropa naftowa. W globalnej transformacji w kierunku odnawialnych systemów energetycznych gaz ziemny jest często postrzegany jako „paliwo przejściowe” ze względu na jego zdolność do utrzymania stabilności i niezawodności dostaw energii1. Poprzez skraplanie gazowy gaz ziemny jest schładzany do postaci cieczy kriogenicznej (LNG), co zmniejsza jego objętość około 600-krotnie, znacznie ułatwiając transport i magazynowanie2. Spiralny wymiennik ciepła (SWHE) jest kluczowym elementem procesu skraplania gazu ziemnego. Ten typ wymiennika ciepła składa się z serii spiralnych rurek zamocowanych w cylindrycznej płaszczu, nawijanych warstwa po warstwie w przeciwnych kierunkach wokół centralnego trzpienia, przy czym dystanse oddzielają poszczególne warstwy, aby zapewnić odpowiedni prześwit dla wymiany ciepła. Dzięki spiralnej konfiguracji SWHE zapewnia dużą powierzchnię wymiany ciepła przy niewielkiej powierzchni zabudowy3. Taka kompaktowa konstrukcja sprawia, że nadaje się on idealnie do integracji w dużych instalacjach, szczególnie na morskich pływających platformach produkcyjnych, gdzie przestrzeń jest ściśle ograniczona. W procesach skraplania z użyciem mieszanego czynnika chłodniczego, powszechnie stosowanych w produkcji LNG, nieazeotropowe węglowodory przepływają w górę wewnątrz rurek, podczas gdy zimny płyn po stronie płaszcza przepływa w dół w sposób przeciwprądowy przez szczeliny w wiązce rurek. W tych warunkach głównym procesem po stronie rurek jest kondensacja nieazeotropowych węglowodorów wewnątrz spiralnych rurek, co wiąże się ze złożonym przepływem dwufazowym gaz-ciecz4,5.

W celu dokładnego przewidzenia charakterystyki przepływu i wymiany ciepła podczas kondensacji wewnątrz rurek przeprowadzono szeroko zakrojone badania. W przypadku alkanów jednoskładnikowych Fries i wsp.6 zmierzyli charakterystykę wymiany ciepła podczas kondensacji propanu w rurkach poziomych, stwierdzając, że spadek ciśnienia wzrastał wraz ze zmniejszaniem się średnicy rurki i ciśnienia nasycenia. Zauważyli również, że grawitacja powodowała, iż współczynnik wymiany ciepła na dnie rurki był niższy niż na jej górnej ściance. Zhuang i wsp.7,8 badali kondensację metanu i etanu w rurkach poziomych, wykazując, że współczynnik wymiany ciepła i frictional pressure drop wzrastały wraz z natężeniem przepływu i jakością pary. Wcześniejsze badanie9 analizowało proces kondensacji propanu w mikrokanałach, potwierdzając, że trendy wymiany ciepła i spadku ciśnienia były podobne do tych w kanałach konwencjonalnych. W odniesieniu do mieszanek chłodniczych, Smit i wsp.10 badali kondensację mieszanin R22/R142b w rurkach poziomych, stwierdzając, że przy niskich strumieniach masy zwiększenie ułamka masowego R142b znacząco obniżało współczynnik wymiany ciepła. Berrada i wsp.11 badali mieszaninę R134a/R23 i stwierdzili, że temperatura glide miała niewielki wpływ na wymianę ciepła przy różnych stosunkach składników. Neeraas przeprowadził eksperymenty na mieszaninach etanu i propanu w rurkach spiralnych, zauważając, że efekt mieszania znacząco wpływał na obliczenie współczynnika wymiany ciepła podczas kondensacji12. W symulacjach numerycznych Li i wsp.13 symulowali proces kondensacji etanu i propanu, wykazując, że współczynnik wymiany ciepła i frictional pressure drop malały wraz ze wzrostem ciśnienia nasycenia. Qiu i wsp.14 wprowadzili efekt porywania cieczy przez parę podczas symulowania kondensacji propanu w rurkach spiralnych; ich wyniki wykazały, że uwzględnienie tego efektu zredukowało rozbieżność między wynikami symulacji a danymi eksperymentalnymi do poziomu poniżej 25%.

Pomimo obszernych badań nad różnymi cieczami roboczymi i konfiguracjami kanałów, wciąż istnieje znacząca luka między cieczami czystymi lub binarnymi, zazwyczaj badanymi w literaturze, a mieszaninami wieloskładnikowymi wykorzystywanymi w przemysłowej produkcji LNG. W szczególności badania numeryczne ukierunkowane na nieazeotropowe mieszaniny węglowodorowe składające się z trzech lub więcej komponentów w złożonych kanałach przepływowych pozostają niezwykle ograniczone15. Ponadto, w odniesieniu do specyficznego zastosowania platform LNG typu offshore, wciąż brakuje kompleksowego zrozumienia tego, w jaki sposób ruch urządzeń wywołany środowiskiem morskim zmienia przepływ kondensacyjny i charakterystykę wymiany ciepła. Aby wypełnić te luki badawcze, niniejsza praca łączy symulacje obliczeniowej mechaniki płynów (CFD) z istniejącymi danymi eksperymentalnymi w celu opracowania szczegółowego trójwymiarowego modelu dwufazowego przepływu kondensacyjnego. Przeprowadzono porównanie i analizę symulowanego współczynnika wymiany ciepła oraz spadku ciśnienia tarcia w oparciu o klasyczne dane eksperymentalne. Na podstawie tego modelu w pracy skupiono się na symulacji procesu kondensacji reprezentatywnych składów z pól gazowych w rurach spiralnych, szczególnie w warunkach przechyłu, aby zbadać podstawowe mechanizmy wpływu złożonego ruchu na wymianę ciepła podczas kondensacji wieloskładnikowej. Badanie to dostarcza wiarygodnych podstaw teoretycznych i wytycznych inżynieryjnych dla projektowania i optymalizacji wydajnych wymienników ciepła w procesach skraplania gazu ziemnego na platformach offshore.

Protokół

Ponieważ niniejsza praca koncentruje się na lokalnych charakterystykach wymiany ciepła i spadku ciśnienia podczas kondensacji wewnątrz rury helikalnej, po pełnym rozwinięciu przepływu można zastosować zredukowaną domenę, co umożliwia dokładne odwzorowanie lokalnego zachowania przepływu i właściwości termicznych. W celu walidacji z danymi eksperymentalnymi skonstruowano model rury helikalnej składający się z trzech sekcji, oparty na modelu fizycznym zaproponowanym przez Neeraasa12, z średnicą rury 14 mm, kątem helisy 10° i średnicą zwoju 2 m. Model składa się z trzech obszarów: sekcji pełnego rozwinięcia (0,6 m), sekcji testowej (0,2 m) oraz sekcji stabilizacji ciśnienia (0,2 m). Sekcja pełnego rozwinięcia zapewnia, że przepływ jest odpowiednio rozwinięty przed wejściem do obszaru zainteresowania. Sekcja testowa służy do porównania z danymi eksperymentalnymi oraz do szczegółowej analizy lokalnych charakterystyk przepływu i wymiany ciepła. Sekcja stabilizacji ciśnienia została zaprojektowana w celu utrzymania stabilności ciśnienia wylotowego i zapobiegania przepływowi zwrotnemu, co pozwala uniknąć zakłóceń wyników uzyskanych w sekcji testowej. Specjalistyczne oprogramowanie do modelowania znajduje się w Tabeli materiałów.

1. Model fizyczny i siatka

  1. Otwórz oprogramowanie do modelowania. W dolnym pasku stanu wybierz Tryb szkicu (Sketch Mode), a następnie kliknij płaszczyznę Z–X, aby wejść do środowiska szkicowania.
  2. W górnym pasku narzędzi wybierz narzędzie Koło (Circle). Narysuj w punkcie początkowym koło o średnicy 14, a następnie naciśnij Enter. Kliknij Powrót do trybu 3D (Return to 3D Mode) w górnym pasku narzędzi. Koło szkicu zostanie wówczas przekonwertowane na powierzchnię.
  3. Wybierz wygenerowaną powierzchnię kołową. Kliknij narzędzie Przesuń (Move) (skrót: M) w górnym pasku narzędzi. Na powierzchni pojawi się manipulator triad (uchwyt trójosiowy). Przeciągnij żółtą sferę w centrum manipulatora do globalnego początku układu współrzędnych (0, 0, 0), który będzie służył jako punkt odniesienia dla obrotu i translacji.
  4. Przeciągnij czerwoną strzałkę wzdłuż osi X, wprowadź 1000 mm i naciśnij Enter. Obrót: Kliknij pierścień obrotu wokół osi X (niebieski lub zielony łuk), wprowadź 10° i naciśnij Enter.
  5. Kliknij narzędzie wyciągania (Pull tool) (skrót: P) w górnym pasku narzędzi i wybierz powierzchnię kołową. W lewym panelu wybierz opcję Obrót (Revolve). Następnie wybierz oś Z globalnego układu współrzędnych jako oś obrotu. Włącz opcję Helisa (Helix) w lewym panelu. Utwórz Objętość 1 (Volume 1): W polu wprowadzania lub lewym panelu wprowadź wysokość 138,87 mm oraz kąt 45,16°, a następnie naciśnij Enter. Pierwsza dziedzina płynu zostanie wygenerowana.
  6. Utwórz Objętość 2 (Volume 2): Wybierz nową powierzchnię końcową Objętości 1. Powtórz operację wyciągania helisowego, ponownie używając osi Z jako osi obrotu. Wprowadź wysokość 69,44 i kąt 22,58°.
  7. Utwórz Objętość 3 (Volume 3): Wybierz nową powierzchnię końcową Objętości 2. Utwórz Objętość 3 z tymi samymi parametrami co Objętość 2, stosując tę samą metodę.
  8. Kliknij kartę Workbench w górnym pasku menu. Kliknij przycisk Udostępnij (Share). Oprogramowanie automatycznie podświetli dwie powierzchnie przecinające się pomiędzy trzema objętościami. Kliknij przycisk Zakończ (Complete) (znacznik wyboru) po prawej stronie.
  9. Kliknij kartę Grupy (Groups) w lewym panelu. Wybierz początkową powierzchnię kołową pierwszej objętości, a następnie kliknij Utwórz nazwany wybór (Create Named Selection) i zdefiniuj ją jako wlot (in).
  10. Wybierz końcową powierzchnię kołową trzeciej objętości, naciśnij Ctrl + G, aby utworzyć grupę, i zdefiniuj ją jako wylot (out).
  11. Wybierz zewnętrzne powierzchnie cylindryczne trzech objętości i zdefiniuj je jako granice ścian: wall1, wall2 i wall3.
  12. W drzewie struktury po lewej stronie, przytrzymując Ctrl, zaznacz trzy bryły. Naciśnij Ctrl + G, aby utworzyć grupę i zmień jej nazwę na fluid.
  13. Połącz wygenerowaną geometrię z modułem Siatka (Mesh) i kliknij dwukrotnie, aby otworzyć oprogramowanie do siatkowania. W lewym drzewie kliknij Siatka (Mesh). W panelu Szczegóły (Details) po lewej stronie na dole rozwiń sekcję Wymiarowanie (Sizing) i ustaw Rozmiar elementu (Element Size) na 3.
  14. Kliknij prawym przyciskiem myszy Siatka (Mesh) w drzewie, wybierz Wstaw (Insert), Wymiarowanie (Sizing). Wybierz powierzchnię wlotową (in) jako geometrię i kliknij Zastosuj (Apply). Ustaw Rozmiar elementu (Element Size) na 0,6.
  15. Kliknij prawym przyciskiem myszy Siatka (Mesh), wybierz Wstaw (Insert), Inflacja (Inflation). Geometria (Geometry): Wybierz wszystkie trzy dziedziny płynu i kliknij Zastosuj (Apply). Granica (Boundary): Wybierz zewnętrzne powierzchnie ścian zdefiniowane jako wall, a następnie kliknij Zastosuj (Apply). Zmień opcję na Grubość pierwszej warstwy (First Layer Thickness). Wysokość pierwszej warstwy (First Layer Height): 0,01 mm. Maksymalna liczba warstw (Maximum Layers): 15. Współczynnik wzrostu (Growth Rate): 1,25.
  16. Kliknij prawym przyciskiem myszy Siatka (Mesh), wybierz Wstaw (Insert), Metoda (Method). Wybierz trzy dziedziny płynu i kliknij Zastosuj (Apply). W rozwijanym menu Metoda wybierz Przesunięcie (Sweep). W sekcji Wybór (Selection) wybierz Źródło ręczne (Manual Source). Wybierz powierzchnię wlotową (in) jako powierzchnię źródłową i kliknij Zastosuj (Apply).
  17. Kliknij prawym przyciskiem myszy Siatka (Mesh) w drzewie i wybierz Generuj siatkę (Generate Mesh). W niniejszym badaniu ściśle kontrolowano jakość siatki. Minimalna Jakość ortogonalna (Orthogonal Quality) wygenerowanej siatki wynosi powyżej 0,90.

2. Obsługa oprogramowania do symulacji

  1. Uruchom oprogramowanie do rozwiązywania zadań. Przejdź do Plik zakładka, a pod Przeczytajwybierz SiatkaNastępnie przejdź do Siatka skali oraz ustawić Siatkę utworzono w milimetrach (mm).
  2. W Ustawienia solverawybierz Solver oparty na ciśnieniu, wybierz Bezwzględny dla sformułowanie prędkościoraz włączone Opcja przejściowa przez określony czas.
    UWAGA: Poprzez nałożenie równania oscylacji na stacjonarny przypadek odniesienia i zaimplementowanie go za pomocą funkcji zdefiniowanej przez użytkownika, ruchomy układ współrzędnych może reprezentować warunki oscylacyjne.
  3. Kliknij Zdefiniowany przez użytkownikanastępnie wybierz FunkcjeW sekcji Interpreted UDFs wczytaj skompilowany plik oscylacji.
    UWAGA: Wynikowy ruch jest wyrażony zgodnie z równaniem (1). Zastosowano i zaimplementowano metodologię statycznej siatki z wykorzystaniem ruchomego układu współrzędnych. Podstawowe zjawiska fizyczne kołysania opierają się na ruchu względnym cieczy względem ścianek zbiornika. Wzbudzenie kołysania jest reprezentowane jako równoważne dynamiczne wyrazu źródłowe przyspieszenia w równaniach pędu, co umożliwia pełne odwzorowanie dynamicznych sił cieczy na stacjonarnej siatce.
    Równanie ruchu harmonicznego X=Xmaxsin(2πt/Tc), wzór, fizyka, analiza fal sinusoidalnych.      (1)
    W równaniu Tc oznacza okres kołysania, a X oznacza przemieszczenie wywołane kołysaniem.
  4. Ustaw Przyspieszenie grawitacyjne w Ykierunek do −9,81 m/s²2Pod Modelewłącz Energia i włącz Równanie energii.
  5. Pod Modelewłącz Lepki i wybrać Model naprężeń Reynoldsa (7 równań). W Ustawienia modelu naprężeń Reynoldsawybierz Liniowa zależność ciśnienie-odkształcenieDla Traktowanie przyściennewybierz Skalowalne funkcje ścianne.
  6. W Fazyzestaw Faza 1 (faza pierwotna) jako gaz i Faza 2 (faza wtórna) jako cieczPod Opcje globalnewłącz Modelowanie siły napięcia powierzchniowegooraz wybierz Model siły powierzchniowej kontinuum.
    UWAGA: Przyjęto równoważne podejście z użyciem pseudopłynu o właściwościach termofizycznych zależnych od temperatury i ciśnienia, co jest powszechnie akceptowaną metodologią w badaniach CFD nad mieszaninami wieloskładnikowymi. Przy ustalonym początkowym składzie mieszaniny obliczono i wygenerowano właściwości termofizyczne zależne od stanu — w tym gęstość, lepkość dynamiczną, przewodność cieplną, właściwą pojemność cieplną oraz charakterystyki nasycenia — korzystając z bazy danych NIST REFPROP w całym zakresie temperatur i ciśnieni roboczych. W niniejszym badaniu mieszanina zachowuje jednorodny makroskład przez cały czas trwania symulacji. Wykorzystanie zmiennych właściwości pochodzących z bazy NIST pozwala na dokładne odwzorowanie nieliniowych charakterystyk termofizycznych płynu wieloskładnikowego, unikając jednocześnie niepotrzebnego obciążenia obliczeniowego.
  7. Przyjmując za przykład mieszaninę etanu i propanu, przy jakości pary wynoszącej 0,56 i ciśnieniu 3,2 MPa, zdefiniuj właściwości fazy ciekłej w Materiały w następujący sposób:
    1. Gęstość: 393,06 kg/m³3
    2. Wymagana ciepło właściwe (Cp): 3866,4 J/(kg·K)
    3. Przewodność cieplna: 0,078798 W/(m·K)
    4. Lepkość: 5,4796 × 10⁻5 Pa·s
    5. Masa cząsteczkowa: 37,115 kg/kmol
    6. Entalpia stanu standardowego: 0
    7. Temperatura odniesienia: 321 K
  8. W sekcji Materiały zdefiniuj właściwości fazy gazowej w następujący sposób:
    1. Gęstość: 67,49 kg/m³3
    2. ciepło właściwe (Cp): 3488,7 J/(kg·K)
    3. Przewodność cieplna: 0,03035 W/(m·K)
    4. Lepkość: 1,129 × 10⁻5 Pa·s
    5. Masa molowa: 34,756 kg/kmol
    6. Entalpia stanu standardowego: 0
    7. Temperatura odniesienia: 321 K
  9. Ustawie warunek brzegowy wlotu Jako Wlot z przepływem masowym(strumień masy 300 kg/(m2·s)), the wylot jako Wylot ciśnieniowy(0 MPa), a warunek brzegowy ścianki jako Strumień ciepła(-10340 W/m2).
  10. Pod Metodywybierz algorytm PISO dla metod rozwiązywania. Dla Ułamek objętościowywybierz Geo-Rekonstrukcja.
    UWAGA: Chociaż metoda objętości płynu (VOF) jest powszechnie akceptowana w śledzeniu ewolucji topologicznej wolnej powierzchni w makroskali w procesach przelewania i termicznej przemiany fazowej, nadal istnieją ograniczenia w zakresie precyzji przechwytywania interfejsu oraz rejestrowania mikrostruktur fluktuacji międzyfazowych. Sformułowanie VOF opiera się zasadniczo na dyskretnych ułamkach objętościowych faz w komórkach. Zastosowany tutaj schemat Geo-Reconstruct znacząco redukuje dyfuzję numeryczną; niemniej jednak rozdzielczość mikrokropel poniżej siatki, formowanie rozpylenia lub mikrostruktury międzyfazowe pozostają ściśle ograniczone przez lokalne zagęszczenie siatki. W przypadku dynamiki przelewania w makroskali, konwekcji termicznej w objętości oraz praw wymiany masy podczas przemiany fazowej, które są priorytetem w niniejszym badaniu, obecny framework VOF z około 1,42 miliona elementów siatki zapewnia optymalny balans między dokładnością topologiczną a kosztem obliczeniowym.
  11. W Monitory, skonfiguruj monitorowanie dla:
    1. Ciśnienie na wlocie i wylocie sekcji testowej.
    2. Temperatura na wlocie i wylocie.
    3. Temperatura ścian.
    4. Ułamek objętości na wlocie i wylocie.
      UWAGA: Kryterium zbieżności dla residuum energii zostało ustawione na 1 × 10⁻8, natomiast dla pozostałych parametrów ustawiono wartość 1 × 10⁻4Kluczowe zmienne globalne, w tym średnia temperatura ważona powierzchniowo oraz całkowity spadek ciśnienia na sekcji pomiarowej, były monitorowane dynamicznie. Obliczenia kontynuowano do momentu, aż zmienne te przestały wykazywać dalsze fluktuacje, co zapewniło osiągnięcie przez pole przepływu stanu w pełni rozwiniętego i stabilnego.
  12. Uruchom standardową inicjalizację wyboru metody, obliczając dane ze wszystkich stref. Po inicjalizacji, w Wykonaj obliczenia panel, zestaw: Wielkość kroku czasowego: 1 × 10⁻4 s, i Liczba kroków czasowych: 1 × 106.

3. Konfiguracja postprocessingu i eksportu danych

  1. W panelu Calculation Activities kliknij Autosave (Every Flow Time), aby otworzyć okno Autosave. W ustawieniach Autosave ustaw wartość Save Data File Every [s] na 0,01 i określ Flow Time jako typ interwału zapisu. Dla parametru Save Associated Case Files Type wybierz Only if Modified, a następnie kliknij OK.
  2. Otwórz okno Contours z panelu Results. W ustawieniach Contours aktywuj opcje Filled, Node Values, Boundary Values, Global Range oraz Auto Range.
  3. Wybierz Phases jako typ contour oraz Volume Fraction jako variable, a następnie określ fazę-1 jako fazę docelową. Na koniec kliknij Save/Display, aby zwizualizować rozkład konturów.
    UWAGA: Współczynnik przenoszenia ciepła jest obliczany jako strumień ciepła przez ściankę podzielony przez siłę napędową temperatury, uzyskaną z różnicy temperatur między wlotem a wylotem sekcji badawczej. W warunkach kołysania przyjęto uśredniony w czasie współczynnik przenoszenia ciepła. Spadek ciśnienia jest określany poprzez monitorowanie różnicy między ciśnieniem wlotowym a wylotowym, a gradient tarcia spadku ciśnienia jest następnie obliczany jako stosunek tego spadku ciśnienia do długości odcinka rury.
  4. Importuj uzyskane dane do programu Excel, takie jak wartości temperatury i ciśnienia na wlocie i wylocie.
  5. Oblicz różnicę temperatur i różnicę ciśnień między wlotem a wylotem zgodnie z metodą obliczeń opisaną w sekcji 3.3.

Wyniki

Wykorzystując zwalidowany model numeryczny, przeprowadzono symulację rzeczywistego procesu skraplania w celu systematycznego zbadania zmian współczynnika przejmowania ciepła oraz spadku ciśnienia tarcia przy różnych parametrach pracy, co zapewniło teoretyczne podstawy do projektowania i optymalizacji wymienników ciepła. Główne wnioski przedstawiają się następująco: w przypadku kondensacji czystego płynu wymiana ciepła ogranicza się przede wszystkim do filmu cieczy przylegającego do ścianki rurki, gdzie temperatura powierzchni międzyfazowej gaz-ciecz jest równa temperaturze rdzenia pary, a obie one odpowiadają temperaturze nasycenia. W przeciwieństwie do tego, kondensacja mieszaniny jest procesem nierównowagowym, charakteryzującym się jednoczesnym transportem ciepła zarówno wewnątrz filmu cieczy, jak i w rdzeniu pary. W konsekwencji temperatura powierzchni międzyfazowej gaz-ciecz odbiega od temperatury nasycenia fazy objętościowej, co towarzyszy przesunięciu stężenia na powierzchni międzyfazowej względem stanu nasycenia w stanie równowagi. W trakcie tego procesu preferencyjnie skrapla się składnik mniej lotny, co powoduje gromadzenie się składnika bardziej lotnego na granicy faz. Akumulacja ta zwiększa lokalne stężenie składnika bardziej lotnego, tworząc gradient stężeń między powierzchnią międzyfazową a parą w objętości. Gradient ten wywołuje znaczny opór transportu masy, który utrudnia kondensację składnika mniej lotnego, prowadząc tym samym do obniżenia współczynnika przejmowania ciepła podczas kondensacji.

Równanie ułamka objętościowego:

Równanie różniczkowe cząstkowe dla dynamiki płynów, obejmujące transport skalarny na schemacie matematycznym.      (2)

Schemat równania dynamiki płynów ∂a/∂t + ∇·(ua) = -S/ρ; zasada zachowania masy.      (3)

Udziały objętościowe fazy gazowej i ciekłej spełniają następujący warunek:

Wzór na równowagę statyczną Σaₗ + aₑ = 1; diagram; edukacyjny koncept fizyczny.    (4)

Równanie energii:

Równanie transportu energii w dynamice płynów; zawiera symbole, operatory różniczkowe, gradient.    (5)

Model Lee dla przejść fazowych:

Równanie termodynamiczne S_al=-r·a_l·ρ_l(T-T_s)/T_s, T≥T_s, powiązane z procesami termicznymi.      (6)

Równanie równowagi statycznej, wzór na rozkład naprężeń, związany z warunkami temperatury.      (7)

gdzie S(αl) oznacza szybkość transportu masy związany ze zmianą fazy na jednostkę objętości i jednostkę czasu; αl oznacza ułamkową zawartość objętościową fazy ciekłej; αg oznacza ułamkową zawartość objętościową fazy gazowej; u⃗ oznacza wspólną prędkość obu faz m/s; ρ to gęstość mieszaniny uzyskana poprzez uśrednianie ważone ułamkiem objętościowym kg/m3; µ oznacza lepkość dynamiczną mieszaniny Pa·s; h jest średnią entalpią fazy gazowej i ciekłej J/kg; λeff to efektywna przewodność cieplna między fazą gazową a ciekłą W/(m·K); r to czynnik relaksacji czasowej 1/s, w niniejszym artykule przyjęto wartość 104; Ts to temperatura nasycenia. Zachowanie mieszaniny czynników roboczych podczas kondensacji różni się od zachowania czystych czynników roboczych, głównie ze względu na lotność składników.

Strumień masy, stopień suchości pary oraz ciśnienie nasycenia mają istotny wpływ na współczynnik przejmowania ciepła podczas kondensacji oraz spadek ciśnienia tarcia. Wraz ze wzrostem strumienia masy rośnie prędkość przepływu, co intensyfikuje zaburzenia filmu pary i tym samym usprawnia wymianę ciepła w obrębie filmu, prowadząc do ogólnego wzrostu współczynnika przejmowania ciepła. Jednocześnie naprężenia ścinające wywierane przez fazę pary na film cieczy wzrastają, co skutkuje większym spadkiem ciśnienia tarcia. Wraz ze wzrostem stopnia suchości pary zwiększają się zarówno stosunek poślizgu między fazami, jak i prędkość mieszaniny, co wzmacnia oddziaływanie ścinające między filmem cieczy a ścianką, a także ścinanie międzyfazowe między fazą pary i cieczy. Poprawia to wydajność wymiany ciepła. W tych warunkach efekty ścinania stają się dominujące, a spadek gęstości mieszaniny dodatkowo przyczynia się do wzrostu spadku ciśnienia tarcia. Ciśnienie nasycenia również odgrywa krytyczną rolę w określaniu charakterystyk przepływu i wymiany ciepła. Przy niskich ciśnieniach nasycenia gęstość pary maleje, podczas gdy prędkość przepływu rośnie, co prowadzi do powstania cieńszego filmu cieczy i zmniejszenia oporu cieplnego, zwiększając tym samym wymianę ciepła. Z kolei przy wyższych ciśnieniach nasycenia temperatura cieczy rośnie, a gęstość i lepkość cieczy maleją, co osłabia oddziaływanie ścinające między filmem cieczy a ścianką, prowadząc do zmniejszenia spadku ciśnienia tarcia. Przy stopniu suchości pary wynoszącym 0,5, gdy strumień masy wzrasta z 450 do 550 kg/(m2·s), współczynnik przejmowania ciepła rośnie z 5118 do 5637 W/(m2·K), co stanowi wzrost o 10%. Jednocześnie spadek ciśnienia tarcia wzrasta z 2523 do 3442 Pa/m, co oznacza znaczny wzrost o 36%.

Wpływy okresu i amplitudy toczenia na proces wymiany ciepła wykazują podobne trendy, w obu przypadkach obserwuje się współistnienie wzmocnienia i pogorszenia wymiany ciepła. Ruch toczący zmienia intensywność turbulencji wewnątrz filmu cieczy, a w konsekwencji wpływa na turbulentną energię kinetyczną filmu. Gdy średnia z cyklu turbulentna energia kinetyczna rośnie, dominuje transport wzmocniony przez turbulencje, co prowadzi do poprawy wymiany ciepła. Przeciwnie, gdy średnia z cyklu turbulentna energia kinetyczna spada, osłabienie turbulencji ogranicza efektywność wymiany ciepła. Jednocześnie ruch toczący intensyfikuje fluktuacje w filmie cieczy i zmienia jego grubość. Zmniejszenie grubości filmu cieczy obniża opór termiczny, a tym samym wzmacnia wymianę ciepła, podczas gdy zwiększenie grubości filmu podnosi opór termiczny i osłabia efektywność wymiany ciepła. Te dwa mechanizmy, mianowicie zmienność turbulentnej energii kinetycznej oraz zmiana grubości filmu cieczy, oddziałują na siebie i wspólnie determinują całkowite zachowanie wymiany ciepła w ciągu cyklu toczenia. W zakresie rozważanym w niniejszym badaniu wpływ okresu toczenia na efektywność wymiany ciepła mieści się w przybliżeniu w granicach ±20%, podczas gdy wpływ amplitudy toczenia wynosi ±10%.

Schemat wymiennika ciepła z oznaczonymi sekcjami stabilizacji ciśnienia i porównania; proces przepływu płynu.
Rycina 1: Schematyczny rysunek symulowanego modelu fizycznego. Ze względu na zbyt wysoki koszt obliczeniowy symulacji pełnowymiarowych rurek helikalnych, przyjęto uproszczony model z zredukowaną domeną, jak pokazano na Rycini 1. W celu walidacji z danymi eksperymentalnymi autorstwa Neeraas12, skonstruowano model trójsekcyjny (średnica rurki: 14 mm, kąt helisy: 10°, średnica zwoju: 2 m). Składa się on z sekcji pełnego rozwinięcia (0,6 m) w celu ustalenia przepływu, sekcji testowej (0,2 m) do porównania danych lokalnych oraz sekcji stabilizacji ciśnienia (0,2 m), aby zapobiec cofaniu się przepływu i utrzymać stabilność ciśnienia na wylocie. Model składa się z trzech części, z których pierwsza została opracowana na podstawie schematu z książki opublikowanej wcześniej przez Cai1. Kliknij tutaj, aby wyświetlić powiększoną wersję tej ryciny.

Wykres wymiany ciepła w stosunku do spadku ciśnienia tarcia; zależność od liczby elementów siatki; analiza wydajności termicznej.
Figura 2: Wyniki badania niezależności siatki. Figura 2 przedstawia wyniki weryfikacji niezależności siatki dla współczynnika przekazywania ciepła i spadku ciśnienia tarcia w funkcji liczby elementów siatki. Jak pokazano na rysunku, zarówno współczynnik przekazywania ciepła, jak i spadek ciśnienia tarcia znacząco maleją wraz ze wzrostem całkowitej liczby komórek od 0,60 mln do 1,33 mln. Powyżej 1,33 mln komórek zmiany obu monitorowanych wartości stabilizują się; dalsze zagęszczanie siatki do 1,85 mln komórek daje odchylenie względne poniżej 0,5%, co wskazuje na osiągnięcie niezależności siatki. Aby zrównoważyć dokładność obliczeń i nakłady zasobów, do wszystkich kolejnych symulacji przyjęto rozdzielczość siatki wynoszącą około 1,42 mln komórek. Ponadto zweryfikowano, że ta rozdzielczość siatki jest odpowiednia zarówno dla warunków stacjonarnych, jak i kołysania. Kliknij tutaj, aby wyświetlić powiększoną wersję tej figury.

Wykres słupkowy porównujący dane symulacyjne z eksperymentalnymi w zależnościności współczynnika przenikania ciepła od jakości pary.
Rysunek 3: Wyniki weryfikacji symulacji numerycznej współczynnika przenikania ciepła oraz eksperymentalnych danych Neeraasa. Przewidziane współczynniki przenikania ciepła są zgodne z danymi eksperymentalnymi w zakresie jakości pary od 0.2 do 0.8. Konkretnie, wyniki symulacji są nieco wyższe niż dane eksperymentalne przy jakości pary 0.2–0.4, podczas gdy wartości eksperymentalne nieznacznie przekraczają przewidywania numeryczne przy jakości pary 0.5–0.8. Na podstawie oceny ilościowej maksymalne odchylenie wynosi 15%. Kliknij tutaj, aby wyświetlić powiększoną wersję tego rysunku.

Wykres słupkowy spadku ciśnienia tarcia w funkcji jakości pary, porównujący dane symulacyjne i eksperymentalne.
Rycina 4: Wyniki weryfikacji symulacji numerycznej spadku ciśnienia tarcia oraz danych eksperymentalnych Neeraasa. Przewidywany spadek ciśnienia tarcia jest ogólnie nieco wyższy niż wyniki eksperymentalne, przy czym maksymalne odchylenie nie przekracza 10%. Kliknij tutaj, aby wyświetlić większą wersję tej ryciny.

Diagram ułamka objętościowego fazy gazowej; natężenia przepływu: G=350, 450, 550 kg/m²·s; przedstawiono skalę kolorów.
Rysunek 5: Ułamek objętościowy fazy gazowej przy różnych strumieniach masy (średnica = 10 mm, stopień nadmiaru pary = 0,5). Rysunek 5 przedstawia rozkłady ułamka objętościowego pary w przekroju wylotowym dla różnych strumieni masy przy tym samym stopniu nadmiaru pary. Jak pokazano na rysunku, minimalny ułamek objętościowy pary wynosi 0, co wskazuje, że ścianka pozostaje całkowicie zwilżona przez film cieczy. Przy niskich strumieniach masy wzorzec przepływu jest determinowany głównie przez grawitację i wykazuje typową strukturę przepływu warstwowego. Wraz ze wzrostem strumienia masy naprężenia ścinające wywierane przez fazę pary na film cieczy stają się progresywnie silniejsze i ostatecznie dominują nad zachowaniem przepływu, powodując stopniowe przejście wzorca przepływu z przepływu warstwowego w przepływ pierścieniowy. Ponadto stopień nadmiaru pary ma również istotny wpływ na ewolucję wzorca przepływu i wraz ze strumieniem masy determinuje zmienność struktury przepływu dwufazowego. Prosimy kliknąć tutaj, aby wyświetlić większą wersję tego rysunku.

Wykres wymiany ciepła w funkcji jakości pary; trzy krzywe dla różnych strumieni masowych (350-550 kg/m²s).
Rysunek 6: Współczynnik wymiany ciepła przy różnych strumieniach masowych. Zmiana współczynnika wymiany ciepła w zależności od różnych strumieni masowych jest przedstawiona na Rysunku 6. Przy stałej jakości pary współczynnik wymiany ciepła wzrasta wraz ze wzrostem strumienia masowego. Podczas procesu kondensacji wzdłuż wewnętrznej ścianki rury tworzy się film parowy. Wraz ze wzrostem strumienia masowego rośnie prędkość przepływu, co intensyfikuje zaburzenia filmu parowego i zwiększa wymianę ciepła wewnątrz filmu, redukując tym samym opór cieplny. W konsekwencji współczynnik wymiany ciepła staje się wyższy przy zwiększonych strumieniach masowych. Jednocześnie wraz ze wzrostem strumienia masowego zwiększa się również liczba Reynoldsa odpowiadająca filmowi cieczy. Ogólnie rzecz biorąc, strumień masowy ma istotny wpływ na współczynnik wymiany ciepła. Aby zobaczyć powiększoną wersję tego rysunku, kliknij tutaj.

Wykres spadku ciśnienia tarcia w funkcji jakości pary. Linie przedstawiają natężenia przepływu G=350, 450, 550 kg/(m²·s).
Rycina 7: Spadek ciśnienia tarcia przy różnych strumieniach masy.Rycina 7 przedstawia zmienność spadku ciśnienia tarcia w zależności od różnych warunków strumienia masy. Wyniki wskazują, że przy tej samej jakości pary spadek ciśnienia tarcia wzrasta znacząco wraz ze wzrostem strumienia masy. Wynika to głównie z faktu, że wyższy strumień masy prowadzi do wyższej prędkości przepływu, co zwiększa siły ścinania wywierane przez fazę par na film ciecz oraz naprężenia ścinające przy ściance, skutkując tym samym większym spadkiem ciśnienia tarcia. Ogólnie rzecz biorąc, strumień masy ma wyraźny wpływ na spadek ciśnienia tarcia. Kliknij tutaj, aby wyświetlić powiększoną wersję tej ryciny.

Wykres ułamka objętościowej zawartości fazy gazowej przedstawiający rozkład zakodowany kolorami dla różnych wartości (0,3, 0,5, 0,7, 0,9).
Rycina 8: Ułamek objętościowy fazy gazowej przy różnych parametrach jakości pary (średnica = 10 mm). Rycina 8 przedstawia rozkłady ułamka objętościowego pary na wylocie dla czterech wartości jakości pary. Ułamek objętościowy gwałtownie rośnie przy niskiej jakości pary, natomiast przy wysokiej jakości pary stabilizuje się w pobliżu wartości 1. Zidentyfikowano cztery odrębne wzorce przepływu: warstwowy, półpierścieniowy, pierścieniowy oraz mgielny. Przy niskiej jakości pary dominuje grawitacja, co skutkuje przepływem warstwowym z parą na górze i cieczą na dole. Wraz ze wzrostem jakości pary ścinanie międzyfazowe zastępuje grawitację jako dominujący mechanizm, prowadząc przepływ przez reżimy półpierścieniowy i pierścieniowy aż do przepływu mgielnego. Prosimy kliknąć tutaj, aby zobaczyć powiększoną wersję tej ryciny.

Diagram udziału objętościowej fazy gazowej; porównanie ciśnień przy 3 MPa i 5 MPa ze skalą kolorystyczną.
Rysunek 9: Udział objętościowy fazy gazowej przy różnych ciśnieniach nasycenia. Wraz ze wzrostem ciśnienia nasycenia gęstość cieczy maleje, natomiast gęstość pary rośnie, co prowadzi do zmiany różnicy gęstości między dwiema fazami i ogólnego wzrostu gęstości mieszaniny. Jednocześnie zmieniają się charakterystyki poślizgu gaz-ciecz, a ścinanie międzyfazowe między dwiema fazami ulega osłabieniu, co skutkuje zmniejszeniem udziału objętościowego pary. Zmiany te są bezpośrednio odzwierciedlone w trendach współczynnika przenikania ciepła i spadku ciśnienia tarcia. Kliknij tutaj, aby wyświetlić powiększoną wersję tego rysunku.

Wykres słupkowy analizujący współczynnik prze conductivity ciepła w funkcji jakości pary przy ciśnieniach 3 MPa i 5 MPa.
Rycina 10: Współczynnik prze conductivity ciepła przy różnych ciśnieniach nasycenia.Rycina 10 przedstawia współczynniki prze conductivity ciepła dla różnych wartości jakości pary i ciśnień nasycenia. Przy stałej jakości pary niższe ciśnienie nasycenia skutkuje wyższym współczynnikiem prze conductivity ciepła. Z punktu widzenia mechanistycznego, wyższe ciśnienie zwiększa gęstość pary, co redukuje prędkość przepływu i naprężenia ścinające na granicy faz. Powoduje to pogrubienie filmu cieczy, zwiększając tym samym opór cieplny i pogarszając wymianę ciepła. Ponadto wpływ ciśnienia nasycenia staje się bardziej wyraźny przy wyższych wartościach jakości pary, gdzie dominuje prędkość pary, a zmiany gęstości indukowane ciśnieniem powodują większe wahania naprężeń ścinających na granicy faz. Prosimy kliknąć tutaj, aby wyświetlić powiększoną wersję tej ryciny.

Wykres słupkowy spadku ciśnienia tarcia w funkcji jakości pary przy 3 Mpa i 5 Mpa, ilustrujący dynamikę przepływu płynu.
Rysunek 11: Spadek ciśnienia tarcia przy różnych ciśnieniach nasycenia. Rysunek 11 przedstawia zmienność spadku ciśnienia tarcia przy różnych ciśnieniach nasycenia. Wyniki wskazują, że przy tej samej jakości pary spadek ciśnienia tarcia maleje wraz ze wzrostem ciśnienia nasycenia. W połączeniu z rozkładem prędkości, polem temperatury przechłodzenia i rozkładem ułamka objętościowego pary przy różnych ciśnieniach nasycenia, wyniki te wskazują, że wyższe ciśnienie nasycenia odpowiada wyższej temperaturze płynu, co wiąże się ze spadkiem zarówno gęstości, jak i lepkości cieczy. W rezultacie oddziaływanie ścinające między filmem cieczy a ścianką ulega osłabieniu, co prowadzi do zmniejszenia spadku ciśnienia tarcia. Aby wyświetlić powiększoną wersję tego rysunku, kliknij tutaj.

Ułamk objętościowy fazy gazowej; wyniki symulacji; różne proporcje czasu; mapowanie kolorów; dynamika płynów.
Rycina 12: Ułamk objętościowy fazy gazowej dla różnych okresów obrotu (jakość pary = 0,5, strumień masy = 550 kg/(m2·s, A = 3 m). Przy stałej amplitudzie obrotu krótszy okres obrotu prowadzi do silniejszego dodatkowego efektu bezwładności indukowanego ruchem oscylacyjnym, co skutkuje intensywniejszymi fluktuacjami prędkości w polu przepływu. Fluktuacje te wykazują również wyraźne zachowanie okresowe z naprzemiennymi fazami przyspieszenia i spowolnienia przepływu. Jednocześnie ruch obrotowy modyfikuje przestrzenny rozkład filmu cieczy i zmienia strukturę przepływu, wpływając tym samym na wymianę ciepła. Wraz ze wzrostem średniej grubości filmu cieczy rośnie opór cieplny filmu, co osłabia efektywność wymiany ciepła. Przeciwnie, gdy średnia grubość filmu cieczy maleje, opór cieplny filmu spada, zwiększając tym samym wymianę ciepła. Klasyfikacja reżimów przepływu opiera się na kryteriach przejścia struktury przepływu zaproponowanych w Literaturze4. Aby wyświetlić powiększoną wersję tej ryciny, kliknij tutaj.

Wykres współczynnika przekazywania ciepła, warunki stacjonarne vs toczenie; wyniki analizy wymiany ciepła.
Rysunek 13: Współczynnik przekazywania ciepła przy różnych okresach toczenia. Rysunek 13 porównuje uśrednione w czasie współczynniki przekazywania ciepła (HTC) w warunkach toczenia w stosunku do stanu stacjonarnego. Toczenie zmienia HTC w zakresie ±20%, wykazując zarówno poprawę, jak i pogorszenie wymiany ciepła. Przy niskich wartościach HTC (niższa jakość pary), toczenie wzmacnia przekazywanie ciepła — szczególnie przy krótszych okresach toczenia — poprzez intensyfikację turbulencji w filmie cieczy i fluktuacji międzyfazowych. I odwrotnie, przy wysokich wartościach HTC (wyższa jakość pary), toczenie upośledza wymianę ciepła poprzez ściskanie rdzenia pary i zwiększanie grubości filmu cieczy (poprzez uśrednione pogrubienie i efekty odśrodkowe w przepływie pierścieniowym), co zwiększa opór cieplny. W związku z tym zaleca się zastosowanie odpowiedniego marginesu projektowego w aplikacjach offshore. Każdy punkt danych na rysunku odpowiada niezależnemu i deterministycznemu przypadkowi symulacji numerycznej. Rozwiązanie CFD równań rządzących nie uwzględnia szumu pomiarowego, pomijając wariancję statystyczną właściwą dla wielokrotnych prób doświadczalnych; dlatego słupki błędów oparte na rozkładach statystycznych nie są właściwe ani konieczne. Kliknij tutaj, aby zobaczyć powiększoną wersję tego rysunku.

Wykres współczynnika przenoszenia ciepła; porównanie w okresach toczenia; zawiera wskaźniki wariancji 10%.
Rysunek 14: Współczynnik przenoszenia ciepła przy różnych amplitudach toczenia. Rysunek 14 przedstawia porównanie uśrednionych w czasie współczynników przenoszenia ciepła (HTC) przy różnych amplitudach toczenia w stosunku do stanu stacjonarnego. Amplituda toczenia zmienia HTC w zakresie ±10%, wykazując zarówno poprawę, jak i pogorszenie parametrów. Przy niskich wartościach HTC (niższa jakość pary), toczenie zwiększa intensywność przenoszenia ciepła — bardziej wyraźnie przy większych amplitudach — poprzez intensyfikację turbulencji filmu cieczy oraz fluktuacji międzypowierzchniowych. I odwrotnie, przy wysokich wartościach HTC (wyższa jakość pary), toczenie pogarsza przenoszenie ciepła poprzez ściskanie rdzenia parowego i pogrubianie filmu cieczy (w wyniku średniego pogrubienia i efektów odśrodkowych w przepływie pierścieniowym), co zwiększa opór termiczny. W konsekwencji, w zastosowaniach morskich zaleca się stosowanie odpowiedniego marginesu projektowego. Kliknij tutaj, aby wyświetlić powiększoną wersję tego rysunku.

Strumień masyParaCiśnienieŚrednica rury (mm)Kąt owinięcia Średnica nawinięcia (m)okres przesunięcia (s)amplituda przechyłu (m)
kg/(m2·s)jakośćMPa°
350–5500.1–0.93–510422–52–3

Tabela 1: Symulowane warunki pracy. Tabela 1 podsumowuje warunki symulacji dla mieszaniny lekkich węglowodorów w sekcji skraplania w rzeczywistym procesie przemysłowym15. Ciecz robocza składa się z metanu, propanu, izopentanu, etylenu i azotu, w stosunku molowym 55,314:1,407:0,04:23,709:19,53. Aby dokładnie odwzorować nieliniowe zachowanie jednorodnej mieszaniny we wszystkich warunkach pracy przy jednoczesnym zminimalizowaniu kosztów obliczeniowych, wykorzystano właściwości pochodne z NIST REFPROP.

Dyskusja

Konfiguracja trójsekcyjna, kluczowa dla zapewnienia wiarygodności symulacji, ustanawia w pełni rozwinięte warunki przepływu przed sekcją badawczą i tłumi przepływ zwrotny na wylocie, optymalizując tym samym dokładność przewidywanych wyników. Podejście to zostało odzwierciedlone również w poprzednich badaniach nad rurami poziomymi1. Podczas procesu generowania siatki w niniejszym badaniu, wysokość pierwszej warstwy siatki, liczba warstw przyściennych oraz minimalny wymagany poziom jakości ortogonalnej są równie istotne, ponieważ bezpośrednio wpływają na dokładność wyników symulacji. Strumień masy, stopień odparowania oraz ciśnienie nasycenia znacząco wpływają na wymianę ciepła podczas kondensacji i spadek ciśnienia. Zwiększenie strumienia masy zwiększa prędkość pary i ścinanie międzyfazowe, co prowadzi do wzrostu zarówno współczynnika wymiany ciepła, jak i tarczowego spadku ciśnienia. Zwiększenie stopnia odparowania również wzmacnia ścinanie międzyfazowe i sprzyja przejściu z przepływu warstwowego w kierunku przepływu pierścieniowego i mglistego. Z kolei zwiększenie ciśnienia nasycenia redukuje zarówno współczynnik wymiany ciepła, jak i tarczowy spadek ciśnienia. Tendencje te są ogólnie zgodne z wcześniejszymi badaniami eksperymentalnymi i numerycznymi dotyczącymi kondensacji węglowodorów6,7,8,9,13. W przypadku mieszanin nieazeotropowych należy również uwzględnić dodatkowy opór transportu masy spowodowany redystrybucją składników w pobliżu granicy faz para–ciecz10,11,12.

Ważnym odkryciem jest to, że ruch toczny może albo zwiększać, albo pogarszać intensywność wymiany ciepła podczas kondensacji. W badanym zakresie okres toczenia zmienia wydajność wymiany ciepła o około ±20%, podczas gdy amplituda toczenia powoduje wahania rzędu około ±10%. Zachowanie to wynika głównie z połączonego wpływu turbulencji filmu cieczy oraz zmienności grubości filmu. Zwiększona turbulencja lub cieńszy film cieczy intensyfikują wymianę ciepła, natomiast zmniejszona turbulencja lub pogrubienie filmu prowadzą do jej pogorszenia. Zatem całkowita odpowiedź wymiany ciepła zależy od konkurencji między tymi dwoma mechanizmami. Przy zastosowaniu tej metody należy wziąć pod uwagę kilka kwestii numerycznych. Ponieważ przewidywana wymiana ciepła i spadek ciśnienia są wrażliwe na grubość filmu cieczy oraz zachowanie na granicy faz, niezbędna jest wystarczająca rozdzielczość siatki przy ściankach oraz odpowiedni krok czasowy. Ponadto zbieżność nie powinna być oceniana wyłącznie na podstawie reszt. Należy również monitorować kluczowe wielkości fizyczne, w tym temperaturę, ciśnienie, ułamkową objętość pary oraz spadek ciśnienia, aby odróżnić oscylacje numeryczne od rzeczywistych fluktuacji wywołanych toczeniem.

W niniejszym badaniu uwzględniono jednak ograniczoną liczbę warunków przechyłu, a w celu uzyskania pełniejszego zrozumienia wpływu dynamicznych warunków pracy na wydajność kondensacji konieczne są szersze badania parametryczne. W praktycznych zastosowaniach morskich LNG wymienniki ciepła mogą być poddane złożonym ruchom sześciu stopni swobody wywołanym ruchem jednostki, w tym połączonymi ruchami przechyłu, pochylenia i odchylenia. Efekty dynamiczne te mogą nieustannie modyfikować pole grawitacyjne, struktury przepływu wtórnego oraz rozkład filmu cieczy wewnątrz rurki helikalnej, wpływając tym samym na lokalne charakterystyki wymiany ciepła i spadku ciśnienia. W związku z tym przyszłe badania powinny dotyczyć sprzężonych efektów różnych amplitud, częstotliwości i kierunków przechyłu, aby opracować bardziej kompletny system oceny wydajności helikalnych wymienników ciepła w środowisku morskim.

Ponadto wymagana jest dalsza walidacja z wykorzystaniem praktycznych danych operacyjnych, szczególnie biorąc pod uwagę różnicę między czynnikiem roboczym przyjętym w niniejszym badaniu a nieazeotropowymi mieszaninami węglowodorów stosowanymi w rzeczywistych przemysłowych procesach LNG. W realnych systemach LNG mieszane chłodziwa zazwyczaj wykazują znaczny poślizg temperaturowy oraz złożoną równowagę fazową ze względu na oddziaływania między wieloma składnikami. Charakterystyki te mogą wpływać na mechanizm kondensacji, międzyfazowy transport masy oraz lokalne właściwości termofizyczne. Choć obecny model skutecznie przewiduje ogólne trendy przepływu i wymiany ciepła, niezbędne są badania eksperymentalne z wykorzystaniem praktycznych pięcioskładnikowych chłodziw mieszanych, takich jak mieszaniny azot/metan/etylen/propan/izopentan, aby dalej zweryfikować wiarygodność modelu i poprawić jego przydatność w warunkach przemysłowych.

Ponadto, stosowalność wybranego modelu turbulencji w warunkach przepływu pierścieniowo-mglistego o wysokiej jakości pary wymaga dalszych badań. W tym reżimie przepływu mogą wystąpić silne deformacje międzyfazowe, porywanie kropel oraz intensywne oddziaływania turbulencyjne, co prowadzi do złożonych mechanizmów wymiany pędu i energii między rdzeniem parowym a fazą ciekłą. Konwencjonalne modele turbulencji mogą wprowadzać niepewności podczas przewidywania tych silnie anizotropowych charakterystyk przepływu dwufazowego. Dlatego w przyszłych badaniach można rozważyć zaawansowane modele turbulencji, udoskonalone korelacje sił międzyfazowych lub numeryczne metody rozdzielające interfejs w celu zwiększenia dokładności przewidywań w ekstremalnych warunkach pracy. Niezawodność wyników numerycznych przy ciśnieniach roboczych znacznie przekraczających zakres zbadany w niniejszej pracy (3–5 MPa) również wymaga dalszej weryfikacji za pomocą dodatkowych danych doświadczalnych. Zmiany ciśnienia mogą silnie wpływać na właściwości termofizyczne czynnika chłodniczego, charakterystykę równowagi fazowej oraz zachowanie kondensacji, co prowadzi do rozbieżności między przewidywaniami numerycznymi a rzeczywistą wydajnością. Podobnie, w niniejszym badaniu analizowano strumienie masy w zakresie 350–550 kg/(m2·s), podczas gdy wymienniki ciepła LNG mogą pracować przy wyższych strumieniach masy. To, czy proponowany model numeryczny zachowuje wystarczającą dokładność i ogólną stosowalność przy wyższych strumieniach masy, pozostaje do potwierdzenia w dalszych badaniach doświadczalnych i numerycznych.

Pomimo tych ograniczeń, niniejsze badanie dostarcza istotnych spostrzeżeń teoretycznych oraz wytycznych ilościowych dla projektowania i optymalizacji spiralnych wymienników ciepła do zastosowań w technologii LNG. W badanym zakresie roboczym zwiększenie marginesu projektowego o około 20% może skutecznie zrekompensować spadek wydajności spowodowany warunkami przechyłu, zapewniając praktyczne podejście inżynieryjne w celu zagwarantowania niezawodnej pracy w dynamicznych środowiskach morskich. Wyniki te przyczyniają się nie tylko do głębszego zrozumienia charakterystyki kondensacji w spiralnych wymiennikach ciepła w warunkach ruchu, ale stanowią również wartościowe odniesienie dla rozwoju bardziej wydajnych i odpornych systemów wymiany ciepła w instalacjach LNG.

Oświadczenia

Autorzy oświadczają, że nie mają żadnych znanych konfliktów interesów finansowych ani relacji osobistych, które mogłyby wpłynąć na pracę opisaną w niniejszej publikacji.

Podziękowania

Badania te są wspierane przez Basic Research Project for Universities of Liaoning Provincial Education Department (LJ212512594008 dla Xianshi Fang) oraz Shenyang Key Laboratory of Industrial Product Testing Technology and Intelligent Testing Equipment (JC2503, JC2512).

Materiały

Lista materiałów użytych w tym artykule
NazwaFirmaNumer katalogowyKomentarze
FluentANSYS2020r1oprogramowanie do symulacji
SpaceClaimANSYS2020r1oprogramowanie do modelowania

Bibliografia

  1. Cai W, Fang X, Chen J. The core liquefaction facility in many floating liquefaction facilities is the spiral-wound heat exchanger. In: Sustainable Liquefied Natural Gas. Elsevier; 2024:85-123.
  2. Yu J, Huo R, Shen H, et al. A simulation study on the condensation flow and thermal control characteristics of mixed refrigerant in a dimpled tube. Appl Therm Eng. 2023;120889.
  3. Li J, Hu H, Wang H. Numerical investigation on flow pattern transformation and heat transfer characteristics of two-phase flow boiling in the shell side of LNG spiral wound heat exchanger. Int J Therm Sci. 2020;152:106289.
  4. Fang X, Qiu G, Chen J, et al. A new frictional pressure drop correlation based on flow patterns for hydrocarbon refrigerants condensation flow. Int J Refrig. 2025;170:214-223.
  5. Fang X, Qiu G, Li Q, et al. A new heat transfer correlation based on flow patterns for hydrocarbon refrigerants condensation flow. Case Stud Therm Eng. 2024;56:104244.
  6. Fries S, Skusa S, Luke A. Heat transfer and pressure drop of condensation of hydrocarbons in tubes. Heat Mass Transfer. 2019;55:33-40.
  7. Zhuang XR, Gong MQ, Zou X, et al. Experimental investigation on flow condensation heat transfer and pressure drop of R170 in a horizontal tube. Int J Refrig. 2016;66:105-120.
  8. Zhuang XR, Chen GF, Zou X, et al. Experimental investigation on flow condensation of methane in a horizontal smooth tube. Int J Refrig. 2017;78:193-214.
  9. López-Belchí A, Illán-Gómez F, García-Cascales JR, et al. Condensing two-phase pressure drop and heat transfer coefficient of propane in a horizontal multiport mini-channel tube: experimental measurements. Int J Refrig. 2016;68:59-75.
  10. Smit FJ, Meyer JP. Condensation heat transfer coefficients of the zeotropic refrigerant mixture R-22/R-142b in smooth horizontal tubes. Int J Therm Sci. 2002;41:625-630.
  11. Berrada N, Marvillet C, Bontemps A, et al. Heat transfer in-tube condensation of a zeotropic mixture of HFC23/HFC134a in a horizontal smooth tube. Int J Refrig. 1996;19:463-472.
  12. Neeraas BO. Condensation of hydrocarbon mixtures in coil-wound LNG heat exchangers: tube-side heat transfer and pressure drop. Trondheim: Norwegian Institute of Technology; 1993.
  13. Li S, Cai W, Chen J, et al. Numerical study on condensation heat transfer and pressure drop characteristics of ethane/propane mixture upward flow in a spiral pipe. Int J Heat Mass Transfer. 2018;121:170-186.
  14. Qiu GD, Cai WH, Wu ZY, et al. Numerical simulation of forced convective condensation of propane in a spiral tube. J Heat Transfer. 2015;137:041502.
  15. Fang X, Guo Z, Tang K, et al. Numerical Study on Condensation Flow and Heat Transfer of Hydrocarbon Mixtures in Inclined Tubes under Static and Swaying Conditions. Front Heat Mass Transf. 2026; 24(2):18.

Przedruki i uprawnienia

Tagi

Wymiennik ciepła z rurkami spiralnymimieszaniny nieazeotropowesymulacja numerycznaprodukcja LNGfrictionalny spadek ciśnieniaamplituda falowaniaokres falowania