การคำนวณเชิงตัวเลข (๒๘) การแก้ปัญหาสมการอนุพันธ์สามัญ ด้วย ODE solvers ของ GNU Octave ตอนที่ ๓ MO Memoir : Saturday 15 February 2563

การเกิดปฏิกิริยาใน fixed-bed catalystic reactor นั้น ถ้า heat of reaction มีค่าต่ำ ก็จะประมาณได้ว่า reactor ตัวนั้นทำงานแบบอุณหภูมิคงที่ (isothermal) และถ้า heat of reaction มีค่าไม่สูงมาก ก็จะออกแบบให้ reactor นั้นทำงานแบบ adiabatic (ไม่มีการระบายความร้อนออกหรือให้ความร้อนแก่ reactor โดยตรง) กล่าวคือปล่อยให้ปฏิกิริยาดำเนินไปข้างหน้าเรื่อย ๆ จนได้ค่า conversion ที่ต้องการ แต่ถ้าหากว่าอุณหภูมิของระบบนั้นสูง/ต่ำเกินไปก่อนจะได้ค่า conversion ที่ต้องการ ก็ต้องมีการนำเอา reactor หลายตัวมาต่ออนุกรมกัน โดยมีเครื่องแลกเปลี่ยนความร้อน (ที่ระบาย/ป้อนความร้อน) อยู่ระหว่าง reactor แต่ละตัว

  

ในกรณีของการทำงานแบบ adiabatic นั้น ปฏิกิริยาดูดความร้อน (endothermic reaction) มีแนวโน้มที่จะหยุดตัวเอง เพราะเมื่อปฏิกิริยาดำเนินไปข้างหน้าเรื่อย ๆ ความเข้มข้นของสารตั้งต้นจะลดต่ำลง และอุณหภูมิของระบบก็จะลดต่ำลง สองปรากฏการณ์นี้ส่งผลให้อัตราการเกิดปฏิกิริยาลดต่ำลง 

    

ในกรณีของปฏิกิริยาคายความร้อน (exothermic reaction) นั้นเมื่อปฏิกิริยาดำเนินไปข้างหน้า ความเข้มข้นสารตั้งต้นที่ลดต่ำลงจะดึงให้อัตราการเกิดปฏิกิริยาลดลง แต่อุณหภูมิที่สูงขึ้นจากความร้อนที่คายออกมาจะทำให้ปฏิกิริยาเกิดเร็วขึ้น สองปรากฏการณ์นี้ให้ผลที่ขัดแย้งกัน ถ้าหากการทำปฏิกิริยานั้นเกิดขึ้นในเครื่องปฏิกรณ์แบบ adiabatic (ไม่มีการระบายความร้อนออก) ผลที่เกิดขึ้นจะเกิดในทิศทางเดียวคืออุณหภูมิของระบบจะเพิ่มสูงขึ้นเรื่อย ๆ จนอาจเหนือการควบคุม (ที่เรียกว่าปฏิกิริยาเกิดการ runaway)

   

เพื่อที่จะ "ลด" โอกาสเกิดเหตุการณ์ดังกล่าว ในกรณีของปฏิกิริยาคายความร้อนสูงจึงมักดำเนินการในรูปแบบที่เรียกว่า non-isothermal non-adiabatic คือมีการระบายความร้อนออกจาก reactor โดยตรง ดังนั้นอุณหภูมิใน reactor จึงขึ้นกับอัตราการคายความร้อนจากปฏิกิริยาและอัตราการระบายความร้อนออก เนื่องจากอัตราการระบายความร้อนออกนั้นขึ้นอยู่กับผลต่างระหว่างอุณหภูมิภายใน reactor กับอุณหภูมิแหล่งรับความร้อน ดังนั้นในช่วงที่ปฏิกิริยาเพิ่งเริ่มเกิด (ด้านขาเข้าของ reactor) อุณหภูมิใน reactor (ที่เพิ่มขึ้นจากความร้อนที่ปฏิกิริยาคายออกมา) ยังไม่สูง ผลต่างอุณหภูมิจึงยังไม่สูง อัตราการระบายความร้อนออกจึงต่ำ ส่งผลให้อุณหภูมิใน reactor เพิ่มสูงขึ้นเรื่อย ๆ เมื่อปฏิกิริยาดำเนินไปข้างหน้า แต่เมื่ออุณหภูมิใน reactor สูงถึงระดับหนึ่งจนทำให้ผลต่างอุณหภูมิสูงมากพอ อุณหภูมิใน reactor จะไม่เพิ่มขึ้นและจะเริ่มลดต่ำลง (ผลจากความเข้มข้นสารตั้งต้นที่ลดต่ำลง)

  

ตัวอย่างของปฏิกิริยาประเภทนี้ได้แก่ปฏิกิริยา partial oxidation ไฮโดรคาร์บอนไปเป็นสารประกอบ oxygenate ที่ใช้อากาศเป็นสารออกซิไดซ์ เนื่องจากออกซิเจนในอากาศมีความว่องไวต่ำในการออกซิไดซ์สารตั้งต้น การเริ่มการเกิดปฏิกิริยาจึงต้องการอุณหภูมิที่สูงมากพอ (ไม่เช่นนั้นปฏิกิริยาจะไม่เกิด หรือเกิดช้ามากจนไม่ได้ conversion ที่ต้องการ) อุณหภูมิของแหล่งรับความร้อนจึงต้องสูงตามไปด้วยเพื่อทำให้สารตั้นต้นที่ป้อนเข้ามานั้นมีอุณหภูมิสูงพอที่จะเกิดปฏิกิริยาได้ แต่ในขณะเดียวกันต้องไม่สูงจนทำให้การดึงความร้อนออกจาก reactor ต่ำเกินไป แหล่งรับความร้อนจากปฏิกิริยาที่ใช้กันจึงเป็นพวก molten salt หรือเกลือหลอมเหลวที่เกิดจากสารผสมในสัดส่วนที่เหมาะสมระหว่างเกลือบางชนิด (ทำให้มีจุดหลอมเหลวต่ำเพื่อง่ายต่อการหลอมให้เป็นของเหลว) ตัวอย่างอุณหภูมิของ molten salt นี้อาจอยู่ที่ระดับประมาณ 250-450ºC ขึ้นอยู่กับปฏิกิริยา สภาวะการทำปฏิกิริยา ความว่องไวของตัวเร่งปฏิกิริยา ความเข้มข้นสารตั้งต้น ฯลฯ

ที่ใช้คำว่า "ลด" ในย่อหน้าข้างบนนั่นก็เพราะว่าแม้ว่าจะมีการดึงความร้อนออกจาก reactor โดยตรงมันก็ยังมีโอกาสที่จะเกิดการ runaway อยู่ถ้าหากว่ามีความร้อนสะสมใน reactor มากเกินไป สิ่งที่ทำให้ความร้อนสะสมใน reactor ได้มีทั้ง ความเข้มข้นสารตั้งต้นที่สูง อัตราการไหลที่ต่ำ อุณหภูมิของแหล่งรับความร้อนที่สูงเกินไป อุณหภูมิขาเข้าที่สูงเกินไป และความว่องไวของตัวเร่งปฏิกิริยาที่สูงเกินไป ดังนั้นช่วงการทำงานของ reactor เหล่านี้จึงต้องอยู่ในช่วงที่เหมาะสมที่สามารถทำให้ได้ค่า conversion ที่สูงมากพอโดยไม่เกิดการ runaway

  

การคาดการณ์การเกิด runaway ทำได้ด้วยการแก้สมการดุลมวลสารและพลังงานของระบบ และทดลองปรับเปลี่ยนค่าพารามิเตอร์ต่าง ๆ เพื่อตรวจสอบดูว่าจะเกิดอะไรขึ้น รูปแบบง่าย ๆ รูปแบบหนึ่งของสมการดุลมวลสารและพลังงานแสดงไว้ข้างล่าง (เรียกว่าแบบจำลอง 1 มิติหรือ 1-dimensional model เพราะสมมุติให้การเปลี่ยนแปลงเกิดขึ้นตามความยาวของ reactor เท่านั้น)

  

ในที่นี้พจน์ C คือความเข้มข้น, T คืออุณหภูมิของ reactor, Tc คืออุณหภูมิของแหล่งรับความร้อน, z คือระยะทางตามความยาวของ reactor, exp(a - b/T) คือ rate constant ที่เปลี่ยนตามอุณหภูมิ, H เป็นผลรวมของความร้อนที่เกิด และ U คือค่าสัมประสิทธิ์การถ่ายเทความร้อนระหว่างภายใน reactor กับแหล่งรับความร้อน สมการดุลมวลสารนั้นประกอบด้วยพจน์ของอัตราการหายไปของสารตั้งต้น ในขณะที่สมการดุลความร้อนนั้นประกอบด้วยความร้อนที่เกิดจากปฏิกิริยาและความร้อนที่มีการดึงออก ระบบสมการทั้งสองเป็นระบบสมการอนุพันธ์สามัญที่มีพจน์ที่เป็น non-linear ร่วมอยู่

   

ในที่นี้การหาคำตอบระบบ 2 สมการข้างต้นจะใช้คำสั่ง lsode, ode45 และ ode15s ของ GNU Octave เปรียบเทียบกันโดยจะหาคำตอบในช่วง z = 0.0 - 0.2 คู่มือของโปรแกรมนั้นไม่ได้ให้รายละเอียดว่าคำสั่งดังกล่าวหาคำตอบด้วยวิธีใด คำสั่ง lsode นั้นให้กำหนดช่วงที่ต้องการหาคำตอบและระบุด้วยว่าต้องการซอยช่วงดังกล่าวออกเป็นช่วงย่อยที่ละเอียดมากน้อยเท่าใด (ยิ่งละเอียดมากก็ยิ่งเข้าใกล้คำตอบที่ถูกต้องมากขึ้น) ในที่นี้ได้ทดลองเลือกแบ่งย่อยเป็น 200 ช่วง (ในโปรแกรมต้องใส่ 201 จุด) ค่าพารามิเตอร์ต่าง ๆ ที่ใช้คือ a = 20, b = 14000, H = 10000, T ขาเข้า = 350ºC, C ขาเข้า = 1.0, และ U = 500 ตัวโปรแกรมที่เขียนขึ้นแสดงไว้ในรูปที่ ๑ ส่วนผลการคำนวณที่ได้ที่ใช้ค่า Tc = 414.0, 414.7 และ 414.8ºC แสดงไว้ในรูปที่ ๒ - ๖ โดยรูปที่ ๒ และ ๓ นั้นได้มาจากคำสั่ง lsode รูปที่ ๔ และ ๕ ได้มาจากคำสั่ง ode45 และรูปที่ ๖ ได้มาจากคำสั่ง ode15s

   

ในกรณีของคำสั่ง lsode และ ode45 นั้นพบว่าสองตัวนี้ให้คำตอบที่เหมือนกันคือการทำงานของ reactor จะเป็นถ้าอุณหภูมิของแหล่งรับความร้อนนั้นไม่เกิน 417.7ºC แต่เมื่อเพิ่มอุณหภูมิของแหล่งรับความร้อนเป็น 417.8ºC พบว่าปฏิกิริยาจะเกิดการ runaway ซึ่งเห็นได้จากผลการคำนวณนั้นให้ค่าอุณหภูมิที่เพิ่มสูงขึ้นกระทันหันโดยให้ค่าอุณหภูมิสูงสุดอยู่ที่หลายพันองศาเซลเซียส การแกว่งของคำตอบที่เกิดขึ้นกับกรณีของ lsode นั้นเกิดจากการที่จำนวนจุดคำนวณนั้นต่ำเกินไป (ใช้เพียงแค่ 200 จุด) แต่ที่ไม่พบในกรณีของ ode45 นั้นเป็นเพราะว่า ode45 ใช้การปรับระยะ step size การคำนวณโดยอัตโนมัติเพื่อให้สอดรับกับอัตราการเปลี่ยนแปลงของคำตอบ แต่สิ่งที่เห็นคือระยะห่างของจุดบริเวณนั้นแคบมาก (อยู่ที่ระดับ 10^-7-10^-6) จำนวนจุดการคำนวณตลอดช่วงอยู่ที่เกือบ 35,000 จุด ทำให้การคำนวณใช้เวลานานมากอย่างเห็นได้ชัด

  

แต่พอเปลี่ยนมาใช้ ode15s ปรากฏว่าคำตอบที่ได้นั้นเปลี่ยนไปเป็นคนละเรื่องเลย คือคำสั่งนี้ไม่สามารถพบตำแหน่งที่ปฏิกิริยาเกิดการ runaway ได้ (ตำแหน่งที่อุณหภูมิและความเข้มข้นเปลี่ยนแปลงอย่างรวดเร็ว) คำตอบที่ได้นั้นเสมือนกับว่าปฏิกิริยาแทบจะไม่เกิดแม้ว่าจะใช้อุณหภูมิ Tc สูงกว่า 415.0ºC ก็ตาม

   

รูปที่ ๑ ชุดคำสั่งที่ใช้ในการคำนวณ ข้างบนเป็นของ lsode ส่วนข้างล่างเป็นของ ode45s

ความแตกต่างของตัวอย่างนี้กับตัวอย่างในตอนที่ ๑ และ ๒ อยู่ตรงที่ในตอนที่ ๑ และ ๒ นั้นการเปลี่ยนแปลงอย่างรวดเร็วหรืออย่างกระทันหันเกิดขึ้นจากจุดเริ่มต้นการคำนวณ แต่ในกรณีนี้เกิดขึ้นระหว่างเส้นการคำนวณไปยังจุดปลายทาง ในกรณีของวิธี multi-step ที่มีการใช้จุดข้อมูลย้อนหลังในการคำนวณ ประกอบกับการเลือกระยะ step size อัตโนมัติ จึงเป็นไปได้ว่าการเปลี่ยนแปลงอย่างช้า ๆ ที่เกิดขึ้นอย่างต่อเนื่องมาก่อนหน้า ทำให้ไม่มีการคาดการณ์ว่าจะมีการเปลี่ยนแปลงอย่างกระทันหันเกิดขึ้นต่อหน้า ทำให้การคำนวณเดินข้ามตำแหน่งดังกล่าวไป 

   

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

   

รูปที่ ๒ กราฟการเปลี่ยนแปลองอุณหภูมิภายใน reactor ที่ได้จากคำสั่ง lsode ที่ค่า Tc = 414.0, 414.7 และ 414.8ºC โดยแบ่งช่วงการคำนวณออกเป็น 200 ช่วงย่อย

   

รูปที่ ๓ กราฟการเปลี่ยนแปลงความเข้มข้นของสารตั้งต้นภายใน reactor ที่ได้จากคำสั่ง lsode ที่ค่า Tc = 414.0, 414.7 และ 414.8ºC โดยแบ่งช่วงการคำนวณออกเป็น 200 ช่วงย่อย

   

รูปที่ ๔ กราฟการเปลี่ยนแปลองอุณหภูมิภายใน reactor ที่ได้จากคำสั่ง ode45 ที่ค่า Tc = 414.0, 414.7 และ 414.8ºC

   

รูปที่ ๕ กราฟการเปลี่ยนแปลงความเข้มข้นของสารตั้งต้นภายใน reactor ที่ได้จากคำสั่ง ode45 ที่ค่า Tc = 414.0, 414.7 และ 414.8ºC

   

รูปที่ ๖ กราฟการเปลี่ยนแปลงอุณหภูมิและความเข้มข้นของสารตั้งต้นภายใน reactor ที่ได้จากคำสั่ง ode15s ที่ค่า Tc = 415.0ºC