Üç farklı ve birbirine karışmayan sıvı veya gaz fazının aynı hesaplama alanında hareket etmesi, hesaplamalı akışkanlar dinamiğinin zor problemlerinden biridir. Arayüzler hareket eder, şekil değiştirir, birleşir, ayrılır ve üç fazın kesiştiği bölgelerde yüzey gerilimi kuvvetleri aynı anda etkili olur. Adam Y. Ghoneim tarafından geliştirilen bu çalışma, söz konusu problemi klasik bir hesaplama ağı oluşturmadan çözmek amacıyla Smoothed Particle Hydrodynamics (SPH) ile Phase-Field yaklaşımını birleştirmektedir.

Yöntemde sıvılar arasındaki arayüz keskin bir sınır olarak değil, belirli kalınlığa sahip geçişli bir Phase-Field bölgesi olarak temsil edilmektedir. Üç fazın evrimi Cahn–Hilliard denklemleriyle, akış alanı Navier–Stokes denklemleriyle çözülmektedir. SPH ağırlık fonksiyonlarının tutarlılığını geliştirmek için Reproducing Kernel (RK) ve Differential Reproducing Kernel (DRK) yaklaşımları kullanılmıştır. Basınç çözümünde ise klasik yoğunluk tabanlı Equation of State yerine sönümleyici bir basınç difüzyonu içeren General Pressure Equation tercih edilmiştir.

Yöntem Couette, Poiseuille ve lid-driven cavity gibi standart benchmarklardan başlayarak üç faz ayrışması, sıvı lensi, yüksek yoğunluk oranlı su–yağ–hava sistemi, düşen damla, yükselen kabarcık ve serbest yüzeyli tank karışımı gibi daha karmaşık problemlerde sınanmıştır. Sıvı lensi testinde analitik sonuçlara göre maksimum bağıl hata %5'in altında kalmış; RK ve DRK yakınsama eğimleri sırasıyla 1,86 ve 1,93 ile yaklaşık ikinci mertebe davranış göstermiştir.

Yöntem yaklaşık 1000:1 yoğunluk oranındaki su–hava sistemini de kararlı biçimde çözebilmiştir. Ancak Phase-Field mobility parametresinin arayüz kalınlığı ve hacim korunumu üzerinde güçlü etkisi vardır: örnek dönen kabarcık testinde M=0,01 için hacim kaybı yalnız %0,3 iken M=0,05'te %20'ye çıkmıştır. Bu nedenle çalışma yalnız yeni bir çözüm yöntemi sunmakla kalmamakta, yöntemin sayısal parametrelerinin dikkatli ayarlanması gerektiğini de ortaya koymaktadır.

Üç fazlı akışı hesaplamak neden daha zor?

İki fazlı akışlarda bile sıvı–gaz veya sıvı–sıvı arayüzünün hareketini doğru biçimde izlemek zordur. Üç fazlı sistemde buna üç ayrı arayüz gerilimi ve üç fazın buluşabildiği triple-junction bölgeleri eklenir. Arayüzün birleşmesi, kopması veya topolojik biçim değiştirmesi klasik mesh tabanlı yöntemlerin yeniden ağ oluşturma veya özel interface-capturing algoritmalarına ihtiyaç duymasına yol açabilir.

Phase-Field Yaklaşımı Üç Fazı Nasıl Temsil Ediyor?

Phase-Field yönteminde her faz için \(\psi_p\) adı verilen sürekli bir order parameter kullanılır. Fazın saf bölgesinde bu değer 0 veya 1'e yaklaşırken iki faz arasındaki arayüzde değer sürekli olarak değişir. Böylece keskin bir sınır yerine sonlu kalınlıkta diffuse interface oluşur.

Üç faz için:

\[ \psi_1+\psi_2+\psi_3=1 \]

koşulu uygulanmaktadır.

Dolayısıyla:

\[ \psi_3=1-\psi_1-\psi_2 \]

olduğundan yalnız iki bağımsız Cahn–Hilliard denklemi çözmek yeterlidir.

Sistemin Helmholtz serbest enerjisi genel olarak:

\[ E(\psi)= \int_{\Omega} \sum_{p=1}^{N} \left[ F(\psi_p)+ \frac{\epsilon^2}{2}|\nabla\psi_p|^2 \right]dx \]

şeklinde tanımlanmaktadır.

İlk terim fazların saf durumlara yönelmesini sağlayan bulk free-energy bileşenidir. Gradient terimi ise arayüz oluşumunun enerji maliyetini temsil eder.

Cahn–Hilliard denklemi ne yapıyor?

Kimyasal potansiyel serbest enerjinin order parameter'a göre varyasyonundan elde edilir ve fazların dağılımı:

\[ \frac{D\psi_p}{Dt} = M\nabla^2\mu_p \]

ile evrilir.

Buradaki \(M\), Phase-Field mobility parametresidir. Bu değer yalnız matematiksel bir sabit değildir; arayüzün difüzyonu, kalınlığının korunması ve hacim kaybı arasında önemli bir denge oluşturur.

SPH Neden Meshfree Olarak Adlandırılıyor?

Smoothed Particle Hydrodynamics'te hesaplama alanı klasik sonlu eleman veya sonlu hacim ağlarına bölünmez. Bunun yerine alan parçacıklardan oluşur ve bir değişkenin belirli bir parçacıktaki değeri komşu parçacıkların kernel ağırlıklı katkısından hesaplanır.

Bir alan değişkeninin SPH yaklaşımı:

\[ \hat{u}_i = \sum_{j=1}^{n} u_j W_{ij} \]

biçimindedir.

Bu yapı sayesinde düzensiz nokta dağılımları, hareket eden parçacıklar ve karmaşık geometriler mesh connectivity oluşturmadan işlenebilir.

RK ve DRK Neden Gerekli?

Klasik SPH kernel fonksiyonlarının özellikle sınır bölgelerinde partition-of-unity ve reproducibility koşullarını tam olarak sağlayamaması doğruluk problemlerine neden olabilir. Reproducing Kernel yaklaşımı, ağırlık fonksiyonlarını moment matrisi üzerinden düzelterek bu koşulları doğrudan uygular.

Differential Reproducing Kernel yaklaşımının temel farkı, yeniden üretilebilirlik koşulunu yalnız ağırlık fonksiyonuna değil onun uzaysal türevlerine de uygulamasıdır. Böylece RK'nin açık türev hesabında gereken moment matrisi türevlerinin hesaplanması önlenebilir.

Bu çalışmada yüksek mertebeli DRK türevleri doğrudan kullanılmamıştır. RK ve DRK yalnız ağırlık fonksiyonu ile birinci mertebe gradyentleri için kullanılmış; Laplasyen daha dayanıklı bir SPH ifadesiyle hesaplanmıştır. Bu nedenle DRK'nin hesaplama avantajı vardır fakat dramatik düzeyde değildir.

General Pressure Equation Neyi Değiştiriyor?

Weakly-compressible SPH'de basınç genellikle yoğunluk değişiminden Equation of State yardımıyla hesaplanır. Çok fazlı sistemlerde yoğunluk arayüzde ani biçimde değiştiğinden bu yaklaşım istenmeyen basınç ve yoğunluk salınımlarına neden olabilir.

Bu çalışma bunun yerine:

\[ \frac{\partial P}{\partial t} = -C^2(\nabla\cdot\rho V) + \nu_{GPE}\nabla^2P \]

biçimindeki General Pressure Equation'ı kullanmıştır.

\(\nu_{GPE}\nabla^2P\) terimi basınç alanındaki sayısal salınımları sönümleyen bir difüzyon mekanizmasıdır. Çözüm yerel ve explicit gerçekleştirilebildiğinden her zaman adımında global Pressure-Poisson lineer sistemi çözmek gerekmez.

Zaman integrasyonu

CHNS sistemi ikinci mertebe Runge–Kutta yaklaşımıyla explicit olarak ilerletilmiştir. Zaman adımı; yapay ses hızı, maksimum viskozite ve Phase-Field mobility kısıtlarının minimumuna göre seçilmektedir:

\[ \Delta t= \min \left( \frac{\gamma_{CFL}h}{C}, \frac{\gamma_{CFL}h^2}{2\nu_{max}}, \frac{\gamma_{CFL}h^2}{2M} \right) \]

ve çalışmada \(\gamma_{CFL}=0,25\) kullanılmıştır.

Particle Shifting Neden Kritik?

Lagrangian SPH'de parçacıklar akışla birlikte hareket ettiği için zamanla birbirlerine yaklaşabilir ve kümelenebilir. RK/DRK'de yerel moment matrisinin tersinin hesaplanması gerektiğinden bu durum yalnız görsel bir düzensizlik değildir; matris kötü koşullu hale gelerek çözümü kararsızlaştırabilir.

Çalışma particle shifting ile yüksek parçacık yoğunluğuna sahip bölgelerden parçacıkları daha seyrek alanlara doğru küçük miktarlarda kaydırmaktadır.

PDF'nin 24. sayfasındaki karşılaştırmada shifting uygulanmadığında belirgin kümelenmeler oluşurken shifting sonrasında parçacıklar yaklaşık homojen dağılıma dönmektedir.

Serbest yüzey parçacıklarında ise shifting miktarı azaltılmaktadır; aksi halde gerçek serbest yüzey geometrisi yapay biçimde bozulabilir.

Çalışmanın Yöntemi ve Bulguları

Couette flow doğrulaması

50 parçacıklı düşük Reynolds sayılı Couette probleminde RK-SPH ve DRK-SPH sonuçları analitik çözüme çok yakın çıkmıştır. Farklı zamanlardaki L2-norm hataları yaklaşık \(10^{-3}\)–\(10^{-2}\) mertebesindedir ve RK–DRK sonuçları arasındaki maksimum fark yaklaşık %0,16'dır.

Poiseuille flow doğrulaması

Üç farklı tepe hızında RK ve DRK çözümleri analitik profile iyi uyum göstermiştir. İki yaklaşımın hata değerlerindeki maksimum farklılık yaklaşık %1,7 olmuştur. RK bazı durumlarda analitik çözüme biraz daha yakınken DRK maksimum hızı hafifçe yüksek hesaplamıştır.

Lid-driven cavity

Re=100 için 50×50 parçacıklı alan çözülmüştür. Eulerian yaklaşımda sabit parçacıkların sonuçları Ghia ve arkadaşlarının klasik referans verileriyle çok iyi uyum göstermiştir. Lagrangian çözümde parçacık düzensizliği nedeniyle hız piklerinde daha fazla hata oluşmuş; çözünürlüğün 10.000 parçacığa yükseltilmesi bu farkı azaltmıştır.

Üç Faz Kendiliğinden Ayrışabiliyor mu?

10.000 parçacıklı keyfi 2B geometride başlangıçta rastgele dağıtılan üç order parameter zamanla ayrı kırmızı, yeşil ve mavi faz bölgelerine ayrılmıştır. Bu süreç sırasında toplam serbest enerji sürekli düşmüştür.

Aynı davranış 64.000 parçacıklı 3B küre ve torus geometrilerinde de elde edilmiştir. Böylece yöntem düzenli bir kare grid veya basit dikdörtgen geometriyle sınırlı kalmamıştır.

Yüzey Gerilimi Arayüz Geometrisini Nasıl Değiştiriyor?

İki fazlı kontrol testinde başlangıçta kare biçimindeki gömülü faz, yüksek eğrilikli köşelerdeki yüzey gerilimi kuvvetleri nedeniyle zamanla dairesel hale gelmiştir. Bu sonuç serbest enerjinin arayüz alanını küçültme eğilimiyle uyumludur.

Üç fazlı durumda yüzey gerilimi oranı değiştirilmiştir:

Sıvı Lensinde Analitik Çözüme Ne Kadar Yaklaşıldı?

Üç fazın buluştuğu sıvı lensi probleminde Neumann yasasıyla teorik temas açıları ve lens genişliği hesaplanmıştır.

Kaynak bütün test edilen yüzey gerilimi oranlarında maksimum bağıl hatanın %5'in altında kaldığını bildirmektedir.

Yakınsama

Bu benchmark için yaklaşık ikinci mertebe doğruluk görülmüştür. 22.500 parçacıktan sonra kararlı durum serbest enerjisinde yalnız küçük değişimler kalmıştır.

1000:1 Yoğunluk Oranı Kararlı Biçimde Çözülebildi mi?

Su–zeytinyağı–hava sistemiyle yüksek yoğunluk oranı testi yapılmıştır.

Su/hava oranı yaklaşık 1000:1 olmasına rağmen RK ve DRK tabanlı PF-SPH çözümleri kararlı kalmıştır. Sonuçlar OpenFOAM multiphaseInterFoam FVM–VOF çözümüyle genel interfacial davranış açısından uyumludur.

Bununla birlikte damla çarpması ve geri tepmesinde PF-SPH damlacık geometrisi OpenFOAM sonucuna göre biraz daha basık görünmektedir. Yazar bunun olası nedenlerinden biri olarak Phase-Field mobility kaynaklı arayüz yumuşamasını ve hacim kaybını tartışmaktadır.

Mobility Parametresi Hacim Korunumunu Nasıl Etkiliyor?

Bu sonuç mobility değerinin büyütülmesinin otomatik olarak daha doğru çözüm üretmediğini göstermektedir. Çok küçük değerlerde arayüz adveksiyon altında yayılırken çok yüksek değerlerde şekil korunmasına rağmen hacim kaybı artabilir.

Üç fazlı düşen damla

Yüzey gerilimi sıfır olduğunda orta yoğunluktaki damla alt sıvının yüzeyine çarptıktan sonra geniş biçimde yayılmıştır. Yüzey gerilimi oranı 1,9:1:1'e çıkarıldığında damlanın yayılması önemli ölçüde engellenmiş ve damla iki faz arasındaki arayüzde tutulmuştur.

Üç fazlı yükselen kabarcık

Yoğunluğu en düşük fazın kabarcığı, daha yoğun alt sıvının içinde yükselmiş ve iki yoğun sıvının arayüzüne ulaştığında yüzey gerilimi nedeniyle burada tutulmuştur. Aynı davranış 10.000 parçacıklı 2B Lagrangian ve 64.000 noktalı 3B Eulerian çözümde gösterilmiştir.

Üç fazlı tank karışımı

Yoğunlukları 100, 50 ve 1 olan üç karışmayan faz başlangıçta yan yana yerleştirilmiştir. Yerçekimi altında sistem yeniden düzenlenmiş ve kararlı durumda en yoğun kırmızı faz tabana, orta yoğunluktaki yeşil faz ortaya, en hafif mavi faz üst bölgeye yerleşmiştir.

Eulerian ve Lagrangian Yaklaşımların Farkı Nedir?

Lagrangian SPH'de parçacıklar fiziksel akışla birlikte hareket eder. Bu yaklaşım serbest yüzeyleri doğal biçimde temsil eder ancak her zaman adımında komşuluk listesinin güncellenmesini ve particle shifting yapılmasını gerektirir.

Eulerian kullanımda parçacıklar uzayda sabit approximation points olarak bırakılmış ve fazlar ile hız alanı denklemler üzerinden advect edilmiştir. Bu yaklaşım özellikle 3B örneklerde hesaplama maliyetini azaltmıştır; çünkü komşu parçacık listesi ve shifting algoritması sürekli yeniden oluşturulmamaktadır.

Çalışmanın desteklediği sonuçlar

• CHNS denklemleri SPH-RK/DRK ile tamamen meshfree biçimde çözülebilmiştir.
• GPE, density-based EOS kullanılmadan basınç çözümü sağlamıştır.
• RK ve DRK benchmarklarda analitik/reference çözümlerle iyi uyum göstermiştir.
• DRK bazı testlerde RK'den hafif daha hızlıdır.
• RK bazı benchmarklarda DRK'den hafif daha doğru sonuç vermiştir.
• Liquid-lens sonuçlarının analitik hatası %5'in altında kalmıştır.
• RK ve DRK liquid-lens yakınsaması yaklaşık ikinci mertebedir.
• 2B ve 3B üç faz ayrışması keyfi geometrilerde çözülebilmiştir.
• Yaklaşık 1000:1 yoğunluk oranında kararlı çözüm elde edilmiştir.
• Particle shifting Lagrangian RK/DRK çözümünün kararlılığı için kritik bulunmuştur.
• Phase-Field yaklaşımı SPH'nin hem Eulerian hem Lagrangian kullanılabilmesini sağlamıştır.

Çalışmanın desteklemediği sonuçlar

• DRK'nin bu uygulamada RK'den 10–25 kat daha hızlı olduğu gösterilmemiştir.
• Bütün benchmarklarda DRK'nin RK'den daha doğru olduğu gösterilmemiştir.
• 1000:1 yoğunluk testi deneysel olarak doğrulanmış değildir.
• OpenFOAM karşılaştırması nicel olarak bütün arayüz geometrilerinin aynı olduğunu göstermemektedir.
• Phase-Field yaklaşımının tam yerel hacim korunumu sağladığı gösterilmemiştir.
• Tek bir mobility değerinin bütün üç fazlı problemler için optimum olduğu gösterilmemiştir.
• Simülasyondaki diffuse-interface kalınlığı gerçek sıvı arayüzünün fiziksel kalınlığı olarak yorumlanamaz.
• Yöntemin bütün endüstriyel yüksek Reynolds sayılı türbülanslı üç fazlı akışlarda doğrulandığı gösterilmemiştir.

Kaynak ve Yöntem Notu

Özgün çalışma: Meshfree Phase-Field Modeling of Three-Phase Flow using Smoothed Particle Hydrodynamics with Differential Reproducing Kernels and Artificial Compressibility.

Yazar: Adam Y. Ghoneim.

Kurum: Department of Mechanical Engineering, University of Manitoba, Winnipeg, Manitoba, Canada.

Yayın: Journal of Computational and Applied Mathematics, Volume 469, Article 116654, 2025.

Yayın tarihi: 1 December 2025.

DOI: 10.1016/j.cam.2025.116654.

Hakemlik durumu: Hakemli dergi makalesidir. Yüklenen dosya yayımlanmış makale sürümünün içeriğini taşımakta ve teşekkür bölümünde anonim hakemlerin değerlendirmelerine teşekkür edilmektedir.

Bibliyografik not: PDF'nin ilk sayfasındaki dipnotta derginin adı “Journal of Computational of Applied Mathematics” şeklinde yazılmıştır. Resmî yayınevi kaydındaki doğru dergi adı “Journal of Computational and Applied Mathematics”tir.

Ana model: Cahn–Hilliard–Navier–Stokes üç fazlı diffuse-interface modeli.

Sayısal yöntem: Smoothed Particle Hydrodynamics; Reproducing Kernels ve Differential Reproducing Kernels.

Basınç çözümü: General Pressure Equation tabanlı artificial compressibility.

Zaman integrasyonu: İkinci mertebe Runge–Kutta.

Ana stabilizasyonlar: GPE pressure diffusion, artificial viscosity ve Lagrangian uygulamada particle shifting.

Başlıca doğrulamalar: Couette flow, Poiseuille flow, lid-driven cavity, surface-tension-driven interface relaxation, liquid lens analytical comparison, OpenFOAM FVM–VOF karşılaştırması.

2B çözünürlük örnekleri: 2.500–40.000 parçacık; birçok üç fazlı testte 10.000 parçacık.

3B çözünürlük: Küre, torus, damla ve kabarcık örneklerinde 64.000 parçacık/approximation point.

Yüksek yoğunluk oranı testi: su–zeytinyağı–hava; yaklaşık 1000:1 maksimum yoğunluk oranı.

Yakınsama: sıvı lensi benchmarkında RK ≈1,86 ve DRK ≈1,93 yakınsama eğimi.

Temel sınırlama: Phase-Field mobility ve interface-thickness parametreleri problem bağımlıdır. Yüksek mobility geometrik arayüz kalınlığını koruyabilir ancak hacim kaybını artırabilir; düşük mobility ise aşırı arayüz difüzyonuna yol açabilir.

Doğrulama sınırı: çalışma sayısal/metodolojik bir CFD makalesidir. Çok fazlı gerçek deney verisine karşı kapsamlı doğrulama yapılmamıştır.

Görsel yeniden çizim notu: Kaynaktaki renkli üç fazlı simülasyon görüntülerini doğrudan kopyalamak yerine Verianla için meshfree parçacık ağı, diffuse interface, üç faz, RK/DRK kernel, GPE basınç çözümü, düşen damla ve yükselen kabarcık mekanizmalarını birleştiren özgün bir bilimsel şema kullanılmalıdır.