จากผลอากาศพลศาสตร์สู่ DroneSim: เตรียมตารางแรง

ส่งปริมาณที่มีความหมายเข้าสู่การจำลอง

จบบทนี้ผู้เรียนจะเตรียมสัญญาข้อมูลแรงและเขียนแผนตรวจส่วนเชื่อมต่อได้ Simulator การเคลื่อนที่ต้องการแรงและโมเมนต์ในกรอบพิกัดที่แน่นอน ส่วน CFD หรือข้อมูลทดลองอาจให้ค่าสัมประสิทธิ์ในคนละนิยาม การแปลงผิดหน่วยหรือกลับเครื่องหมายทำให้ลำเร่งผิดทิศได้ แม้ข้อมูลต้นทางคำนวณถูก

บทนี้มี ส่วนเชื่อมตารางแรงต้านที่ติดตั้งและทดสอบใน DroneSim แล้ว พร้อมตัวอย่างสังเคราะห์เพื่อพิสูจน์เส้นทางข้อมูล ส่วนสูตรใบพัดและแบบบันทึกยังใช้เป็นพื้นฐานแยกจาก adapter นี้ ไม่มีชุดแรงที่ผ่าน validation ให้ใช้งานบิน และไม่ได้ใช้ผลจานกำหนดแรงจากบท CFD เป็นแรงต้านลำ

แผนภาพนี้เป็นกระบวนการทั่วไปสำหรับข้อมูลอากาศพลศาสตร์ในอนาคต รวมถึงข้อมูลใบพัดที่อาจต้องใช้ RPM ส่วน adapter ที่ทำงานจริงในบทนี้รับเฉพาะ airframe drag แบบ fixed-direction และ fixed-density โดยไม่ใช้ RPM

ห้องทดลองต่อยอด: จากตารางแรงสู่การเคลื่อนที่

เริ่มส่วนนี้หลังอ่านสามกรณี ใกล้พื้น, รอบอาคาร และ สิ่งบรรทุก คุณจะได้ฝึกตัดสินใจว่าแรงที่มีใช้กับโมเดลใดได้ ตรวจการเคลื่อนที่ต่อเนื่อง และเตรียมข้อมูลสำหรับเทียบการวัดจริง

ดาวน์โหลด สมุดงานตรวจและส่งต่อแรง มีเครื่องมือตรวจหลักฐาน Python ตัวอย่างเทียบค่าพร้อมความไม่แน่นอน และโค้ดสาธิตวิถีต่อเนื่องสำหรับ DroneSim พร้อมผลที่รันไว้ Python ใช้ได้บน Windows/Linux โดยไม่ติดตั้งไลบรารีเพิ่ม ส่วนการรัน DroneSim ต้องใช้โครงการและส่วนเชื่อมตาม README

แรงเป็นของใคร: ตรวจให้ถูกก่อนแปลงแกน

ทั้ง 14 กรณีผ่านเกณฑ์จบเชิงตัวเลขของชุดฝึก แต่ ยังไม่มีข้อมูลวัดที่ยืนยันว่าตรงกับรูปทรง สมบัติของไหล และขอบเขตของกรณีเหล่านี้ เครื่องมือในสมุดงานจึงเก็บผลไว้พร้อมเหตุผลที่ยังส่งต่อไม่ได้ ไม่เปลี่ยนป้ายให้เป็นแรงโดรน นี่เป็นการใช้ข้อมูลอย่างถูกความหมายตามหลัก NASA Verification and Validation

การเคลื่อนที่ต่อเนื่องที่ตรวจคำตอบได้

ตัวอย่างใหม่ใช้ ข้อมูลสังเคราะห์ กับ stepDynamics ตัวจริงของ DroneSim ต่อสถานะจากก้าวก่อนหน้าโดยไม่ย้ายลำกลับไปจุดเริ่มต้น ตั้งลำระดับที่ความสูง 10 m ความเร็วแนวนอนเริ่มต้น 5 m/s และตั้งรอบมอเตอร์เริ่มต้นให้แรงยกสมดุลน้ำหนัก คงแรงดันแบตเตอรี่ไว้ แบบฝึกนี้จึงศึกษาแรงต้านตามแนวนอนในสภาพที่จำกัด ไม่ใช่ภารกิจบินที่มีผู้ควบคุมหรือลมกระโชก

ตารางตัวอย่างในช่วง 0–5 m/s กำหนด Fₓ = −kvₓ โดย k = 0.06 N·s/m เครื่องหมายลบหมายถึงแรงต้านทิศการเคลื่อนที่ เมื่อนำไปแทน m dvₓ/dt = Fₓ จะได้คำตอบวิเคราะห์:

vₓ(t) = v₀ exp(−kt/m) และ x(t) − x₀ = (mv₀/k)[1 − exp(−kt/m)]

m คือมวล kg, t คือเวลา s, v₀ = 5 m/s และ x คือระยะทาง m ค่า kt/m ไม่มีหน่วย สูตรนี้ใช้ได้ในช่วงตารางที่เป็นเส้นตรงดังกล่าวเท่านั้น ตารางช่วง 5–10 m/s มีความชันต่างกันจึงห้ามใช้สูตรเดิมข้ามช่วงโดยไม่พิจารณาใหม่

ลองเปรียบเทียบก้าวเวลา 0.002, 0.001 และ 0.0005 s กับคำตอบนี้ เพื่อแยกความคลาดเคลื่อนของการเดินเวลาออกจากความถูกต้องของฟิสิกส์ ดูค่า x และ v ที่เวลาเดียวกัน ไม่เปรียบเทียบแถวลำดับเดียวกันซึ่งอาจหมายถึงคนละเวลา ส่วนแรงต้านเดิมของ DroneSim เป็นอีกโมเดลหนึ่ง จึงไม่ควรคาดหวังวิถีเท่ากับตารางสังเคราะห์

กราฟซ้ายแสดงผลรันต่อเนื่องจริงจากโค้ด แต่แรงต้านที่ป้อนเป็นค่าตั้งขึ้นเพื่อทดสอบ กราฟขวาแสดงผลต่างจากคำตอบวิเคราะห์ ไม่ใช่ error จากการบินจริง เปิดกราฟขนาดเต็ม

ผลที่เวลา 2 s ใช้มวลจากลำตัวอย่าง 0.5043 kg และก้าวเวลา 0.0005 s:

ความต่างนี้เกิดจากการเลือกโมเดลแรงต้าน ไม่ได้บอกว่าแบบใดตรงกับโดรนจริงมากกว่า ในไฟล์ dronesim/result/trajectory.csv มีทุกก้าวเวลาให้ตรวจซ้ำ ส่วน JSON บันทึกมวล การตั้ง trim และผลตรวจเงื่อนไข ข้อจำกัดของตารางยังเหมือนเดิม: ทิศความเร็วสัมพัทธ์ใน body frame คงที่ ความหนาแน่นคงที่ และอัตราหมุนลำเป็นศูนย์ หากหลุดเงื่อนไข ส่วนเชื่อมจะหยุดพร้อมข้อผิดพลาด

เตรียมเทียบค่าที่วัดโดยไม่เติมความแม่นยำที่ยังไม่มี

ก่อนเปรียบเทียบ ต้องตรงกันทั้งวัตถุ รูปทรง ปริมาณที่วัด แกน จุดอ้างอิงโมเมนต์ ความหนาแน่น ความหนืด ลมเข้า ขอบเขต และวิธีเฉลี่ยเวลา เครื่องมือจะปฏิเสธคู่ที่เงื่อนไขต่างกันหรือยังไม่ทราบความไม่แน่นอน การกำหนด ID ให้ตรงกันเป็นเพียงสัญญาข้อมูล ผู้ใช้ยังต้องตรวจหลักฐานจริงประกอบ

ให้ S เป็นค่าจำลองและ D เป็นค่าที่วัดในหน่วยเดียวกัน ผลต่างคือ E = S − D หาก uₛ และ uᴅ เป็น ความไม่แน่นอนมาตรฐาน ที่ประเมินแล้ว และแหล่งความไม่แน่นอนเป็นอิสระต่อกัน จะรวมได้เป็น u꜀ = √(uₛ² + uᴅ²) ตาม NIST TN 1297 หากมีความสัมพันธ์ต้องพิจารณาพจน์ covariance เพิ่ม เครื่องมือตัวอย่างไม่รับกรณีนั้น

ตัวอย่างฝึกที่ สมมติทั้งหมด: S = 0.30 N, D = 0.28 N, uₛ = 0.01 N และ uᴅ = 0.02 N จะได้ E = 0.02 N, u꜀ ≈ 0.02236 N และ |E|/u꜀ ≈ 0.894 ตัวเลขสุดท้ายบอกขนาดผลต่างเทียบกับความไม่แน่นอนรวม ไม่ใช่คะแนนผ่านการรับรอง และไม่ได้พิสูจน์ว่าโมเดลใช้ได้กับทุกสภาพ

ลองทำ: เปลี่ยนความหนืดของไฟล์ตัวอย่างเพียงด้านเดียว แล้วตรวจว่าระบบปฏิเสธการเปรียบเทียบ จากนั้นคืนค่าและลบความไม่แน่นอนออก ระบบควรให้แก้ข้อมูล ไม่แทนช่องว่างด้วยศูนย์ อย่านำผลต่างสองกริดจากเฟส 2 มาใส่เป็น uₛ โดยตรง เพราะยังไม่มีการประเมิน error bound ตาม NASA Spatial Convergence

ในไฟล์เปรียบเทียบ flow_velocity_m_s หมายถึงลมสัมพัทธ์ที่ไหลผ่านวัตถุ ส่วน velocity_direction_frd ของ DroneSim หมายถึงตัวลำเคลื่อนผ่านอากาศ สองนิยามนี้กลับเครื่องหมายกัน เช่น ลม +x ผลักวัตถุหยุดนิ่งไป +x แต่ลำเคลื่อน +x ในอากาศนิ่งจะถูกแรงต้านไป −x

ส่วนถัดไปอธิบายสัญญาข้อมูล แกน และโค้ดส่วนเชื่อมเดิมอย่างละเอียด ส่วนงานใบพัดท้ายบทเป็นอีกกรณีศึกษาหนึ่ง มีสถานะหลักฐานแยกจากชุดทดลองนี้

นิยามที่ต้องตรงกันก่อนสร้างตาราง

หากต้องการฝึกกับข้อมูลวัดจริงแทนตัวเลขสมมติ อ่าน ข้อมูลทดลอง NASA MTB2 เพื่อเลือก Run/Point ตรวจหน่วยและสภาพอุโมงค์ลม ชุดนั้นเป็นแรง hub ของโรเตอร์ จึงยังไม่ใช่ตารางแรงต้านลำสำหรับ adapter นี้

ให้ n เป็นรอบต่อวินาที (ไม่ใช่ rpm), D เป็นเส้นผ่านศูนย์กลาง m, V เป็นความเร็วตามแกนอ้างอิง m/s, T เป็นแรงขับ N และ P เป็นกำลัง W ใช้นิยามมาตรฐาน J = V/(nD), C_T = T/(ρn²D⁴), C_P = P/(ρn³D⁵) ดังนั้น T = C_Tρn²D⁴ และ P = C_Pρn³D⁵

เมื่อ n > 0 โมเมนต์บิดเพลา Q = P/(2πn) หน่วย N·m สัญลักษณ์ Q ตรงนี้คือ torque ต่างจากอัตราไหล Q ในบทช่องไหล ต้องระบุชื่อคอลัมน์ให้ชัด สูตรที่หารด้วย n ไม่ใช้ตรง ๆ ที่รอบศูนย์

UIUC Propeller Database เป็นแหล่งข้อมูลอุโมงค์ลมของใบพัดขนาดเล็ก โดยนิยาม Reynolds ในกราฟ static ใช้คอร์ดและความเร็วหมุนที่สถานี 75% ของรัศมี จึงไม่ควรนำไปเทียบ Re ที่ใช้เส้นผ่านศูนย์กลางโดยไม่มีการแปลง บทนี้เชื่อมไปยังข้อมูลเจ้าของโดยตรง ไม่คัดลอกฐานข้อมูลหรืออ้างสิทธิ์เผยแพร่ต่อโดยอัตโนมัติ

ตัวอย่างแปลงหน่วยที่ติดป้ายว่าสมมติ

กำหนด n = 100 รอบ/s = 6,000 rpm, D = 0.2 m, ρ = 1.225 kg/m³, C_T = 0.10 และ C_P = 0.05 จะได้ T = 1.96 N, P = 19.6 W และ Q ≈ 0.03119 N·m ค่าเหล่านี้ตั้งขึ้นเพื่อฝึกสูตร ไม่ใช่ข้อมูลใบพัดรุ่นใดหรือผลจาก UIUC

หาก V = 5 m/s จะได้ J = 0.25 แต่ไม่ควรใช้ C_T เดิมโดยอัตโนมัติ เพราะสัมประสิทธิ์อาจเปลี่ยนตาม J และ Reynolds กำลัง P ในสูตรเป็นกำลังเพลา ไม่ใช่กำลังไฟจากแบตเตอรี่ ต้องมีโมเดลมอเตอร์/ESC และประสิทธิภาพเพิ่มเติมจึงเชื่อมกับการใช้พลังงานได้

สัญญาข้อมูลที่ต้องตกลงกับปลายทาง

ตารางใน สมุดงานช่องไหล ยังเป็นแบบบันทึกสำหรับวางแผน ส่วน adapter ที่นำไปใช้จริงด้านล่างกำหนดสัญญาเฉพาะของแรงต้านลำ พร้อมแยก provenance เป็น synthetic, simulation หรือ experiment; ป้ายที่มาเหล่านี้ไม่ใช่ใบรับรองความแม่นยำ

ส่วนเชื่อมที่ทำงานจริงและขอบเขตข้อมูล

บน baseline DroneSim commit ต่อไปนี้:

b1e26271a08090c278e6c6a45ed64e9f7f7a0cec

เพิ่ม AirframeForceProvider เป็น อาร์กิวเมนต์ลำดับที่หกแบบ optional ของ stepDynamics ผู้เรียกเดิมยังใช้ drag model เดิม เมื่อส่ง provider เข้ามา จะใช้แรงและโมเมนต์จากตารางแทน params.dragCoef เฉพาะ airframe drag เท่านั้น แรงขับ rotor, reaction torque, แรงโน้มถ่วง และระบบมอเตอร์ยังอยู่ในโมเดลเดิม ไม่มี UI toggle หรือการโหลดตารางอัตโนมัติ

สร้าง provider ด้วยคำสั่ง:

createAerodynamicProvider(data, { densityKgM3, windWorldMps })

ฟังก์ชันตรวจ schema, หน่วย SI, ที่มา และ force_scope ที่เป็น airframe_drag_only ก่อนสร้าง provider ตารางใช้กรอบ FRD: x ไปหน้า, y ไปขวา, z ลง ส่วน body frame ของ DroneSim เป็น FLU: x ไปหน้า, y ไปซ้าย, z ขึ้น จึงหมุนแรงและโมเมนต์ด้วย (x, −y, −z) ก่อนแปลงแรงไป world frame

reference_from_cg_m เป็นเวกเตอร์จาก CG ไปยังจุดอ้างอิงโมเมนต์ของตาราง adapter คำนวณ M_CG = M_reference + r × F หลังแปลงเข้ากรอบ FLU ต้องตรวจจุดอ้างอิงกับต้นทางจริงเพื่อไม่ให้เกิดโมเมนต์ซ้ำ

ตารางรุ่นนี้รองรับความเร็วหนึ่งมิติใน ทิศการเคลื่อนที่ผ่านอากาศคงที่ใน body frame, ความหนาแน่นคงที่ และ body angular rates เป็นศูนย์ โดยหัก wind world ออกจากความเร็วลำก่อนแปลงแกน ทิศนี้ตรงข้ามกับทิศลมที่พุ่งเข้าหาลำ ใช้ linear interpolation ภายในช่วงความเร็วเท่านั้น ไม่ extrapolate และไม่ปรับ density หรือมุมให้อัตโนมัติ อินพุตนอกช่วง/ผิดเงื่อนไขจะ throw แทน silent fallback

สาธิตหนึ่งก้าวเวลาและดาวน์โหลดโค้ด

ดาวน์โหลด สมุดงาน DroneSim adapter พร้อมโค้ด test, patch และผล JSON ภายในมีสำเนาไฟล์ original adapter ห้าไฟล์, patch เฉพาะ dynamics.ts, baseline และ SHA256 manifest ไม่รวมทั้ง simulator หรือการแก้ไขอื่นในโครงการ อ่าน README.md เพื่อ apply กับ checkout ใหม่ที่ baseline ตรงกัน แล้วจาก app รัน:

pnpm exec vitest run --maxWorkers=2
pnpm typecheck
pnpm build
New-Item -ItemType Directory -Force output
pnpm exec tsx scripts/aerodynamic-demo.ts output/aerodynamic-demo-reproduced.json

ตั้งชื่อผลลัพธ์ใหม่ทุกครั้ง เพราะ demo ไม่เขียนทับไฟล์เดิม ผลที่บันทึกใช้ตาราง synthetic ที่ 0, 5 และ 10 m/s ให้ Fx = 0, −0.3 และ −1.2 N ตามลำดับ แต่ละการเปรียบเทียบเริ่มใหม่ที่ความสูง 10 m ไม่ได้ต่อกันเป็น trajectory:

กำหนด vx เริ่มต้น 5 m/s และ Δt = 0.001 s เท่ากันทั้งสองกรณี:

การรัน integration ผ่าน 650 tests ด้วย --maxWorkers=2, typecheck, build และ lint เฉพาะ TypeScript ที่แก้ ผลเหล่านี้ตรวจ software behavior รวมถึงแกน โมเมนต์ interpolation และการปฏิเสธข้อมูลผิดเงื่อนไข ไม่ใช่การ validate การบิน แรงโน้มถ่วงทำให้ทิศความเร็วเปลี่ยนได้ทันที จึงไม่ควรใช้ตาราง fixed-direction นี้กับ trajectory ทั่วไป

CFD กับ adapter เป็นสองการทดลองแยกกัน: source 0.08 N ของบทก่อนเป็นแรง rotor-like ที่กำหนดให้ของไหล ไม่ใช่แรงต้าน airframe ที่ solver ทำนาย การส่งเข้า adapter จะทำให้ ownership ของแรงผิด ก่อนเชื่อมข้อมูลจริงต้องมีแรงต้านลำที่ geometry, density, ทิศ, จุดอ้างอิง และช่วงเงื่อนไขตรงกันโดยไม่มี rotor หรือ gravity รวมอยู่

แผนตรวจการเชื่อมต่อก่อนเชื่อผล

เริ่มทดสอบการแปลง rpm เป็น rad/s และ n แล้วป้อนแรงสังเคราะห์ที่ทราบทิศเพื่อเช็กแกน ทดสอบแรงที่เยื้องศูนย์มวลเพื่อเช็ก r × F และแยกโมเมนต์ปฏิกิริยาจาก torque เพลา ถัดมาทดสอบ interpolation ภายในช่วง และตรวจว่าอินพุตนอกช่วงถูกจัดการตามสัญญา รวมถึงกรณีรอบศูนย์

เมื่อมีข้อมูลจริงจึงเปรียบเทียบแรงกับการทดลองที่เงื่อนไขตรงกัน ประเมินความไวต่อความหนาแน่น รอบและวิธีเฉลี่ย แล้วค่อยทดสอบการเคลื่อนที่รวม การที่ FLOWUnsteady รองรับงานอากาศพลศาสตร์ตามเวลาไม่ได้หมายความว่าผลใด ๆ ใช้กับ DroneSim ได้ทันที หรือว่าได้จำลอง rotor–rotor interaction ครบในตารางใบพัดเดี่ยว

กรณีศึกษาอุปกรณ์จริง: APC Thin Electric 10×5E

เลือกใบพัดรุ่นนี้เป็นจุดเริ่มต้นเพราะมีชุดวัดอุโมงค์ลมของ UIUC และเอกสารผู้ผลิตให้ตรวจข้ามแหล่ง ขนาดระบุ 10 นิ้ว × pitch 5 นิ้ว; เส้นผ่านศูนย์กลางที่ใช้แปลงหน่วยคือ D = 0.254 m การเลือกรุ่นนี้เป็นกรณีศึกษาไม่ได้หมายความว่าเหมาะกับมอเตอร์หรือโดรนทุกลำ

ตัวอย่างแปลงสัมประสิทธิ์: แถว UIUC ที่ 5,869 rpm ให้ CT = 0.0977 และ CP = 0.0373 เมื่อสมมติความหนาแน่นสำหรับการคำนวณ ρ = 1.225 kg/m³ และ n = 5,869/60 รอบ/s สูตรด้านบนให้ T ≈ 4.7664 N, P ≈ 45.2118 W, Q ≈ 0.073563 N·m นี่คือค่าที่คำนวณกลับภายใต้ ρ ที่กำหนด ไม่ใช่การอ้างว่าความหนาแน่นระหว่างทดลองเท่ากัน หรือเป็นกำลังไฟแบตเตอรี่

ลำดับงานตรวจเทียบคือเลือก static หรือ dynamic ให้ตรงโจทย์ → บันทึกไฟล์ รุ่น และเงื่อนไข → คำนวณ CT/CP จากแบบจำลองด้วยนิยามเดียวกัน → เปรียบเทียบที่รอบและ J ตรงกัน → วิเคราะห์ความต่างและความไม่แน่นอน การประมาณระหว่างจุดสำหรับการเรียนรู้ต้องอยู่ในช่วงข้อมูล และไม่แทนการตรวจเทียบกับข้อมูลอิสระ

ชุดจานแรง 0.08 N เดิมยังไม่ใช่แบบจำลองของใบพัดรุ่นนี้ ต้องสร้าง geometry/โมเดลใบพัดและเงื่อนไขที่ตรงก่อนใช้ผล UIUC ตรวจเทียบ ส่วน adapter ที่แทนแรงต้านลำตัวก็ไม่ใช่ทางนำเข้าตารางแรงขับใบพัด

ไฟล์ข้อมูลและภาพเชื่อมจากเจ้าของโดยตรง ไม่รวมสำเนาฐานข้อมูลใน ZIP อ่านวิธีอ้างอิงที่ UIUC Propeller Database และตรวจสิทธิ์เพิ่มเติมก่อนเผยแพร่ข้อมูลหรือ geometry ต่อ

ตรวจรูปทรงก่อนสร้าง CFD ของใบพัดรุ่นนี้

เปิดแผนภาพขนาดเต็ม

เราเปิด geometry archive ของ APC เดือนกุมภาพันธ์ 2026 และตรวจไฟล์ 10x5E-PERF.PE0 แล้ว พบหัวไฟล์ v2025-1001 วันที่ 24 กุมภาพันธ์ 2026 ให้ station/chord เป็นนิ้วและ twist เป็นองศา รวมถึงความหนา sweep และ rake ไฟล์ระบุ E63 และ APC12 พร้อมจุดเริ่ม/สิ้นสุด transition; APC อธิบายว่า APC12 เทียบเท่า NACA 4412 แต่ข้อความนี้เพียงอย่างเดียวยังไม่ระบุวิธี blend และผิวสามมิติทุกจุด

เทียบกับ ตาราง geometry UIUC ซึ่งให้ r/R, c/R และ β โดยข้อมูลสองชุดไม่ได้ยืนยันว่าเป็น specimen หรือ revision เดียวกัน และ APC แยกมุมตาม datum ของ leading/trailing edge ออกจาก pitch ที่วัดด้วยเครื่องมือ จึงไม่ควรเรียกส่วนต่าง β กับ twist ว่า error การวัดทันที

ตัวอย่างตรวจหน่วย: ใช้รัศมี nominal R = 0.127 m ที่ r/R = 0.75 จะได้ r = 0.09525 m; แถว UIUC ให้ c/R = 0.128 จึงได้ chord c = 0.016256 m หรือ 16.256 mm และ β = 13.39° นี่เป็นการแปลงข้อมูล ไม่ใช่การสร้าง airfoil จากความกว้างใบเพียงค่าเดียว

ดาวน์โหลด ชุดตรวจและเปรียบเทียบ geometry ด้วย Python แล้วดาวน์โหลดตาราง UIUC และ archive APC จากเจ้าของแยกต่างหากตาม README เครื่องมือจับคู่สถานีรัศมีสามจุดด้วย linear interpolation บันทึก hash และแสดงส่วนต่างแบบอธิบายข้อมูล ไม่ทำนายแรงขับและไม่ขยายข้อมูลนอกช่วง

ผลจากไฟล์ที่ตรวจจริง: ที่ r/R = 0.75 ตาราง APC หลัง interpolation ให้ chord ≈ 17.137 mm และ twist ≈ 11.983° ส่วน UIUC ให้ chord 16.256 mm และ β 13.39° ความกว้างต่างกันประมาณ 0.881 mm แต่ยังสรุปไม่ได้ว่าชุดใดผิด เพราะยังไม่ยืนยัน revision ของใบพัดและ datum ของมุมว่าตรงกัน เก็บผลนี้เป็นหลักฐานตรวจข้อมูลก่อนเลือกแบบจำลอง

ก่อนรันใบพัดจริงต้องตรวจเพิ่ม: พิกัดหน้าตัดและวิธี transition, นิยามแกน/ทิศหมุน, hub/root/tip, revision ที่ตรงกับใบพัดทดลอง และผิวปิดที่ mesh ได้ จากนั้นจึงเลือกรอบกับ J ให้ตรงข้อมูล UIUC และศึกษากริด/ขอบเขต ส่วนไฟล์ที่มีอยู่เพียงพอสำหรับฝึกตรวจข้อมูลรูปทรงแล้ว แต่ยังไม่ยืนยัน CAD ที่ตรง specimen เดิม

ชุดข้อมูลพร้อมศึกษาต่อ: MS1101 จาก ENOLA

ภาพสร้างจากข้อมูลรูปทรงของ UPV / ENOLA และเผยแพร่ภาพดัดแปลงภายใต้ CC BY-SA 4.0 เป็นภาพอ้างอิง geometry ไม่ใช่ภาพ CFD หรือภาพสินค้า เปิดภาพขนาดเต็ม

เพื่อให้เดินหน้าต่อได้โดยไม่ต้องเดาหน้าตัด APC เราเลือก T-Motor MS1101 ในชุด ENOLA เป็นกรณีศึกษาข้อมูลจริงเพิ่มเติม ดาวน์โหลดและเปิดตรวจ archive แล้ว พบทั้ง OBJ, ตารางใบพัด, ผลวัดแรงขับและแรงบิด และผล OpenFOAM ของผู้วิจัย ข้อมูลเหล่านี้เพียงพอสำหรับเริ่มตรวจ geometry และจัดกรณีเปรียบเทียบ แต่ยังไม่ใช่ผล CFD ที่เว็บนี้รันเองหรือผลยืนยันความแม่นยำของแบบจำลอง ข้อมูลต้นทางและดาวน์โหลด

เริ่มจากไฟล์ใด

ดาวน์โหลด Propeller_Database.zip จากเจ้าของ แล้วเปิดโฟลเดอร์ BBDD/T-Motor/MS1101/ บน Windows หรือ Linux:

ดาวน์โหลดคู่มือของเรา: บันทึกตรวจข้อมูลและ case sheet มีรายการไฟล์, hash, ตัวอย่างคำนวณ และขั้นตอนเตรียมกรณีศึกษา โดยไม่รวมฐานข้อมูลต้นทาง ผู้จัดทำอนุญาตข้อมูลภายใต้ CC BY-SA 4.0 และกำหนดการใช้เพื่อศึกษา/วิจัย ไม่ใช้ตัดสินใจการบินหรือภารกิจที่ความปลอดภัยขึ้นกับผลนี้

ตัวอย่างคำนวณจากผลวัดจริง

เลือกแถว 5088 RPM, V = 0 m/s ซึ่งบันทึกแรงขับ 4.6688 N และแรงบิด 0.05954 N·m:

n = 5088 ÷ 60 = 84.8 รอบ/s

P = 2πnQ = 2π × 84.8 × 0.05954 ≈ 31.724 W

นี่คือ กำลังกลที่เพลา ไม่ใช่กำลังไฟจากแบตเตอรี่ การหากำลังไฟต้องรู้ความสูญเสียของมอเตอร์และ ESC เพิ่มเติม ขณะ hover สูตร TV/P ให้ศูนย์ จึงไม่ควรนำไปสรุปว่าใบพัดไม่มีประสิทธิภาพ ใช้ตัวชี้วัด hover ที่เหมาะสมแยกต่างหาก

สิ่งที่ตรวจแล้ว และจุดที่ต้องรักษาไว้ในแบบจำลอง

• OBJ มี 335,812 vertices และ 671,620 สามเหลี่ยม ตรวจ indexed edges แล้วแต่ละเส้นติดกับสองหน้า จึงไม่พบขอบเปิดด้วยวิธีนี้ แต่ยังไม่รับรองว่าไม่มี self-intersection หรือสร้าง volume mesh ได้ดี
• Excel ระบุ D nominal = 0.2794 m ขณะที่ระยะ X ใน OBJ เท่ากับ 0.280288 ตรวจเพิ่มกับ ผู้ผลิต MS1101 ซึ่งระบุขนาดวัดจริง 280.3 mm จึงมีหลักฐานรองรับหน่วยเมตรของ OBJ โดยต่างเพียง 0.012 mm เก็บผิวเดิมและแยกขนาด nominal สำหรับ normalization ไม่ย่อผิวเพื่อให้ตรงตัวเลข 11 นิ้ว
• ตาราง twist ที่ r/R = 0.9 ให้ 8.75° แต่ช่อง radians ให้ 0.16 ซึ่งไม่ตรงกัน จึงใช้ OBJ เป็นรูปทรงหลัก ไม่สร้างผิวใหม่จากตารางโดยแก้ข้อมูลเงียบ ๆ
• CQ ที่เป็น NaN หมายถึงไม่มีข้อมูล ไม่ใช่ศูนย์ ส่วนสภาพอากาศและความไม่แน่นอนต้องตรวจตามกรณีทดลองก่อนกล่าวว่า validation ผ่าน

ลำดับทดลองที่ทำซ้ำได้

1. เก็บไฟล์ต้นฉบับและตรวจ hash จากชุดคู่มือ ดู geometry เพื่อกำหนดแกนหมุนและเครื่องหมายแรงบิด
2. เริ่ม hover ที่ 5088 RPM ตรงกับแถวทดลอง ใช้ Ω ≈ 532.814 rad/s บันทึกความหนาแน่นและความหนืดที่เลือก หากยังไม่ตรงการทดลองให้ระบุว่าเป็น preliminary comparison
3. เริ่ม MRF แล้วศึกษาขนาดโดเมนและกริดอย่างน้อยสามระดับ ตรวจแรงขับ แรงบิด y+ และ residual ร่วมกัน ก่อนขยายไป URANS
4. เปรียบเทียบ T และ Q แยกกัน ไม่ตั้งเกณฑ์ผ่านเป็นเปอร์เซ็นต์เองโดยไม่มีข้อมูลความไม่แน่นอน หากใช้ผล 3000/6000 RPM จากไฟล์ OpenFOAM ต้องระบุการ interpolation ของแถวทดลองด้วย
5. เมื่อมีผลที่ตรวจแล้วจึงออกแบบส่วนเชื่อมแรงใบพัดใน DroneSim แยกจาก airframe-drag adapter ที่มีอยู่

งาน Aular et al. (2026) สนับสนุนการตรวจข้อจำกัด MRF เมื่อการไหลเอียง แต่ใช้ใบพัดอ้างอิงคนละรูปทรง จึงนำค่าความคลาดเคลื่อนของงานนั้นมาเป็นเกณฑ์ผ่านของ MS1101 โดยตรงไม่ได้

เปรียบเทียบผลต้นทางก่อนใช้เป็นเกณฑ์

กราฟดัดแปลงจาก UPV / ENOLA v1 ภายใต้ CC BY-SA 4.0 เป็นผลของผู้วิจัยต้นทาง ไม่ใช่ผล CFD ที่เรารันเอง ไม่มี error bar เพราะยังไม่มีข้อมูลความไม่แน่นอนของชุดทดลองนี้ เปิดกราฟขนาดเต็ม

ที่ 6000 RPM ใช้ผลวัด 5815 และ 6068 RPM ทำ linear interpolation ได้ T ≈ 6.46244 N, Q ≈ 0.0814705 N·m ส่วนผล MRF ต้นทางให้ 5.5896 N และ 0.1134 N·m จึงต่างประมาณ −13.51% และ +39.19% ตามลำดับ ตัวเลขนี้เป็นส่วนต่างเชิงพรรณนา ไม่ใช่คะแนนผ่าน/ไม่ผ่าน และจุดที่ประมาณขึ้นมาไม่ใช่การวัดใหม่

สูตรตรวจซ้ำ: ส่วนต่าง (%) = 100 × (ค่าจำลอง − ค่าทดลองที่ประมาณ) ÷ ค่าทดลองที่ประมาณ หากตัวหารเป็นศูนย์ต้องใช้ส่วนต่างสัมบูรณ์แทน เครื่องมือจะไม่เทียบจุดที่อยู่นอกช่วง RPM ที่วัด

ดาวน์โหลดเครื่องมือเปรียบเทียบผล ENOLA ด้วย Python ใช้ได้บน Windows/Linux โดยไม่ติดตั้งไลบรารีเพิ่ม ดาวน์โหลด archive ต้นทางแยกตาม README เครื่องมือตรวจ hash ก่อนอ่านไฟล์และไม่เขียนทับรายงานเดิม

การตรวจ normalization ที่แถว 5088 RPM พบว่าเมื่อใช้นิยามมาตรฐานและ D = 0.2794 m การคำนวณกลับจาก T/CT ให้ค่าความหนาแน่น 1.20984 แต่จาก Q/CP ให้ 1.19584 kg/m³ ตัวเลขเหล่านี้เป็นเพียงเครื่องมือตรวจความสอดคล้อง ไม่ใช่ความหนาแน่นที่วัดได้ การเปลี่ยนเส้นผ่านศูนย์กลางเพียงค่าเดียวก็ไม่แก้ความต่างได้ทุกแถว จึงต้องขอนิยามการเฉลี่ยและ normalization จากต้นทางก่อนสรุปสาเหตุ

สิ่งที่ยังยืนยันแทนผู้วิจัยไม่ได้: สภาพอากาศจริงระหว่างวัด การสอบเทียบเซนเซอร์ และ uncertainty budget ไม่อยู่ใน README/CSV ที่ตรวจ จึงไม่ควรนำค่าอุณหภูมิหรือความหนาแน่นที่บทความใช้ตั้งค่า CFD มาอ้างว่าเป็นสภาพทดลอง ค่าที่เราตั้งในการรันต้องติดป้ายว่าเป็น simulation assumptions

ผลรันใบพัดจริง: ตรวจความพร้อมก่อนเชื่อค่าจำลอง

ผลล่าสุด: ปิดปัญหาการลู่เข้าได้หนึ่งกรณี

กรณี MS1101 กริด 156,104 เซลล์ผ่านเกณฑ์เชิงตัวเลขแล้ว หลังเปลี่ยนวิธีเชื่อมความดันกับความเร็วจาก SIMPLE เป็น SIMPLEC ตาม แนวทาง OpenCFD เราใช้ consistent=yes, pressure relaxation=1 และ U/k/omega relaxation=0.7 เป็นชุดการตั้งค่าที่ทดลองร่วมกัน รูปทรง ของไหล mesh และเกณฑ์ผ่านคงเดิม จึงยังไม่แยกว่าการตั้งค่าใดเพียงตัวเดียวเป็นสาเหตุทั้งหมด

solver ผ่านเกณฑ์ residual 10⁻⁵ ที่รอบ 3,554 แล้วเราปิดเฉพาะการหยุดอัตโนมัติเพื่อรันตรวจต่ออีก 200 รอบ โดยประเมินด้วยเกณฑ์เดิม ผลที่รอบ 3,754 คือแรงขับ 3.38252 N และโมเมนต์ต้าน 0.0848474 N·m residual สูงสุด 1.25 × 10⁻⁶ ช่วงแกว่งแรงขับ 100 รอบท้ายเพียง 0.000649% และ 200 รอบท้าย 0.00515% จึงผ่านทั้ง residual และความนิ่งของแรงในช่วงที่ตรวจ

นี่คือการปิดปัญหา iterative convergence ของกรณีนี้ ยังไม่ใช่การยืนยันว่าแรงตรงกับของจริง: ด้านล่างมีผลคัดกรองความไวต่อกริด/ขอบเขตและการกระจายชั้นใกล้ผิว ซึ่งยังไม่พอประกาศ grid independence ค่า yPlus ล่าสุดของการตรวจต่ออยู่ที่ checkpoint 3,750 ไม่ใช่แรงรอบ 3,754 และยังมีช่วงประมาณ 2.61–90.40 หลักฐานแยกไว้ใน simplec กับ simplec-hold ของ ZIP MS1101

อ่านค่า y⁺ เพื่อวางกริดใกล้ผิว

บริเวณติดผิวใบพัด ความเร็วเปลี่ยนเร็วในระยะสั้น กริดที่ดูละเอียดเมื่อเทียบกับขนาดใบพัดจึงอาจยังละเอียดไม่พอสำหรับชั้นนี้ ระยะไร้มิติ y⁺ = y₁uτ/ν ช่วยเชื่อมระยะจากผิวถึงศูนย์กลางเซลล์แรก y₁ กับความหนืดจลน์ ν โดยนิยามความเร็วเสียดทาน uτ = √(|τw|/ρ); τw คือความเค้นเฉือนที่ผิว หน่วย Pa ส่วน uτ มีหน่วย m/s

ตัวอย่างฝึกที่สมมติ uτ: ใช้ ν = 1.5 × 10⁻⁵ m²/s และ uτ = 1 m/s หากตั้งเป้า y⁺ = 1 จะได้ y₁ = 15 µm; หากตั้งเป้า 30 จะได้ y₁ = 0.45 mm ตัวเลขนี้แสดงผลของการเลือกเป้า ไม่ใช่ข้อกำหนดกริด MS1101 เพราะ uτ จริงเปลี่ยนตามตำแหน่งและต้องเลือก wall treatment ให้เข้ากันก่อน

หากใช้เซลล์ชั้นแรกที่มีศูนย์กลางใกล้กึ่งกลางความหนา อาจเริ่มประมาณ Δy₁ ≈ 2y₁ แล้วตรวจระยะจริงหลังสร้างกริด สำหรับ N ชั้นที่ขยายความหนาด้วยอัตรา g คงที่ ความหนารวมคือ H = Δy₁(gᴺ − 1)/(g − 1) เมื่อ g ≠ 1 เช่น Δy₁ = 30 µm, g = 1.2 และ N = 10 ได้ H ≈ 0.779 mm สูตรนี้ยังไม่ยืนยันว่าชั้นกริดครอบคลุม boundary layer หรือว่าสร้างได้ครบตามขอบใบพัด ดูขั้นตอน OpenCFD Layer addition

สำหรับกริดอ้างอิง 156,104 เซลล์ ที่ checkpoint 3,750 ตรวจค่ารายหน้าผิว 15,536 หน้าแล้วพบว่า 76.70% ของจำนวนหน้า อยู่ในช่วง 5 ≤ y⁺ < 30 แม้ค่าเฉลี่ยอยู่ที่ 22.99 จึงควรตรวจการกระจายประกอบเสมอ สัดส่วนนี้ นับจำนวนหน้า ไม่ได้ถ่วงพื้นที่ เพราะแต่ละหน้ามีพื้นที่ไม่เท่ากัน ช่วงที่แบ่งใช้เพื่ออ่านการกระจาย ไม่ใช่เกณฑ์รับรอง wall function

ค่าในไฟล์ใช้ ฟังก์ชัน yPlus ของ OpenFOAM 2512 ตามค่าเริ่มต้น useWallFunction=true ซึ่งอาศัยนิพจน์ของ nut wall function ที่เลือก จึงต้องระบุวิธีคำนวณเมื่อเทียบกับค่าจาก wall shear โดยตรง ชุดปัจจุบันยังไม่มี prism layers และยังไม่ปิดการศึกษาความไวใกล้ผิว

ฝึกตรวจซ้ำ: เปิด wall-yplus.csv และ wall-resolution.json ใน ZIP แต่ละกรณี รวมจำนวนหน้าในทุกช่วงให้เท่าจำนวนหน้าทั้งหมด แล้วตรวจ field_iteration ก่อนเทียบกับแรง ค่าผิวอาจอยู่ใน checkpoint ก่อนรอบแรงสุดท้ายตามรอบการบันทึกไฟล์

ความไวต่อกริดและขอบเขต: หลักฐานที่มีตอนนี้

เราเทียบค่าเฉลี่ย 100 รอบท้ายของการรัน SIMPLEC หลังตรวจว่าค่าฟิสิกส์และไฟล์ตั้งค่าที่ควรตรึงตรงกัน กรณี L4 และกริดอ้างอิง L5 ผ่านเกณฑ์ residual และแรง ส่วน L6 รันต่อจาก checkpoint 1,300 ถึง 1,800 รอบแล้ว ค่า residual รอบสุดท้ายผ่าน แต่มีเพียง 48 จาก 100 รอบท้ายที่ผ่านเกณฑ์ จึงยังไม่คำนวณร้อยละต่างเทียบกริดอ้างอิง

ผลนี้ชี้ว่าการเพิ่มความละเอียดเฉพาะผิวทำให้แรงที่คำนวณต่างกันระหว่าง L4 กับ L5 และการขยายโดเมนครั้งเดียวก็ยังเห็นผลต่าง จึงยังตัดผลจากกริดและขอบเขตทิ้งไม่ได้ แต่ยัง คำนวณ GCI หรือประกาศ grid independence ไม่ได้: การ refine เป็น local ไม่ได้ลดขนาดเซลล์อย่างเป็นระบบทั่วใบพัดและ wake, กรณี L6 ยังไม่ผ่าน residual window แม้รันเพิ่ม, ชุด relaxation ที่ทดลองแยกก็ไม่ผ่านเกณฑ์แรง, ยังไม่มี prism layers และใช้ first-order upwind ส่วนโดเมนที่ขยายหนึ่งระดับยังไม่พอทดสอบการเป็นอิสระจากขอบเขต. NASA แนะนำให้ใช้คำตอบที่ลู่เข้าแล้วบนกริดที่ละเอียดขึ้นเป็นลำดับและตรวจช่วง asymptotic ก่อนประมาณ Richardson/GCI (NASA Spatial Convergence).

ค่า y⁺ ของกรณี L4 มีช่วง 5.84–151.78 และ L5 reference 2.61–90.40; L6 checkpoint 1,800 อยู่ที่ 1.03–53.51 โดย 96.00% ของจำนวนหน้าผิวอยู่ในช่วง 5–30 ค่าในวงเล็บเหล่านี้นับจำนวนหน้า ไม่ได้ถ่วงพื้นที่ และไม่ใช่หลักฐานความแม่นของแรงเฉือน หน้าคู่มือ OpenFOAM อธิบายว่า y⁺ ต้องอ่านให้สอดคล้องกับ wall function ที่เลือก (OpenFOAM Wall Functions); ค่า L6 ที่ดูต่ำลงจึงไม่ทำให้ผล L6 ใช้เปรียบเทียบได้เมื่อ residual window ยังไม่ผ่าน.

ตรวจต่อพบ residual ของ Ux และ k ใน L6 แกว่งซ้ำใกล้เกณฑ์ แม้ค่าแรงและทอร์กท้ายชุดจะนิ่งมาก การทดลองวินิจฉัยแยกจาก checkpoint 1,800 โดยลด under-relaxation ของ U/k/omega จาก 0.7 เป็น 0.35 ทำให้แรงขับใน 100 รอบท้ายแกว่ง 14.19% และ residual สูงสุดในช่วงนั้นเป็น 2.28×10⁻³ (รอบสุดท้าย 1.49×10⁻³) จึงไม่รับผลชุดนี้และไม่ใช้เทียบกริด

อีกการทดลองหนึ่งลดเฉพาะ field relaxation ของความดัน p จาก 1.0 เป็น 0.3 และคงค่า U/k/omega ที่ 0.7 จาก checkpoint เดิม ผลช่วง 100 รอบ (1,901–2,000) ผ่าน residual เพียง 1/100 รอบ; Ux ผ่าน 1, Uy ผ่าน 4 และ Uz ผ่าน 9 รอบ ขณะที่ p/k/omega ผ่านครบ แรงขับและทอร์กนิ่งมาก (ช่วง peak-to-peak 0.00344% และ 0.00157%) แต่ไม่ผ่าน residual gate จึงยังใช้ไม่ได้ การทดสอบสองแบบชี้ว่าการลด relaxation อย่างเดียวไม่แก้ปัญหานี้ ต้องตรวจสมการความเร็ว/mesh/discretization ต่อ และใช้ protocol เดียวกันกับทุกกริดหากเปลี่ยนค่าตัวแก้สมการ ผลทั้งหมดเป็นการวินิจฉัยเชิงตัวเลข ไม่ใช่ผล validation. หลักการของ OpenFOAM SIMPLE และ under-relaxation อธิบายว่า relaxation ช่วยจำกัดการเปลี่ยนตัวแปร/สมการ แต่ไม่ได้เป็นเกณฑ์ยืนยัน convergence ด้วยตัวเอง.

ตรวจสมมติฐานว่า “เพิ่มจำนวนรอบแล้วจะผ่าน”

เราใช้สำเนาแยกจาก checkpoint L6 ที่รอบ 1,800 แล้วต่อถึง 2,000 โดยเปลี่ยนเฉพาะ endTime; hash ของ mesh, geometry, physics, schemes, fvSolution และสนาม checkpoint ตรงกับต้นฉบับ รอบ 1,901–2,000 ผ่าน residual รวม 46/100 เทียบกับ 48/100 ก่อนต่อ โดย Ux และ k ผ่านตัวแปรละ 68/100 ส่วน Uy ผ่าน 95/100 และค่าสูงสุดเป็น 1.36 × 10⁻⁵

แม้ residual ยังไม่ผ่าน ค่าเฉลี่ยแรงขับและโมเมนต์ต้านระหว่างหน้าต่าง 1,701–1,800 กับ 1,901–2,000 เปลี่ยนเพียง 0.0000073% และ 0.0000040% ตามลำดับ จึงสรุปได้เพียงว่าแรงเกือบคงที่ แต่เกณฑ์ residual ที่ตั้งไว้ยังไม่ผ่าน เครื่องมือตรวจรูปแบบ residual พบการเกิดซ้ำของ Ux/k ที่ lag ประมาณ 18 SIMPLE iterations; นี่เป็นรูปแบบในรอบแก้สมการ ไม่ใช่คาบเวลาทางกายภาพ และยังไม่ระบุสาเหตุ

เราเขียน quality fields ของ checkMesh และ cell-centre coordinates จากสำเนา mesh แยก แล้วจัดอันดับเซลล์ตาม non-orthogonality, skewness, aspect ratio, volume และ volume ratio พบจุด non-orthogonality สูงสุดบริเวณรัศมีราว 35 mm และ skewness สูงสุดราว 140 mm ภายในขอบเขต rotor ตามแบบจำลอง (r ≤ 150 mm, −25 ≤ z ≤ 25 mm) ตำแหน่งเหล่านี้เป็น cell centres สำหรับชี้พื้นที่ตรวจต่อ ไม่ได้ยืนยันว่าจุดนั้นตัดผิวใบพัดหรือเป็นสาเหตุของ residual cycle; เกณฑ์ Mesh OK และสถิติ/ลำดับนี้เป็นข้อมูลวินิจฉัย ไม่ได้สร้าง pass/fail gate ใหม่

CSV/JSON รายเซลล์อยู่ใน ชุด MS1101 โดยไม่บรรจุ mesh 87 MB หรือ raw quality fields เพิ่มเข้า ZIP วิธีทำซ้ำใช้ OpenFOAM checkMesh -writeAllFields และ postProcess -func writeCellCentres ตาม เอกสาร checkMesh ของ OpenFOAM 2512 และ writeCellCentres. เชื่อมตำแหน่งเหล่านี้กับแผนที่ residual และ gradient แล้ว; ผล diagnostic นี้ไม่ใช้เทียบกริดหรือยืนยันผลทดลอง

ตำแหน่งที่สนาม Ux และ k เปลี่ยนระหว่าง checkpoint

เราเปรียบเทียบค่ารายเซลล์ของ Ux และ k ที่บันทึกไว้รอบ 1,800→1,950 และ 1,950→2,000 บน L6 จำนวน 399,881 เซลล์ โดย mesh ในสำเนาวินิจฉัยมี hash ของ constant/polyMesh ตรงกับกรณี solver จึงจับคู่ดัชนีเซลล์ได้ ผลต่อไปนี้เป็น การเปลี่ยนของคำตอบที่บันทึกไว้ ไม่ใช่ residual หรือ gradient:

การถ่วงน้ำหนักใช้ Σ(V·|Δfield|²) ในแต่ละช่วงรัศมี หารด้วยผลรวมทั้งโดเมน และ r = √(x²+y²); แถบรัศมีจึงรวมทุกค่า z จุดที่เปลี่ยนมากที่สุดทั้งสี่กรณีอยู่ที่ r ≈ 0.153–0.167 m, z ≈ 0.0093 m และมี non-orthogonality ราว 25.2° ส่วนซองเรขาคณิตของ rotor (r ≤ 0.15 m, |z| ≤ 0.025 m) มีส่วนแบ่งการเปลี่ยน Δ² เพียง 0.10–0.42% ตามสนามและช่วงที่เลือก ค่ามากจึงกระจุกอยู่นอกซองนั้นในแถบรัศมีรอบนอก

นี่ช่วยชี้บริเวณให้ตรวจต่อ แต่ ไม่บอกสาเหตุ ของ residual cycle: ความต่างของ checkpoint ไม่ใช่สมการ residual หรือเวลาทางกายภาพ และไม่ใช่ข้อมูล validation สำหรับผลเทียบกริด ดู JSON/CSV และสคริปต์ทำซ้ำได้จาก ชุด MS1101 รุ่นล่าสุด

ตรวจ gradient ของความเร็วและพลังงานจลน์ปั่นป่วน

คำนวณ grad(U) และ grad(k) เพิ่มจาก checkpoint 1,800, 1,950 และ 2,000 ด้วย postProcess ของ OpenFOAM 2512 โดยใช้ fvSchemes ชุดเดิม ตัวดำเนินการ grad ให้เวกเตอร์เมื่อรับสนามสเกลาร์ และให้เทนเซอร์เมื่อรับสนามเวกเตอร์ ตาม OpenFOAM 2512 API เราเทียบองค์ประกอบทั้งหมดของสนาม gradient ระหว่าง checkpoint ไม่ได้เทียบเฉพาะขนาด เพื่อไม่ทิ้งการเปลี่ยนทิศทาง:

ΔgU = || grad(U)b − grad(U)a ||F      [1/s]
Δgk = | grad(k)b − grad(k)a |         [m/s²]
ส่วนแบ่งช่วงรัศมี = Σ(V · Δg²)ช่วง / Σ(V · Δg²)ทั้งโดเมน

การเปลี่ยน gradient ส่วนใหญ่เกิดในแถบรัศมี 0.15–0.18 m เช่นเดียวกับความต่างของ Ux/k จุดเปลี่ยน gradient ที่สูงสุดซ้ำในสองช่วงอยู่ใกล้ cell 375687 และ 375631 ที่ r ≈ 0.15365 m, z ≈ 0.00929 m และ non-orthogonality ≈ 25.3° อย่างไรก็ดี เซลล์ที่อยู่ใน 1% สูงสุดทั่วทั้งกริดเมื่อจัดอันดับด้วย non-orthogonality หรือ skewness รวมกันมีส่วนต่อผลต่าง gradient กำลังสองต่ำกว่า 0.0018% ผลนี้ระบุเพียงตำแหน่งที่ gradient เปลี่ยนระหว่าง checkpoint; ไม่ได้พิสูจน์ว่า mesh quality เป็นหรือไม่เป็นสาเหตุของ residual cycle

ค่า RMS ถ่วงปริมาตรกับ p95 แบบนับเซลล์ใช้การถ่วงน้ำหนักต่างกัน จึงอาจเห็น RMS ของ grad(k) สูงกว่า p95 มากเมื่อการเปลี่ยนกระจุกในเซลล์ปริมาตรเล็ก ช่วง checkpoint ห่างกัน 50 หรือ 150 รอบแก้สมการ steady SIMPLE ซึ่งไม่ใช่เวลาที่อากาศไหลจริง ที่สำคัญ grad(U/k) ไม่ใช่ local equation-residual: residual ที่บันทึกเป็นค่ารวมของ linear solver เราจึงยังต้อง instrument สมการอย่างตรงนิยาม หากต้องการแผนที่ residual เชิงตำแหน่งและหาสาเหตุของรูปแบบ Ux/k

ทำซ้ำด้วย extract_checkpoint_gradients.sh, วิเคราะห์ด้วย analyze_checkpoint_gradients.py และสร้างภาพด้วย plot_checkpoint_gradient_changes.py; ZIP มี logs, SHA256 provenance, CSV/JSON, tests และโค้ด โดยไม่รวมสนาม gradient ขนาดใหญ่หรือ mesh

แผนที่ local equation residual และ probe แนว wake

ในสำเนาวินิจฉัยแยก เราเปิด solverInfo ของ OpenFOAM 2512 พร้อม writeResidualFields โดยคง fvSolution และ fvSchemes เดิม แล้วต่อจาก checkpoint 2,000 ถึง 2,100 เพื่อเก็บ local initial residual ของสมการ Ux, Uy, Uz, p, k, omega ทุก 10 รอบ รวม 60 แผนที่ บน mesh 399,881 เซลล์ แผนที่คือ residual รายแถวของระบบสมการเชิงเส้นก่อนการแก้แต่ละสมการ ส่วนประวัติ residual ปกติเป็น norm รวมที่ solver ทำ normalization ตามเมทริกซ์และวิธีแก้สมการ ตาม OpenFOAM solverInfo และ นิยาม residual ของ OpenFOAM

สำหรับแผนที่ที่บันทึก แถบ r = 0.15–0.18 m มีสัดส่วนผลรวมค่าสัมบูรณ์ของ residual row ประมาณ 60.9% ใน Ux และ 81.3% ใน omega; ส่วน k อยู่ประมาณ 75.2% แผนที่ Ux และ k ที่ปรับให้ผลรวม L1 เท่ากันมีรูปร่างเกือบซ้ำระหว่างรอบ 2010 กับ 2100 (ระยะ L1 ประมาณ 6.8 × 10⁻⁷ และ 1.2 × 10⁻⁷ ตามลำดับ) ขณะที่ global initial residual สองสมการแทบคงค่าเดิมที่ 9.30 × 10⁻⁶ และ 9.60 × 10⁻⁶ ตามลำดับ นี่ช่วยระบุรูปแบบเชิงพื้นที่และการเกิดซ้ำในช่วงที่วัด แต่ไม่ได้พิสูจน์ว่าแถบรัศมีหรือคุณภาพ mesh เป็นสาเหตุของ cycle เพราะเป็น unweighted algebraic rows ไม่ได้ถ่วงด้วยปริมาตร และไม่ใช่ residual ของ nonlinear SIMPLE outer iteration ที่แจกแจงสาเหตุได้

เราวาง 24 จุด probe จากช่วงรัศมีและแนว wake ที่สนใจ บวก upstream references และจุด hotspot จาก gradient แล้วสั่ง postProcess ด้วย function object probes ที่รอบ 2000, 2050 และ 2100 ตาม คู่มือ OpenFOAM probes ที่ 24/24 จุดด้าน downstream ในแนว −z มี Uz < 0 ทุก checkpoint ซึ่งสอดคล้องกับอากาศที่ถูกเหนี่ยวนำให้เคลื่อนสวนทางแรงบนใบพัดซึ่งเป็น +z ความเร็วแกนกลางที่ z = −0.16 m อยู่ราว −6.20 m/s ณ รอบ 2100 แต่เป็นค่าจากสนาม steady solver ไม่ใช่การวัดจริงหรือตัวแทนความเร็ว wake ที่ผ่านการตรวจสอบแล้ว จุด probe ใช้ cell interpolation จึงไม่ใช่เส้น streamline ต่อเนื่อง

แผนที่ residual, จุด probe, ผลสรุป, metadata และสคริปต์วิเคราะห์อยู่ใน workbook; ไม่รวม raw mesh หรือ field ขนาดใหญ่ ค่า iteration เป็นรอบแก้สมการ steady SIMPLE ไม่ใช่เวลาไหลจริง ทั้งสองการตรวจจึงช่วยชี้บริเวณและตรวจเครื่องหมายทิศการไหลในแบบจำลองเท่านั้น ไม่ได้ปิด residual cycle หรือยืนยันความแม่นยำทางกายภาพ

NASA แนะนำให้ติดตามทั้ง residual และปริมาณวิศวกรรม เพราะอาจลู่เข้าต่างอัตรา ขณะที่ OpenFOAM ระบุว่า residual เป็นค่าที่ขึ้นกับตัว solver และวิธี normalization แหล่งเหล่านี้ช่วยตีความประวัติการคำนวณ แต่เราไม่เปลี่ยนเกณฑ์คัดกรองของชุดฝึกจากหลักฐานนี้ (NASA: iterative convergence, OpenFOAM: residuals).

ดาวน์โหลด ชุด MS1101 รุ่นล่าสุด เพื่อดู sensitivity-simplec-study.csv/json, checkpoint field และ gradient diagnostics, log/provenance, wall-resolution.json และ wall-yplus.csv รวมถึงกรณีที่ไม่ผ่านเกณฑ์ โค้ดวิเคราะห์แสดงค่า candidate lag เพื่อช่วยอ่านรูปแบบเท่านั้น ไม่สรุปสาเหตุ และตัวตรวจผลต่างยังปฏิเสธกรณีที่ residual window ไม่ครบหรือค่าตั้งทางฟิสิกส์ไม่ตรงกัน

ตรวจกรอบหมุนกับคำตอบสมการก่อนใช้กับใบพัด

เพื่อแยกปัญหาการตั้ง MRF ออกจากความซับซ้อนของใบพัด เราสร้างโจทย์ของไหลระหว่างทรงกระบอกสองวงร่วมแกน: วงในรัศมี 1 m หมุน 1 rad/s วงนอกรัศมี 2 m อยู่นิ่ง กำหนด ν = 0.1 m²/s และสนามสองมิติ ไม่มีการไหลตามแกนหรือปลายทรงกระบอก แนวโจทย์อ้างอิง OpenCFD Rotating cylinders โดยไฟล์และผลรันในชุดนี้เราสร้างเอง

เมื่อเป็น steady, laminar และสมมาตรรอบแกน สมการความเร็วตามแนวสัมผัสวงกลมลดเหลือ:

d²uθ/dr² + (1/r) duθ/dr − uθ/r² = 0
uθ = A r + B/r
uθ(1 m) = 1 m/s, uθ(2 m) = 0
A = −1/3 s⁻¹, B = 4/3 m²/s

ตัวอย่างที่ r = 1.5 m ได้ uθ = 7/18 ≈ 0.388889 m/s เราคำนวณจากสนาม CFD ด้วย uθ = (−yUx+xUy)/r เพื่อรักษาเครื่องหมายทิศหมุน แล้วเทียบกับคำตอบนี้ที่ตำแหน่งเดิม 19 จุดภายในโดเมน

เมื่อเพิ่มความละเอียดตามแนวรัศมีและรอบวงสองเท่า error ลดลงเกือบสี่เท่า คำนวณอันดับที่สังเกตได้ด้วย p = log(eหยาบ/eละเอียด)/log(2) ได้ 1.966 และ 1.979 ค่านี้รวมผลการดึงค่าที่จุดตัวอย่าง จึงไม่ใช่ error ของทุกเซลล์ในโดเมน กรณีละเอียดมี error สูงสุดประมาณ 0.0995% ของความเร็ววงใน 1 m/s

ดาวน์โหลดชุดตรวจ MRF สามกริด พร้อมโค้ดและผลจริง เปิดผลและทดสอบ Python ได้บน Windows/Linux; รัน OpenFOAM ซ้ำบน Ubuntu หรือ Windows ผ่าน WSL ตาม README การผ่านชุดนี้ตรวจได้เฉพาะปัญหากรอบหมุนแบบ laminar ที่กำหนด ไม่รับรองโมเดล turbulence หรือสมรรถนะใบพัด

แบบฝึก: ทำไม error ลดสี่เท่าจึงได้อันดับประมาณสอง? เพราะ log(4)/log(2)=2 ลองคำนวณด้วยค่าจริงในตาราง แล้วอธิบายว่าเหตุใดไม่เท่ากับสองพอดี

หลักฐานการทดลองก่อนพบวิธีที่ลู่เข้า

กราฟนี้เป็น ผลที่เรารันเอง ต่างจากกราฟเปรียบเทียบผลต้นทางด้านบน ใช้ผิว MS1101 จาก ENOLA และคำนวณแรงจากความดันรวมกับความหนืด ไม่ได้ป้อนแรงขับทดลองเข้าไปบังคับผล เปิดกราฟขนาดเต็ม

สิ่งที่ชุดทดลองทำได้แล้ว

สร้าง mesh รอบใบพัดจริงและรัน OpenFOAM 2512 แบบ steady MRF / k-omega SST ที่ 5088 RPM ได้ พร้อมบันทึกแรงขับ แรงบิด residual และ yPlus เลือก ρ = 1.225 kg/m³ และ ν = 1.5 × 10⁻⁵ m²/s เป็น สมมติฐานของการจำลอง ซึ่งยังไม่ยืนยันว่าตรงกับสภาพทดลอง ENOLA

ผิวที่ใช้คงพิกัดเดิม ตัดออกเฉพาะสามเหลี่ยมพื้นที่ศูนย์ 2 หน้า จากทั้งหมด 671,620 หน้า แล้วตรวจผิวปิดและทิศ normal ก่อนสร้าง mesh ชุดทดลองตรวจ hash ของ archive และ OBJ ก่อนอ่าน ไม่รวมสำเนารูปทรงต้นทางหรือ OpenFOAM ไว้ในไฟล์ดาวน์โหลด

ดาวน์โหลดชุดรัน MS1101 พร้อมหลักฐาน CFD สำหรับ Ubuntu 24.04 amd64 หรือ Windows ผ่าน WSL Ubuntu 24.04 มี README ขั้นตอนรันและผลตรวจแต่ละกรณี ต้องดาวน์โหลด ENOLA จากเจ้าของแยกต่างหาก และใช้ชื่อโฟลเดอร์ผลลัพธ์ใหม่ทุกครั้ง

อ่านผลอย่างไรโดยไม่สรุปเกินหลักฐาน

รันจบตามจำนวนรอบไม่เท่ากับลู่เข้า เราตรวจ 100 รอบท้าย: ค่าเฉลี่ยระหว่างครึ่งแรกกับครึ่งหลังเปลี่ยนต่ำกว่า 1% และช่วงสูงสุด–ต่ำสุดต่ำกว่า 2% ของค่าเฉลี่ย ทั้งแรงขับและแรงบิด พร้อม initial residual ของตัวแปรทั้งหกต่ำกว่า 10⁻⁵ เกณฑ์นี้เป็นเกณฑ์คัดกรองเชิงตัวเลขของแบบฝึก ไม่ใช่ใบรับรองความถูกต้องทางกายภาพ

กริด 85,289 cells ที่รอบ 800 ให้แรงขับรอบสุดท้ายประมาณ 2.5165 N แต่แรงขับ 100 รอบท้ายยังแกว่งเป็นช่วงประมาณ 14.12% ของค่าเฉลี่ย ส่วนกริด 156,104 cells ให้ประมาณ 2.7336 N และช่วงแกว่ง 8.83% จึงยังไม่ผ่านการลู่เข้า ไม่ควรเลือกตัวเลขรอบสุดท้ายไปใส่ DroneSim เป็นสมรรถนะที่ยืนยันแล้ว

เมื่อลด relaxation และรันกริด 156,104 cells ต่อถึงรอบ 1600 แรงขับ 100 รอบท้ายแกว่งเหลือประมาณ 1.24% และผ่านเกณฑ์ plateau ของแรงทั้งสองตัว แต่ residual สูงสุดยังประมาณ 8.42 × 10⁻⁴ จึง ยังไม่ผ่านเกณฑ์รวม ตัวอย่างนี้แสดงว่ากราฟแรงดูนิ่งขึ้นเพียงอย่างเดียวยังไม่เพียงพอ

การศึกษากริดที่สาม โดเมนใหญ่ขึ้น และการลด relaxation มีบันทึกใน results/*/summary.json และกราฟด้านบน:

กรณี finer มีค่า yPlus จาก checkpoint รอบ 350 จึงไม่อ้างว่าเป็นค่ารอบ 358 และยังไม่สรุปสาเหตุการหยุดจากหลักฐานที่มี การปรับ relaxation เปลี่ยนวิธีแก้สมการ ไม่ได้ปรับความหนาแน่นหรือแรงให้เข้ากับข้อมูลทดลอง

อัปเดต 1.12.6: รันกริด finer ต่อจาก checkpoint 350 ในสำเนาแยกจนถึง 800 สำเร็จแล้ว โดยไม่เปลี่ยน mesh หรือเงื่อนไข แรงขับรอบสุดท้าย 2.8932 N, แรงบิด 0.08292 N·m แต่แรงขับ 100 รอบท้ายยังมีช่วงแกว่ง 12.33% และ residual สูงสุด 1.37 × 10⁻³ จึงยังไม่ลู่เข้า ค่า yPlus รอบ 800 อยู่ประมาณ 1.10–64.48 (เฉลี่ย 16.76) การรันครบช่วยปิดช่องว่างการประมวลผล แต่ไม่ทำให้การตรวจเทียบผ่านโดยอัตโนมัติ

ผลรันต่อที่ตรวจเพิ่ม: เราต่อกรณี relaxed จากรอบ 1,600 ถึง 3,200 โดยตรวจว่า geometry และค่าตั้ง 24 ไฟล์คงเดิม เปลี่ยนเพียง endTime โปรแกรมจบปกติ ได้ T ≈ 3.02979 N, Q ≈ 0.0948042 N·m แต่แรงขับ 100 รอบท้ายมีช่วงแกว่ง 5.13% และ residual สูงสุด 7.76 × 10⁻⁴ จึงยังไม่ผ่านทั้งเกณฑ์แรงและ residual ช่วงที่เคยดูนิ่งในรอบ 1,600 ไม่คงอยู่เมื่อรันนานขึ้น ผลชุด relaxed-3200 และ provenance รวมไว้ใน ZIP แล้ว

ข้อจำกัดสำคัญ: mesh ยังไม่มี prism layers และ yPlus กระจายหลายช่วง; ใช้ convection แบบ first-order upwind; การแกว่งตาม SIMPLE iteration ไม่ใช่การแกว่งตามเวลาจริงของใบพัด แม้มีหลายกริดก็ยังอ้าง grid independence, Richardson extrapolation หรือ GCI ไม่ได้เมื่อผลยังไม่ลู่เข้า

ทางไปสู่แบบจำลองที่ตรวจเทียบได้

1. ตรวจชั้นใกล้ผิวและความไวต่อ mesh/ขอบเขต พร้อมความสอดคล้องของแกนหมุนและแรง หาก steady MRF ยังไม่ลู่เข้า ให้ศึกษา transient พร้อม timestep และช่วงเฉลี่ยที่เพียงพอ
2. ขอข้อมูลสภาพทดลอง นิยาม normalization การติดตั้งและความไม่แน่นอนจากเจ้าของข้อมูล แล้วเก็บแยกจากค่าที่ตั้งสมมติใน CFD
3. เมื่อผ่านการตรวจเชิงตัวเลขแล้วจึงเทียบ T/Q และสัมประสิทธิ์ ณ รอบและสภาวะเดียวกัน พร้อม uncertainty ก่อนส่งค่าผ่าน adapter เข้า DroneSim

ชุดนี้จึง พร้อมใช้ฝึกรันและวิเคราะห์สาเหตุที่ผลยังไม่น่าเชื่อถือ ส่วนการยืนยันสมรรถนะใบพัดยังไม่ผ่าน ไม่ใช้ผลนี้ตัดสินใจการบินจริง

ข้อมูลเพิ่มจากการตรวจงานวิจัยและวิธีรายงานผล

เราเปิดตรวจ งาน ICAS 2024 ของทีม ENOLA เพิ่มแล้ว พบว่าใช้ใบพัดคนละรุ่นกับ MS1101 จึงไม่ยกข้อมูลการสอบเทียบหรือวิธีเก็บตัวอย่างของงานนั้นมาเป็น metadata ของชุดนี้ ส่วน บทความ AST ที่เกี่ยวข้อง ระบุว่าขอข้อมูลจากผู้วิจัยได้ การค้นรอบนี้ยังไม่ยืนยันสภาพอากาศและ uncertainty ของแถว MS1101 ที่ใช้ตรวจเทียบ

NASA แนะนำให้ดูทั้ง residual และปริมาณวิศวกรรมที่สนใจ ค่าแรงที่นิ่งและสมการที่ลู่เข้าอาจเกิดคนละช่วงเวลา เกณฑ์ 10⁻⁵ ในชุดฝึกเป็นเกณฑ์ที่เราประกาศใช้ ไม่ใช่ค่ามาตรฐานสากลที่ NASA รับรองให้ทุกกรณี

ต่อมาพบ ฉบับเต็ม AST ที่ผู้เขียนเผยแพร่ ซึ่งมีตารางความไม่แน่นอนของเครื่องชั่งแรง แต่ยังเป็นการทดลองใบพัดอีกชุด จึงเพิ่มคำถามเรื่องเครื่องชั่ง การสอบเทียบ แกน และ coverage factor ในคำขอข้อมูล แทนการนำค่าจากตารางนั้นมาใช้กับ MS1101 ทันที

ค่า yPlus ในชุดนี้ใช้ค่าเริ่มต้น useWallFunction = true ตาม OpenFOAM 2512 จึงต้องอ่านร่วมกับ wall function ที่เลือก ไม่ถือว่าเป็นการตรวจ shear จากสนามไหลโดยตรงเสมอไป ส่วน kOmegaSSTLM รองรับการเปลี่ยนจาก laminar เป็น turbulent แต่ต้องเตรียมตัวแปรและเงื่อนไขเพิ่ม สูตร turbulence intensity ที่หารด้วยความเร็วขาเข้าใช้ตรง ๆ กับความเร็วศูนย์ของ hover ไม่ได้ จึงยังไม่สลับโมเดลแล้วอ้างว่าปัญหาหาย

เตรียมความไม่แน่นอนก่อนเทียบ CT และ CP

ภาพใช้ค่าจาก uncertainty-example.json เพื่อแสดงว่าแต่ละอินพุตมีส่วนต่อความแปรปรวนมากเพียงใด เมื่อเปลี่ยนค่าหรือ uncertainty สัดส่วนก็เปลี่ยน ไม่มีผลการวัด ENOLA อยู่ในตัวอย่างนี้ เปิดกราฟขนาดเต็ม

เมื่อ CT = T/(ρn²D⁴) และ CP = 2πQ/(ρn²D⁵) โดย n = RPM/60 หากอินพุตเป็นอิสระและความไม่แน่นอนเล็กพอสำหรับการประมาณอันดับหนึ่ง จะได้:

u(CT)/CT = √[(u(T)/T)² + (u(ρ)/ρ)² + (2u(n)/n)² + (4u(D)/D)²]

u(CP)/CP = √[(u(Q)/Q)² + (u(ρ)/ρ)² + (2u(n)/n)² + (5u(D)/D)²]

u คือ standard uncertainty ไม่ใช่ค่าความคลาดเคลื่อนสูงสุดจากสเปกโดยอัตโนมัติ หากตัวแปรสัมพันธ์กันต้องรวม covariance ตาม แนวทาง NIST TN 1297 จุดที่พลาดได้ง่ายคือคำนวณ P = 2πnQ แล้วนำ P กับ n ไปถือว่าเป็นอิสระอีกครั้ง จึงควรลดสูตรให้เหลือ Q และ n ก่อนอย่างที่แสดงด้านบน

ตัวอย่างสมมติเพื่อเรียนรู้: หากมีเพียง u(D)/D = 1% และอินพุตอื่นไม่มีความไม่แน่นอน สูตรให้ u(CT)/CT = 4% และ u(CP)/CP = 5% นี่แสดงความไวต่อเส้นผ่านศูนย์กลาง ไม่ใช่ uncertainty ที่วัดได้ของ ENOLA

ชุด Python เปรียบเทียบ ENOLA เพิ่ม uncertainty.py และ uncertainty-example.json แล้ว รัน python uncertainty.py uncertainty-example.json new-result.json โดยตัวอย่างทั้งหมดติดป้ายสมมติ เครื่องมือไม่เติมค่าที่ขาดเป็นศูนย์ ไม่ประเมิน model bias และไม่ประกาศผล validation อัตโนมัติ

เทียบกับแถวทดลองจริงที่ 5,088 RPM

เราเปิดตรวจไฟล์ทดลอง 0° ของ MS1101 ในฐานข้อมูล UPV ENOLA v1 ซึ่งบันทึกไว้ตรงจุด hover 5,088 RPM, V = 0 m/s จึงไม่ต้องประมาณค่าระหว่างรอบ ผลวัดแถวนี้คือ T = 4.6688 N, Q = 0.05954 N·m และกำลังเพลา P = ΩQ = 31.72 W; ค่า CQ ในต้นฉบับเป็น NaN จึงคงสถานะ missing ไม่แทนด้วยศูนย์ (ข้อมูลต้นทางและดาวน์โหลด, เงื่อนไขสิทธิ์ใช้ข้อมูล).

จุด CFD ที่นำมาเทียบผ่านเกณฑ์ตรวจเชิงตัวเลขของชุดงานแล้ว แต่สมมติ ρ = 1.225 kg/m³ และ ν = 1.5 × 10⁻⁵ m²/s เพราะไฟล์ทดลองที่ตรวจไม่ระบุสภาวะอากาศรายจุดให้จับคู่ เรารายงานผลต่างแบบมีเครื่องหมายด้วย 100 × (CFD − experiment) / experiment และคำนวณกำลังจาก P = ΩQ, Ω = 2π(RPM/60):

ความต่าง L4/L5 และโดเมนขยายเป็น ผล sensitivity คนละคำถามกับการเทียบข้อมูลทดลอง: คู่ L4→L5 เปลี่ยนแรงขับ 5.82% และแรงบิด 1.51%; การขยายโดเมน L5 เปลี่ยนแรงขับ 1.65% และแรงบิด 0.125% เมื่อเทียบ L5 อ้างอิง ชุด L6 ที่มีอยู่เดิมไม่อยู่ในตาราง เพราะ residual ผ่านเพียง 48/100 รอบท้าย แม้แรงจะนิ่ง จึงไม่ผ่านเกณฑ์รับผลของโครงการ

การศึกษา co-refinement เพิ่มระดับผิวใบพัด/บริเวณ rotor/wake พร้อมกันจาก 4/2/1 เป็น 5/3/2 และ 6/4/3 โดยคง domain, ฟิสิกส์และ scheme เดิม ทั้งสามกริดผ่าน checkMesh, force plateau, solver-end และ residual window 100/100; fine case รัน OpenFOAM 2512 แบบ Open MPI 4.1.6 จำนวน 16 ranks ใช้หน่วยความจำราว 3.7 GB บน WSL จุดสำคัญคือคำตอบ ไม่เปลี่ยนไปทางเดียวกัน: จากหยาบไปกลาง T เปลี่ยน −2.39%, Q −5.54% แต่จากกลางไปละเอียด T กลับ +42.55%, Q +11.24% แม้ fine grid แรงทรงตัวในหน้าต่างท้าย 100 รอบอย่างแคบมาก ผลนี้เป็นหลักฐานของ grid sensitivity ที่ไม่ monotonic ไม่ใช่ grid convergence

ดังนั้นเครื่องมือจงใจไม่คำนวณหรือรายงาน Richardson extrapolation/GCI สำหรับทั้งแรงขับและแรงบิด การมีสามกริดที่ solver ผ่านเพียงอย่างเดียวยังไม่พอ หากคำตอบไม่เป็นลำดับและยังระบุช่วง asymptotic ไม่ได้ การใช้สูตร GCI กับข้อมูลนี้จะให้ตัวเลขที่ชวนเข้าใจผิด NASA แนะนำอย่างน้อยสามระดับเพื่อประมาณอันดับและตรวจ asymptotic range ก่อนใช้ extrapolation/GCI (NASA NPARC: spatial convergence); ขั้นตอนมาตรฐานสำหรับรายงานความไม่แน่นอนเชิง discretization อธิบายโดย Celik และคณะ (ASME J. Fluids Engineering, 2008).

ข้อสรุปที่ใช้ได้ตอนนี้: solver fine case ผ่าน numerical gates และแรงขับเข้าใกล้จุดทดลอง แต่แรงบิดยังสูงกว่าการวัด 52.01% และผล co-refinement ไม่ monotonic จึงยังสรุป grid uncertainty ไม่ได้ และทำได้เพียง quantitative preliminary comparison — ยังไม่ใช่ formal physical validation และห้ามนำตัวเลขไปคาดการณ์สมรรถนะหรือกำหนดภารกิจ การผ่าน iterative convergence เป็นคนละข้อกับ spatial convergence; ต้องตรวจลำดับกริดและช่วง asymptotic แยกกัน (NASA: iterative convergence, NASA: spatial/grid convergence).

ไฟล์ทดลอง 0° จำนวน 21 แถวและผลคำนวณที่ตรวจย้อน hash ได้อยู่ใน ชุดข้อมูล MS1101 เปรียบเทียบผลวัด; มี CSV, JSON, สคริปต์และคำอธิบายแหล่งที่มา โดยอ้างอิง ENOLA CC BY-SA 4.0 ชุดจำลองทำซ้ำแยกได้จาก MS1101 CFD workbook ใช้ Python และ OpenFOAM 2512 ตาม README

ข้อมูล ENOLA อธิบายว่ามีการทดสอบที่อุโมงค์ลม Francisco Payri และแท่นวัดแรง/โมเมนต์ custom balance พร้อมอ่าน load cells ผ่านระบบ National Instruments DAQ; หน้าโครงการระบุขนาดอุโมงค์ 2.8 × 2.8 × 22 m และความเร็วลมสูงสุด 40 m/s (แท่นทดลองของ ENOLA, สิ่งอำนวยความสะดวกโครงการ). อย่างไรก็ตาม แหล่งที่ตรวจยังไม่จับคู่ข้อมูลนี้กับแถว MS1101 ที่ 5,088 RPM โดยเฉพาะ และไม่ได้ให้ uncertainty/calibration รายจุด สภาพอากาศในเวลาวัด หรือหลักฐานยืนยันว่า OBJ เป็น specimen เดียวกับชิ้นงานทดลอง จึงห้ามนำรายละเอียดแท่นทดลองทั่วไปมาแทน metadata รายจุด เราไม่สร้าง error bar หรือยืม uncertainty จากใบพัดคนละรุ่น

ตรวจผิวใบพัด: ความต่างของแรงมาจากส่วนไหน

เมื่อผลสามกริดไม่เรียงเข้าหาค่าเดียวกัน เราควรย้อนดูสิ่งที่กริดแทนอยู่ การตรวจรอบนี้อ่านพิกัดหน้าผิวและสนามที่บันทึกไว้ รอบ 4,000 เหมือนกันทั้งสามกริด แล้วรวมแรงจากความดันอีกครั้ง ผลตรงกับแรงดันใน log ของ OpenFOAM ภายในเกณฑ์ตรวจสัมพัทธ์ 10⁻⁸ จึงยืนยันได้ว่าเครื่องมืออ่านข้อมูลและทิศแรงสอดคล้องกัน แต่ยังไม่บอกว่าแบบจำลองตรงกับการทดลอง

อ่านกราฟ: แท่งด้านบนแสดงพื้นที่จริงที่แต่ละช่วง y⁺ ครอบคลุม จึงต่างจากสถิติที่นับจำนวนหน้าในกรณีศึกษาเก่า เส้นด้านล่างรวมแรงขับจากความดันในวงแหวนกว้าง 0.1R โดย R = 0.140144 m ตามรูปทรงที่นำเข้า ความดันรวมต่างระหว่างกริดละเอียดกับกริดกลาง 1.32606 N; ในครึ่งนอกของใบพัด (r/R ≥ 0.5) ต่าง 1.02623 N หรือประมาณ 77.4% ของผลต่างนี้ แรงขับรวมที่เพิ่มจึงมาจากความดันเป็นหลัก มิใช่การเพิ่มแรงหนืดโดยตรง อย่างไรก็ตาม wall treatment สามารถเปลี่ยนสนามความเร็วและความดันผ่านการแก้สมการร่วมกันได้ จึงยังตัดบทบาทของมันไม่ได้

สูตรที่ใช้ตรวจซ้ำมีดังนี้:

• ค่าเฉลี่ยถ่วงพื้นที่ ȳ⁺ = Σ(Aᵢyᵢ⁺) / ΣAᵢ โดย Aᵢ คือขนาดเวกเตอร์พื้นที่หน้าที่ i หน่วย m²
• แรงจากความดัน Fₚ = ρ Σ(pₖᵢSᵢ) โดย pₖ คือความดันหารด้วยความหนาแน่น หน่วย m²/s² ซึ่งเป็นค่าที่เก็บในไฟล์ p ของกรณี incompressible นี้; Sᵢ ชี้ออกจากของไหลเข้าสู่ตัวใบพัด และ ρ = 1.225 kg/m³
• โมเมนต์จากความดัน Mₚ = Σ[(cᵢ − c₀) × (ρpₖᵢSᵢ)] หน่วย N·m; cᵢ คือจุดศูนย์กลางหน้า และ c₀ = (0,0,0) คือศูนย์อ้างอิงเดิม

ตัวอย่างฝึก: ผิวสองหน้ามีพื้นที่ 1 และ 3 mm² และ y⁺ เท่ากับ 10 และ 30 ตามลำดับ การนับหน้าจะเฉลี่ยได้ 20 แต่การถ่วงพื้นที่ได้ (1×10 + 3×30)/(1+3) = 25 ดังนั้นต้องระบุวิธีเฉลี่ยก่อนเทียบตัวเลข

ปริมาตรที่ใกล้ OBJ มากขึ้นช่วยตรวจรูปทรงโดยรวม แต่ไม่ได้ยืนยันความละเอียดที่ขอบนำ ขอบตาม หรือชั้นใกล้ผิว กริดชุดนี้ยังมีเซลล์ระดับ 0 ใน far field ประมาณ 29,000–31,000 เซลล์ แม้บริเวณใบพัดและ wake จะละเอียดขึ้น จึงเป็น local co-refinement ไม่ใช่การย่อระยะกริดทุกส่วนด้วยอัตราส่วนเดียว NASA อธิบายว่าการศึกษากริดต้องพิจารณาระยะใกล้ผิว ขอบเขต และบริเวณต่อเชื่อมให้เป็นระบบด้วย ไม่ใช่อาศัยจำนวนเซลล์รวมอย่างเดียว (NASA: grid considerations)

ช่วง 5–30 ในตารางเป็นช่วงแบ่งเพื่อวิเคราะห์ ไม่ใช่เกณฑ์รับรองหรือเกณฑ์ล้มเหลวอัตโนมัติ เราใช้ค่า yPlus ที่คำนวณผ่าน wall function ตามที่ระบุในหัวข้อก่อนหน้า หลักฐานนี้ชี้ให้ตรวจ wall treatment ร่วมกับชั้นกริดใกล้ผิวและการแทนรูปทรง ก่อนสร้างกริดชุดใหม่ คู่มือ OpenFOAM อธิบายการเพิ่มเซลล์เป็นชั้นตามผิวใน addLayersControls และการเลือกความหนาสัมบูรณ์/สัมพัทธ์ ซึ่งต้องตรวจจำนวนชั้นที่สร้างได้จริงและคุณภาพหลังสร้างด้วย (OpenFOAM: layer addition)

เปิด surface-faces.csv เพื่อดูรายหน้า, surface-audit.json เพื่อดู hash ของ mesh/สนาม และ surface-family-summary.csv เพื่อดูตารางสรุปใน MS1101 CFD workbook การอ่านสนามละเอียดที่รันแบบขนานใช้ reconstructPar -time 4000 -fields '(p yPlus nut)' ก่อนวิเคราะห์ ไม่ได้แก้ผล solver ต้นฉบับ

ทดสอบ wall function โดยมีกรณีควบคุม

เพื่อทดสอบสมมติฐานเรื่องชั้นใกล้ผิว เราแยกสำเนากริดกลางจาก checkpoint 4,000 เป็นสองกรณี และกำหนดรันแบบ serial ด้วย OpenFOAM 2512 ต่อถึง 4,500 เหมือนกัน กรณีควบคุมใช้ค่าทั้งหมดเดิม ส่วนกรณีทดสอบเปลี่ยนเฉพาะ nutkWallFunction เป็น nutUSpaldingWallFunction ในสนาม nut; กริด ความเร็วรอบ ของไหล omega wall function และวิธีแก้สมการคงเดิม มี hash ของไฟล์ก่อน/หลังเพื่อยืนยันขอบเขตการเปลี่ยน

Spalding wall function ใช้ความเร็วสร้างความสัมพันธ์ของความหนืดปั่นป่วนที่ต่อเนื่องถึงผิว แทนการประเมินผ่าน k แบบเดิม การเลือกนี้เป็น การทดลองความไวหนึ่งตัวแปร ไม่ใช่คำยืนยันว่าเหมาะสมกว่ากับใบพัดนี้ (OpenFOAM: nutUSpaldingWallFunction) ผลต้องผ่าน residual ทุกสมการใน 100 รอบท้ายและความนิ่งของแรงด้วยเกณฑ์เดิมก่อนรายงานเป็นผลเปรียบเทียบที่ลู่เข้า และค่า yPlus หลังเปลี่ยน wall function ต้องระบุวิธีคำนวณ เพราะนิพจน์ที่ใช้ประเมินเปลี่ยนตามแบบจำลอง

บทเรียนจากผลจริง: การเปลี่ยน wall function ทำให้ solver ต้องปรับสนามใหม่ แม้แรงดูนิ่ง แต่กรณีทดสอบยังไม่ผ่าน residual ในกรอบ 500 รอบที่กำหนด เราจึงไม่รายงานเปอร์เซ็นต์ความต่างว่าเป็นผลของ wall model ที่ลู่เข้า ไม่สรุปว่า Spalding แก้ non-monotonic response และไม่เลือกใช้เพราะแรงเข้าใกล้ข้อมูลทดลอง

เข้าใจ tolerance สามชั้นก่อนปรับ solver

คำว่า “คำนวณลู่เข้าแล้ว” ต้องระบุว่าหมายถึงระดับไหน ในตัวอย่างนี้มีอย่างน้อยสามระดับที่ไม่ควรนำค่ามาเทียบกันตรง ๆ:

1. ภายใน wall function: หาความเร็วเสียดทาน uτ ที่ทำให้ความสัมพันธ์ใกล้ผิวสอดคล้องกัน โดยวนวิธี Newton–Raphson
2. ภายใน linear solver: แก้ระบบสมการที่ประกอบขึ้นในรอบนั้น เช่น Aφ = b; การแก้ระบบนี้ได้ดีไม่ได้แปลว่าสนามทุกตัวสอดคล้องกันแล้ว
3. ระหว่างรอบ SIMPLE: ปรับความเร็ว ความดัน และตัวแปรความปั่นป่วนร่วมกัน โครงการนี้ตรวจ initial residual ของทั้งหกสมการทุกครั้งใน 100 รอบท้าย และตรวจความคงที่ของแรงอีกเงื่อนไขหนึ่ง

สำหรับ Spalding ความสัมพันธ์พื้นฐานคือ

y⁺ = u⁺ + [exp(κu⁺) − 1 − κu⁺ − (κu⁺)²/2 − (κu⁺)³/6] / E
u⁺ = Urel / uτ       y⁺ = y uτ / ν       uτ = √(|τw| / ρ)

Urel คือขนาดความเร็วสัมพัทธ์ที่ใช้ในแบบจำลองผนัง, y คือระยะจากผิว, ν คือความหนืดจลน์, τw คือความเค้นเฉือนที่ผิว และ ρ คือความหนาแน่น ส่วน κ และ E เป็นค่าคงที่ของ wall law; y⁺ และ u⁺ ไม่มีหน่วย สมการนี้ช่วยปิดแบบจำลองบริเวณใกล้ผิว แต่ไม่แทนสมการนาเวียร์–สโตกส์ทั้งสนาม ดูรายละเอียดและข้อจำกัดใน เอกสาร Spalding ของ OpenCFD

ตัวอย่างหน่วย: ถ้า uτ = 0.6 m/s, y = 0.0005 m และ ν = 1.5 × 10⁻⁵ m²/s จะได้ y⁺ = 20 หาก Urel = 6 m/s จะได้ u⁺ = 10 ตัวเลขนี้ใช้ฝึกแทนค่าเท่านั้น ยังไม่ได้ยืนยันว่าคู่นี้เป็นคำตอบของสมการ Spalding หรือเหมาะกับกริดของเรา

ต้องตรวจรุ่นที่ใช้จริง: ตารางในคู่มือรุ่น 2312 ระบุค่าเริ่มต้น tolerance 0.0001 แต่ ซอร์ส constructor รุ่น 2512 ระบุ 0.01 และ maxIter 10 จึงไม่ควรคัดลอกตัวเลขจากคู่มือข้ามรุ่นโดยไม่ตรวจ เอกสารยังอธิบายว่าการ restart ของเงื่อนไขนี้อาจให้ค่า nut ต่างเล็กน้อย และยกตัวอย่าง tolerance 10⁻⁷ กับ maxIter 100 สำหรับการคำนวณภายในที่เข้มขึ้น นี่เป็นแนวทางตั้งการทดลองแยก ไม่ใช่หลักฐานว่าปัญหาปัจจุบันเกิดจาก tolerance หรือว่าปรับแล้วจะผ่านแน่นอน

ฝึกแยกประเด็น: หากแรงขับเปลี่ยนเพียง 0.1% แต่ initial residual ของ Uz ยังสูงกว่าเกณฑ์ จะรับผลได้หรือไม่? สำหรับ protocol นี้ยังรับไม่ได้ และการเพิ่ม maxIter ของ wall function ไม่ใช่การเพิ่มจำนวนรอบ SIMPLE ต้องบันทึกให้ชัดว่ากำลังเปลี่ยนตัวแปรใด

เพิ่มจำนวนรอบแล้วผ่านหรือไม่: ผลทดสอบถึงรอบ 6,000

เราเริ่มทั้งสองกรณีใหม่จาก checkpoint 4,000 เดียวกัน กำหนดหยุดที่ 6,000 ก่อนเริ่มรัน และใช้ OpenFOAM 2512 แบบขนาน 8 ranks เท่ากันตาม แนวทาง OpenCFD คงกริด สมการ เงื่อนไขขอบเขตอื่นและเกณฑ์รับผลไว้เดิม จึงต้องเทียบภายในคู่การทดลองใหม่นี้ ไม่ใช้ความต่างระหว่าง serial กับ parallel เป็นผลทางฟิสิกส์

ผลที่สรุปได้: คู่นี้ผ่านเกณฑ์เชิงตัวเลขของโครงการทั้งคู่แล้ว บนกริดกลางเดิม การเปลี่ยนเป็น Spalding ทำให้แรงขับลด 0.991% และแรงบิดเพิ่ม 2.220% โดยใช้สูตร Δ% = 100 × (ค่าทดสอบ / ค่าควบคุม − 1) ค่านี้เป็น wall-model sensitivity บนกริดเดียว ยังไม่พิสูจน์ว่าแก้ความไม่สอดคล้องของสามกริด ไม่ใช่ค่า uncertainty และยังไม่ใช่การยืนยันว่าตรงกับการทดลองจริง

อ่านกราฟให้ถูก: เส้นบนคือ initial residual ที่สูงที่สุดจากทั้งหกสมการในแต่ละรอบ เส้นประคือเกณฑ์ 10⁻⁵ และพื้นที่แรเงาคือ 100 รอบท้ายที่กำหนดไว้ก่อนรัน แกนนอนเป็นรอบ SIMPLE ไม่ใช่วินาทีบินจริง การเห็นเส้นต่ำกว่าเกณฑ์เพียงบางช่วงยังไม่พอ ส่วนกราฟแรงทั้งสองรูปช่วยตรวจคำตอบที่เราจะนำไปใช้ร่วมด้วย

ตรวจย้อนกลับได้: solver เขียน End ครบทั้งคู่ที่ 6,000 และประวัติแรงกับ residual ตรวจซ้ำผ่าน สคริปต์ครอบงานเดิมเกิดข้อผิดพลาดหลัง solver จบจากการแก้ไฟล์ขณะยังทำงาน จึงบันทึก wrapper exit 1 ไว้ตามจริง ไม่แทนด้วยสถานะสำเร็จ สคริปต์แก้ให้ส่งต่อการทำงานด้วย exec และรันทดสอบสั้นแยกใหม่แล้วได้ exit 0 เหตุการณ์นี้ไม่ถูกนำไปปะปนกับค่าฟิสิกส์หรือเกณฑ์รับผล

ดาวน์โหลด หลักฐานและกราฟ wall-function study สำหรับ CSV ประวัติทั้งหมด log, input hashes, สถานะการรัน และโค้ดสร้างกราฟ ชุดผล serial เดิมที่ 4,500 ยังอยู่ใน results/wall-trial/ และยังมี qualified delta เป็น null; ชุดใหม่แยกอยู่ใน results/wall6000/ เพื่อรักษาประวัติทั้งผลที่ผ่านและไม่ผ่าน

ชั้นกริดใกล้ผิว: ตั้งไว้ 8 ชั้น ได้จริงเท่าไร

ชิดผิวใบพัด ความเร็วเปลี่ยนเร็วในระยะสั้น เราจึงทดลองเพิ่มเซลล์บาง ๆ เรียงตามผิว เรียกว่า layer mesh ก่อนนำไปแก้สมการการไหล ชุดนี้ใช้ OpenFOAM 2512 และสำเนากริดเดิมสามระดับ เก็บกริดต้นฉบับไว้ ตั้งเป้า 8 ชั้น อัตราขยายความหนา 1.2 และความหนาชั้นแรก 100, 50, 25 ไมโครเมตร ตามลำดับ

เราเปรียบเทียบวิธีขยับกริดเพื่อเปิดที่ให้ชั้นใหม่สองแบบ ได้แก่ displacementMedialAxis ซึ่งเป็นค่าเริ่มต้น และ displacementMotionSolver ที่แก้สมการการขยับกริด วิธีหลังต้องตั้งค่าใน fvSolution และ fvSchemes เพิ่มด้วย ตาม คู่มือ OpenCFD เรื่อง Layer addition ผลต่อไปนี้มาจาก log การรันจริง ไม่ใช่ค่าที่คาดจากคู่มือ

เกณฑ์ที่ประกาศก่อนรัน: checkMesh ต้องผ่าน เซลล์ไม่เกิน 3 ล้าน และอย่างน้อย 95% ของหน้าผิวที่ร้องขอต้องสร้างชั้นได้ เกณฑ์ 95% เป็นเกณฑ์คัดกรองของโครงการนี้ ทั้งหกกรณีผ่าน checkMesh แต่สองกรณีไม่ผ่านความครอบคลุม เราเก็บผลเหล่านั้นไว้ด้วย

คำว่า 95.31% ไม่ได้หมายถึงผิว 95.31% มีครบ 8 ชั้น เพราะนับจำนวนหน้าที่มีการสร้างชั้นอย่างน้อยบางส่วน ไม่ได้ถ่วงด้วยพื้นที่ และค่าเฉลี่ยจริงของกริดหยาบมีเพียง 2.82 ชั้น การดูแค่ข้อความ Mesh OK หรือค่า nSurfaceLayers 8 จึงยังไม่เพียงพอ ชุดนี้ยังไม่ได้แก้สมการการไหลบนกริดที่มีชั้น จึงยังไม่มีผลแรงหรือ yPlus ใหม่

คำนวณความหนาก่อนสร้าง แล้วตรวจค่าที่ได้จริง

กำหนดความหนาชั้นแรกเป็น t₁ อัตราขยายเป็น r และจำนวนชั้นเป็น n:

tᵢ = t₁rⁱ⁻¹ และ H = t₁(rⁿ − 1)/(r − 1) เมื่อ r ≠ 1; ถ้า r = 1 จะได้ H = nt₁

ตัวอย่างกริดกลาง: t₁ = 50 µm, r = 1.2, n = 8 ให้ความหนารวมตามแบบ H ≈ 0.825 mm แต่ log ของ Motion solver รายงานความหนารวมประมาณ 0.679 mm พร้อมจำนวนชั้นเฉลี่ย 4.81 ค่ารายงานทั้งสองใช้บอกว่าการสร้างจริงถูกตัดทอน ไม่ควรแทนค่าจำนวนชั้นเฉลี่ยกลับเข้าอนุกรมแล้วคาดว่าจะได้ความหนารายงานตรงกันทุกจุด

สำหรับเซลล์ที่เกือบตั้งฉากกับผิว อาจประมาณระยะศูนย์กลางชั้นแรก y₁ ≈ t₁/2 เพื่อเริ่มวางแผน y⁺ = y₁uτ/ν แต่ uτ ต้องได้จากแรงเฉือนที่ผิวหลังแก้สมการ จึงห้ามเรียก yPlus ที่คาดไว้ว่าเป็นค่าที่ตรวจวัดได้แล้ว

ทำให้ละเอียดทั้งโดเมน: เหตุใดต้องตรวจกริดที่สร้างจริง

เราตรึงรูปทรง โดเมน 0.84 m และระดับ refinement ผิว/rotor/wake 5/3/2 เปลี่ยนเฉพาะจำนวนเซลล์พื้นฐาน N แล้วใช้ h = 0.84/N กับความสูงชั้นแรก t₁ = 50 µm × 32/N เกณฑ์เดิมทุก protocol คือ checkMesh ผ่าน, เซลล์ไม่เกิน 3 ล้าน และสร้างชั้นบนอย่างน้อย 95% ของหน้าผิวที่ร้องขอ

N=28 ผ่านสัดส่วน coverage แต่ไม่ผ่าน checkMesh จึงต้องตัดออกตาม protocol เดิม จากนั้นได้ลงทะเบียน protocol v3 ก่อนรันชุดใหม่ โดยไม่ปรับ threshold; ทั้ง N=30/32/36/44 ผ่าน construction gate จาก log และมีหลักฐาน mesh hash ครบ ค่า coverage นับจำนวนหน้า ไม่ได้ถ่วงพื้นที่ และหมายถึงมีชั้นบางส่วน ไม่ได้หมายความว่าแต่ละหน้ามีครบ 8 ชั้น

ขั้น solve จบครบทั้งสี่กริดแล้ว: ทั้งสี่กริดเริ่มจาก fields ตั้งต้นเดียวกัน ใช้ OpenFOAM 2512 simpleFoam/steady SIMPLE, Open MPI 4.1.6 จำนวน 8 ranks ถึง iteration 6,000 และเกณฑ์ residual/force/y⁺ เดียวกัน หาก launcher ชนเพดานเวลา จะต่อจาก checkpoint ล่าสุดโดยไม่เปลี่ยนสมการหรือเกณฑ์ และเก็บ log/hash แยกช่วงตาม protocol addendum ผลตรวจพบว่ามีเพียง N32 ผ่านทุกเกณฑ์ จึงไม่รับเป็น qualified grid sensitivity และงด Richardson/GCI การคัดกรองกริดอย่างเดียวไม่ใช่ผลแรงและไม่ใช่การยืนยันทางกายภาพกับ ENOLA

Checkpoint คืออะไร: เป็นสนามการไหลที่บันทึกไว้เพื่อเริ่มคำนวณต่อ ตัวอย่างเช่น ถ้างานหยุดระหว่างรอบ 1,167 แต่เขียนครบล่าสุดที่ 1,150 จะเริ่มต่อจาก 1,150 โดยต้องตรวจว่าไฟล์ทั้ง 8 ส่วนอยู่ในรอบเดียวกันและครบถ้วน การรันต่อยังต้องผ่านเกณฑ์เดิม และไม่รับประกันผลเหมือนรันรวดเดียวทุกบิต เพราะข้อมูล ASCII มีจำนวนหลักที่บันทึกจำกัด อ่าน บันทึกวิธีตรวจ checkpoint และรับผล

ผลที่ตรวจแล้ว: N30 และ N32 รันครบ 6,000 รอบ มี checkpoint ครบทั้ง 8 ranks และผ่าน mesh/ความนิ่งของแรงเหมือนกัน แต่ได้ผล residual ต่างกัน:

นี่เป็นตัวอย่างว่า แรงที่นิ่งไม่ได้แทนการตรวจ residual ของสนามทั้งหมด N30 ยังไม่ผ่านจึงไม่ใช้ส่วนต่างแรง N30–N32 เป็นผลยืนยันความไวต่อกริด แม้ N32 ผ่านทุกเกณฑ์ของกรณีเดียว ก็ยังไม่ยืนยันว่าไม่ขึ้นกับกริดหรือถูกต้องตรงกับการทดลองจริง

ผลสุดท้ายครบสี่กริด

เปิดกราฟขนาดเต็ม · รายงานสรุปผลและข้อจำกัด · คู่มือเริ่มใช้ชุดดาวน์โหลด

ปิดการรันครบตามแผน แต่ยังไม่ยืนยัน grid independence: N32 ผ่านเกณฑ์รายกรณี ส่วน N30/N36/N44 ไม่ผ่าน residual และ N36 ไม่ผ่านความนิ่งของแรงด้วย ตรวจแล้ว inputs และ protocol ร่วมตรงกัน แต่เงื่อนไขนี้ไม่ได้ชดเชยกรณีที่ไม่ผ่าน จึงไม่รายงาน qualified grid sensitivity หรือ GCI

N44 มีช่วงแกว่งแรงเพียง 0.00378% และแรงบิด 0.01582% แต่ residual สูงสุด 2.449×10⁻⁵ เกินเกณฑ์10⁻⁵ จึงเป็นอีกตัวอย่างว่าต้องตรวจทั้งสนามและปริมาณที่ต้องการใช้ ผลนี้ไม่ใช่ formal physical validation กับ ENOLA

ระวังการเทียบกริดที่ใกล้กันเกินไป: คู่ N30→N32 มีอัตราส่วนระยะ r = h₃₀/h₃₂ = 32/30 ≈ 1.067 ขณะที่ NASA แนะนำ r อย่างน้อย 1.1 เพื่อช่วยแยกผลของกริดจากความคลาดเคลื่อนอื่น ดังนั้นความต่างของแรงที่น้อยในคู่นี้เพียงคู่เดียว ยังไม่แปลว่าแบบจำลองแม่นยำ เราเก็บทุกกริดตามแผนและพิจารณาลำดับทั้งหมดก่อนสรุป

ทดลองด้วยตนเอง: ดาวน์โหลด ชุด Layer mesh: โค้ดและหลักฐานการรัน เพื่อเปรียบเทียบ summary กับ log รายกรณี อ่าน protocol v2, protocol v3 และ รายงานผลคัดกรอง v3 ข้อมูลรูปทรงต้นทางต้องจัดหาตามสิทธิ์ของ ENOLA แยกต่างหาก ชุดดาวน์โหลดไม่รวมรูปทรงหรือผลทดลองที่สร้างขึ้นเอง

แบบฝึกหัด: ถ้าลด t₁ ครึ่งหนึ่งโดยคง r และ n เดิม H จะเปลี่ยนอย่างไร? คำตอบคือครึ่งหนึ่ง แต่ยังสรุปไม่ได้ว่า yPlus หรือความคลาดเคลื่อนแรงจะลดครึ่งหนึ่ง เพราะต้องแก้สมการและตรวจสนามการไหลก่อน

อย่าสับสนการฟิตกับการตรวจยืนยัน: ในเครื่องมือ GCI การฟิตอันดับ p จากสามจุด แล้วนำสามจุดเดิมกลับไปทดสอบความสอดคล้อง เป็นเพียงการตรวจสมการคำนวณ ไม่ใช่หลักฐานอิสระว่าถึง asymptotic range โค้ดจึงแยก richardson_fit_consistency_ratio ออกจาก asymptotic_range_verified ซึ่งยังเป็น false และต้องพิจารณากริด/หลักฐานเพิ่มเติม ผลสามกริด MS1101 ที่มีอยู่ยังถูกปฏิเสธก่อนขั้นฟิตเพราะไม่ monotonic

แบบฝึกทบทวนและแนวคำตอบ

โจทย์: ในตัวอย่างสมมติ หากเพิ่ม n เป็น 200 รอบ/s โดยตรึง ρ, D และสมมติ C_T, C_P คงเดิม จะได้ T และ P เท่าไร? สมมติฐานใดต้องตรวจใหม่?

แนวคำตอบ: T เพิ่มสี่เท่าเป็น 7.84 N และ P เพิ่มแปดเท่าเป็น 156.8 W แต่ต้องตรวจ Reynolds, tip Mach, J และช่วงข้อมูล เพราะสัมประสิทธิ์อาจไม่คงที่ การคูณตามสูตรไม่ได้ขยายช่วง validation

สิ่งที่พร้อมใช้คือสูตร แบบบันทึก โค้ด adapter ข้อมูล geometry อ้างอิง ชุดตรวจ MRF กับคำตอบสมการ และกรณี MS1101 ที่เทียบข้อมูลทดลองตรงแถว 5,088 RPM ได้ พร้อมหลักฐานการลู่เข้าและ sensitivity แยกตามกริด/โดเมน ดูไฟล์ทำซ้ำได้ใน ชุดข้อมูลเปรียบเทียบ ENOLA และ ชุดทดลอง CFD; ผลต่างที่พบยังไม่ผ่านการยืนยันทางกายภาพแบบมี uncertainty เพราะ metadata จุดทดลองและ specimen traceability ยังขาด โมเดลหลายมุมและ angular rates สำหรับการบินทั่วไปก็ยังเป็นงานต่อยอด กลับไปตรวจ กระบวนการ CFD หรือทบทวน ภาพรวมบทเรียน

อ่านผล CFD ให้เป็น: คำนวณครบไม่ได้หมายความว่าผ่าน

กรณี N36 ต่อไปนี้เป็นผลจริงที่รันครบ 6,000 รอบ แต่ไม่ผ่านเกณฑ์ ใช้เรียนรู้การตรวจผลก่อนนำแรงไปใช้ใน DroneSim ขณะนี้ N44 จบแล้ว แต่ไม่ผ่าน residual เช่นกัน อ่านผลครบสี่กริดและข้อจำกัดด้านล่าง

คำศัพท์สำหรับเริ่มอ่านกราฟ

• Mesh (กริด): เซลล์เล็ก ๆ ที่แบ่งพื้นที่ของไหลเพื่อคำนวณ การผ่าน checkMesh ช่วยตรวจคุณภาพเชิงเรขาคณิต แต่ไม่ยืนยันว่ากริดละเอียดพอ
• Residual: ตัวชี้วัดความไม่สมดุลของสมการที่ solver กำลังแก้ ค่าที่แสดงเป็น initial residual ก่อนการแก้ระบบเชิงเส้นแต่ละครั้ง ไม่ใช่เปอร์เซ็นต์ความผิดพลาดของแรงเมื่อเทียบของจริง
• Convergence (การลู่เข้า): ต้องระบุว่าพูดถึงการวนรอบ การลดขนาดกริด หรือการเทียบผลทดลอง แต่ละอย่างใช้หลักฐานต่างกัน
• y⁺: ระยะจากผนังในรูปไร้มิติ ใช้พิจารณาความละเอียดใกล้ผนังร่วมกับ wall treatment ไม่มีค่าเดียวที่เหมาะกับทุกแบบจำลอง
• Verification / validation: อย่างแรกตรวจว่าคำนวณตามแบบจำลองได้ถูกต้องเพียงใด อย่างหลังตรวจว่าแบบจำลองแทนสิ่งที่วัดจริงได้ดีเพียงใด

เริ่มจากภาพ

กราฟบนแสดงแรง กราฟกลางแสดงแรงบิด กราฟล่างแสดง initial residual แยกหกตัวแปร สีและชื่อใน legend ใช้บอกตัวแปร ไม่ใช่คะแนนคุณภาพ เส้นประสีทองคือค่าเฉลี่ย และเส้นจุดแนวตั้งแบ่งข้อมูลออกเป็นสองครึ่ง ตัวเลขบนแกนนอนคือรอบ SIMPLE ไม่ใช่วินาที จึงนำช่วงการแกว่งไปแปลงเป็นความถี่การสั่นจริงของโดรนไม่ได้

เส้นทุกเส้นมาจากข้อมูลดิบ ไม่มีการทำ smoothing สำหรับ pressure ซึ่งแก้สมการสองครั้งต่อ iteration ใช้ residual ค่าสูงสุดของรอบนั้น กราฟนี้เป็นหลักฐานทางตัวเลข ไม่ใช่ภาพการวัดจากโดรนจริง

จบโปรแกรม กับ ลู่เข้า ต่างกันอย่างไร

คำว่า End หมายความว่า solver จบการทำงานตามการตั้งค่า ส่วนการยอมรับผลต้องตรวจเกณฑ์เพิ่มเติม แนวทาง NASA แนะนำให้ติดตามทั้ง residual และผลที่ต้องการใช้งาน เช่น แรงยกหรือแรงต้าน เพราะแต่ละปริมาณอาจลู่เข้าด้วยอัตราต่างกัน (NASA: Iterative Convergence)

เกณฑ์ของ การศึกษาชุดนี้ กำหนดไว้ล่วงหน้า: residual ทั้งหกตัวแปรต้องต่ำกว่า 10⁻⁵ ทุก iteration ใน 100 รอบสุดท้าย และแรงกับแรงบิดต้องผ่านตัวชี้วัดทั้งสองด้านล่าง ตัวเลขเหล่านี้ไม่ใช่ค่ามาตรฐานสากลที่นำไปใช้กับทุก solver ได้โดยอัตโนมัติ

สมการที่ใช้

ให้ qᵢ เป็นแรงหรือแรงบิดในแต่ละรอบ และมีข้อมูล 100 ค่า:

ค่าเฉลี่ย: q̄ = (q₁ + q₂ + … + q₁₀₀) / 100

ช่วงแกว่ง: A_pp = (q_max − q_min) / |q̄| × 100%

ความต่างสองครึ่ง: D_half = |q̄₅₁:₁₀₀ − q̄₁:₅₀| / |q̄| × 100%

A_pp วัดช่วงแกว่งสูงสุดถึงต่ำสุดเทียบกับค่าเฉลี่ย ส่วน D_half วัดการเปลี่ยนของค่าเฉลี่ยระหว่างครึ่งแรกและครึ่งหลัง เกณฑ์ที่ใช้คือ A_pp < 2% และ D_half < 1% ทั้งแรงและแรงบิด ถ้าค่าเฉลี่ยใกล้ศูนย์ สัดส่วนนี้จะตีความยาก ต้องกำหนดสเกลอ้างอิงที่เหมาะสมใน protocol ก่อนใช้ ไม่ควรเปลี่ยนสเกลหลังเห็นผล

ตัวอย่างแรง: (2.8154121841 − 2.4917468178) / 2.702727095717 × 100 = 11.9755%

แม้ค่าเฉลี่ยสองครึ่งของแรงจะใกล้กัน แต่กราฟขึ้นแล้วลง ทำให้ความแตกต่างระหว่างค่าเฉลี่ยเล็กได้ จึงต้องตรวจช่วงแกว่งด้วย ในกรณีนี้ residual ยังผ่าน 0/100 รอบ โดยค่าสูงสุดเท่ากับ 0.007300501907 จึงมีหลักฐานไม่ผ่านมากกว่าหนึ่งด้าน

ผังการตัดสินใจสำหรับนักศึกษา

เปิดแผนภาพขนาดเต็ม

การมี y+ รอบสุดท้ายเป็นเพียงเกณฑ์ความครบของข้อมูลในชุดนี้ ไม่ได้ยืนยันว่า wall treatment เหมาะสมทุกตำแหน่ง และการผ่านรายกรณีก็ยังไม่ใช่การยืนยันความแม่นยำทางกายภาพ

แบบฝึกพร้อมแนวตอบ

1. ถ้าพิจารณาแรงเฉพาะ D_half จะตัดสินผิดอย่างไร? อาจคิดว่าแรงนิ่งแล้ว ทั้งที่ A_pp เกินเกณฑ์เกือบหกเท่า
2. N36 จบพร้อม exit 0 แล้วนำค่าเฉลี่ยไปคำนวณ GCI ได้หรือไม่? ยังไม่ได้ เพราะรายกรณีไม่ผ่าน อีกทั้ง GCI ต้องมีเงื่อนไขด้านกริดและพฤติกรรมการลู่เข้าเพิ่มเติม (NASA: Grid Convergence)
3. กราฟขึ้นลงนี้พิสูจน์ว่าเกิด vortex shedding จริงหรือไม่? ไม่พิสูจน์ รอบ SIMPLE ไม่ใช่ physical time และข้อมูลนี้ยังแยกสาเหตุทางตัวเลขจากพฤติกรรมทางฟิสิกส์ไม่ได้
4. หากต้องศึกษาเพิ่ม ควรแก้ค่า solver แล้วแทนผลเดิมเลยหรือไม่? ไม่ควร ต้องเก็บผลเดิม ตั้งสมมติฐานและ protocol ใหม่ก่อนทดลอง เพื่อเปรียบเทียบอย่างตรวจสอบย้อนกลับได้

ตารางตรวจรับก่อนนำแรงไปใช้

ลงมือทำ: ดาวน์โหลดชุด Layer mesh แล้วเปิดโฟลเดอร์ results/layers/v3-global36-solver ตรวจตัวเลขกับตารางนี้ ชุดปัจจุบันมีผลครบ N30/N32/N36/N44 พร้อม QUICKSTART-TH.md เปิดกราฟขนาดเต็ม หรือ อ่านบทเรียนพร้อมวิธีสร้างกราฟซ้ำ

บรรณานุกรม

• OpenFOAM — nutUSpaldingWallFunction

• OpenFOAM — Layer addition in snappyHexMesh

• NIST TN 1297 — Combined Standard Uncertainty

• NASA NPARC — Examining Iterative Convergence

• NASA NPARC — Examining Spatial (Grid) Convergence

• Celik et al. — Procedure for Estimation and Reporting of Uncertainty Due to Discretization in CFD Applications (2008)

• ENOLA — improved propeller test bench

• OpenFOAM — Residuals and normalization

• OpenFOAM 2512 — solverInfo source and residual fields

• OpenFOAM — probes function object

• OpenFOAM 2512 — yPlus implementation

• García-Tíscar et al. — On the effect of popular additive manufacturing technologies (AST 169, 2026)

• UPV — ENOLA numerical and experimental propeller database, v1 (2026)

• ENOLA — Data licensing and reuse

• Aular et al. — Assessing the Fidelity of Steady-State MRF Modeling (2026)

• Brandt, Deters, Ananda, Dantsker & Selig — UIUC Propeller Database, Vols 1–4

• BYU FLOW Lab — FLOWUnsteady

• NASA NPARC — Overview of CFD Verification and Validation

• UIUC — APC Thin Electric 10×5 geometry table

• APC — Propeller Geometry Data และ archive 202602, ตรวจไฟล์ 10x5E-PERF.PE0 วันที่ 26 กันยายน 2569

ตรวจแหล่งต้นทาง 26 กันยายน 2569