Channel-flow lab: analytic solution versus numerical methods
แก้ plane Poiseuille flow ระหว่างแผ่นนิ่ง เปรียบเทียบ finite differences และฝึกอ่านความคลาดเคลื่อนของอัตราไหล
เป้าหมายและขอบเขต
จบบทนี้ผู้เรียนจะตรวจหน่วย แก้ความเร็วระหว่างแผ่น และอธิบายได้ว่าความเร็วตรงที่จุดกริดยังไม่ได้ทำให้อัตราไหลจากการอินทิเกรตตรงทั้งหมด แบบฝึกเป็น ปัญหา steady laminar 1D ที่ลดรูปจาก Navier–Stokes ไม่ใช่ solver สนามไหลทั่วไปหรือ CFD ใบพัด
สมมติของไหล Newtonian ไม่อัดตัว ρ และ μ คงที่ แผ่นนิ่งห่าง H กว้างมากเมื่อเทียบกับช่องว่าง ไม่มีผลผนังด้านข้าง สนามพัฒนาเต็มที่ u = u(y) และไม่มีแรงปริมาตรตามแนวไหล พิจารณาบริเวณห่างจากทางเข้าออก ไม่จำลองช่วงพัฒนาของโปรไฟล์ เงื่อนไข no-slip คือ u(0) = u(H) = 0 ตรงกับโจทย์ MIT plane Poiseuille flow
จากแรงดันสู่ความเร็ว
นิยาม Δp = pin − pout ≥ 0 และ L เป็นระยะที่ใช้วัดความดันตก จึงได้ dp/dx = −Δp/L สมการ μd²u/dy² = −Δp/L อินทิเกรตสองครั้งแล้วใส่เงื่อนไขผนังจะได้
u(y) = Δp y(H − y)/(2μL)
umax = ΔpH²/(8μL), Ū = ΔpH²/(12μL), Q′ = ŪH = ΔpH³/(12μL)
Q′ เป็นปริมาตรไหลต่อหน่วยความกว้าง หน่วย m²/s ไม่ใช่ m³/s หากความกว้างจริง W มากพอให้ละเลยขอบข้างได้จึงประมาณ Q ≈ WQ′ ส่วน umax/Ū = 1.5 เป็นคุณสมบัติของโปรไฟล์นี้ ความหนาแน่นไม่ปรากฏในสูตรความเร็วเมื่อกำหนด Δp แต่ยังใช้ตรวจ Reynolds ได้
ตัวอย่างอ้างอิงแม้ปิด JavaScript
ใช้ ρ = 1.225 kg/m³, μ = 0.0000181 Pa·s, H = 0.01 m, L = 1 m และ Δp = 0.01 Pa
| ปริมาณ | ค่าคำตอบวิเคราะห์ |
|---|---|
| ความเร็วสูงสุดที่ y = H/2 | 0.00690608 m/s |
| ความเร็วเฉลี่ย Ū | 0.00460405 m/s |
| Q′ ต่อหน่วยความกว้าง | 0.0000460405 m²/s |
| Dh ≈ 2H | 0.02 m |
| ReDh = ρŪ(2H)/μ | ประมาณ 6.232 |
ค่า Re ต่ำนี้สอดคล้องกับการเลือกตัวอย่างที่ความหนืดเด่น แต่การเพิ่ม Δp จนสูงมากไม่ทำให้คำตอบพาราโบลาเป็นคำตอบทางกายภาพที่เชื่อได้เสมอ เครื่องคำนวณยังแสดงสูตรเดิม จึงต้องประเมินสมมติฐานใหม่เอง
วิธีเชิงตัวเลขที่แก้อย่างอิสระ
แบ่ง H เป็น N ช่วงเท่ากัน มี N + 1 จุด ใช้ h = H/N และจุดภายใน i = 1,…,N − 1 ประมาณอนุพันธ์อันดับสองด้วย central difference:
ui−1 − 2ui + ui+1 = −Δp h²/(μL)
นำค่าผนังศูนย์เข้าสู่ระบบสมการสามแนวทแยง แล้วแก้ค่าภายในโดยไม่แทนสูตรพาราโบลาเป็นคำตอบเชิงตัวเลข จากนั้นเปรียบเทียบแต่ละจุดกับสูตรวิเคราะห์ เครื่องคำนวณและสมุดงานใช้วิธีนี้เพื่อให้ตรวจทั้งการตั้งระบบและการแก้ระบบได้
มีจุดน่าสังเกต: อนุพันธ์อันดับสองแบบ central difference ให้ค่าตรงสำหรับพหุนามกำลังสองบนกริดสม่ำเสมอ กรณีนี้จึงอาจมี error ที่จุดเพียงระดับปัดเศษ แม้กริดหยาบ ไม่ใช่หลักฐานว่า solver ทั่วไปแม่นเท่านี้ การหาปริมาตรไหลด้วยกฎคางหมูยังมี error โดย Q′trap = Q′(1 − 1/N²) สำหรับโปรไฟล์และกริดนี้
| N ช่วง | ความคลาดเคลื่อนสัมพัทธ์ของ Q′ จากคางหมู |
|---|---|
| 10 | 1% ต่ำกว่าค่าจริง |
| 20 | 0.25% ต่ำกว่าค่าจริง |
| 40 | 0.0625% ต่ำกว่าค่าจริง |
ที่ N = 20 จะได้ Q′trap ≈ 0.0000459254 m²/s การเพิ่ม N สองเท่าลด error ของการอินทิเกรตประมาณสี่เท่า นี่เป็นผลคำนวณของแบบฝึกนี้ ไม่ใช่ผลทดลอง NASA หรือ MIT และยังไม่ใช่ validation กับของไหลจริง ตามความแตกต่างที่ NASA อธิบายใน verification
ลงมือทดลองและส่งชิ้นงาน
- เปิดเครื่องคำนวณท้ายบท คำนวณค่าเริ่มต้นและตรวจหน่วยตามตาราง
- เพิ่ม Δp สองเท่า แล้วตรวจว่า u และ Q′ เพิ่มสองเท่า
- เปลี่ยน N จาก 10 เป็น 20 และ 40 เก็บ CSV และเปรียบเทียบ error ของ Q′
- ดาวน์โหลด สมุดงาน Python และไฟล์อ้างอิง อ่าน README แล้วรัน solver และ tests ตามคำสั่งในนั้น เก็บเวอร์ชัน Python พร้อมผลที่ได้
โจทย์: หากเพิ่ม H สองเท่า แต่ตรึง μ, L และ Δp ความเร็วเฉลี่ยกับ Q′ เปลี่ยนอย่างไร?
แนวคำตอบ: Ū เพิ่มสี่เท่า และ Q′ เพิ่มแปดเท่า ส่วน ReDh เพิ่มแปดเท่า จึงต้องกลับไปตรวจขอบเขต laminar และเงื่อนไขช่องกว้างอีกครั้ง ก่อนตีความว่าเป็นผลของช่องจริง
ต่อไป: จากแบบฝึกหนึ่งมิติสู่งาน CFD ที่ตรวจสอบได้
บรรณานุกรม
ตรวจแหล่งต้นทาง 26 กันยายน 2569
ห้องทดลองการไหลระหว่างแผ่นขนาน
แบบจำลองคงตัว 1 มิติ: แผ่นอยู่นิ่ง การไหลพัฒนาเต็มที่ ลามินาร์ ไม่อัดตัว ไม่มีแรงปริมาตรในแนวไหล ไม่ใช่ CFD ของใบพัด ข้อมูลคำนวณในเบราว์เซอร์ ไม่ส่งออกจากเครื่อง