Konservasi kuantitas fisik ketika menggunakan kondisi batas Neumann diterapkan pada persamaan adveksi-difusi

25

Saya tidak mengerti perilaku berbeda dari persamaan advection-difusion ketika saya menerapkan kondisi batas yang berbeda. Motivasi saya adalah simulasi kuantitas fisik nyata (kepadatan partikel) di bawah difusi dan adveksi. Kepadatan partikel harus dilestarikan di interior kecuali mengalir keluar dari tepi. Dengan logika ini, jika saya menerapkan syarat batas Neumann pada ujung sistem seperti (di sisi kiri dan kanan) maka sistem harus "ditutup" yaitu jika fluks pada batas adalah nol maka tidak ada partikel yang bisa lolos.ϕx=0

Untuk semua simulasi di bawah ini, saya telah menerapkan diskritisasi Crank-Nicolson ke persamaan advection-difusion dan semua simulasi memiliki kondisi batas. Namun, untuk baris pertama dan terakhir dari matriks (baris kondisi batas) saya izinkan diubah secara independen dari nilai interior. Ini memungkinkan titik akhir sepenuhnya tersirat.βϕx=0β

Di bawah ini saya membahas 4 konfigurasi berbeda, hanya satu di antaranya yang saya harapkan. Pada akhirnya saya membahas implementasi saya.

Batas difusi saja

Di sini istilah advection dimatikan dengan mengatur kecepatan ke nol.

Hanya difusi, dengan β = 0,5 (Crank-Niscolson) di semua titik

Hanya difusi (Batas Neumann dengan beta = 0,5)

Kuantitas tidak dikonservasi seperti yang terlihat oleh pengurangan area pulsa.

Hanya difusi, dengan = 0,5 (Crank-Niscolson) pada titik interior, dan = 1 (implisit penuh) pada batasβββ

Hanya difusi (batas Neumann dengan beta = 0,5 untuk interior, beta = 1 sepenuhnya tersirat) batas

Dengan menggunakan persamaan sepenuhnya implisit pada batas saya mencapai apa yang saya harapkan: tidak ada partikel yang lolos . Anda dapat melihat ini dengan area yang dikonservasi karena partikelnya berdifusi. Mengapa pilihan pada titik batas mempengaruhi fisika situasi? Apakah ini bug atau yang diharapkan?β

Difusi dan adveksi

Ketika istilah advection dimasukkan, nilai pada batas tampaknya tidak mempengaruhi solusi. Namun, untuk semua kasus ketika batas tampaknya "terbuka" yaitu partikel bisa lolos dari batas. Mengapa demikian?β

Adveksi dan Difusi dengan = 0,5 (Crank-Niscolson) di semua titikβ

Advection-Difusion (Batas Neumann dengan beta = 0,5)

Adveksi dan Difusi dengan = 0,5 (Crank-Niscolson) pada titik interior, dan = 1 (implisit penuh) pada batasβββ

Adveksi dan Difusi (batas Neumann dengan beta = 0,5 untuk interior, beta = 1 sepenuhnya implisit) batas

Penerapan persamaan advection-difusion

Dimulai dengan persamaan advection-difusion,

ϕt=D2ϕx2+vϕx

Menulis menggunakan Crank-Nicolson memberi,

ϕjn+1ϕjnΔt=D[1β(Δx)2(ϕj1n2ϕjn+ϕj+1n)+β(Δx)2(ϕj1n+12ϕjn+1+ϕj+1n+1)]+v[1β2Δx(ϕj+1nϕj1n)+β2Δx(ϕj+1n+1ϕj1n+1)]

Perhatikan bahwa = 0,5 untuk Crank-Nicolson, = 1 untuk sepenuhnya implisit, dan, = 0 untuk sepenuhnya eksplisit.βββ

Untuk menyederhanakan notasi, mari kita lakukan substitusi,

s=DΔt(Δx)2r=vΔt2Δx

dan pindahkan nilai yang diketahui dari turunan waktu ke sisi kanan,ϕjn

ϕjn+1=ϕjn+s(1β)(ϕj1n2ϕjn+ϕj+1n)+sβ(ϕj1n+12ϕjn+1+ϕj+1n+1)+r(1β)(ϕj+1nϕj1n)+rβ(ϕj+1n+1ϕj1n+1)

Anjak istilah memberi,ϕ

β(rs)ϕj1n+1+(1+2sβ)ϕjn+1β(s+r)ϕj+1n+1Aϕn+1=(1β)(sr)ϕj1n+(12s[1β])ϕjn+(1β)(s+r)ϕj+1nMϕn

yang dapat kita tulis dalam bentuk matriks sebagai mana,Aϕn+1=Mϕn

A=(1+2sββ(s+r)0β(rs)1+2sββ(s+r)β(rs)1+2sββ(s+r)0β(rs)1+2sβ)

M=(12s(1β)(1β)(s+r)0(1β)(sr)12s(1β)(1β)(s+r)(1β)(sr)12s(1β)(1β)(s+r)0(1β)(sr)12s(1β))

Menerapkan ketentuan batas Neumann

NB bekerja melalui derivasi lagi saya pikir saya telah melihat kesalahan. Saya mengasumsikan skema sepenuhnya implisit ( = 1) ketika menulis perbedaan kondisi batas yang terbatas. Jika Anda menganggap skema Crank-Niscolson di sini kompleksitasnya menjadi terlalu besar dan saya tidak bisa menyelesaikan persamaan yang dihasilkan untuk menghilangkan node yang berada di luar domain. Namun, tampaknya mungkin, ada dua persamaan dengan dua yang tidak diketahui, tetapi saya tidak bisa mengelolanya. Ini mungkin menjelaskan perbedaan antara plot pertama dan kedua di atas. Saya pikir kita dapat menyimpulkan bahwa hanya plot dengan = 0,5 pada titik batas yang valid.ββ

Dengan asumsi fluks di sisi kiri diketahui (dengan asumsi bentuk sepenuhnya implisit),

ϕ1n+1x=σL

Menulis ini sebagai perbedaan terpusat memberi,

ϕ1n+1xϕ2n+1ϕ0n+12Δx=σL

oleh karena itu, ϕ0n+1=ϕ2n+12ΔxσL

Perhatikan bahwa ini memperkenalkan simpul yang berada di luar domain masalah. Node ini dapat dihilangkan dengan menggunakan persamaan kedua. Kita dapat menulis simpul sebagai,ϕ0n+1j=1

β(rs)ϕ0n+1+(1+2sβ)ϕ1n+1β(s+r)ϕ2n+1=(1β)(sr)ϕj1n+(12s[1β])ϕjn+(1β)(s+r)ϕj+1n

Mengganti nilai ditemukan dari kondisi batas memberikan hasil sebagai berikut untuk = 1 baris,ϕ0n+1j

(1+2sβ)ϕ1n+12sβϕ2n+1=(1β)(sr)ϕj1n+(12s[1β])ϕjn+(1β)(s+r)ϕj+1n+2β(rs)ΔxσL

Melakukan prosedur yang sama untuk baris terakhir (at = ) menghasilkan,jJ

2sβϕJ1n+1+(1+2sβ)ϕJn+1=(1β)(sr)ϕJ1n+(12s(1β))ϕJn+2β(s+r)ΔxσR

Akhirnya membuat baris batas tersirat (pengaturan = 1) memberi,β

(1+2s)ϕ1n+12sϕ2n+1=ϕj1n+1ϕjn+2(rs)ΔxσL

2sϕJ1n+1+(1+2s)ϕJn+1=ϕJn+2(s+r)ΔxσR

Karenanya dengan syarat batas Neumann kita dapat menulis persamaan matriks, ,Aϕn+1=Mϕn+bN

dimana,

A=(1+2s2s0β(rs)1+2sββ(s+r)β(rs)1+2sββ(s+r)02s1+2s)

M=(100(1β)(sr)12s(1β)(1β)(s+r)(1β)(sr)12s(1β)(1β)(s+r)001)

bN=(2(rs)ΔxσL002(s+r)ΔxσR)T

Pemahaman saya saat ini

  • Saya pikir perbedaan antara plot pertama dan kedua dijelaskan dengan memperhatikan kesalahan yang diuraikan di atas.

  • Mengenai konservasi kuantitas fisik. Saya percaya penyebabnya adalah bahwa, seperti yang ditunjukkan di sini , persamaan adveksi dalam bentuk yang saya tulis tidak memungkinkan propagasi dalam arah sebaliknya sehingga gelombang hanya melewati bahkan dengan kondisi batas nol-fluks . Intuisi awal saya tentang konservasi hanya diterapkan ketika jangka waktu adveksi adalah nol (ini adalah solusi dalam plot nomor 2 di mana area tersebut dilestarikan).

  • Bahkan dengan kondisi batas nol-fluks Neumann massa masih dapat meninggalkan sistem, ini karena kondisi batas yang benar dalam kasus ini adalah kondisi batas Robin di mana total fluks ditentukan . Selain itu, kondisi Neunmann menentukan bahwa massa tidak dapat meninggalkan domain melalui difusi , ia tidak mengatakan apa pun tentang adveksi. Pada dasarnya apa yang kita dengar adalah kondisi batas tertutup untuk difusi dan kondisi batas terbuka untuk kemajuan. Untuk informasi lebih lanjut lihat jawabannya di sini, Penerapan kondisi batas gradien nol dalam persamaan adveksi-difusiϕx=0j=Dϕx+vϕ=0.

Apakah kamu setuju?

Boyfarrell
sumber
Tampaknya syarat batas tidak diterapkan dengan benar. Bisakah Anda menunjukkan kepada kami bagaimana Anda telah memberlakukan persyaratan batas?
David Ketcheson
OK saya memperbarui dengan implementasi dan saya pikir saya melihat kesalahan mengenai penerapan = 0,5 pada baris batas saja. Saya telah memperbarui "pemahaman saya saat ini" di bagian bawah pertanyaan. Apakah Anda punya komentar? β
boyfarrell
Jadi ... bagaimana diskritisasi batas-batas itu terlihat dalam kasus batas Robin? Anda telah menunjukkannya untuk batas Neumann, tetapi bukan batas Robin.

Jawaban:

15

Saya pikir salah satu masalah Anda adalah (seperti yang Anda amati dalam komentar Anda) Kondisi Neumann bukanlah kondisi yang Anda cari , dalam arti bahwa itu tidak menyiratkan konservasi kuantitas Anda. Untuk menemukan kondisi yang benar, tulis ulang PDE Anda sebagai

ϕt=x(Dϕx+vϕ)+S(x,t).

Sekarang, istilah yang muncul dalam tanda kurung, adalah fluks total dan ini adalah jumlah yang harus Anda masukkan ke nol pada batas yang harus dilestarikan . (Saya telah menambahkan untuk kepentingan umum dan untuk komentar Anda.) Maka syarat batas yang harus Anda terapkan adalah (seandainya domain ruang angkasa Anda adalah )Dϕx+vϕ=0ϕS(x,t)(10,10)

Dϕx(10)+vϕ(10)=0

untuk sisi kiri dan

Dϕx(10)+vϕ(10)=0

untuk sisi kanan. Ini adalah kondisi batas Robin (perhatikan bahwa Wikipedia secara eksplisit mengatakan ini adalah kondisi isolasi untuk persamaan advection-difusion).

Jika Anda menetapkan kondisi batas ini, Anda mendapatkan properti konservasi yang Anda cari. Memang, berintegrasi dengan domain luar angkasa, kami miliki

ϕtdx=x(Dϕx+vϕ)dx+S(x,t)dx

Menggunakan integrasi oleh bagian-bagian di sisi kanan, kami punya

ϕtdx=(Dϕx+vϕ)(10)(Dϕx+vϕ)(10)+S(x,t)dx

Sekarang, dua istilah sentral lenyap berkat kondisi batas. Integrasi dalam waktu, kami dapatkan

0Tϕtdxdt=0TS(x,t)dxdt

dan jika kita diizinkan untuk beralih dua integral pertama ,

ϕ(x,T)dxϕ(x,0)dx=0TS(x,t)dx

Ini menunjukkan bahwa domain terisolasi dari eksterior. Secara khusus, jika , kita mendapatkan konservasi .S=0ϕ

Dr_Sam
sumber
Saya menyadari sekarang mengapa ini hanya bekerja ketika = 0; karena ini akan menyiratkan konservasi mengikuti pendekatan Anda di atas. Apa yang akan menjadi konsekuensi dari penggunaan kondisi batas ini di atas, apakah gelombang akan mencerminkan? Saya pikir ini tidak mungkin karena tidak ada dalam persamaan yang akan memberi saya kecepatan negatif? v
boyfarrell
Cara terbaik untuk mengetahui mungkin untuk mencoba! Tetapi jika ini berperilaku benar (dan IMO itu), Anda akan melihat sejumlah yang mulai terakumulasi di sisi kiri domain: advection mendorong ke arah itu tetapi batas ditutup. Akumulasi berhenti ketika difusi cukup besar untuk menyeimbangkannya. Jadi tidak, seharusnya tidak ada gelombang yang dipantulkan. ϕϕ
Dr_Sam
@DrSam Hanya pertanyaan tentang implementasi. Saya mengerti maksud Anda tentang bagaimana membuat kuantitas nol di sisi kiri. Tetapi ketika Anda mengatakan "di sebelah kanan hanya istilah batas" apa artinya itu? Saya pikir syarat batas haruslah Neumann atau Dirichlet (atau gabungan keduanya)?
boyfarrell
@boyfarrell Kiri / kanan dalam jawaban mengacu pada derivasi kondisi batas yang benar, bukan cara penerapannya (diedit untuk kejelasan). Kondisi Robin adalah kondisi klasik, meskipun ada yang kurang dikenal daripada Dirichlet dan Neumann.
Dr_Sam
Jadi mengenai implementasi, apakah Anda pikir saya harus mendapatkan syarat batas Robin untuk kedua batas? Juga, jika persamaan memiliki istilah reaksi (misalnya apakah syarat batas juga perlu mencakup istilah ini?
ϕt=x(Dϕx+vϕ)+S(x,t)
boyfarrell