Heat Modeling
Heat Modeling
This is a special issue published in “Mathematical Problems in Engineering.” All articles are open access articles distributed under the
Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the
original work is properly cited.
Editorial Board
Mohamed Abd El Aziz, Egypt Abdel-Ouahab Boudraa, France Horst Ecker, Austria
Farid Abed-Meraim, France Taha Boukhobza, France Karen Egiazarian, Finland
Silvia Abrahão, Spain Francesco Braghin, Italy Ahmed El Hajjaji, France
Paolo Addesso, Italy Michael J. Brennan, UK Mohsen Elhafsi, USA
Claudia Adduce, Italy Gunther Brenner, Germany Fouad Erchiqui, Canada
Ramesh Agarwal, USA Maurizio Brocchini, Italy Anders Eriksson, Sweden
Juan C. Agüero, Australia Julien Bruchon, France Giovanni Falsone, Italy
Ricardo Aguilar-Lpez, Mexico Javier Buldufh, Spain Hua Fan, China
Tarek Ahmed-Ali, France Tito Busani, USA Yann Favennec, France
Hamid Akbarzadeh, Canada Pierfrancesco Cacciola, UK Roberto Fedele, Italy
Muhammad N. Akram, Norway Salvatore Caddemi, Italy Giuseppe Fedele, Italy
Salvatore Alfonzetti, Italy Jose E. Capilla, Spain Jacques Ferland, Canada
Francisco Alhama, Spain Ana Carpio, Spain Jose R. Fernandez, Spain
Tofigh Allahviranloo, Iran Miguel E. Cerrolaza, Spain Simme Douwe Flapper, Netherlands
Juan A. Almendral, Spain Mohammed Chadli, France Thierry Floquet, France
Saiied Aminossadati, Australia Gregory Chagnon, France Eric Florentin, France
Lionel Amodeo, France Ching-Ter Chang, Taiwan Francesco Franco, Italy
Igor Andrianov, Germany Michael J. Chappell, UK Tomonari Furukawa, USA
Sebastian Anita, Romania Kacem Chehdi, France Mohamed Gadala, Canada
Renata Archetti, Italy Xinkai Chen, Japan Matteo Gaeta, Italy
Felice Arena, Italy Chunlin Chen, China Zoran Gajic, USA
Sabri Arik, Turkey Francisco Chicano, Spain Ciprian G. Gal, USA
Fumihiro Ashida, Japan Hung-Yuan Chung, Taiwan Rafael Gallego, Spain
Hassan Askari, Canada Joaquim Ciurana, Spain Ugo Galvanetto, Italy
Mohsen Asle Zaeem, USA John D. Clayton, USA Akemi Glvez, Spain
Francesco Aymerich, Italy Carlo Cosentino, Italy Rita Gamberini, Italy
Seungik Baek, USA Paolo Crippa, Italy Maria Gandarias, Spain
Khaled Bahlali, France Erik Cuevas, Mexico Arman Ganji, Canada
Laurent Bako, France Peter Dabnichki, Australia Zhong-Ke Gao, China
Stefan Balint, Romania Luca D’Acierno, Italy Xin-Lin Gao, USA
Alfonso Banos, Spain Weizhong Dai, USA Giovanni Garcea, Italy
Roberto Baratti, Italy Purushothaman Damodaran, USA Fernando Garca, Spain
Martino Bardi, Italy Farhang Daneshmand, Canada Laura Gardini, Italy
Azeddine Beghdadi, France Fabio De Angelis, Italy Alessandro Gasparetto, Italy
Tarak Ben Zineb, France Stefano de Miranda, Italy Vincenzo Gattulli, Italy
Abdel-Hakim Bendada, Canada Filippo de Monte, Italy Jürgen Geiser, Germany
Ivano Benedetti, Italy Xavier Delorme, France Oleg V. Gendelman, Israel
Elena Benvenuti, Italy Luca Deseri, USA Mergen H. Ghayesh, Australia
Jamal Berakdar, Germany Yannis Dimakopoulos, Greece Anna M. Gil-Lafuente, Spain
Enrique Berjano, Spain Zhengtao Ding, UK Hector Gómez, Spain
Jean-Charles Beugnot, France Ralph B. Dinwiddie, USA Rama S. R. Gorla, USA
Simone Bianco, Italy Mohamed Djemai, France Oded Gottlieb, Israel
David Bigaud, France Alexandre B. Dolgui, France Antoine Grall, France
Jonathan N. Blakely, USA George S. Dulikravich, USA Jason Gu, Canada
Daniela Boso, Italy Bogdan Dumitrescu, Finland Quang Phuc Ha, Australia
Ofer Hadar, Israel Yaguo Lei, China Ben T. Nohara, Japan
Masoud Hajarian, Iran Thibault Lemaire, France Mohammed Nouari, France
Frédéric Hamelin, France Stefano Lenci, Italy Mustapha Nourelfath, Canada
Zhen-Lai Han, China Roman Lewandowski, Poland Sotiris K. Ntouyas, Greece
Thomas Hanne, Switzerland Qing Q. Liang, Australia Roger Ohayon, France
Xiao-Qiao He, China Panos Liatsis, UK Mitsuhiro Okayasu, Japan
Marı́a I. Herreros, Spain Wanquan Liu, Australia Javier Ortega-Garcia, Spain
Vincent Hilaire, France Yan-Jun Liu, China Alejandro Ortega-Moux, Spain
Eckhard Hitzer, Japan Peide Liu, China Naohisa Otsuka, Japan
Jaromir Horacek, Czech Republic Peter Liu, Taiwan Erika Ottaviano, Italy
Muneo Hori, Japan Jean J. Loiseau, France Alkiviadis Paipetis, Greece
Andrs Horvth, Italy Paolo Lonetti, Italy Alessandro Palmeri, UK
Gordon Huang, Canada Luis M. López-Ochoa, Spain Anna Pandolfi, Italy
Sajid Hussain, Canada Vassilios C. Loukopoulos, Greece Elena Panteley, France
Asier Ibeas, Spain Valentin Lychagin, Norway Manuel Pastor, Spain
Giacomo Innocenti, Italy Fazal M. Mahomed, South Africa Pubudu N. Pathirana, Australia
Emilio Insfran, Spain Yassir T. Makkawi, UK Francesco Pellicano, Italy
Nazrul Islam, USA Noureddine Manamanni, France Haipeng Peng, China
Payman Jalali, Finland Didier Maquin, France Mingshu Peng, China
Reza Jazar, Australia Paolo Maria Mariano, Italy Zhike Peng, China
Khalide Jbilou, France Benoit Marx, France Marzio Pennisi, Italy
Linni Jian, China Gefhrard A. Maugin, France Matjaz Perc, Slovenia
Bin Jiang, China Driss Mehdi, France Claudio Pernechele, Italy
Zhongping Jiang, USA Roderick Melnik, Canada Francesco Pesavento, Italy
Ningde Jin, China Pasquale Memmolo, Italy Maria do Rosrio Pinho, Portugal
Grand R. Joldes, Australia Xiangyu Meng, Canada Antonina Pirrotta, Italy
Joaquim Joao Judice, Portugal Jose Merodio, Spain Vicent Pla, Spain
Tadeusz Kaczorek, Poland Luciano Mescia, Italy Javier Plaza, Spain
Tamas Kalmar-Nagy, Hungary Laurent Mevel, France Jean-Christophe Ponsart, France
Tomasz Kapitaniak, Poland Philippe Micheau, Canada Mauro Pontani, Italy
Haranath Kar, India Y. Vladimirovich Mikhlin, Ukraine Stanislav Potapenko, Canada
Konstantinos Karamanos, Belgium Aki Mikkola, Finland Sergio Preidikman, USA
C. Masood Khalique, South Africa Hiroyuki Mino, Japan Christopher Pretty, New Zealand
Nam-Il Kim, Korea Pablo Mira, Spain Carsten Proppe, Germany
Do Wan Kim, Korea Vito Mocella, Italy Luca Pugi, Italy
Oleg Kirillov, Germany Roberto Montanini, Italy Yuming Qin, China
Alexander Klimenko, Australia Gisele Mophou, France Dane Quinn, USA
Manfred Krafczyk, Germany Rafael Morales, Spain Jose Ragot, France
Frederic Kratz, France Aziz Moukrim, France K. Ramamani Rajagopal, USA
Jurgen Kurths, Germany Emiliano Mucchi, Italy Gianluca Ranzi, Australia
Kyandoghere Kyamakya, Austria Domenico Mundo, Italy Sivaguru Ravindran, USA
Davide La Torre, Italy Jose J. Muoz, Spain Alessandro Reali, Italy
Risto Lahdelma, Finland Giuseppe Muscolino, Italy Giuseppe Rega, Italy
Hak-Keung Lam, UK Marco Mussetta, Italy Oscar Reinoso, Spain
Antonino Laudani, Italy Hakim Naceur, France Nidhal Rezg, France
Aime’ Lay-Ekuakille, Italy Hassane Naji, France Ricardo Riaza, Spain
Mark Leeson, UK Dong Ngoduy, UK Gerasimos Rigatos, Greece
Marek Lefik, Poland Tatsushi Nishi, Japan José Rodellar, Spain
Rosana Rodriguez-Lopez, Spain Luciano Simoni, Italy Raoul van Loon, UK
Ignacio Rojas, Spain Christos H. Skiadas, Greece Alain Vande Wouwer, Belgium
Carla Roque, Portugal Michael Small, Australia Pandian Vasant, Malaysia
Aline Roumy, France Francesco Soldovieri, Italy M. E. Vázquez-Méndez, Spain
Debasish Roy, India Raffaele Solimene, Italy Josep Vehi, Spain
R. Ruiz Garcı́a, Spain Ruben Specogna, Italy Kalyana C. Veluvolu, Korea
Antonio Ruiz-Cortes, Spain Victor Sreeram, Australia Fons J. Verbeek, Netherlands
Ivan D. Rukhlenko, Australia Sri Sridharan, USA Franck J. Vernerey, USA
Mazen Saad, France Ivanka Stamova, USA Georgios Veronis, USA
Kishin Sadarangani, Spain Rolf Stenberg, Finland Anna Vila, Spain
Mehrdad Saif, Canada Yakov Strelniker, Israel Rafael J. Villanueva, Spain
Miguel A. Salido, Spain Sergey A. Suslov, Australia U. E. Vincent, UK
Roque J. Saltarén, Spain Thomas Svensson, Sweden Mirko Viroli, Italy
Alessandro Salvini, Italy Andrzej Swierniak, Poland Michael Vynnycky, Sweden
Angel Snchez, Spain Yang Tang, Germany Junwu Wang, China
Maura Sandri, Italy Sergio Teggi, Italy Shuming Wang, Singapore
Miguel A. F. Sanjuan, Spain Roger Temam, USA Yan-Wu Wang, China
Juan F. San-Juan, Spain Alexander Timokha, Norway Yongqi Wang, Germany
Roberta Santoro, Italy Rafael Toledo-Moreo, Spain Jeroen A. S. Witteveen, Netherlands
Ilmar Ferreira Santos, Denmark Gisella Tomasini, Italy Yuqiang Wu, China
José A. Sanz-Herrera, Spain Francesco Tornabene, Italy Dash Desheng Wu, Canada
Nickolas S. Sapidis, Greece Antonio Tornambe, Italy Xuejun Xie, China
E. J. Sapountzakis, Greece Fernando Torres, Spain Guangming Xie, China
Themistoklis P. Sapsis, USA Fabio Tramontana, Italy Gen Qi Xu, China
Andrey V. Savkin, Australia Sébastien Tremblay, Canada Hang Xu, China
Valery Sbitnev, Russia Irina N. Trendafilova, UK Xinggang Yan, UK
Thomas Schuster, Germany George Tsiatas, Greece Luis J. Yebra, Spain
Mohammed Seaid, UK Antonios Tsourdos, UK Peng-Yeng Yin, Taiwan
Lotfi Senhadji, France Vladimir Turetsky, Israel Ibrahim Zeid, USA
Joan Serra-Sagrista, Spain Mustafa Tutar, Spain Qingling Zhang, China
Leonid Shaikhet, Ukraine Efstratios Tzirtzilakis, Greece Huaguang Zhang, China
Hassan M. Shanechi, USA Filippo Ubertini, Italy Jian Guo Zhou, UK
Sanjay K. Sharma, India Francesco Ubertini, Italy Quanxin Zhu, China
Bo Shen, Germany Hassan Ugail, UK Mustapha Zidi, France
Babak Shotorban, USA Giuseppe Vairo, Italy Alessandro Zona, Italy
Zhan Shu, UK Kuppalapalle Vajravelu, USA
Dan Simon, USA Robertt A. Valente, Portugal
Contents
Mathematical Modeling of Heat and Mass Transfer in Energy Science and Engineering 2014,
Zhijun Zhang, Hua-Shu Dou, Ireneusz Zbicinski, Zhonghua Wu, and Jun Liu
Volume 2015, Article ID 609382, 3 pages
Numerical Investigations of the Effect of Nonlinear Quadratic Pressure Gradient Term on a Moving
Boundary Problem of Radial Flow in Low-Permeable Reservoirs with Threshold Pressure Gradient,
Wenchao Liu and Jun Yao
Volume 2015, Article ID 275057, 12 pages
A Numerical Study of Natural Convection Heat Transfer in Fin Ribbed Radiator, Hua-Shu Dou,
Gang Jiang, and Lite Zhang
Volume 2015, Article ID 989260, 13 pages
Numerical Study of Buoyancy Convection of Air under Permanent Magnetic Field and Comparison
with That under Gravity Field, Kewei Song, Wenkai Li, Yang Zhou, and Yuanru Lu
Volume 2014, Article ID 494585, 13 pages
System Model of Heat and Mass Transfer Process for Mobile Solvent Vapor Phase Drying Equipment,
Shiwei Zhang, Yufang Zhu, Baozhen Qiao, and Zhijun Zhang
Volume 2014, Article ID 267276, 11 pages
Research on the Impact of Wind Angles on the Residential Building Energy Consumption, Zou Huifen,
Yang Fuhua, and Zhang Qian
Volume 2014, Article ID 794650, 15 pages
Local Fractional Laplace Variational Iteration Method for Nonhomogeneous Heat Equations Arising in
Fractal Heat Flow, Shu Xu, Xiang Ling, Carlo Cattani, Gong-Nan Xie, Xiao-Jun Yang, and Yang Zhao
Volume 2014, Article ID 914725, 7 pages
Numerical Simulation and Experimental Research on Coal Ash Collecting and Grading System,
Yuanhua Xie, Xianjin Li, Liwei Wang, Hong Yu, Bing Bai, Zhizhou Xu, and Tong Zhu
Volume 2014, Article ID 373967, 12 pages
Simulation of Microstructure during Laser Rapid Forming Solidification Based on Cellular Automaton,
Zhi-jian Wang, Shuai Luo, Hong-wu Song, Wei-dong Deng, and Wen-yi Li
Volume 2014, Article ID 627528, 9 pages
Application of CFD, Taguchi Method, and ANOVA Technique to Optimize Combustion and Emissions
in a Light Duty Diesel Engine, Senlin Xiao, Wanchen Sun, Jiakun Du, and Guoliang Li
Volume 2014, Article ID 502902, 9 pages
Mathematical Simulation of Heat and Mass Transfer Processes at the Ignition of Liquid Fuel by
Concentrated Flux of Radiation, Olga V. Vysokomornaya, Genii V. Kuznetsov, and Pavel A. Strizhak
Volume 2014, Article ID 156150, 7 pages
Heat and Mass Transfer of Droplet Vacuum Freezing Process Based on Dynamic Mesh, Lili Zhao,
Yuekai Zhang, Zhijun Zhang, Xun Li, and Wenhui Zhang
Volume 2014, Article ID 798040, 6 pages
Numerical and Experimental Research of Heat and Mass Transfer at the Heterogeneous System Ignition
by Local Energy Source with Limited Heat Content, Dmitrii O. Glushkov, Genii V. Kuznetsov,
and Pavel A. Strizhak
Volume 2014, Article ID 281527, 9 pages
An Improved Dispatch Strategy of a Grid-Connected Hybrid Energy System with High Penetration
Level of Renewable Energy, Yan Zhang, Jie Meng, Bo Guo, and Tao Zhang
Volume 2014, Article ID 602063, 18 pages
Analysis of Combined Power and Refrigeration Generation Using the Carbon Dioxide Thermodynamic
Cycle to Recover the Waste Heat of an Internal Combustion Engine, Shunsen Wang, Kunlun Bai,
Yonghui Xie, Juan Di, and Shangfang Cheng
Volume 2014, Article ID 689398, 12 pages
Mathematical Modeling of Eddy-Current Loss for a New Induction Heating Device, Hai Du, Junyuan Li,
and Yanbin Qu
Volume 2014, Article ID 923745, 7 pages
Energy Loss in Pulse Detonation Engine due to Fuel Viscosity, Weipeng Hu, Zichen Deng,
and Gongnan Xie
Volume 2014, Article ID 735926, 5 pages
Pore Network Analysis of Zone Model for Porous Media Drying, Yuan Yuejin, Zhao Zhe, Nie Junnan,
and Xu Yingying
Volume 2014, Article ID 624145, 8 pages
Model-Based Water Wall Fault Detection and Diagnosis of FBC Boiler Using Strong Tracking Filter,
Li Sun, Junyi Dong, Donghai Li, and Yuqiong Zhang
Volume 2014, Article ID 504086, 8 pages
Application of D-S Evidence Fusion Method in the Fault Detection of Temperature Sensor, Zheng Dou,
Xiaochun Xu, Yun Lin, and Ruolin Zhou
Volume 2014, Article ID 395057, 6 pages
Mixed Convection Unsteady Stagnation-Point Flow towards a Stretching Sheet with Slip Effects,
Hui Chen
Volume 2014, Article ID 435697, 7 pages
Effects of Wall Shear Stress on MHD Conjugate Flow over an Inclined Plate in a Porous Medium with
Ramped Wall Temperature, Arshad Khan, Ilyas Khan, Farhad Ali, and Sharidan Shafie
Volume 2014, Article ID 861708, 15 pages
An Optimization Model Based on Electric Power Generation in Steel Industry, Jing-yu Liu and Jiu-ju Cai
Volume 2014, Article ID 924960, 10 pages
Identification of Shaft Centerline Orbit for Wind Power Units Based on Hopfield Neural Network
Improved by Simulated Annealing, Kun Ren and Jihong Qu
Volume 2014, Article ID 571354, 6 pages
Hindawi Publishing Corporation
Mathematical Problems in Engineering
Volume 2015, Article ID 609382, 3 pages
[Link]
Editorial
Mathematical Modeling of Heat and Mass Transfer in Energy
Science and Engineering 2014
Zhijun Zhang,1 Hua-Shu Dou,2 Ireneusz Zbicinski,3 Zhonghua Wu,4 and Jun Liu5
1
School of Mechanical Engineering and Automation, Northeastern University, Shenyang 110004, China
2
Faculty of Mechanical Engineering and Automation, Zhejiang Sci-Tech University, Hangzhou, Zhejiang 310018, China
3
Faculty of Process and Environmental Engineering, Technical University of Lodz, 90-924 Lodz, Poland
4
Institute of Drying and Dewatering, College of Mechanical Engineering, Tianjin University of Science and Technology,
1038 Daguan Road, Hexi District, Tianjin 300222, China
5
Institute of Biology and Chemistry, Shenyang University, Shenyang 110044, China
Copyright © 2015 Zhijun Zhang et al. This is an open access article distributed under the Creative Commons Attribution License,
which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Heat and mass transfer process is the basis of energy research pressure gradient” by W. Liu and J. Yao, based on these
from the boiler, gas turbine to fuel cell, solar power. It is concerns, in consideration of the QPGT, a moving boundary
the need of energy conversion and transfer process. The model of radial flow in low-permeable reservoirs with the
heat and mass transfer process mainly include the following TPG for the case of a constant flow rate at the inner boundary
aspects: thermodynamics, power and fluid machinery, heat is constructed.
and mass transfer, combustion, and multiphase flow. The idea In “Numerical study of buoyancy convection of air under
of this special issue is to consider the study and applications permanent magnetic field and xomparison with that under
of mathematical modeling method on energy science and gravity field” by K. Song et al., magneto thermal free con-
technology. vection of air in a square enclosure under a nonuniform
This special issue contains 27 papers, the contents of magnetic field provided by a permanent neodymium-iron-
which are summarized as follows. boron magnet is numerically studied.
The paper entitled “A numerical study of natural con- In “System model of heat and mass transfer process for
vection heat transfer in fin ribbed radiator” by H.-S. Dou mobile solvent vapor phase drying equipment” by S. Zhang
et al. numerically investigates the physical mechanism of flow et al., on the basis of necessary simplification and assumption
instability and heat transfer of natural convection in a cavity for MVPD equipment and process, a heat and mass transfer
fixed with fin arrays. mathematical model including 40 mathematical equations
In “A differential-algebraic model for the once-through is established, which represents completely thermodynamics
steam generator of MHTGR-based multi-modular nuclear laws of phase change and transport process of solvent, water,
plants” by Z. Dong, based upon the conservation laws of mass, and air in MVPD technological processes and describes in
energy, and momentum, a new differential-algebraic model detail the quantitative relationship among important physical
for the OTSGs of the MHTGR based multimodular nuclear quantities such as temperature, pressure, and flux in key
plants is given. equipment units and process.
In “Numerical investigations of the effect of nonlinear The paper entitled “Research on the impact of wind angles
quadratic pressure gradient term on a moving boundary prob- on the residential building energy consumption” by Z. Huifen
lem of radial flow in low-permeable reservoirs with threshold et al., combined with natural ventilation under various wind
2 Mathematical Problems in Engineering
angles, gives the best recommended solution of building numerical and experimental investigations were executed
direction in Shenyang. for determination of macroscopic regularities of heat and
In “Local fractional Laplace variational iteration method mass transfer processes under the conditions of the phase
for nonhomogeneous heat equations arising in fractal heat transformations and chemical reaction at the ignition of
flow” by S. Xu et al., the approximate solutions are nondiffer- vapors coming from fabrics impregnated by combustible
entiable functions and their plots are also given to show the liquid into oxidant area at the local power supply.
accuracy and efficiency to implement the previous method. The paper entitled “An improved dispatch strategy of a
In “Numerical simulation and experimental research on grid-connected hybrid energy system with high penetration
coal ash collecting and grading system” by Y. Xie at al., focused level of renewable energy” by Y. Zhang et al. focuses on the
on single coal ash particle, Matlab software was used to fluctuation alleviation and power quality improvement of
simulate the force conditions and separation parameters of grid-connected HES with high penetration level of RES.
various diameter coal ash particles in airflow. Fluent software In “Analysis of combined power and refrigeration gen-
was used to simulate the nozzle fluidization domain shape eration using the carbon dioxide thermodynamic cycle to
and to determine optimal jet flux. recover the waste heat of an internal combustion engine” by
In “Simulation of microstructure during laser rapid form- S. Wang et al., a novel thermodynamic system is proposed
ing solidification based on cellular automaton” by Z. Wang to recover the waste heat of an internal combustion engine
et al., the grain microstructure of molten pool during the (ICE) by integrating the transcritical carbon dioxide (CO2 )
solidification of TC4 titanium alloy in the single point laser refrigeration cycle with the supercritical CO2 power cycle,
cladding was investigated based on the CAFE model which is and eight kinds of integration schemes are developed.
the cellular automaton (CA) coupled with the finite element In “Mathematical modeling of eddy-current loss for a new
(FE) method. induction heating device” by H. Du et al., this device can
In “Application of CFD, Taguchi method, and ANOVA convert mechanical energy into heat energy by utilizing eddy
technique to optimize combustion and emissions in a light currents, which are induced by rotating permanent magnets.
duty diesel engine” by S. Xiao et al., in order to understand A mathematical model is established for estimating eddy-
the combined effect of EGR rate, pilot fuel quantity, and current loss of the device.
main injection timing on the NO𝑥 (oxides of nitrogen), In “Energy loss in pulse detonation engine due to fuel
soot, and ISFC (indicated specific fuel consumption), CFD viscosity” by W. Hu et al., to analyze the energy loss in
(computational fluid dynamics) simulation together with the pulse detonation engine (PDE) due to the viscosity of
the Taguchi method and the ANOVA (analysis of variance) the fuel, the energy loss in the Burgers model excited by
technique was applied as an effective research tool. periodic impulses is investigated based on the generalized
In “Mathematical simulation of heat and mass transfer multisymplectic method in this paper.
processes at the ignition of liquid fuel by concentrated flux of The paper entitled “Pore network analysis of zone model
radiation” by O. V. Vysokomornaya et al., the physical and for porous media drying” by Y. Yuejin et al. fused the physical
forecasting mathematical models of heat and mass transfer parameters of porous media, such as porosity, pore mean
with phase transformations and chemical reactions under diameter, and pore size distribution into the model param-
heating and following ignition of typical liquid fuel by using eters, and a sand bed drying experiment was conducted to
concentrated flow of radiation were developed. verify the validity of this model.
In “Heat and mass transfer of droplet vacuum freezing In “Model-based water wall fault detection and diagnosis
process based on dynamic mesh” by L. Zhao et al., the initial of FBC boiler using strong tracking filter” by L. Sun et al., a
droplet diameter, initial droplet temperature, and vacuum model-based approach is presented to estimate internal states
chamber pressure effect are studied. and heat transfer coefficient dually from the noisy measurable
In “Entropy generation analysis of power-law non- outputs.
Newtonian fluid flow caused by micropatterned moving In “Application of D-S evidence fusion method in the fault
surface” by M. H. Yazdi et al., the first and second law analyses detection of temperature sensor” by Z. Dou et al., based on
of power-law non-Newtonian flow over embedded open the idea of information fusion and the requirements of D-S
parallel microchannels within micropatterned permeable evidence method, a D-S evidence fusion structure with two
continuous moving surface are examined at prescribed layers was introduced to detect the temperature sensor fault
surface temperature. in drying process.
In “Numerical investigation on a prototype centrifugal The paper entitled “Mixed convection unsteady stagna-
pump subjected to fluctuating rotational speed” by Y.-L. Zhang tion-point flow towards a stretching sheet with slip effects” by
et al., in order to study the transient response characteristic H. Chen studies the unsteady mixed convection flow of an
of a prototype centrifugal pump subjected to fluctuating incompressible viscous fluid about a stagnation point on a
rotational speed, a closed-loop pipe system including the stretching sheet in presence of velocity and thermal slips.
pump is built to accomplish unsteady flow calculations in The paper entitled “Effects of wall shear stress on MHD
which the boundary conditions at the inlet and the outlet of conjugate flow over an inclined plate in a porous medium
the pump are not required to be set. with ramped wall Temperature” by A. Khan et al. investigates
In “Numerical and experimental research of heat and mass the effects of an arbitrary wall shear stress on unsteady
transfer at the heterogeneous system ignition by local energy magnetohydrodynamic (MHD) flow of a Newtonian fluid
source with limited heat content” by D. O. Glushkov et al., with conjugate effects of heat and mass transfer.
Mathematical Problems in Engineering 3
Research Article
Numerical Investigations of the Effect of Nonlinear Quadratic
Pressure Gradient Term on a Moving Boundary Problem of
Radial Flow in Low-Permeable Reservoirs with Threshold
Pressure Gradient
Copyright © 2015 W. Liu and J. Yao. This is an open access article distributed under the Creative Commons Attribution License,
which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
The existence of a TPG can generate a relatively high pressure gradient in the process of fluid flow in porous media in low-permeable
reservoirs, and neglecting the QPGTs in the governing equations, by assuming a small pressure gradient for such a problem, can
cause a significant error in predicting the formation pressure. Based on these concerns, in consideration of the QPGT, a moving
boundary model of radial flow in low-permeable reservoirs with the TPG for the case of a constant flow rate at the inner boundary
is constructed. Due to strong nonlinearity of the mathematical model, a numerical method is presented: the system of partial
differential equations for the moving boundary problem is first transformed equivalently into a closed system of partial differential
equations with fixed boundary conditions by a spatial coordinate transformation method; and then a stable, fully implicit finite
difference method is used to obtain its numerical solution. Numerical result analysis shows that the mathematical models of radial
flow in low-permeable reservoirs with TPG must take the QPGT into account in their governing equations, which is more important
than those of Darcy’s flow; the sensitive effects of the QPGT for the radial flow model do not change with an increase of the
dimensionless TPG.
pressure profiles due to a relatively high formation pressure where 𝑘 is the permeability of the porous medium, 10−3 ⋅ 𝜇m2 ;
gradient, generated in fluid flow in porous media with TPG. 𝜇 is the fluid viscosity, mPa⋅s; r is the radial distance, m; 𝜐 is
What is more, in modern times, with the development of the seepage velocity, m⋅s−1 ; and 𝜆 is the TPG, MPa⋅m−1 .
advanced analysis methods and improved resolutions of The continuous equation for the radial flow in the porous
pressure measurement machines [29], it is necessary to study medium is as follows [15]:
quantitatively the effect of QPGT on these moving boundary
problems. 1 𝜕 𝜕 (𝜌𝜙)
Hence, based on these concerns, in this paper, QPGT is − (𝑟𝜌𝜐) = , 0 ≤ 𝑟 ≤ 𝑠 (𝑡) , (4)
𝑟 𝜕𝑟 𝜕𝑡
incorporated in the governing equation for a basic moving
boundary problem of radial flow in an infinite reservoir where 𝑡 is time, h and 𝑠 is the moving boundary, m.
with TPG for the case of a constant flow rate at the inner Substituting (1)–(3) into (4), the governing (mass bal-
boundary. Owing to the existence of the moving boundary ance) equation for the radial flow in the low-permeable
and the nonlinearity of partial differential equations, it is reservoir, considering the nonlinear QPGT, can be deduced
really difficult to obtain its exact analytical solution. Here, we as follows (see Appendix):
adopted a spatial coordinate transform based finite difference
method to solve the nonlinear moving boundary model for 𝜕2 𝑝 1 𝜕𝑝 𝜆 𝜕𝑝 2 𝜇𝜙𝑖 𝐶𝑡 𝜕𝑝
the radial flow case. The numerical method has been strictly + ⋅ − + 𝐶𝑓 ⋅ ( ) = ⋅ , (5)
verified for solving such moving boundary problems with 𝜕𝑟2 𝑟 𝜕𝑟 𝑟 𝜕𝑟 𝑘 𝜕𝑡
good accuracy by a one-dimensional flow case in the recently
published literature [21]. Furthermore, the effect of this where 𝐶𝑡 is the total compression coefficient, MPa−1 .
QPGT on the numerical solutions of the moving boundary The initial condition is as follows:
problem can be discussed and analyzed quantitatively, by
𝑠 (0) = 0,
using the numerical results with respect to different values of
(6)
a dimensionless TPG. 𝑝𝑡=0 = 𝑝𝑖 .
2. Mathematical Modeling by The inner boundary conditions with constant flow rate
Considering QPGT are
2𝜋𝑘ℎ𝑟 𝜕𝑝
The problem considered involves radial flow in an infinite
( − 𝜆) = 𝑞 ⋅ 𝐵, (7)
reservoir with TPG for the case of a constant flow rate at 𝜇 𝜕𝑟 𝑟=𝑟
𝑤
the inner boundary; the porous medium is homogeneous,
isotropic, and isothermal; the single-phase horizontal flow
where 𝑞 is the constant flow rate, m3 ⋅d−1 ; h is the reservoir
does not have any gravity effect; for a basic problem, skin
thickness, m; 𝑟𝑤 is the wellbore radius, m; and B is the volume
effect and wellbore storage are not considered here; finally,
factor, dimensionless.
the fluid and porous medium are slightly compressible.
The moving boundary conditions are:
The fluid density is as follows [20]:
𝜌 = 𝜌𝑖 exp (−𝐶𝑓 (𝑝𝑖 − 𝑝)) , (1) 𝑝𝑟=𝑠(𝑡) = 𝑝𝑖 ,
𝜕𝑝 (8)
where 𝜌 is the fluid density, kg⋅m−3 ; 𝜌𝑖 is the initial fluid = 𝜆.
𝜕𝑟 𝑟=𝑠(𝑡)
density, kg⋅m−3 ; 𝑝𝑖 is the initial pressure, MPa; 𝑝 is the
pressure, MPa; and 𝐶𝑓 is the compression coefficient of the
The moving boundary conditions (8) physically mean
fluid, MPa−1 . that the seepage flow only happens at the area near the
The porosity of the porous medium is as follows [20]: wellbore inside the moving boundary, where the formation
pressure gradient is larger than TPG; however, outside the
𝜙 = 𝜙𝑖 exp (−𝐶𝜙 (𝑝𝑖 − 𝑝)) , (2) moving boundary, the formation pressure gradient is smaller
than TPG (equal to zero), so there is no seepage flow behavior,
where 𝜙 is the porosity of the porous medium, fraction; 𝜙𝑖
and the formation pressure also keeps the initial pressure;
is the initial porosity, fraction; and 𝐶𝜙 is the compression
in addition, on the moving boundary, the pressure gradient
coefficient of the porosity, MPa−1 . is just equal to TPG; with the time increasing, the moving
The modified Darcy’s law for the fluid flow in the porous boundary also moves outside gradually. The existence of
medium with TPG is as follows [8]: moving boundary is the main difference between the models
d𝑝 of fluid flow in porous media with TPG and the classical
{
{ 0, 0≤ ≤ 𝜆, Darcy’s flow models.
{
{ d𝑟
{ Equations (5)–(8) together form the moving boundary
𝜐={ (3)
{
{ problem of radial flow in the infinite reservoir with TPG for
{− 𝑘 ⋅ ( d𝑝 − 𝜆) ,
{ d𝑝
> 𝜆, the case of a constant flow rate at the inner boundary, in
{ 𝜇 d𝑟 d𝑟 consideration of QPGT.
4 Mathematical Problems in Engineering
Define the following dimensionless variables: Letting 𝑥𝐷 = 𝛿(𝑡𝐷) on both sides of (10) and then
substituting (14) yield
𝑟
𝑟𝐷 = ,
𝑟𝑤
𝜕𝑃𝐷 𝜕2 𝑃𝐷
= 2
− 𝛼𝐷 ⋅ 𝜆2𝐷. (19)
3.6 × 10−3 𝑘𝑡 𝜕𝑡𝐷 𝑟𝐷 =𝛿(𝑡𝐷 ) 𝜕𝑟𝐷 𝑟𝐷 =𝛿(𝑡𝐷 )
𝑡𝐷 = ,
𝜇𝜙𝑖 𝐶𝑡 𝑟𝑤2
𝑠 By (18) and (19), the velocity of the moving boundary can
𝛿= , be written as follows:
𝑟𝑤
(9)
𝑘ℎ [𝑝𝑖 − 𝑝]
𝑃𝐷 = , 𝜕𝛿 1 𝜕2 𝑃𝐷
1.842𝑞𝜇𝐵 = 2
− 𝛼𝐷 ⋅ 𝜆 𝐷. (20)
𝜕𝑡𝐷 𝜆 𝐷 𝜕𝑟𝐷 𝑟 𝐷 =𝛿(𝑡𝐷 )
𝑘ℎ𝑟𝑤 𝜆
𝜆𝐷 = ,
1.842𝑞𝜇𝐵 From (20), it can be seen that the existence of the QPGT
(the dimensionless compressibility 𝛼𝐷 is not equal to zero)
1.842𝑞𝜇𝐵𝐶𝑓 can slow down the velocity of the moving boundary of the
𝛼𝐷 = ,
𝑘ℎ radial flow in low-permeable reservoir with TPG.
where 𝑟𝐷 is the dimensionless radial distance, 𝑡𝐷 is the
dimensionless time, 𝑃𝐷 is the dimensionless pressure, 𝛼𝐷 is 3. Numerical Method
the dimensionless compressibility, 𝜆 𝐷 is the dimensionless
TPG, and 𝛿 is the dimensionless moving boundary. Consider Due to the existence of the moving boundary in the dimen-
sionless mathematical model, the flow region is not fixed and
1 𝜕𝑃𝐷 𝜕2 𝑃𝐷 𝜆 𝐷 𝜕𝑃𝐷 2 𝜕𝑃𝐷 (𝑟𝐷, 𝑡𝐷) expands outward continuously with time increasing. In order
+ 2
+ − 𝛼𝐷 ( ) = , to overcome this difficulty in the space discretization for the
𝑟𝐷 𝜕𝑟𝐷 𝜕𝑟𝐷 𝑟𝐷 𝜕𝑟𝐷 𝜕𝑡𝐷 (10) transient flow region with the moving boundary, a spatial
1 ≤ 𝑟𝐷 ≤ 𝛿 (𝑡𝐷) , coordinate transformation method is used as follows [21, 38]:
𝑃𝐷𝑡𝐷 =0 = 0, (11) 𝑟𝐷 − 1
𝑦 (𝑟𝐷, 𝑡𝐷) = , 1 ≤ 𝑟𝐷 ≤ 𝛿 (𝑡𝐷) . (21)
𝛿 (𝑡𝐷) − 1
𝛿 (0) = 0, (12)
𝜕𝑃𝐷 Through (21), the dynamic flow region for the moving
𝑟𝐷 = − (1 + 𝜆 𝐷) ,
𝜕𝑟𝐷 𝑟𝐷 =1
(13) boundary problem [0, 𝛿(𝑡𝐷)] can be transformed to a fixed
region [0, 1], and the dimensionless pressure 𝑃𝐷(𝑟𝐷, 𝑡𝐷) can
𝜕𝑃𝐷 be transformed as a new unknown function 𝜂 of two variables
= −𝜆 𝐷,
𝜕𝑟𝐷 𝑟𝐷 =𝛿(𝑡𝐷 )
(14) 𝑦 and 𝑡𝐷, that is, 𝜂(𝑦, 𝑡𝐷) equivalently; correspondingly,
the differential variables can be transformed, respectively, as
𝑃𝐷𝑟𝐷 =𝛿(𝑡𝐷 ) = 0. (15) follows:
Substituting (22)–(24) into (13)–(15) and (20), respec- From (26) and (27), the following equation can be
tively, yields obtained [21]:
Substituting (29) and (30) into (25) to cancel the variables 𝑗+1 𝑗+1 𝑗+1 𝑗+1 2
𝜂𝑖+1 − 𝜂𝑖 1 𝜂 − 𝜂0
𝜕𝛿/𝜕𝑡𝐷 and 𝛿 yields [− ⋅( 1 )
Δ𝑦 1 + 𝜆𝐷 Δ𝑦
[
2
𝜕𝜂 [ 1 𝜕𝜂 𝑦 𝑗+1 𝑗+1
⋅ − ⋅ ( ) − 𝑖Δ𝑦 𝑖Δ𝑦 𝜂 − 𝜂0
𝜕𝑦 1 + 𝜆𝐷
𝜕𝑦 𝑦=0 1 + 𝜆𝐷 − ⋅ (− ⋅ 1 + 1)
[ 1 + 𝜆𝐷 1 + 𝜆𝐷 Δ𝑦
𝑗+1
from the well; 𝑗 denotes the index of a time step; and Δ𝑡𝐷 solve these nonlinear difference equations; when 𝑃𝐷𝑖 (𝑖 =
denotes the time step size [21]. 0, 1, 2, . . . , 𝑁 − 1) are numerically solved for, let 𝑗 be equal
From (28), we have to 𝑗 + 1, and, in the same manner, the numerical solutions of
𝑁 difference equations with respect to 𝑁 unknown variables
𝑗+1
𝜂𝑁 = 0. (35) 𝑃𝐷𝑖 𝑗+2 (𝑖 = 0, 1, 2, . . . , 𝑁 − 1) at the (𝑗 + 2)th time step can also
be numerically solved for; the rest can be deduced by analogy
Then, from (34) and (35), the difference equation corre- [21, 39]. Eventually, numerical solutions for the transformed
sponding to the (𝑁 − 1)th spatial grid can be expressed as nonlinear partial differential equation system with respect to
follows: 𝜂(𝑦, 𝑡𝐷) can be obtained.
𝑗+1 2 The difference equation of (21) is
𝑗+1 𝑗+1
−𝜂𝑁−1 1 𝜂 − 𝜂0 (𝑁 − 1) Δ𝑦
⋅ [− ⋅( 1 ) − 𝑗+1
𝑟𝐷𝑖 = 𝑖Δ𝑦 ⋅ (𝛿𝑗+1 − 1) + 1. (39)
Δ𝑦 1 + 𝜆𝐷 Δ𝑦 1 + 𝜆𝐷
[
𝑗+1 𝑗+1 The difference equation of (30) is
(𝑁 − 1) Δ𝑦 𝜂1 − 𝜂0
⋅ (− ⋅ + 1) 𝑗+1 𝑗+1
1 + 𝜆𝐷 Δ𝑦 1 𝜂 − 𝜂0
𝛿𝑗+1 − 1 = − ⋅ 1 . (40)
𝑗+1 𝑗+1
1 + 𝜆𝐷 Δ𝑦
1 𝜂𝑁−2 − 2𝜂𝑁−1 2
⋅( 2
⋅ (1 + 𝜆 𝐷) Substituting (40) into (39) yields
𝜆𝐷 (Δ𝑦)
𝑗+1 𝑗+1
𝑗+1 𝑗+1 2 𝑗+1 𝜂1 − 𝜂0
𝜂1 − 𝜂0 𝑟𝐷𝑖 = −𝑖 ⋅ + 1. (41)
− 𝛼𝐷 ⋅ 𝜆 𝐷 ⋅ ( ) )] 1 + 𝜆𝐷
Δ𝑦
] By (41), numerical solutions of 𝜂(𝑦, 𝑡𝐷) can be trans-
𝑗+1
𝜂𝑁−2 −
𝑗+1
2𝜂𝑁−1 (𝑁 − 1) Δ𝑦
𝑗+1
𝜂1 −
𝑗+1
𝜂0 formed as the ones of 𝑃𝐷(𝑟𝐷, 𝑡𝐷) in the process of numerical
+ 2
⋅ (− ⋅ + 1) solutions at every time step [21]. From (41), it can be further
(Δ𝑦) 1 + 𝜆𝐷 Δ𝑦 concluded that the presented spatial coordinate transforma-
3 tion lets the time-dependent space discretization in the spa-
𝑗+1 𝑗+1 2 𝑗+1 𝑗+1
𝜂1 − 𝜂0 1 𝜂 − 𝜂0 tial coordinate of 𝑟𝐷 be transformed as a time-independent
×( ) + 𝜆 𝐷( )( 1 )
Δ𝑦 1 + 𝜆𝐷 Δ𝑦 space discretization in the new spatial coordinate of y and
then makes the numerical solutions by the finite difference
𝑗+1 2 𝑗+1 𝑗+1 method more applicable and simpler [21].
𝜂𝑁−1 (𝑁 − 1) Δ𝑦 𝜂1 − 𝜂0
− 𝛼𝐷( ) ⋅ (− ⋅ + 1)
Δ𝑦 1 + 𝜆𝐷 Δ𝑦
4. Results and Discussions
𝑗+1 𝑗+1 𝑗+1 𝑗
𝜂1 − 𝜂0 𝜂𝑁−1 − 𝜂𝑁−1 4.1. Effect of QPGT under Different Values of TPG. Figures
×( )−
Δ𝑦 Δ𝑡𝐷 2–4 show the effect of the QPGT on numerical solutions of
the model, with respect to the dimensionless formation pres-
2 𝑗+1 𝑗+1 3
1 𝜂 − 𝜂0 sure distribution, dimensionless transient wellbore pressure,
⋅( ) ⋅( 1 ) and dimensionless transient distance of moving boundary,
1 + 𝜆𝐷 Δ𝑦
respectively, under different values of the dimensionless
𝑗+1 𝑗+1 TPG. From Figures 2–4, it can be clearly seen that with
(𝑁 − 1) Δ𝑦 𝜂1 − 𝜂0
⋅ (− ⋅ + 1) = 0. an increase of the dimensionless TPG, the effect of the
1 + 𝜆𝐷 Δ𝑦 QPGT on the mathematical model solutions become more
(36) and more obvious; the reason can be explained as follows:
the bigger the TPG, the steeper the formation pressure
The difference equation of (33) is as follows: gradient, which generates the solution deviations, resulting
𝑗+1 𝑗+1 𝑗+1
from neglecting the QPGT in the governing equation, more
𝜂1 − 𝜂0 1 + 𝜆 𝐷 𝜂𝑁−1 seriously. Moreover, from Figures 2–4, it can also be indicated
+ ⋅ = 0. (37)
Δ𝑦 𝜆𝐷 Δ𝑦 that the dimensionless transient wellbore pressure for the
case of neglecting the QPGT is larger than that for the case
From (32), the initial conditions are obtained as follows: of considering the QPGT; and the dimensionless transient
distance of moving boundary for the case of neglecting the
𝜂𝑖0 = 0, 𝑖 = 0, 1, 2, . . . , 𝑁 − 1. (38) QPGT is also larger than that for the case of considering the
QPGT. The QPGT can slow down the velocity of the moving
Equations (34), (36), and (37) together form 𝑁 difference boundary in low-permeable reservoir with TPG.
equations at the (𝑗 + 1)th time step and also contain 𝑁 Tables 1 and 2 show the selected data of calculated relative
𝑗+1
unknown variables 𝑃𝐷𝑖 (𝑖 = 0, 1, 2, . . . , 𝑁 − 1). The Newton- errors 𝜀𝑟 from the already computed dimensionless formation
Raphson iterative method [21, 39] is used to numerically pressure distribution and dimensionless transient wellbore
Mathematical Problems in Engineering 7
35 350
30 300
25 250
PD (rD , 1 × 105 )
20 200
𝛿(tD )
15 150
10 100
5 tD = 1 × 105 50
0 0
0 50 100 150 200 250 300 350 101 102 103 104 105
rD tD
𝜆D = 0.034, 𝛼D = 0 𝜆D = 0.067, 𝛼D = 0.005 𝜆D = 0.034, 𝛼D = 0 𝜆D = 0.067, 𝛼D = 0.005
𝜆D = 0.034, 𝛼D = 0.005 𝜆D = 0.1, 𝛼D = 0 𝜆D = 0.034, 𝛼D = 0.005 𝜆D = 0.1, 𝛼D = 0
𝜆D = 0.067, 𝛼D = 0 𝜆D = 0.1, 𝛼D = 0.005 𝜆D = 0.067, 𝛼D = 0 𝜆D = 0.1, 𝛼D = 0.005
Figure 2: The effect of QPGT on dimensionless formation pressure Figure 4: The effect of QPGT on dimensionless transient distance
distribution under different values of dimensionless TPG. of moving boundary under different values of dimensionless TPG.
160
30
140
25 120
20 100
PD (0, tD )
𝜀r (%)
80
15
60
10
40
5
20
5% error line
0 0
0.0 2.0 × 104 4.0 × 104 6.0 × 104 8.0 × 104 1.0 × 105 0 50 100 150 200 250 300 350
tD rD
Table 1: Comparison data for the computed dimensionless formation pressure distribution and relative errors.
Table 2: Comparison data for the computed dimensionless transient wellbore pressure and relative errors.
10.0
Darcy’s flow, stay a lower level for the change of relative errors
7.5 of the dimensionless transient wellbore pressure with the
5% error line increase of dimensionless time, whereas the change of relative
5.0
errors of the dimensionless formation pressure distribution
2.5 with the increase of dimensionless distance is still very
considerable.
0.0 In conclusion, the QPGT has a great effect on the math-
2.0 × 104 4.0 × 104 6.0 × 104 8.0 × 104 1.0 × 105
ematical model solutions of radial flow in low-permeable
tD
reservoirs with TPG, especially for a large value of the
𝛼D = 0.005 dimensionless TPG; the larger the dimensionless distance
𝜆D = 0.034 and dimensionless time, the bigger the effect of the QPGT on
𝜆D = 0.067 the dimensionless formation pressure and the dimensionless
𝜆D = 0.1
transient wellbore pressure. In summary, the mathematical
models of radial flow in low-permeable reservoirs with TPG
Figure 6: Relative error curves for computing dimensionless tran- should take the QPGT into account in their governing
sient wellbore pressure by neglecting quadratic pressure gradient equations, which is more important than those of Darcy’s
term.
flow.
Mathematical Problems in Engineering 9
250
25
200
20
PD (rD , 1 × 105 )
15 150
𝛿(tD )
10 100
tD = 1 × 105
5 50
𝜆D = 0.067
𝜆D = 0.067
0
0
0 50 100 150 200 250
101 102 103 104 105
rD tD
𝛼D = 0.003 𝛼D = 0.007 𝛼D = 0.003 𝛼D = 0.007
𝛼D = 0.005 𝛼D = 0.009 𝛼D = 0.005 𝛼D = 0.009
Figure 7: The sensitive effect of QPGT on dimensionless formation Figure 9: The sensitive effect of QPGT on dimensionless transient
pressure distribution. distance of moving boundary.
25 5. Conclusions
The existence of a TPG can lead to relatively steep formation
20 pressure gradients, and, thus, it is not appropriate to neglect
the QPGT for the fluid flow in porous media with the TPG.
PD (0, tD )
Equation (A.2) can be rewritten as follows: Expanding the right-hand side of (4) yields
𝜕𝜌 𝜕𝑝 𝜕 (𝜌𝜙) 𝜕𝜙 𝜕𝜌
= 𝐶𝑓 ⋅ 𝜌 . (A.3) =𝜌 +𝜙 . (A.7)
𝜕𝑟 𝜕𝑟 𝜕𝑡 𝜕𝑡 𝜕𝑡
In the same manner as above, from (1) and (2), the Substituting (A.4) into the right-hand side of (A.7) yields
following equations can also be deduced as
𝜕 (𝜌𝜙) 𝜕𝑝 𝜕𝑝
𝜕𝜌 𝜕𝑝 = 𝜌 ⋅ 𝐶𝜙 ⋅ 𝜙 ⋅ + 𝜙 ⋅ 𝐶𝑓 ⋅ 𝜌 ⋅
= 𝐶𝑓 ⋅ 𝜌 , 𝜕𝑡 𝜕𝑡 𝜕𝑡
𝜕𝑡 𝜕𝑡 (A.8)
(A.4) 𝜕𝑝 𝜕𝑝
𝜕𝜙 𝜕𝑝 =𝜌⋅𝜙⋅ ⋅ (𝐶𝜙 + 𝐶𝑓 ) = 𝜌 ⋅ 𝜙 ⋅ ⋅ 𝐶𝑡 .
= 𝐶𝜙 ⋅ 𝜙 . 𝜕𝑡 𝜕𝑡
𝜕𝑡 𝜕𝑡
The left-hand side of (4) can be expanded as follows: Substituting (A.6) and (A.8) into (4), the governing
equation in consideration of the QPGT can be obtained as
1 𝜕 follows:
− ⋅ (𝑟 ⋅ 𝜌 ⋅ 𝜐)
𝑟 𝜕𝑟
𝑘 1 𝜕 𝜕𝑝 𝑘𝜌 [ 2
[ 𝜕 𝑝 + 1 ⋅ 𝜕𝑝 − 𝜆 + 𝜕𝑝 2 ] ]
= ⋅ ⋅ (𝑟 ⋅ 𝜌 ⋅ ( − 𝜆)) [ 𝐶 𝑓 ⋅ ( ) ]
𝜇 𝑟 𝜕𝑟 𝜕𝑟 𝜇 ⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟
𝜕𝑟2 𝑟 𝜕𝑟 𝑟 𝜕𝑟
⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟
(A.9)
[ Main Term Quadratic Gradient Term]
𝑘 1 𝜕𝑝 𝑘 1 𝜕𝜌 𝜕𝑝
= ⋅ ⋅ (𝜌 ⋅ ( − 𝜆)) + ⋅ ⋅ (𝑟 ⋅ ⋅ ( − 𝜆)) 𝜕𝑝
𝜇 𝑟 𝜕𝑟 𝜇 𝑟 𝜕𝑟 𝜕𝑟 =𝜌⋅𝜙⋅ ⋅ 𝐶𝑡 .
𝜕𝑡
𝑘 1 𝜕2 𝑝
+ ⋅ ⋅ (𝑟 ⋅ 𝜌 ⋅ 2 ) Equation (A.9) can be equivalently simplified, by cancel-
𝜇 𝑟 𝜕𝑟 ing the variable 𝜌 in its both sides, as follows:
𝑘 1 𝜕𝑝 𝑘 1
= ⋅ ⋅ (𝜌 ⋅ ( − 𝜆)) + ⋅ 𝜕2 𝑝 1 𝜕𝑝 𝜆 𝜕𝑝 2 𝜇𝜙𝐶𝑡 𝜕𝑝
𝜇 𝑟 𝜕𝑟 𝜇 𝑟 + ⋅ − + 𝐶𝑓 ⋅ ( ) = . (A.10)
𝜕𝑟2 𝑟 𝜕𝑟 𝑟 𝜕𝑟 𝑘 𝜕𝑡
𝜕𝑝 𝜕𝑝
⋅ (𝑟 ⋅ 𝐶𝑓 ⋅ 𝜌 ⋅ ⋅ ( − 𝜆))
𝜕𝑟 𝜕𝑟 Conflict of Interests
2
𝑘 1 𝜕𝑝 Wenchao Liu and Jun Yao declare that there is no conflict of
+ ⋅ ⋅ (𝑟 ⋅ 𝜌 ⋅ 2 )
𝜇 𝑟 𝜕𝑟 interests regarding the publication of this paper.
𝑘𝜌 𝜕2 𝑝 1 𝜕𝑝 𝜆 𝜕𝑝 𝜕𝑝
= [ 2 + ⋅ − + 𝐶𝑓 ⋅ ⋅ ( − 𝜆)] Acknowledgments
𝜇 𝜕𝑟 𝑟 𝜕𝑟 𝑟 𝜕𝑟 𝜕𝑟
The authors would like to acknowledge the funding by the
𝑘𝜌 Projects Grants nos. 11102237 and 51404232, both sponsored
=
𝜇 by the Natural Science Foundation of China (NSFC), and the
Program for Changjiang Scholars and Innovative Research
Team in University (Grant no. IRT1294).
[ 𝜕2 𝑝 1 𝜕𝑝 𝜆 𝜕𝑝 2 𝜕𝑝 ]
×[ + ⋅ −
[⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟ + 𝐶 ⋅ ( ) − 𝐶𝑓 ⋅ 𝜆 ⋅ ] .
𝜕𝑟 2 𝑟 𝜕𝑟 𝑟 𝑓
𝜕𝑟
⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟ 𝜕𝑟 ]
⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟
References
[ Main Term Quadratic Gradient Term Small Term
]
(A.5) [1] Y. Z. Huang, Z. M. Yang, Y. He, and X. W. Wang, “An overview
on nonlinear porous flow in low permeability porous Media,”
Because 𝜆 ≪ 1 and 𝐶𝑓 ≪ 1, the small term in the right- Theoretical & Applied Mechanics Letters, vol. 3, no. 2, Article ID
hand side of (A.5) can be neglected, and then (A.5) can be 022001, 2013.
rewritten as follows: [2] P. J. M. Monteiroa, C. H. Rycroftc, and G. I. Barenblatt, “A
1 𝜕 mathematical model of fluid and gas flow in nanoporous
− ⋅ (𝑟 ⋅ 𝜌 ⋅ 𝜐) media,” Proceedings of the National Academy of Sciences of the
𝑟 𝜕𝑟
United States of America, vol. 109, no. 50, pp. 20309–20313, 2012.
(A.6) [3] R. Z. Yu, Y. N. Bian, Y. Li et al., “Non-Darcy flow numerical sim-
𝑘𝜌 [ 2
[ 𝜕 𝑝 + 1 ⋅ 𝜕𝑝 − 𝜆 + 𝜕𝑝 2 ] ]. ulation of XPJ low permeability reservoir,” Journal of Petroleum
= [ 𝐶𝑓 ⋅ ( ) ]
𝜇 ⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟
𝜕𝑟 2 𝑟 𝜕𝑟 𝑟 𝜕𝑟
⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟ Science and Engineering, vol. 92-93, pp. 40–47, 2012.
[ Main Term Quadratic Gradient Term] [4] J. Guo, S. Zhang, L. Zhang, H. Qing, and Q. Liu, “Well testing
analysis for horizontal well with consideration of threshold
In (A.6), the QPGT is retained for the deduction of the pressure gradient in tight gas reservoirs,” Journal of Hydrody-
governing equation. namics, vol. 24, no. 4, pp. 561–568, 2012.
Mathematical Problems in Engineering 11
[5] B. Q. Zeng, L. S. Cheng, and C. L. Li, “Low velocity non-linear [22] W. Liu, J. Yao, and Z. Chen, “Analytical solution of a dou-
flow in ultra-low permeability reservoir,” Journal of Petroleum ble moving boundary problem for nonlinear flows in one-
Science and Engineering, vol. 80, no. 1, pp. 1–6, 2011. dimensional semi-infinite long porous media with low perme-
[6] J. C. Cai and B. M. Yu, “A discussion of the effect of tortuosity on ability,” Acta Mechanica Sinica, vol. 30, no. 1, pp. 50–58, 2014.
the capillary imbibition in porous media,” Transport in Porous [23] S. L. Marshall, “Nonlinear pressure diffusion in flow of com-
Media, vol. 89, no. 2, pp. 251–263, 2011. pressible liquids through porous media,” Transport in Porous
[7] Y. D. Yao and J. L. Ge, “Characteristics of non-Darcy flow in Media, vol. 77, no. 3, pp. 431–446, 2009.
low-permeability reservoirs,” Petroleum Science, vol. 8, no. 1, pp. [24] M. Bai, Q. Ma, and J. Roegiers, “A nonlinear dual-porosity
55–62, 2011. model,” Applied Mathematical Modelling, vol. 18, no. 11, pp. 602–
[8] F. Civan, Porous Media Transport Phenomena, John Wiley & 610, 1994.
Sons, Hoboken, NJ, USA, 2011. [25] X. L. Cao, D. K. Tong, and R. H. Wang, “Exact solutions for
[9] Y. Xiang’an, W. Haoguang, Z. Lijuan, Z. Renbao, and Z. nonlinear transient flow model including a quadratic gradient
Yongpan, “Low pressure gas percolation characteristic in ultra- term,” Applied Mathematics and Mechanics, vol. 25, no. 1, pp.
low permeability porous media,” Transport in Porous Media, vol. 102–109, 2004.
85, no. 1, pp. 333–345, 2010. [26] D. K. Tong, H. Q. Zhang, and R. H. Wang, “Exact solution and
[10] F. Hao, L. S. Cheng, O. Hassan, J. Hou, C. Z. Liu, and J. D. its behavior characteristic of nonlinear dual-porosity model,”
Feng, “Threshold pressure gradient in ultra-low permeability Applied Mathematics and Mechanics, vol. 26, no. 10, pp. 1277–
reservoirs,” Petroleum Science and Technology, vol. 26, no. 9, pp. 1283, 2005.
1024–1035, 2008. [27] A. S. Odeh and D. K. Babu, “Comparison of solutions of the
[11] X. Wei, L. Qun, G. Shusheng, H. Zhiming, and X. Hui, “Pseudo nonlinear and linearized diffusion equations,” SPE Reservoir
threshold pressure gradient to flow for low permeability reser- Engineering, vol. 3, no. 4, pp. 1202–1206, 1988.
voirs,” Petroleum Exploration and Development, vol. 36, no. 2,
[28] Y. Wang and M. B. Dusseault, “The effect of quadratic gradient
pp. 232–236, 2009.
terms on the borehole solution in poroelastic Media,” Water
[12] J. Cai, X. Hu, D. C. Standnes, and L. You, “An analytical Resources Research, vol. 27, no. 12, pp. 3215–3223, 1991.
model for spontaneous imbibition in fractal porous media
[29] C. Chakrabarty, S. M. Farouq Ali, and W. S. Tortike, “Effects of
including gravity,” Colloids and Surfaces A Physicochemical and
the nonlinear gradient term on the transient pressure solution
Engineering Aspects, vol. 414, pp. 228–233, 2012.
for a radial flow system,” Journal of Petroleum Science and
[13] J. C. Cai, “A fractal approach to low velocity non-Darcy flow in Engineering, vol. 8, no. 4, pp. 241–256, 1993.
a low permeability porous medium,” Chinese Physics B, vol. 23,
no. 4, Article ID 044701, 2014. [30] C. Chakrabarty, S. M. F. Ali, and W. S. Tortike, “Analytical
solutions for radial pressure distribution including the effects
[14] H. Pascal, “Nonsteady flow through porous media in the
of the quadratic-gradient term,” Water Resources Research, vol.
presence of a threshold gradient,” Acta Mechanica, vol. 39, no.
29, no. 4, pp. 1171–1177, 1993.
3-4, pp. 207–224, 1981.
[31] S. Braeuning, T. A. Jelmert, and S. A. Vik, “The effect of the
[15] Y. S. Wu, K. Pruess, and P. A. Witherspoon, “Flow and displace-
quadratic gradient term on variable-rate well-tests,” Journal of
ment of Bingham non-Newtonian fluids in porous media,” SPE
Petroleum Science and Engineering, vol. 21, no. 3-4, pp. 203–222,
Reservoir Engineering, vol. 7, no. 3, pp. 369–376, 1992.
1998.
[16] S. Fuquan, L. Ciqun, and L. Fanhua, “Transient pressure of per-
colation through one dimension porous media with threshold [32] W. Li, X. Li, S. Li, and Q. Li, “The similar structure of solutions
pressure gradient,” Applied Mathematics and Mechanics (English in fractal multilayer reservoir including a quadratic gradient
Edition), vol. 20, no. 1, pp. 27–35, 1999. term,” Journal of Hydrodynamics, vol. 24, no. 3, pp. 332–338,
2012.
[17] G. Q. Feng, Q. G. Liu, G. Z. Shi, and Z. H. Lin, “An unsteady
seepage flow model considering kickoff pressure gradient for [33] M. Dewei, J. Ailin, J. Chengye, Z. Qian, and H. Dongbo,
low-permeability gas reservoirs,” Petroleum Exploration and “Research on transient flow regulation with the effect of
Development, vol. 35, no. 4, pp. 457–461, 2008. quadratic pressure gradient,” Petroleum Science and Technology,
vol. 31, no. 4, pp. 408–417, 2013.
[18] J. Lu and S. Ghedan, “Pressure behavior of vertical wells in
low-permeability reservoirs with threshold pressure gradient,” [34] R. S. Nie, F. Ge, and Y. L. Liu, “The researches on the nonlinear
Special Topics and Reviews in Porous Media, vol. 2, no. 3, pp. flow model with quadratic pressure gradient and its application
157–169, 2011. for double porosity reservoir,” in Proceedings of the International
[19] J. Lu, “Pressure behavior of uniform-flux hydraulically fractured Forum on Porous Flow and Applications: Flow in Porous Media:
wells in low-permeability reservoirs with threshold pressure From Phenomena to Engineering and Beyond, Wuhan, China,
gradient,” Special Topics & Reviews in Porous Media, vol. 3, no. 2009.
4, pp. 307–320, 2012. [35] R. S. Nie, Y. L. Jia, J. Yu, and Y. L. Liu, “The transient well test
[20] W. Liu, J. Yao, and Y. Wang, “Exact analytical solutions of analysis of fractured- vuggy triple-porosity reservoir with the
moving boundary problems of one-dimensional flow in semi- quadratic pressure gradient term,” in Proceedings of the Latin
infinite long porous media with threshold pressure gradient,” American and Caribbean Petroleum Engineering Conference,
International Journal of Heat and Mass Transfer, vol. 55, no. 21- Cartagena de Indias, Cartagena, Colombia, 2009.
22, pp. 6017–6022, 2012. [36] Y. D. Yao, Y. S. Wu, and R. L. Zhang, “The transient flow analysis
[21] J. Yao, W. C. Liu, and Z. X. Chen, “Numerical solution of a of fluid in a fractal, double-porosity reservoir,” Transport in
moving boundary problem of one-dimensional flow in semi- Porous Media, vol. 94, no. 1, pp. 175–187, 2012.
infinite long porous media with threshold pressure gradient,” [37] J. C. Guo and R. S. Nie, “Nonlinear flow model for well
Mathematical Problems in Engineering, vol. 2013, Article ID production in an underground formation,” Nonlinear Processes
384246, 7 pages, 2013. in Geophysics, vol. 20, no. 3, pp. 311–327, 2013.
12 Mathematical Problems in Engineering
Research Article
A Differential-Algebraic Model for the Once-Through Steam
Generator of MHTGR-Based Multimodular Nuclear Plants
Zhe Dong1,2
1
Institute of Nuclear and New Energy Technology, Tsinghua University, Beijing 100084, China
2
Key Laboratory of Advanced Reactor Engineering and Safety, Ministry of Education of China, Beijing 100084, China
Copyright © 2015 Zhe Dong. This is an open access article distributed under the Creative Commons Attribution License, which
permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Small modular reactors (SMRs) are those fission reactors whose electrical output power is no more than 300 MWe . SMRs usually
have the inherent safety feature that can be applicable to power plants of any desired power rating by applying the multimodular
operation scheme. Due to its strong inherent safety feature, the modular high temperature gas-cooled reactor (MHTGR), which uses
helium as coolant and graphite as moderator and structural material, is a typical SMR for building the next generation of nuclear
plants (NGNPs). The once-through steam generator (OTSG) is the basis of realizing the multimodular scheme, and modeling of the
OTSG is meaningful to study the dynamic behavior of the multimodular plants and to design the operation and control strategy. In
this paper, based upon the conservation laws of mass, energy, and momentum, a new differential-algebraic model for the OTSGs of
the MHTGR-based multimodular nuclear plants is given. This newly-built model can describe the dynamic behavior of the OTSG
in both the cases of providing superheated steam and generating saturated steam. Numerical simulation results show the feasibility
and satisfactory performance of this model. Moreover, this model has been applied to develop the real-time simulation software for
the operation and regulation features of the world first underconstructed MHTGR-based commercial nuclear plant—HTR-PM.
Shutdown rod
Reflector
prevents SMRs from the hazards of core-melting, radiological SMRs usually have the safety features of strong natural-
release, and LOCA (Loss of Coolant Accident). Furthermore, circulation and self-pressurized and are usually tightly cou-
SMRs can offer simpler, safer, and standardized modular pled with the once-through steam generators (OTSGs) to
design by being factory built, requiring smaller initial capital provide superheated steam [3, 4]. The second group are
investment, and having shorter construction period, and have those SMRs cooled either by liquid metal or liquid salt
been viewed by International Atomic Energy Agency (IAEA) such as the 4S (Super-Safe, Small and Simple) fast reactor
as one of the active developing trends of nuclear energy. designed by Toshiba, fueled with enriched uranium or
There are three major groups of SMR designs that are plutonium and cooled by sodium. The schematic view of
actively being developed in the US, China, Japan, Korea, the 4S nuclear energy system is shown in Figure 2 from
and other countries. As shown in Figure 1, the first group of which we can see that the OTSG is also utilized here for
SMRs is based on the design concepts of proven and widely providing super-heated steam [4]. The third group consists of
used light water reactors (LWRs) such as the IRIS (interna- gas-cooled SMRs. The modular high temperature gas-cooled
tional reactor innovative and secure), NuScale and mPower reactor (MHTGR), which uses helium as coolant and graphite
designed by US, and the SMART designed by Korea. These as both moderator and structural materials, is a crucial
Mathematical Problems in Engineering 3
Steam turbine
Generator
Hot chamber
OTSG
Coaxial duct
Steam valve
Condenser
Feedwater pump
gas-cooled SMR and has been seen as one the best candidates temperature, the OTSG is more proper for building multi-
for building the next generation of nuclear plant (NGNP). SMR based nuclear plants [10]. Therefore, it is clear that the
The inherent safety feature of the MHTGR is given by its low dynamic model of the OTSG is very crucial for studying
power density, strong negative temperature feedback effect, behavior of the SMR-based plant and also for verifying the
and slim reactor shape [5]. The first MHTGR of China, that operation and control strategies. Up to now, there have
is, 10 MWth pebble-bed high temperature gas-cooled reactor already been some works in building the dynamic models of
HTR-10, which was developed by the institute of nuclear OTSGs. Through dividing the OTSG secondary side into the
and new energy technology (INET) of Tsinghua University, subcooled, evaporating, and superheated sections, Ray gave a
achieved its full power-level in 2003 [6]. Based upon HTR-10, nonlinear dynamic model to the OTSG of a solar plant [11, 12].
the high temperature gas-cooled reactor pebble-bed module By further dividing the aforementioned evaporating section
(HTR-PM) plant has already been designed by INET and is into the two regions of nucleate boiling and film boiling,
now under construction. From Figure 3, HTR-PM is a two- Tzanos and Abdalla independently proposed the four region
modular nuclear plant that is the first MHTGR-based under- moving boundary models for the helical-coil OTSGs of the
constructed commercial plant. Each module of HTR-PM is a advanced liquid metal reactor (ALMR) [13, 14]. Recently, Li
nuclear steam supply system (NSSS) constituted by a pebble- et al. gave a three-region dynamic model for the OTSG of the
bed one-zone MHTGR, a side-by-side arranged helical-coil HTR-10 plant [15].
OTSG and some connecting pipes [7, 8]. Similar to the first The above OTSG models are all developed under the
and second groups of SMRs, the outlet steam of HTR-PM’s assumption that the outlet steam is superheated. However,
NSSS is also superheated. in some actual cases such as the plant startup or operation
From the above introduction, it is clear to see that the in very low power, the OTSG might also generate saturated
OTSG is key equipment in the SMR-based nuclear energy steam. The OTSG dynamics in these cases is quite significant
systems. Actually, this is given by the necessity of building for designing the plant operation strategy, but it cannot be
large nuclear plants based on parallel-operated multi-SMRs. described by the models presented in [11–15]. Therefore, it
The precondition of applying this parallel operation scheme is very necessary to give OTSG models for describing the
is that the pressure of the steam generated by each NSSS dynamic behavior in the cases of generating saturated steam.
must be equal to each other. Since the widely-utilized U-tube In this paper, a moving boundary region dynamic model
steam generators (UTSG) can only provide saturated steam of the OTSG of the MHTGR-based multimodular nuclear
whose pressure and temperature must obey the one-to-one plant is proposed. This model can be used to describe the
relationship, parallel-operating SMRs based on the UTSGs OTSG dynamics in the cases of providing superheated steam
may lead to the drift in the steady values of the coolant tem- and generating saturated steam. Then, numerical simulation
peratures [9]. Since the OTSG can provide superheated steam results show the feasibility of this new model. This model
which does not satisfy the one-to-one map from pressure to has already been adopted to develop the real-time simulation
4 Mathematical Problems in Engineering
Since the water near point 3 in Figure 4 is saturate, from From Lemma 1 and by integrating (1), (2), and (7) along
the thermodynamic features of saturate water, we can derive the direction of 𝑥-axis in the boiling section, we can obtain
that that
𝑑ℎ𝑠3 𝜕ℎ𝑠3 𝑑𝑃𝑠3
= .
𝜕𝑃𝑠3 𝑃𝑠3 𝑑𝑡
(9) 𝑑 (𝜌𝑠4 𝑙35 )
𝑑𝑡 𝑑𝑙
+ 𝜌𝑠3 13 + 𝐷𝑠5 − 𝐷𝑠3 = 0, (16)
𝑑𝑡 𝑑𝑡
Moreover, it is not loss generality to assume that
𝑑 (𝜌𝑠4 ℎ𝑠4 𝑙35 ) 𝑑𝑙
𝐺𝑠2 = 𝐷𝑠2 𝐴 𝑠2 , + 𝜌𝑠3 ℎ𝑠3 13 + 𝐷𝑠5 ℎ𝑠5 − 𝐷𝑠3 ℎ𝑠3
𝑑𝑡 𝑑𝑡
(17)
𝐺 + 𝐺𝑠3 𝑄 𝑑 (𝑃𝑠4 𝑙35 ) 𝑑𝑙
𝐺𝑠2 = 𝑠1 , = 𝑠4 + + 𝑃𝑠3 13 ,
2 𝐴 𝑠4 𝑑𝑡 𝑑𝑡
ℎ𝑠1 + ℎ𝑠3 (10)
2
ℎ𝑠2 = , 𝑃𝑠5 = 𝑃𝑠3 − 𝜌𝑠4 𝑔𝑙35 sin 𝜃 − 𝐹𝑠4 𝑙35 𝐷𝑠4 . (18)
2
𝑃𝑠1 + 𝑃𝑠3
𝑃𝑠2 = , Based upon relationship between (14)–(15) and (18) and
2
by subtracting the multiplication of (16) and ℎ4 from (17), we
where 𝐺𝑠 is the flowrate of the secondary side. have
Based on (8)–(10) and by subtracting the multiplication
of (5) and ℎ2 from (6), we have
1 𝑑ℎ
𝜌𝑠4 𝑙35 𝑠5
1 2 𝑑𝑡
[𝜌 (ℎ − ℎ𝑠3 ) + (𝑃𝑠3 − 𝑃𝑠1 )
2 𝑠3 𝑠1 1
+ {𝜌 (ℎ − ℎ𝑠5 ) + (𝑃𝑠5 − 𝑃𝑠3 )
𝜕ℎ𝑠3 𝐺2 𝑑𝑙13 2 𝑠3 𝑠3
+ 𝑙13 (1 − 𝜌𝑠2 ) (𝜌𝑠2 𝑔 sin 𝜃 + 𝐹𝑠2 2𝑠2 )]
𝜕𝑃𝑠3 𝑃𝑠3 𝐴 𝑠2 𝑑𝑡 𝜕ℎ𝑠3 𝐺2
+ 𝑙35 [(2 − 𝜌𝑠4 ) (𝜌𝑠2 𝑔 sin 𝜃 + 𝐹𝑠2 2𝑠2 )
1 𝑑ℎ 𝐹 𝐺2 𝑙2 𝜕ℎ 𝑑𝐺𝑠3 𝜕𝑃𝑠3 𝑃𝑠3 𝐴 𝑠2
+ 𝜌𝑠2 𝑙13 𝑠1 + 𝑠2 𝑠22 13 (1 − 𝜌𝑠2 𝑠3 )
2 𝑑𝑡 2𝐴 𝑠2 𝜕𝑃𝑠3 𝑃𝑠3 𝑑𝑡 2
𝐺𝑠4 𝑑𝑙
− (𝜌𝑠4 𝑔 sin 𝜃 + 𝐹𝑠4 2
)]} 13
𝑄𝑠2 1 𝐺𝑠1 𝐺𝑠3 𝐴 𝑠4 𝑑𝑡
= + ( + ) (ℎ𝑠1 − ℎ𝑠3 ) .
𝐴 𝑠2 2 𝐴 𝑠1 𝐴 𝑠3 1 𝜕ℎ 𝐹 𝐺 𝑙
(11) + 𝑙35 [(1 − 𝜌𝑠4 𝑠3 ) 𝑠2 2𝑠2 13
2 𝜕𝑃𝑠3 𝑃𝑠3 𝐴 𝑠2
By considering the transport inertia in the subcooled
section, the dynamics of the secondary-side flowrate and 2𝐹𝑠4 𝐺𝑠4 𝑙13 𝑑𝐺𝑠3
+ ]
water enthalpy at point 1 can be described by 𝐴2𝑠2 𝑑𝑡
𝑑𝐺𝑠3 𝐺𝑠1 − 𝐺𝑠3 𝑄𝑠4 + 𝐺𝑠4 (ℎ𝑠3 − ℎ𝑠5 )
= , (12) = .
𝑑𝑡 𝜏𝐺 𝐴 𝑠4
(19)
𝑑ℎ𝑠1 ℎ𝑓𝑤 − ℎ𝑠1
= , (13)
𝑑𝑡 𝜏𝐻
3.1.3. State-Space Model of the Secondary Side. Define the
respectively, where both 𝜏𝐺 and 𝜏𝐻 are given positive con- state-vector and input-vector of the OTSG secondary side as
stants.
𝑇
3.1.2. Boiling Section. Due to the high velocity and relative x𝑠 = [𝑙13 ℎ𝑠5 ℎ𝑠1 𝐺𝑠3 ] , (20)
low density of the fluid in the boiling section, we assume that
𝑇
u𝑠 = [𝐺𝑠1 ℎ𝑓𝑤 𝑄𝑠2 𝑄𝑠4 ] , (21)
𝐺𝑠3 = 𝐺𝑠4 = 𝐺𝑠5 , (14)
𝑄𝑠2 1 𝐺𝑠1 𝐺𝑠3 the state-space model of the tube wall can be written as
+ ( + ) (ℎ𝑠1 − ℎ𝑠3 )
[ 𝐴 𝑠2 2 𝐴 𝑠1 𝐴 𝑠3 ]
[ ]
[ 𝑄𝑠4 + 𝐺𝑠4 (ℎ𝑠3 − ℎ𝑠5 ) ] ẋ𝑚 = A𝑚 (x)̇ x𝑚 + B𝑚 (x) u𝑚 ,
[ ] (30)
[ ]
[ 𝐴 𝑠4 ]
f𝑠 (x𝑠 , u𝑠 ) = [
[
].
]
[ ℎ𝑓𝑤 − ℎ𝑠1 ] where
[ ]
[ 𝜏 ]
[ 𝐻 ]
[ 𝐺 −𝐺 ] ̇ −1 1
𝑙13
𝑠1 𝑠3
A𝑚 (x)̇ = [ ],
[ 𝜏𝐺 ] 𝑙15 −1 1
(23) (31)
𝑇
1 1 1
B𝑚 (x) = diag ([ ] ).
3.2. Differential Equations of the Tube-Wall. Since there is 𝐶𝑚 𝜌𝑚 𝐴 𝑚2 𝑙13 𝐴 𝑚4 (𝑙15 − 𝑙13 )
no flow inside the metal tube wall between the primary and
secondary sides, it is clear that partial differential equation (2) 3.3. Algebraic Equations of the Primary Side. According to the
describing the energy conservation law can be simplified as nodalization scheme illustrated in Figure 4, the primary side
𝜕 (𝜌ℎ) 𝑞 is also divided into two sections. Since the primary coolant
= . (24) is single-phase helium whose flowing velocity is much faster
𝜕𝑡 𝐴
than that of the two-phase flow inside the secondary side,
Here, for the simplicity of the model, we assume that the the algebraic equations are used to give the energy and
tube temperature is linearly distributed along 𝑥-axis defined temperature relationship of the first side.
Mathematical Problems in Engineering 7
From the conservation law of energy, it is clear that Based on (32)–(34), the algebraic equation for describing
the operation characteristics of the primary side can be
𝑄𝑝2 = 𝐶𝑝 𝐺𝑝 (𝑇𝑝3 − 𝑇𝑝1 ) , written as
(32)
𝑄𝑝4 = 𝐶𝑝 𝐺𝑝 (𝑇𝑝5 − 𝑇𝑝3 ) , A𝑝 (x𝑠 ) x𝑝 = b𝑝 (x𝑚 , x𝑝 ) , (35)
where 𝑇𝑠𝑖 is the temperature of the two-phase flow inside the 3.4. Differential-Algebraic State-Space Model. Based upon the
secondary side at point 𝑖 (𝑖 = 1, . . . , 5), 𝐾𝑝2 and 𝐾𝑝4 are the above analysis and derivation, define the state-vector x and
heat transfer coefficients between the primary side and tube- input-vector u as
wall of the subcooled and boiling sections, respectively, and
𝐾𝑠2 and 𝐾𝑠4 are, respectively, the heat transfer coefficients 𝑇
x = [x𝑠𝑇 x𝑚
𝑇
x𝑝𝑇 ] ,
between the tube-wall and secondary side of the subcooled (39)
and boiling sections. 𝑇
u = [𝐺𝑠1 ℎ𝑓𝑤 𝐺𝑝 𝐺𝑝 𝑇𝑝5 ] ,
Moreover, assume that the helium temperature distribu-
tion along the 𝑥-axis defined in Figure 4 is linear, and then we where x𝑠 , x𝑚 , and x𝑝 are given by (20), (28), and (36),
can easy obtain that respectively. Then, the differential-algebraic model of the
OTSG for MHTGR-based multimodular nuclear plants can
𝑇𝑝1 = 2𝑇𝑝2 − 𝑇𝑝3 , be summarized as
𝑇𝑝2 𝑙35 + 𝑇𝑝4 𝑙13 (34) E (x) ẋ= f (x) + G (x) u, (40)
𝑇𝑝3 = .
𝑙15 where
Ξ (x) O6×4
E (x) = [ ],
O4×6 O4
𝑇
G (x) = [G𝑇1 (x) O4 G𝑇2 (x)] ,
−1
𝐺𝑠3 (2𝐴 𝑠3 ) (ℎ𝑠1 − ℎ𝑠3 ) + 𝐾𝑠2 𝑙13 𝐴−1 𝑠2 (𝑇𝑚2 − 𝑇𝑠2 )
[ ]
[ 𝐴−1 [𝐺 (ℎ − ℎ ) + 𝐾 (𝑙 − 𝑙 ) (𝑇 − 𝑇 )] ]
[ 𝑠4 𝑠4 𝑠3 𝑠5 𝑠4 15 13 𝑚4 𝑠4 ]
[ ]
[ −ℎ 𝜏 −1 ]
[ 𝑠1 𝐻 ]
[ ]
[ −𝐺𝑠3 𝜏𝐺 −1 ]
[ ]
[ ]
[(𝐶 𝜌 𝐴 )−1 [𝐾 𝑇 + 𝐾 𝑇 − (𝐾 + 𝐾 ) 𝑇 ]]
[ 𝑚 𝑚 𝑚2 𝑝2 𝑝2 𝑠2 𝑠2 𝑝2 𝑠2 𝑚2 ]
f (x) = [
[ −1
],
]
[(𝐶𝑚 𝜌𝑚 𝐴 𝑚4 ) [𝐾𝑝4 𝑇𝑝4 + 𝐾𝑠4 𝑇𝑠4 − (𝐾𝑝4 + 𝐾𝑠4 ) 𝑇𝑚4 ]]
[ ]
[ ]
[ 2𝑇𝑝2 − (𝑇𝑝1 + 𝑇𝑝3 ) ]
[ ]
[ ]
[ 𝑙15 𝑇𝑝3 − (𝑙15 − 𝑙13 ) 𝑇𝑝2 − 𝑙13 𝑇𝑝4 ]
[ ]
[ ]
[ 𝐾𝑝2 𝑙13 (𝑇𝑝2 − 𝑇𝑚2 ) ]
[ ]
[ 𝐾𝑝4 (𝑙15 − 𝑙13 ) (𝑇𝑝4 − 𝑇𝑚4 ) ]
8 Mathematical Problems in Engineering
𝜌𝑠2 𝑙13
[ 𝐸𝑠11 (x𝑚 ) 0
2
𝐸𝑠14 (x𝑚 ) 0 0]
[ 𝜌 ]
[ 𝑠4 𝑙35 ]
[ 𝐸𝑠21 (x𝑚 ) 0 𝐸𝑠24 (x𝑚 ) 0 0]
[ 2 ]
[ ]
Ξ (x) = [
[ 0 0 𝜏𝐻 0 0 0]
],
[ 0 0 0 𝜏𝐺 0 0]
[ ]
[ ]
[(𝑇𝑚2 − 𝑇𝑚4 ) 𝑙 −1
0 0 0 1 0]
[ 15 ]
−1
[(𝑇𝑚2 − 𝑇𝑚4 ) 𝑙15 0 0 0 0 1]
−1
(2𝐴 𝑠1 ) (ℎ𝑠1 − ℎ𝑠3 ) 0 0 0
[ ]
[ 0 0 0 0]
G1 (x) = [
[
],
[ 0 −1
𝜏𝐻 0 0]
]
[ 𝜏𝐺−1 0 0 0]
0 0 𝐶𝑝 (𝑇𝑝1 − 𝑇𝑝3 ) 0
G2 (x) = [ ].
0 0 𝐶𝑝 𝑇𝑝3 −𝐶𝑝
(41)
In the process of developing dynamic model (40), there temperature 𝑇𝑝1 are shown in Figure 5. Here, steam quality
is no assumption on the outlet steam to be superheated, 𝑥𝑀 is defined by
which means that this model can describe the dynamic
behavior of the OTSG in the case of generating saturate ℎ𝑠5 − ℎ𝑙
𝑥𝑀 = , (42)
steam and also can also roughly describe the dynamics in ℎ𝑔𝑙
the case of providing super-heated steam. Moreover, due to
the simplicity of model (40), it can be used to design steam where ℎ𝑙 is specific enthalpy of saturated steam corresponding
temperature controller for the OTSG. to outlet steam pressure 𝑃𝑠5 , and ℎ𝑔𝑙 is the latent heat of
vaporization corresponding to 𝑃𝑠5 . If 𝑥𝑀 > 1, then the outlet
steam is superheated. If 𝑥𝑀 ≤ 1, then the outlet steam is
4. Numerical Simulation Results saturate.
In order to verify the feasibility of differential-algebraic OTSG
model (40), we applied it to simulate the dynamic behavior of Case B: Step Decrease in the Feedwater Flowrate at 75% RFP.
the OTSG in the HTR-PM plant. The numerical simulation The system operates at the steady state of 75% RFP for 2000 s,
is done in MATLAB/ SIMULINK environment. The model and then a step decrease in the feed-water flowrate with the
parameters are given by physical and thermal design of value of 5 kg/s is added. The transient responses of outlet
this OTSG. In the following, both the steady and transient steam temperature 𝑇𝑠5 , length of the subcooed section 𝑙13 ,
simulations are performed, and some necessary discussions outlet steam quality 𝑥𝑀, and outlet cold helium temperature
are also given. 𝑇𝑝1 are all illustrated in Figure 6.
Case A: Step Decease in the Helium Flowrate at 75% RFP. After 4.3. Discussions. By comparing the simulated and designed
the system operates at the steady state of 75% RFP for 2000 s, a values given in Table 1, we can clearly see that the maxi-
step decrease of 5 kg/s in the primary helium flowrate is added mal relative error is no more than 2%, which shows that
to the system, and the corresponding transient responses of differential-algebraic model (40) has a high steady precision.
outlet steam temperature 𝑇𝑠5 , the length of the subcooed In case A, since the inlet hot helium temperature is con-
section 𝑙13 , the outlet steam quality 𝑥𝑀, and oulet cold helium stant, the step decrease in the helium flowrate results in
Mathematical Problems in Engineering 9
580 34
560
32
Ts5 (∘ C)
l13 (m)
540
30
520
500 28
2000 2200 2400 2000 2200 2400
Time (s) Time (s)
1.85 230
1.8
Tp1 (∘ C)
xM
1.75 225
1.7
1.65 220
2000 2200 2400 2000 2200 2400
Time (s) Time (s)
the decrease of the heat transferred from the primary to Figure 5, the above physical analysis well copes with the
the secondary side, which then induce the decreases in the numerical results.
outlet steam temperature, steam quality and length of the In the case of feedwater flowrate decrease at 75% RFP,
boiling section. Moreover, thermal power transferred from since both the inlet helium temperature and primary helium
the primary to the secondary sides is also reduced, which flowrate remain constant when the step decrease in the
further causes the decrease of the primary average helium feedwater flowrate occurs, the thermal power transferred
temperature. Since the inlet hot temperature is not changed, from the primary helium flow to the tube-wall is not changed
the decrease in average helium temperature certainly results at the beginning. Then, the decrease in the feedwater flowrate
in the decrease of outlet cold helium temperature. From must lead to the increases of the outlet steam temperature and
10 Mathematical Problems in Engineering
640 30
620
28
Ts5 (∘ C)
l13 (m)
600
26
580
560 24
2000 2200 2400 2000 2200 2400
Time (s) Time (s)
2 240
1.95
235
Tp1 (∘ C)
xM
1.9
230
1.85
1.8 225
2000 2200 2400 2000 2200 2400
Time (s) Time (s)
steam quality and lengthening of the boiling section. Since Finally, from the above discussion, we can easily see that
the specific enthalpy of the feedwater is constant, the increase the steady precision of dynamic model (40) is very high,
of the steam temperature reduces the temperature difference and the transient responses of this model can well cope
between the two sides of the OTSG, which leads to the with the physical trend of the OTSG. Moreover, this model
decrease of the thermal power transferred from the primary can well describe the dynamic behavior of the OTSG in
to the secondary sides, and certainly further results in the case of generating saturate steam. Actually, this model has
increase of the temperature of the primary outlet helium flow. already been adopted to develop the real-time software for
From Figure 6, the numerical phenomenon is in accordance the operation and control features of the HTR-PM plant, as
with the above physical analysis in this case. shown in Figure 8 [16].
In case C, after the occurrence of the large step increase
in the feedwater flowrate, the outlet steam temperature 5. Conclusions
is quickly and largely decreased since the thermal power
transferred from the primary side is nearly not changed at Due to the inherent safety feature of the SMR, SMR-based
the initial stage. The decrease in the steam temperature is nuclear plants are an important developing trend of the
certainly equivalent to the decreases of the boiling section nuclear energy systems. Based upon the multimodular oper-
length and steam quality, and results in a larger temperature ation strategy, SMRs can be used to build nuclear plants
difference between the two sides of the OTSG. This difference with any desired power rating and inherent safety. The
certainly enlarges thermal power transferred to the secondary OTSG is key equipment of any SMR-based multimodular
side, and then results in the temperature decrease of the pri- nuclear plants, and developing the dynamic model for the
mary outlet cold helium. The above analysis well accords with OTSG is very meaningful to study the dynamic behavior
the numerical results in Figure 7. Moreover, from Figure 7, of the multimodular nuclear plants. In this paper, based on
the step increase of the feedwater flowrate is so large that the conservation laws of mass, energy, and momentum, a
the OTSG outlet steam becomes saturated quickly, which differential-algebraic model for the OTSG of those MHTGR-
means that differential-algebraic model (40) can be used for based multimodular nuclear plants is presented. This model
transient simulation in case of generating saturate steam. can describe the dynamic behavior of the OTSG in both
Mathematical Problems in Engineering 11
700 60
600
50
Ts5 (∘ C)
l13 (m)
500
40
400
300 30
2000 2200 2400 2000 2200 2400
Time (s) Time (s)
2 250
1.5 200
Tp1 (∘ C)
xM
1 150
0.5 100
2000 2200 2400 2000 2200 2400
Time (s) Time (s)
Conflict of Interests
The author declares that there is no conflict of interests
regarding the publication of this paper.
References
Figure 8: Real-time simulation software for the HTR-PM operation [1] D. T. Ingersoll, “Deliberately small reactors and the second
by using model (40). nuclear era,” Progress in Nuclear Energy, vol. 51, no. 4-5, pp. 589–
603, 2009.
[2] J. Vujić, R. M. Bergmann, R. Škoda, and M. Miletić, “Small
modular reactors: simpler, safer, cheaper?” Energy, vol. 45, no.
the cases of providing superheated steam and generating 1, pp. 288–295, 2012.
saturated steam. Moreover, since the relative simplicity of this [3] M. D. Carelli, L. E. Conway, L. Oriani et al., “The design and
model, it can also be applied to design steam temperature safety features of the IRIS reactor,” Nuclear Engineering and
control laws for the OTSGs. Numerical simulations results Design, vol. 230, no. 1–3, pp. 151–167, 2004.
show that this model has a satisfactory steady precision, and [4] IAEA, Status of Small and Medium Sized Reactor Designs, Inter-
its transient responses well accord with the thermodynamic national Atomic Energy Agency, Vienna, Austria, 2012.
behavior of the OTSG. This model has already been applied [5] G. H. Lohnert, “Technical design features and essential safety-
to develop a real-time simulation software for the operation related properties of the HTR-module,” Nuclear Engineering and
and control strategy design and verification of the HTR-PM Design, vol. 121, no. 2, pp. 259–275, 1990.
12 Mathematical Problems in Engineering
[6] Z. Wu, D. Lin, and D. Zhong, “The design features of the HTR-
10,” Nuclear Engineering and Design, vol. 218, pp. 25–32, 2002.
[7] Z. Zhang and Y. Sun, “Economic potential of modular reactor
nuclear power plants based on the Chinese HTR-PM project,”
Nuclear Engineering and Design, vol. 237, no. 23, pp. 2265–2274,
2007.
[8] Z. Zhang, Z. Wu, D. Wang et al., “Current status and technical
description of Chinese 2 × 250 MW𝑡ℎ HTR-PM demonstration
plant,” Nuclear Engineering and Design, vol. 239, no. 7, pp. 1212–
1219, 2009.
[9] K. K. Kim and J. A. Bernard, “Considerations in the control of
PWR-type multimodular reactor plants,” IEEE Transactions on
Nuclear Science, vol. 41, no. 6, pp. 2686–2697, 1994.
[10] S. R. P. Perillo, B. R. Upadhyaya, and F. Li, “Control and instru-
mentation strategies for multi-modular integral nuclear reactor
systems,” IEEE Transactions on Nuclear Science, vol. 58, no. 5,
pp. 2442–2451, 2011.
[11] A. Ray, “Dynamic modelling of once-through subcritical steam
generator for solar applications,” Applied Mathematical Mod-
elling, vol. 4, no. 6, pp. 417–423, 1980.
[12] A. Ray, “Nonlinear dynamic model of a solar steam generator,”
Solar Energy, vol. 26, no. 4, pp. 297–306, 1981.
[13] C. P. Tzanos, “Movable boundary model for once-through
steam generator analysis,” Nuclear Technology, vol. 82, no. 1, pp.
5–17, 1988.
[14] M. A. Abdalla, “A four-region, moving-boundary model of a
once-through, helical-coil steam generator,” Annals of Nuclear
Energy, vol. 21, no. 9, pp. 541–562, 1994.
[15] H. Li, X. Huang, and L. Zhang, “A lumped parameter dynamic
model of the helical coiled once-through steam generator with
movable boundaries,” Nuclear Engineering and Design, vol. 238,
no. 7, pp. 1657–1663, 2008.
[16] Z. Dong and X. Huang, “Real-time simulation platform for
the design and verification of the operation strategy of the
HTR-PM,” in Proceedings of the 21st International Conference on
Nuclear Engineering (ICONE ’13), Chengdu, China, August 2013.
Hindawi Publishing Corporation
Mathematical Problems in Engineering
Volume 2015, Article ID 989260, 13 pages
[Link]
Research Article
A Numerical Study of Natural Convection Heat
Transfer in Fin Ribbed Radiator
Copyright © 2015 Hua-Shu Dou et al. This is an open access article distributed under the Creative Commons Attribution License,
which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This paper numerically investigates the thermal flow and heat transfer by natural convection in a cavity fixed with a fin array. The
computational domain consists of both solid (copper) and fluid (air) areas. The finite volume method and the SIMPLE scheme
are used to simulate the steady flow in the domain. Based on the numerical results, the energy gradient function 𝐾 of the energy
gradient theory is calculated. It is observed from contours of the temperature and energy gradient function that the position where
thermal instability takes place correlates well with the region of large 𝐾 values, which demonstrates that the energy gradient method
reveals the physical mechanism of the flow instability. Furthermore, the effects of the fin height, the fin number, and the fin shape on
the heat transfer rate are also investigated. It is found that the thermal performance of the fin array is determined by the combined
effect of the fin space and fin height. It is also observed that the effect of fin shape on heat transfer is insignificant.
with a fixed fin height and fin thickness. They found that numerical simulations to investigate the effect of fin parame-
there was an optimum fin spacing which was not related to ters on dynamic natural convection from long horizontal fin
the temperature difference. However, optimum fin spacing arrays. They observed that the optimum fin spacing decreases
was inversely proportional with the fin height. Mobedi and significantly with the fin height and increases slightly with the
Yüncü [6] numerically investigated the steady state natural fin length. Furthermore, their observations agreed well with
convection heat transfer in longitudinally short rectangular the results reported in the literature.
fin arrays on a horizontal base. They observed two types of All the above experimental and numerical research
flow patterns. For the fin arrays with narrow fin spacing, focused on optimum design of the fin arrays in order to
air could only enter into the channel from the end regions. improve heat transfer rate. Almost all the related factors
However, for the fin arrays with wide fin spacing, the air was result in laminar convection involving low values of heat
also entrained into the channel from the region between the transfer coefficients. However, it is well known that the heat
fins, turned 180 degrees at the base, and then moved up along transfer rate will be prompted when the base flow loses its
the fin surface, while it flowed into the central part of the stability and flows in a turbulent manner [15]. Although
channel. some optimum fin arrays design can obtain a satisfying heat
More recently, lots of researchers tried to enhance heat transfer rate which resulted from flow instability, the physical
transfer rate of fin arrays using fin arrays. Arquis and Rady mechanism of flow instability of natural convection is still
[7] investigated natural convection heat transfer and fluid not fully understood. Recently, Dou et al. [16–24] suggested
flow characteristics from a horizontal fluid layer with finned an energy gradient method which can reasonably reveal the
bottom surface, and observed that the number of convection physical mechanism of flow instability. Dou and Phan-Thien
cells between two adjacent fins is a function of the values of [16] described the rules of fluid material stability from the
the fin height and Rayleigh number. Liu [8] considered an viewpoint of energy field. They claimed that the instability of
optimum design problem for the longitudinal fin arrays with natural convection could not be resolved by Newton’s three
a constant heat transfer coefficient in a fuzzy environment, laws, for the reason that a material system moving in some
where the grid requirements to strictly satisfy the total fin cases is not simply due to the role of forces. This approach
volume and array width and maximize the heat dissipation explains the mechanism of flow instability from physics and
rate are softened. Harahap et al. [9] conducted some exper- derives the criteria of turbulence transition. Accordingly,
iments to investigate the effects of miniaturizing the base this method does not attribute Rayleigh-Benard problem to
plate dimensions of vertically based straight rectangular fin forces, but to energy gradient. It postulates that when the fluid
arrays on the steady state heat dissipation performance under is placed on a horizontal plate and is heated from below, the
dominant natural convection conditions. They found that fluid density in the bottom becomes low which leads to an
the relevant correlations proposed for large fin arrays were energy gradient 𝜕𝐸/𝜕𝑦 > 0 along 𝑦-coordinate. Only when
not applicable to the experimental data obtained from the 𝜕𝐸/𝜕𝑦 is larger than a critical value will the flow become
miniaturized vertical rectangular fin arrays. Subsequently, unstable and then fluid cells of vorticities will be formed.
Harahap et al. [10] conducted concurrent calorimetric and More recently, Dou et al. [25] applied the energy gradient
interferometric measurements to investigate the effect that method to natural convection and the results from numerical
the reduction of the base plate dimensions has on the steady simulations accord well with those predicted based on the
state performance of the rate of natural convection heat criteria originated from energy gradient method.
transfer from miniaturized horizontal single plate-fin systems This study is focused on the research of effects of fin arrays
and plate-fin arrays. Their conclusions suggested that the fin parameters on convection heat transfer coefficient. Then, the
height 𝐿 and the fin number are the prime geometric variables energy gradient method is employed to reveal the physical
for generalization. Dogan and Sivrioglu [11] experimentally mechanism of flow instability and explain the reason why the
investigated mixed convection heat transfer, and the results optimum fin arrays can result in better heat transfer rate.
obtained showed that the optimum fin spacing which yielded
the maximum heat transfer is 𝑆 = 8-9 mm and the optimum 2. Computational Geometry and
fin spacing depends on the value of Ra. Azarkish et al. Numerical Procedures
[12] used a modified genetic algorithm to maximize the
objective function which is defined as the net heat transfer 2.1. Computational Geometry. The computational geometry
rate from the fin surface for a given height. Their results is shown in Figure 1. Here, the geometry is simplified from the
show that the number of the fins is not affected by the fin 3D solid of GH-4 ribbed radiator model [26]. The simplified
profile, but the heat transfer enhancement for the arrays cavity in this study is a two-dimensional (2D) square, in
with optimum fin profile is about 1–3 percent more than which the length of the square cavity is 250 mm, and the
that for the arrays with conventional fin profiles. Giri and origin of the coordinates is at the lower left corner of the
Das [13] numerically performed laminar mixed convection cavity. The fin arrays are fixed at the hot bottom of the cavity
over shrouded vertical rectangular fin arrays attached to a with an equal distance. Here, 𝑆 is the fin space, 𝐻 is the fin
vertical base. They found that the drop in pressure defect for height, and 𝐻 and 𝑆 are variable. In addition, 𝑓𝑡 means the
forced convection, induced velocity for mixed convection and thickness of the fin arrays which is fixed at 2 mm in this study
overall Nusselt number for mixed convection are correlated and can be neglected by comparing to the fin height and
well with governing parameters of the considered problem. the length of the cavity, while 𝑏𝑡 means the thickness of the
Wong and Huang [14] made some three-dimensional (3D) bottom plate which is equal to 3 mm.
Mathematical Problems in Engineering 3
Cold top
are the velocity components in the 𝑥 and 𝑦 directions,
respectively, 𝑔 is the acceleration due to gravity, 𝛽 is the
coefficient of thermal expansion, 𝜌 is the fluid density, 𝑘 is
the thermal diffusivity, and 𝜐 is the kinematic viscosity.
The computational geometry is divided into two areas, Here, 𝜆 is the thermal conductivity, Δ𝑇 is the temperature
which contains a solid area and a fluid area. The solid area difference, and 𝑐𝑝 is the specific heat capacity.
includes the fin arrays and the base plate, in which the
material is copper, while the rest of the geometry is the fluid
area where the air fills in it. All the walls of the cavity are solid 3. Grid Independent Test
and no-slip. The wall style of the left and the right walls of the
cavity is symmetric. The top wall of the cavity is given a fixed In order to examine the independence of the gird size
temperature of 300 K, and the base of the plate is given a fixed to the computing result, three different sized meshes are
temperature of 360 K. Furthermore, the fluid area is given an constructed. Here, grid (a) is uniformly meshed with Δ𝑥 =
initial temperature of 300 K. 0.004 mm, grid (b) is uniformly meshed with Δ𝑥 =
0.002 mm, and grid (c) is uniformly meshed with Δ𝑥 =
2.2. Governing Equations. For the problem described in the 0.001 mm.
previous section, the development of natural convection Figure 2 shows the temperature contours in three dif-
in the fluid area is governed by the following continuity ferent mesh sizes. It can be quantitatively observed that the
equations, two-dimensional steady Navier-Stokes equations temperature contours are almost the same in three different
and an energy equation. The heat transfer in the solid area is
mesh sizes by making a comparison. Figure 3 shows the
governed by a heat conduction equation. All these equations
are based on Boussinesq approximation and listed as follows: temperature versus iterations at the same monitor point
with the three different mesh sizes. It can be found that
𝜕𝑢 𝜕V the temperature difference is getting smaller as the dif-
+ = 0, ference of the mesh size is getting smaller. Furthermore,
𝜕𝑥 𝜕𝑦
for the two finer meshes the difference of the steady-state
𝜕𝑢 𝜕𝑢 𝜕𝑢 1 𝜕𝑝 𝜕 2 𝑢 𝜕2 𝑢 temperature is negligible. In what follows, the medium
+𝑢 +V =− + 𝜐( 2 + 2),
𝜕𝑡 𝜕𝑥 𝜕𝑦 𝜌 𝜕𝑥 𝜕𝑥 𝜕𝑦 (or fine) mesh resolution is adopted for all the calcula-
tions.
𝜕V 𝜕V 𝜕V 1 𝜕𝑝 𝜕2 V 𝜕 2 V
+𝑢 +V =− + 𝜐 ( 2 + 2 ) + 𝑔𝛽 (𝑇 − 𝑇0 ) ,
𝜕𝑡 𝜕𝑥 𝜕𝑦 𝜌 𝜕𝑦 𝜕𝑥 𝜕𝑦
4. Application of Energy Gradient Method
𝜕𝑇 𝜕𝑇 𝜕𝑇 𝜕 2 𝑇 𝜕2 𝑇
+𝑢 +V = 𝑘( 2 + 2 ), 4.1. Energy Gradient Method. It is observed in Figure 2 that
𝜕𝑡 𝜕𝑥 𝜕𝑦 𝜕𝑥 𝜕𝑦
the base flow loses its stability by forming vorticities and
𝜕2 𝑇 𝜕2 𝑇 by moving with a wave form, which is helpful to enhance
+ = 0, heat transfer. However, the physical mechanism of the flow
𝜕𝑥2 𝜕𝑦2
(1) instability of natural convection is still not fully understood.
Energy gradient method will be applied in natural convection
where 𝑥 and 𝑦 are the horizontal and vertical coordinates, 𝑡 to explain the physical mechanism of flow instability in this
is the time, 𝑇 is the temperature, 𝑝 is the pressure, 𝑢 and V study.
4 Mathematical Problems in Engineering
(a) (b)
Grid (c)
(c)
Figure 2: Temperature contours with three different meshes: grid (a) Δ𝑥 = 0.004 mm, grid (b) Δ𝑥 = 0.002 mm, and grid (c) Δ𝑥 = 0.001 mm.
Dou et al. [16–24] suggested an energy gradient method field variable (function) and means the ratio of the transversal
which can be treated as the criteria of flow instability. And the energy gradient and the rate of the energy loss along the
criteria of flow instability can be expressed as follows: streamline. 𝐻𝑠 is the loss of the total mechanical energy per
unit volumetric fluid along the streamline for finite height.
Δ𝐸 (𝜕𝐸/𝜕𝑛) (2𝐴/𝜋) Here, 𝐸 = 𝑝 + 1/2 𝜌 𝑉2 expresses the total mechanical energy
𝐹= = per unit volumetric fluid, 𝑠 is along the streamwise direction,
Δ𝐻 (𝜕𝐻𝑠 /𝜕𝑠) (𝜋/𝜔𝑑 ) 𝑢
(3) and 𝑛 is along the transverse direction.
2 𝐴𝜔𝑑 2 V𝑚
= 𝐾 = 2𝐾 < Const,
𝜋2 𝑢 𝜋 𝑢 4.2. Application of Energy Gradient Method in Natural Con-
vection. In the present study, we will use the energy gradient
where function 𝐾 to reveal the physical mechanism of flow instabil-
𝜕𝐸/𝜕𝑛 ity in natural convection. The fluid fills in a 2D simplified fin
𝐾= . (4) radiator model and the base flow is initially stationary. Based
𝜕𝐻𝑠 /𝜕𝑠
on the energy gradient method and the particular condition
of the base flow, the energy gradient function 𝐾 of natural
Here, 𝑢 is the streamwise velocity of main flow, 𝐴 is the
convection can be expressed as follows:
amplitude of the disturbance distance, 𝜔𝑑 is the frequency of
the disturbance, V𝑚 = 𝐴𝜔𝑑 is the amplitude of the disturbance 𝜕𝐸 2 𝜕𝐸 2
of velocity, and 𝜋 is the circumference ratio. Furthermore, 𝐹 𝐾 = √( ) +( ) . (5)
𝜕𝑥 𝜕𝑦
is a function of coordinates which expresses the ratio of the
energy gained in a half-period by the particle and the energy It should be noted that the influence of the gravitational
loss due to viscosity in the half-period. 𝐾 is a dimensionless energy is neglected in terms of the total mechanical energy in
Mathematical Problems in Engineering 5
320
in the center of the cavity is large is that the heat flux is
focused in these areas and moves downward driven by the
315 buoyancy. However, the heat flux between the adjacent fins
is rare and the buoyancy can be neglected attributed to the
310 small temperature difference between adjacent fins; hence the
velocity between adjacent fins is very weak. Thirdly, there
305 are two blue areas where the velocity is very weak above
the fin arrays. This is due to the convection of heat along
300 the symmetrical wall, and these blue areas are not affected
0 500 1000 1500 2000 2500 3000 3500 4000 evidently by convection.
Iterations
Similarly, the findings from Figure 5(b) can be summa-
Grid (a) rized as follows. Firstly, the total pressure decreases gradually
Grid (b) from the top to the bottom. The fluid flow and the heat
Grid (c) flux along the symmetrical wall move upward driven by the
buoyancy and concentrates near the top wall of the cavity;
Figure 3: Temperature profile calculated with three different grid
thus the total pressure difference will be produced. Secondly,
meshes: grid (a) Δ𝑥 = 0.004 mm, grid (b) Δ𝑥 = 0.002 mm, and grid
two circulating regions occur symmetrically above the fin
(c) Δ𝑥 = 0.001 mm.
arrays. This is due to the circulation flow of the heat flux along
the symmetrical walls. Thirdly, the total pressure gradient is
y irregular above the fin arrays, while the total pressure gradient
between the adjacent fins is regular. The circulation flow of
𝜕E/𝜕y heat flux above the fin arrays which results in the irregular
K distribution of total pressure gradient. The regular gradient
distribution between the fins is attributed to the symmetrical
movement of heat flux along the fin arrays surfaces and the
weak buoyancy between the fins.
Figures 6(a) and 6(b) show the contours of temperature
and the contours of the value of 𝐾 at the iterations of
o
𝜕E/𝜕x x 4000, respectively. It can be seen in Figure 6(a) that the base
flow downstream of the fin loses its stability forming small
Figure 4: Calculation of 𝐾. vorticities and by moving in a wave form. It can be seen in
Figure 6(b) that the value of 𝐾 in the red area and the yellow
area is much higher than that in other areas. In the meantime,
this study since the fluid is air; then 𝐸 ∼ 𝑝0 is obtained; con- there exist two symmetrical green areas with high value of
sequently the energy gradient function 𝐾 of transient natural 𝐾 on the trajectory of the heat flux circulation. Based on
convection can be written as 𝐾 = √(𝜕𝑝0 /𝜕𝑥)2 + (𝜕𝑝0 /𝜕𝑦)2 , the criteria of flow instability of the energy gradient method,
which is shown in Figure 4. Here, 𝑝0 represents the total there is a critical value 𝐾𝑐 , above which the flow will lose its
stability. Thus, the flow instability could most likely occur in
pressure.
the red, yellow, and green areas. Making a further comparison
between Figures 6(a) and 6(b), it is easy to observe that
5. Results and Discussions the positions where instabilities take place are in accordance
with the area with the higher value of 𝐾. This phenomenon
5.1. Physical Mechanism of Flow Instability. Patterson and indicates that the application of the energy gradient method
Imberger [15] announced that the flow instability could in the simplified mode of ribbed radiator is reliable, and the
enhance heat transfer rate of natural convection, and the heat energy gradient method can reveal the physical mechanism
transfer rate increases with the addition of the intensity of of flow instability of natural convection.
flow instability. In the following, the physical mechanism of
flow instability will be discussed using the energy gradient Moreover, it can be found that the regions with high
method. velocity magnitude accord well with the areas with large value
Figures 5(a) and 5(b) show the contours of velocity of 𝐾. Consequently, the flow instability of natural convection
and total pressure respectively with iterations of 4000. The has an instinct affection on the velocity of the fluid flow.
following observations can be made in Figure 5(a). Firstly,
the velocity in the whole flow field distributes symmetrically 5.2. Effect of Fin Height. In order to investigate the effect
along the vertical center line. Secondly, the magnitude of of fin height on heat transfer rate, here two groups of fin
velocity near the left and right walls and in the center of arrays are chosen for comparison. In the first group, all the
6 Mathematical Problems in Engineering
fin spaces are fixed at 61 mm, and the fin height changes from method, which in turn demonstrates that the application of
25 mm to 105 mm with an equal increment of 20 mm. In the energy gradient method to natural convection is reasonable.
second group, all the fin spaces are fixed at 26 mm and the fin Secondly, it is found that the intensity of flow instability
height changes over the same range in the first group. In the changes dramatically in Figure 7 when the fin height changes
following section, the Nusselt number is selected to represent from 45 mm to 65 mm with a fixed fin space of 61 mm. It also
the heat transfer rate. means that the heat transfer rate will be increased fiercely.
Figure 7 shows the temperature contours and contours However, it is found in Figure 8 that the flow instability varies
of the value of 𝐾 with different fin height where the fin in a very small range, which means that the variation of
space is fixed at 61 mm. Figure 8 shows the temperature heat transfer rate is limited. At the same time, it is found
contours and contours of the value of 𝐾 with different fin from the second row of Table 1 that the variation range of
height where the fin space is fixed at 26 mm. Table 1 shows Nusselt number is notable, while the Nusselt number alters
the Nusselt number in different numerical cases. It is easy within a limited extent detected from the third row of Table 1.
to observe the following results by comparing Figures 7 This accordance between flow instability and heat transfer
and 8 and Table 1. Firstly, Figures 7 and 8 show that the rate validates the hypothesis made by Patterson and Imberger
areas where flow instabilities occur in temperature contours [15]. The heat transfer rate increases with the addition of the
accord well with the region where the value of 𝐾 is very intensity of flow instability. Thirdly, it is found in Figure 7 that
large in contours of the value of 𝐾. These results validate the heat flux moves right by jumping over the fins when the
the criteria of flow instability based on the energy gradient fin height is 25 mm and 45 mm, respectively, whereas the flux
Mathematical Problems in Engineering 7
Table 1: Effect of fin height, 𝐻 and 𝑆 mean fin height and fin space respectively.
Nu 𝐻 = 25 mm 𝐻 = 45 mm 𝐻 = 65 mm 𝐻 = 85 mm 𝐻 = 105 mm
𝑆 = 61 mm 1.92 1.97 2.35 2.35 2.34
𝑆 = 26 mm 2.37 2.38 2.38 2.39 2.39
Figure 7: Temperature contours and contours of the value of 𝐾 with different fin heights where fin space is fixed at 61 mm: (a) 𝐻 = 25 mm,
(b) 𝐻 = 45 mm, (c) 𝐻 = 65 mm, (d) 𝐻 = 85 mm, and (e) 𝐻 = 105 mm.
Figure 8: Temperature contours and contours of the value of 𝐾 with different fin heights where fin space is fixed at 26 mm: (a) 𝐻 = 25 mm,
(b) 𝐻 = 45 mm, (c) 𝐻 = 65 mm, (d) 𝐻 = 85 mm, and (e) 𝐻 = 105 mm.
8 Mathematical Problems in Engineering
Table 2: Effect of fin space, 𝐻 and 𝑆 mean fin height and fin space, respectively.
Figure 9: Temperature contours with a fixed fin height of 45 mm and the fin space changes from 82 mm to 23.2 mm characterized by the
increase of fin number (𝑛): (a) 𝑛 = 2, (b) 𝑛 = 3, (c) 𝑛 = 4, (d) 𝑛 = 5, (e) 𝑛 = 6, (f) 𝑛 = 7, (g) 𝑛 = 8, and (h) 𝑛 = 9.
moves upward along the symmetric wall when the fin height second set of numerical cases, the fin height is fixed at
is 65 mm, 85 mm, and 105 mm, respectively. It is observed 85 mm and the fin space decreases characterized by the same
in Figure 8 that the heat flux moves upward symmetrically pattern in the first set. Figures 9 and 10 show the temperature
in all the numerical cases. The difference of flow behavior is contours and contours of the value of 𝐾 for the first set of
jointly affected by the fin height and the fin space. When the numerical cases described above, in which all the fin heights
fin height is small, the heat flux jumps over the fins driven are given as 45 mm. Figures 11 and 12 show the temperature
by the buoyancy. When the fin height is small, the buoyancy contours and contours of the value of 𝐾 for numerical cases,
between the fins is weak, and it is not very easy for the heat in which all the fin heights are fixed at 85 mm. Table 2 lists the
flux to jump over the fins. In addition, when the fin space Nusselt number of different numerical cases.
reduces characterized by the increase of the fin number, the Comparing Figures 9–12 and Table 2, some observations
difficulty for the heat flux jumps over the fins will increase. can be obtained. Firstly, when the fin number is 2, the
Moreover, the temperature difference between adjacent fins strongest intensity of flow instability occurs, which possesses
is very weak which results in a negligible buoyancy between the most efficient heat transfer rate. The heat flux focuses
adjacent fins. Thus, the decrease of the fin space reduces the in the center of the cavity, and it moves upward along the
driven force, which is needed for the heat flux to jump over vertical center line of the cavity. Once the largely extracted
the fins. heat flux loses its stability and spreads randomly, the intensity
It is summarized from the above discussion that the of flow instability is much stronger than any other patterns of
flow behavior is dominated by the combined influences of flow instability. Secondly, Figures 9 and 10 indicate that the
fin space and fin height, which affects the heat transfer intensity of flow instability is very weak when the fin number
rate ultimately. The influences on heat transfer rate can be is 3. However, the intensity difference of flow instability can
explained with the energy gradient theory. be neglected when the fin number is larger than 3. This
result can be validated by the data in Table 2. When the
5.3. Effect of Fin Space. Two sets of numerical cases are fin number is 3, the heat flux moves right by jumping over
selected to study the effect of fin space on flow instability the fins and focuses in the adjacent fins near the right wall.
or heat transfer rate. In the first set of numerical cases, the Hence, little heat flux moves upward, which results in a weak
fin height is fixed at 45 mm, and the fin space decreases intensity of flow instability. However, the heat fluxes in all
characterized by the increasing of the fin number. In the the cavities move upward symmetrically, which results in
Mathematical Problems in Engineering 9
Figure 10: Contours of the value of 𝐾 with a fixed fin height of 45 mm and the fin space changes from 82 mm to 23.2 mm characterized by
the increase of fin number (𝑛): (a) 𝑛 = 2, (b) 𝑛 = 3, (c) 𝑛 = 4, (d) 𝑛 = 5, (e) 𝑛 = 6, (f) 𝑛 = 7, (g) 𝑛 = 8, and (h) 𝑛 = 9.
Figure 11: Temperature contours with a fixed fin height of 85 mm and the fin space changes from 82 mm to 23.2 mm characterized by the
increase of fin number (𝑛): (a) 𝑛 = 2, (b) 𝑛 = 3, (c) 𝑛 = 4, (d) 𝑛 = 5, (e) 𝑛 = 6, (f) 𝑛 = 7, (g) 𝑛 = 8, and (h) 𝑛 = 9.
a negligible heat transfer rate when the fin number is more the fin space is small, the heat conduction between adjacent
than 3. Thirdly, it is easy to find that the intensity of flow fins plays a leading role and it restricts flow instability.
instability changes slightly when the fin number is more than At the same time, more heat flux moves upward which
2 by comparing Figures 11 and 12 and Table 2. This can be results in stronger flow instability. The change of the fin
attributed to the effect of fin space. When the fin space is space results in the change of specific weight of the heat
large, the convection between adjacent fins plays a leading conduction and heat convection between adjacent fins, as
role and it can trigger flow instability. Meanwhile, more heat well as the change of the accumulation of heat flux between
flux focuses between adjacent fins and little heat flux moves adjacent fins. At last, the intensity of flow instability changes
upward; thus the flow instability will be restricted. When slightly.
10 Mathematical Problems in Engineering
Figure 12: Contours of the value of 𝐾 with a fixed fin height of 85 mm and the fin space changes from 82 mm to 23.2 mm characterized by
the increase of fin number (𝑛): (a) 𝑛 = 2, (b) 𝑛 = 3, (c) 𝑛 = 4, (d) 𝑛 = 5, (e) 𝑛 = 6, (f) 𝑛 = 7, (g) 𝑛 = 8, and (h) 𝑛 = 9.
Figure 13: Temperature contours and contours of the value of 𝐾 in different fin shape with a fixed fin space of 82 mm.
Consequently, in addition to the validation of criteria of straight fin arrays are chosen as the base group as model-(A)
flow instability based on the energy gradient method, another shown in Figure 13. The shapes of models-(B–D) shown in
two conclusions can be summarized as follows. (1) The cavity Figure 13 are all different from each other and a fixed fin space
possesses the most efficient heat transfer rate when the fin of 26 mm is chosen for all cases. The fin arrays in model-(B)
number is 2. (2) The change of the fin space results in the are curved for one piece in the same direction with a radius of
change of the specific weight of heat conduction and heat 70 mm. The fin arrays in model-(C) are curved for one piece
convection as well as the accumulation of heat flux between with a radius of 70 mm, which are fixed symmetrically. The fin
adjacent fins; hence the heat transfer rate changes slightly. arrays in model-(D) are curved for four pieces with a radius
of 10 mm.
5.4. Effect of Fin Shape. In this part, the effect of the fin Figures 13–15 show the temperature contours and con-
shape on flow instability or heat transfer rate is discussed. The tours of the value of 𝐾 when the fin space is 82 mm, 48.4 mm,
Mathematical Problems in Engineering 11
Figure 14: Temperature contours and contours of the value of 𝐾 in different fin shape with a fixed fin space of 48.4 mm.
Figure 15: Temperature contours and contours of the value of 𝐾 in different fin shape with a fixed fin space of 29.5 mm.
and 29.5 mm, respectively. Figure 16 shows the heat transfer in the center of the cavity and spreads around and produces
rate obtained with four different fin shapes, and fin space is a large intensity of flow instability. However, the heat fluxes
characterized by the fin number. in models-(B–D) move upward symmetrically and produce
It is easy to get the following results by comparing Figures relatively small flow instability. Thus, this result is validated
13–15 to Figure 16. Firstly, it is found in Figure 13 that the in Figure 16, which shows that the most efficient heat transfer
largest intensity of flow instability occurs in model-(A), and rate occurs in model-(A) when the fin space is 82 mm.
the intensity difference in other models can be neglected. Secondly, it is found in Figure 14 that the flow instability
These results accord well with the data in Figure 16. This in model-(C) is very weak, and the intensity variation of
is attributed to the effect of fin shape, which results in the flow instability in other models is very small. This heat
different form of flow instability. Convection occurs mainly transfer rate difference is due to the bent direction of the fins.
12 Mathematical Problems in Engineering
(2) The fin height and the fin space jointly dominate the
2.40 flow behavior of the base flow and thus affect the heat
transfer rate.
2.35
(3) For a given fin height, the change of the fin space
2.30 results in the change of the accumulation and move-
ment of the heat flux between adjacent fins, as well
2.25 as the relative magnitude of heat conduction and
heat convection between adjacent fins; hence the heat
2 3 4 5 6 7 8 9 transfer rate changes slightly.
Fin number
(4) The effect of the fin shape on heat transfer is distinct
Model-(A) Model-(C) when the fin space is large.
Model-(B) Model-(D)
Figure 16: Nusselt number versus fin space characterized with fin
Conflict of Interests
number in different fin shape. The authors declare that there is no conflict of interests
regarding the publication of this paper.
[8] F. Liu, “A fuzzy approach to the convective longitudinal fin array [25] H. Dou, G. Jiang, and C. Lei, “Numerical simulation and
design,” International Journal of Thermal Sciences, vol. 44, no. 3, stability study of natural convection in an inclined rectangular
pp. 211–217, 2005. cavity,” Mathematical Problems in Engineering, vol. 2013, Article
[9] F. Harahap, H. Lesmana, and A. S. Dirgayasa, “Measurements ID 198695, 12 pages, 2013.
of heat dissipation from miniaturized vertical rectangular fin [26] Y. Gao, A heat transfer study on the horizontal rectangular fin
arrays under dominant natural convection conditions,” Heat arrays for the LED street lamp GH-4 [M.S. thesis], Chongqing
and Mass Transfer, vol. 42, no. 11, pp. 1025–1036, 2006. University, 2011.
[10] F. Harahap, H. Lesmana, and P. L. Sambegoro, “Concurrent
calorimetric and interferometric studies of steady-state natural
convection from miniaturized horizontal single plate-fin sys-
tems and plate-fin arrays,” Heat and Mass Transfer, vol. 46, no.
8-9, pp. 929–942, 2010.
[11] M. Dogan and M. Sivrioglu, “Experimental investigation of
mixed convection heat transfer from longitudinal fins in a
horizontal rectangular channel,” International Journal of Heat
and Mass Transfer, vol. 53, no. 9-10, pp. 2149–2158, 2010.
[12] H. Azarkish, S. M. H. Sarvari, and A. Behzadmehr, “Optimum
design of a longitudinal fin array with convection and radiation
heat transfer using a genetic algorithm,” International Journal of
Thermal Sciences, vol. 49, no. 11, pp. 2222–2229, 2010.
[13] A. Giri and B. Das, “A numerical study of entry region
laminar mixed convection over shrouded vertical fin arrays,”
International Journal of Thermal Sciences, vol. 60, pp. 212–224,
2012.
[14] S. Wong and G. Huang, “Parametric study on the dynamic
behavior of natural convection from horizontal rectangular fin
arrays,” International Journal of Heat and Mass Transfer, vol. 60,
no. 1, pp. 334–342, 2013.
[15] J. Patterson and J. Imberger, “Unsteady natural convection in a
rectangular cavity,” Journal of Fluid Mechanics, vol. 100, no. 1,
pp. 65–86, 1980.
[16] H. S. Dou and N. Phan-Thien, “Instability of fluid material
systems,” in Proceedings of the 15th Australasian Fluid Mechan-
ics Conference, The University of Sydney, Sydney, Australia,
December 2004.
[17] H. S. Dou, “Viscous instability of inflectional velocity profile,”
in Proceedings of the 4th International Conference on Fluid
Mechanics, Tsinghua University Press & Springer, 2004.
[18] H. Dou, “Mechanism of flow instability and transition to
turbulence,” International Journal of Non-Linear Mechanics, vol.
41, no. 4, pp. 512–517, 2006.
[19] H. S. Dou, “Three important theorems for flow stability,”
in Proceedings of the 5th International Conference on Fluid
Mechanics, Tsinghua University Press & Springer, 2007.
[20] H. Dou and B. C. Khoo, “Mechanism of wall turbulence in
boundary layer flow,” Modern Physics Letters B, vol. 23, no. 3,
pp. 457–460, 2009.
[21] H. Dou and B. C. Khoo, “Criteria of turbulent transition in
parallel flows,” Modern Physics Letters B, vol. 24, no. 13, pp. 1437–
1440, 2010.
[22] H. Dou and B. C. Khoo, “Investigation of turbulent transition in
plane couette flows using energy gradient method,” Advances in
Applied Mathematics and Mechanics, vol. 3, no. 2, pp. 165–180,
2011.
[23] H. Dou, B. C. Khoo, and K. S. Yeo, “Energy loss distribution
in the plane Couette flow and the Taylor-Couette flow between
concentric rotating cylinders,” International Journal of Thermal
Sciences, vol. 46, no. 3, pp. 262–275, 2007.
[24] H. Dou, B. C. Khoo, and K. S. Yeo, “Instability of Taylor-Couette
flow between concentric rotating cylinders,” International Jour-
nal of Thermal Sciences, vol. 47, no. 11, pp. 1422–1435, 2008.
Hindawi Publishing Corporation
Mathematical Problems in Engineering
Volume 2014, Article ID 494585, 13 pages
[Link]
Research Article
Numerical Study of Buoyancy Convection of Air under
Permanent Magnetic Field and Comparison with That under
Gravity Field
Received 9 May 2014; Revised 12 August 2014; Accepted 14 August 2014; Published 3 September 2014
Copyright © 2014 Kewei Song et al. This is an open access article distributed under the Creative Commons Attribution License,
which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Magnetothermal free convection of air in a square enclosure under a nonuniform magnetic field provided by a permanent
neodymium-iron-boron magnet is numerically studied. The natural convection under the gravity field alone is also studied for
comparison. The physical fields of magnetizing force, velocity, and temperature as well as the local distribution characteristic of
Nusselt number are all presented in this paper. The results show that the buoyancy convection of air in the square enclosure under
magnetic field is quite different from that under the gravity field. The local value of Nusselt number under the magnetic field supplied
by a permanent magnet with a residual magnetic flux density of about 4.5 Tesla can reach a high value of about three times larger
than the maximum local value of Nusselt number under the gravity field. Relatively uniform distributions of temperature gradient
and Nusselt number can be obtained along the cold wall of the enclosure under the magnetic field. A permanent magnet with high
magnetic energy product with 𝐵𝑟 reaching to 3.5 Tesla can play a comparative role on the averaged Nusselt number compared with
that under the gravity environment.
Fx− Fx+
Cold
Hot
z
Fz−
Cold
Hot
x
g
z
S
x
Magnet
(a) (b)
Figure 1: Configuration of simulation models. (a) Permanent magnetic field and (b) gravity field.
The momentum equation including the magnetizing Magnetic susceptibility of air is a function of temperature;
force alone is as follows: this follows the Curie law:
𝜕u 1 𝐶
𝜌[ + (u ⋅ ∇) u] = −∇𝑝 + 𝜇∇2 u + 𝜌𝜒∇B2 . (2) 𝜒= . (10)
𝜕𝑡 2𝜇𝑚 𝑇
And the momentum equation including the gravity buoy- Here, 𝐶 is a constant; then
ancy force alone is 𝜕𝜒 𝜕 𝐶 𝜒
= ( )=− . (11)
𝜕u 𝜕𝑇 𝜕𝑇 𝑇 𝑇
𝜌[ + (u ⋅ ∇) u] = −∇𝑝 + 𝜇∇2 u − 𝜌𝑔. (3)
𝜕𝑡 The volumetric coefficient of expansion for an ideal gas
can be given from the ideal gas law:
We assume that the air is an incompressible Newtonian
fluid and the Boussinesq approximation is employed. When 𝑝 = 𝜌𝑅𝑔 𝑇. (12)
we consider the isothermal state, the fluid density 𝜌 and
susceptibility 𝜒 would be constant, and convection does not Then, from (7)
arise in both the gravity and magnetic fields. Parameters 1 𝜕𝜌 1 𝜕 𝑝 𝑝 1
under this state are 𝜌0 , 𝑝0 , and 𝜒0 at reference temperature 𝛽=− =− ( )= 2
= ,
𝜌 𝜕𝑇 𝜌 𝜕𝑇 𝑅𝑔 𝑇 𝜌𝑅𝑔 𝑇 𝑇
𝑇0 . Hence we get (13)
1 𝜕𝜌
0 = −∇𝑝0 + 𝜌 𝜒 ∇B2 = −𝜌𝛽
2𝜇𝑚 0 0 (4)
𝜕𝑇
so that
0 = −∇𝑝0 − 𝜌0 𝑔.
𝜌 𝜒
When there is a temperature difference, the magnetic 𝜌𝜒 − (𝜌𝜒)0 = (−𝜒 − 𝜌 ) (𝑇 − 𝑇0 )
𝑇 𝑇 0 (14)
susceptibility and density also change with temperature.
Subtracting (4) from (2) and (3) gives = −2𝜌0 𝛽0 𝜒0 (𝑇 − 𝑇0 ) .
𝜕u Given 𝑝 = 𝑝0 + 𝑝 , the momentum equation (5) can be
𝜌[ + (u ⋅ ∇) u]
𝜕𝑡 written as follows:
𝜌𝜒 − 𝜌0 𝜒0 2 𝜕u
= −∇ (𝑝 − 𝑝0 ) + 𝜇∇2 u + ∇B [ + (u ⋅ ∇) u]
2𝜇𝑚 𝜕𝑡
𝜕u 1 𝜇 𝛽 𝜒 (𝑇 − 𝑇0 ) 2
𝜌[ + (u ⋅ ∇) u] =− ∇𝑝 + ∇2 u − 0 0 ∇B
𝜕𝑡 𝜌0 𝜌0 𝜇𝑚
(15)
2
= −∇ (𝑝 − 𝑝0 ) + 𝜇∇ u − (𝜌 − 𝜌0 ) 𝑔. 𝜕u
[ + (u ⋅ ∇) u]
(5) 𝜕𝑡
1 𝜇
The density can be expressed by a Taylor expansion =− ∇𝑝 + ∇2 u + 𝑔𝛽0 (𝑇 − 𝑇0 ) .
(keeping two terms only) around a reference state 𝑇0 as 𝜌0 𝜌0
follows: Thus, the government equations for the natural con-
𝜕𝜌 vection caused by the magnetic buoyancy force or gravity
𝜌 = 𝜌0 + ( ) (𝑇 − 𝑇0 ) + ⋅ ⋅ ⋅ . (6) buoyancy force are summarized as follows.
𝜕𝑇 0
Continuity equation is
On the other hand, the thermal volume coefficient of
expansion 𝛽 is defined for fluid as follows [23]: ∇ ⋅ u = 0. (16)
1 𝜕𝑉 1 𝜕1/𝜌 1 𝜕𝜌 Momentum equation, for magnetic field alone is
𝛽= ( ) = ( ) =− ( ) . (7)
𝑉 𝜕𝑇 𝑝 1/𝜌 𝜕𝑇 𝑝 𝜌 𝜕𝑇 𝑝 𝜕u
[ + (u ⋅ ∇) u]
Then, (6) can be written as follows: 𝜕𝑡
(17)
𝜕𝜌 1 𝜇 𝛽 𝜒 (𝑇 − 𝑇0 ) 2
𝜌 − 𝜌0 = ( ) (𝑇 − 𝑇0 ) = −𝜌0 𝛽0 (𝑇 − 𝑇0 ) . (8) = − ∇𝑝 + ∇2 u − 0 0 ∇B .
𝜕𝑇 0 𝜌0 𝜌0 𝜇𝑚
The difference of 𝜌𝜒 can be written as And for gravity field alone
𝜕𝜌𝜒 𝜕u
𝜌𝜒 − (𝜌𝜒)0 = ( ) (𝑇 − 𝑇0 ) [ + (u ⋅ ∇) u]
𝜕𝑇 0 𝜕𝑡
(9) (18)
𝜕𝜌 𝜕𝜒 1 𝜇
= (𝜒 + 𝜌 ) (𝑇 − 𝑇0 ) . =− ∇𝑝 + ∇2 u + 𝑔𝛽0 (𝑇 − 𝑇0 ) .
𝜕𝑇 𝜕𝑇 0 𝜌0 𝜌0
4 Mathematical Problems in Engineering
Grid system Nu
62 × 62 × 62 4.540
82 × 82 × 82 4.523
102 × 102 × 102 4.520
4. Results and Discussion 4.2. Magnetizing Force Fields. Figure 4 shows the gravity
4 buoyancy force and magnetizing buoyancy force fields in
All the computations are carried out from Ra = 10 to Ra =
106 , and the residual flux densities of permanent magnet the enclosure. From these figures we can find that the
changes from 0.5 T to 4.5 T. The flow field, temperature field, magnetizing buoyancy force is different from the gravity
and the distribution of local Nusselt number along the walls buoyancy force. The gravity buoyancy force parallels along
are obtained under the magnetic field with 𝐵𝑟 = 4.5 T; the 𝑧-direction, but the magnetizing buoyancy force’s direction
difference of the effect of magnetic field and gravity field on changes in the enclosure. The magnetizing buoyancy force
natural convection is also studied. increased with increase of magnet strength. There is large
magnetizing force near the bottom wall of the enclosure
and the force near the hot wall is larger than that near the
4.1. Temperature Fields. Figure 3 shows the temperature fields cold wall. The magnetizing force decreases quickly along
of the studied physical models at Ra = 104 and 105 . For 𝑧-direction. But for the gravity field, the distribution of
6 Mathematical Problems in Engineering
300.0 300.2 300.3 300.5 300.7 300 301 303 304 305 306 308
Ra = 104 Ra = 105
(a)
300.0 300.2 300.3 300.5 300.7 300 301 303 304 305 306 308
Ra = 104 Ra = 105
(b)
Figure 3: Temperature fields at Ra = 104 and 105 . (a) Under gravity field and (b) under magnetic field, 𝐵𝑟 = 4.5 T.
gravity buoyancy force is quite different. The larger gravity vortex in the center of the enclosure at Ra = 104 . When
buoyancy force locates in the region near the top and hot Ra increases to 105 , the core of the vortex is divided into
wall, and the gravity buoyancy force decreases from hot wall two cores, one is located near the top and hot walls, and the
to cold wall and increases from bottom wall to top wall. The other one locates close to the cold and bottom walls. When
maximum value of magnetizing buoyancy force is larger than there is magnetic field only, the velocity field is different
the maximum gravity buoyancy force. from that under the gravity field. A vortex is also formed
in the enclosure, but the core of the vortex locates close to
4.3. Velocity Fields. Figure 5 shows the velocity fields under the bottom wall at Ra = 104 . When Ra increases to 105 ,
gravity field and magnetic field. When there is gravity field there is still only one core of the vortex and the core moves
only, as shown in Figure 5(a), the large value of velocity much closer to the hot wall. Figures 6(a) and 6(b) show the
locates near the walls and the value of velocity decreases streamlines in the enclosure under gravity field and magnetic
from the walls to the center of the enclosure. There is one field, respectively. The streamlines related to gravity field look
Mathematical Problems in Engineering 7
Ra = 104 Ra = 105
0.1 N 0.5 N
(a)
Ra = 104 Ra = 105
0.1 N 0.5 N
(b)
Figure 4: Gravity and magnetizing buoyancy force fields. (a) 𝑔 and (b) 𝐵𝑟 = 4.5 T.
rotationally symmetric about the center of the enclosure; z3 are at the positions 𝑧/𝐿 = 0.1, 0.5, and 0.8, respectively.
large velocity gradient normal to the wall locates near the In the region near the hot wall, the temperature under the
bottom segment of the hot wall and the top segment of magnetic field along the line z1 is lower than that under the
the cold wall. But under the magnetic field, the streamlines gravity field, but along the lines z2 and z3, the temperature
are not symmetric and there are dense streamlines near the under the magnetic field is higher than the temperature
corner of the hot and bottom walls where there is large value under the gravity field. In the region near the cold wall,
of magnetizing force as shown in Figure 4(b). In the region the temperature related to magnetic field is higher than the
near the top wall, the velocity gradient normal to the wall is temperature related to gravity field along the lines z1 and z2,
significantly less than that in other regions. and the temperature related to magnetic field is lower than the
temperature related to the gravity field along line z3. The high
temperature near the cold wall means a high temperature
4.4. Distribution of Temperature along Lines. The distribution gradient and high temperature near the hot wall means a
curves of temperature at Ra = 105 along the lines parallel low temperature gradient. For the magnetic field, the largest
to 𝑥-axis are plotted in Figure 7; these three lines z1, z2, and temperature gradient near the cold wall is along the line z2
8 Mathematical Problems in Engineering
Ra = 104 Ra = 105
Ra = 104 Ra = 105
Figure 5: Comparison of velocity vectors of physical models. (a) Under gravity field and (b) under magnetic field.
and the largest temperature gradient near the hot wall is along Nusselt number along the hot wall. The value of the Nusselt
the line z1. We can obtain the same conclusion as in Figure 7 number increases nearly linearly from bottom wall and gets
that the distribution of the temperature under the magnetic the peak value near the top wall; then it decreases again till
field is quite different from that under the gravity field. to the top wall. Figure 8(b) shows the distribution curves of
local Nusselt number under the magnetic field with 𝐵𝑟 =
4.5. Heat Transfer on the Walls. The distribution curves of 4.5 T. We can see that the distribution of the local Nusselt
the local Nusselt number along the hot and cold walls are number for the model under magnetic field is different from
plotted in Figure 8. The distribution of the local Nusselt that shown in Figure 8(a). The local Nusselt number gets a
number under gravity field is shown in Figure 8(a). The peak high peak value at the bottom of the hot wall and it decreases
value along the hot wall locates near the bottom wall and quickly for a short distance and then continues decreasing at
decreases nearly linearly to the top wall. The distribution of a relatively slow rate till it reaches the top wall. The change of
local Nusselt number along the cold wall is symmetric about the local Nusselt number along the cold wall is much gentler
the central line of the enclosure with the distribution of local compared with that along the hot wall. The local Nusselt
Mathematical Problems in Engineering 9
Ra = 104 Ra = 105
(a)
Ra = 104 Ra = 105
(b)
Figure 6: Comparison of streamlines. (a) Under gravity field and (b) under magnetic field.
number increases slowly from the top wall to the position magnetic flux density of a permanent magnet ranging from
near the bottom wall and then it decreases quickly till it 1.5 T to 4.5 T. The average Nu under the magnetic field has
reaches the bottom wall. a similar distribution curve compared with that under the
The comparison of the distribution curves of the local gravity field. The averaged Nu increases with increase of 𝐵𝑟 .
Nusselt number under gravity field and magnetic field are The averaged Nusselt numbers on the hot wall with Ra =
presented in Figure 9. The value of local Nusselt number 105 are summarized in Table 3. When 𝐵𝑟 is about 3.5T,
along the hot wall is larger than the value under the gravity the averaged Nu is nearly the same as that under gravity
filed in a short line segment from the bottom wall; the peak field; this means that the magnetizing buoyancy force has a
value of local Nusselt number related to magnetic field is comparative effect on Nu in the enclosure compared with the
more than three times larger than the peak value related to gravity buoyancy force.
the gravity field. On the cold wall, the local Nusselt number
related to magnetic field changes much more gently and has 5. Conclusion
larger value on the lower part of the wall compared with that
related to gravity field. Natural convection of air in a two-dimensional square enclo-
Figure 10 shows the distribution of the average Nusselt sure under a nonuniform magnetic field provided by a per-
number with Ra ranging from 104 to 106 and the residual manent magnet is carried out. The results are compared with
10 Mathematical Problems in Engineering
Nomenclature
306 A: Vector potential
B: Magnetic flux density (T)
T 𝐵𝑟 : Residual magnetic flux density of a
304
permanent magnet (T)
f𝑚 : Magnetizing force (N/m3 )
𝑔: Acceleration coefficient of gravity (m/s2 )
H: Magnetic field strength (A/m)
302 M: Magnetization
Nu: Nusselt number
𝑝: Pressure (N/m2 )
300 𝑝0 : Reference pressure without convection of
0 0.2 0.4 0.6 0.8 1 gas (N/m2 )
x (L)
𝑝 : Perturbed pressure due to convection
z1 Gravity field Pr: Prandtl number
z2 Magnetic field r: Unit vector
z3 Ra: Rayleigh number
𝑡: Time (s)
Figure 7: Temperature distribution along three lines parallel to 𝑥- 𝑇: Temperature (K)
axis.
u: Velocity vector (𝑢, 𝑤)(m/s)
𝑥: Coordinate (m)
Table 3: Summary of computed results. 𝑧: Coordinate (m).
Model 𝐵𝑟 Nu
0.5 T 1.289 Greek Symbols
1.5 T 2.632
Magnetic field 2.5 T 3.617
𝛼: Thermal diffusivity of gas (m2 /s)
3.5 T 4.434
𝛽: Volume coefficient of expansion of gas due
4.5 T 5.149
to temperature difference (1/K)
Gravity field 𝑔 4.523 𝜇: Viscosity of gas (Pa⋅s)
𝜇𝑚 : Magnetic permeability of free space (H/m)
]: Kinematic viscosity of gas (m2 /s)
that under the gravity field. As the magnetizing buoyancy 𝜌: Density of gas (kg/m3 )
force of air in the magnetic field is quite different from the 𝜌0 : Density of gas at a reference state of no
gravitational buoyancy force, the natural convection in the convection (kg/m3 )
enclosure under a magnetic field provided by a permanent 𝜒: Mass magnetic susceptibility of air
magnet has quite different characteristics. For the model (m3 /kg)
studied in this paper, the velocity and temperature fields are 𝜆: Conduction coefficient (W/(m⋅K)).
rotationally symmetric about the center of the enclosure. The
heat transfer Nusselt numbers on the hot wall and the cold
wall are symmetrical about the center line of the enclosure. Subscripts
But the velocity and temperature fields under the magnetic
field do not have such characteristics. Heat transfer Nusselt
number is different on the hot and cold walls. The value of 0: Reference state
local Nusselt number has a peak value on the bottom of the 𝐶: That of a cold wall
hot wall and decreases quickly along the hot wall, but the local ℎ: That of a hot wall.
Nusselt number changes with a relatively gentle rate on the
cold wall. The peak value on the hot wall is three times larger
than that under the gravity field when 𝐵𝑟 = 4.5 T and 2 times Conflict of Interests
larger when 𝐵𝑟 = 2.5 T. When 𝐵𝑟 is about 3.5 T, the value of
averaged Nu is nearly the same with that under gravity field. The authors declare that there is no conflict of interests
This means that in the zero gravity environment, a permanent regarding the publication of this paper.
Mathematical Problems in Engineering 11
20 70
56
15
42
Nu local
Nu local
10
28
5
14
0 0
0 0.25 0.5 0.75 1 0 0.25 0.5 0.75 1
z (L) z (L)
Figure 8: Distribution curves of Nulocal along the walls. (a) Under gravity field and (b) under magnetic field.
70 21
56
10
14
42
Nu local-h
Nu local-c
28
7
0
0.25 0.5 0.75 1
14
0 0
0 0.25 0.5 0.75 1 0 0.25 0.5 0.75 1
z (L) z (L)
Figure 9: Comparison of Nulocal along the hot wall and cold wall. (a) Along hot wall and (b) along cold wall.
12 Mathematical Problems in Engineering
Research Article
System Model of Heat and Mass Transfer Process for Mobile
Solvent Vapor Phase Drying Equipment
Copyright © 2014 Shiwei Zhang et al. This is an open access article distributed under the Creative Commons Attribution License,
which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
The solvent vapor phase drying process is one of the most important processes during the production and maintenance for large oil-
immersed power transformer. In this paper, the working principle, system composition, and technological process of mobile solvent
vapor phase drying (MVPD) equipment for transformer are introduced in detail. On the basis of necessary simplification and
assumption for MVPD equipment and process, a heat and mass transfer mathematical model including 40 mathematical equations
is established, which represents completely thermodynamics laws of phase change and transport process of solvent, water, and air
in MVPD technological processes and describes in detail the quantitative relationship among important physical quantities such
as temperature, pressure, and flux in key equipment units and process. Taking a practical field drying process of 500 KV/750 MVA
power transformer as an example, the simulation calculation of a complete technological process is carried out by programming with
MATLAB software and some relation curves of key process parameters changing with time are obtained such as body temperature,
tank pressure, and water yield. The change trend of theoretical simulation results is very consistent with the actual production
record data which verifies the correctness of mathematical model established.
Figure 2: MVPD8.2 system principle diagram of MVPD equipment with two evaporators.
MVPD equipment is made up of the following cell modules vapor enters into the transformer tank, condenses on the cold
[6]: surface of transformer body, and releases the latent heat of
phase change. The condensed solvent liquid is discharged
(a) main big module 1: vacuum unit and condensation from the bottom of transformer tank and is then delivered
module; back to the evaporator to reheat and vaporize. In this case,
(b) main big module 2: evaporator and electric heater the transformer body is heated by the circular solvent vapor,
module; so the body temperature is gradually getting higher. With the
(c) small equipment module: coarse filter and solvent temperature increasing, moisture in insulation material starts
circulating pump module; to evaporate and diffuse gradually.
At the heating stage, the transformer tank is consistently
(d) auxiliary equipment module 1: solvent tank and evacuated by the leakage pump in order to ensure pressure
wasted oil tank module; difference between the transformer tank and the evaporator.
(e) auxiliary equipment module 2: water tank and water The evacuated gas in transformer tank is the mixture includ-
chiller unit, air compressor, and air tank; ing uncondensed solvent vapor, released water vapor, and
(f) auxiliary component module: box loading connection leaking air. The mixed gas enters into the main condenser.
pipes (two wide pipes and three thin tubes); Solvent vapor and water vapor condense into liquids and
then flow into the collecting tank. Due to the different
(g) thermal insulation blanket and tank heating device density, solvent and water are naturally separated by gravity
module. settling. Water deposits at the bottom of collecting tank and
is discharged while the solvent at upper layer is collected and
2.3. Working Processes of MVPD Equipment. According to delivered to evaporator, thus forming another cycle.
the time order, there are mainly the following processes when The temperature of solvent vapor should be controlled at
the MVPD equipment works [1, 14, 15]. 130∘ C to 135∘ C which is determined according to the max-
imum allowable temperature of class A insulation materials
(1) Preparation Stage. At first, condensation system and under the anaerobic condition. In the several times of heating
evaporator are evacuated by the leakage pump. Then enough stage, about 90% of the water content in insulation material
solvent liquid is injected into the evaporator from the sol- may be drawn out.
vent storage tank. The vacuum system is started and the
transformer tank is evacuated until its pressure drops below (3) Intermediate Pressure Lower Stage. The condensation of
700 Pa. In order to save time, the evaporator is preheated at solvent vapor will be restrained after the insulation material
the preparation stage and its temperature is controlled at 80∘ C absorbing enough solvent, which will hinder both subsequent
to 90∘ C. solvent vapor from condensation and internal moisture
from migrating outwards. In order to increase the solvent
(2) Heating Stage. Close the main vacuum valve and stop the condensation and the moisture migration, it is needed to stop
vacuum system. Open the solvent vapor valve. Solvent in the heating stage and enter the intermediate pressure lower (IPL)
evaporator is heated and changed into solvent vapor. Solvent stage.
4 Mathematical Problems in Engineering
Table 1: Symbol definition in the simplified calculation process model of MVPD system.
Process units Process logistics
(1) Liquid solvent enters into the evaporator
(A) Solvent storage tank
(2) Solvent vapor enters into transformer tank and starts to heat transformer
(B) Evaporator (3) Mixed gas enters into condensation system and condenses
(C) Condenser (4) Leakage pump system evacuates the transformer tank through condensation system
(D) Transformer (5) Solvent condensed by condensation system returns evaporator and is reused
(E) High vacuum system (6) Condensed liquid solvent flows into evaporator from the bottom of transformer tank and is reused
(F) Leakage pump system (7) High vacuum system vacuumizes the transformer tank to drain out condensed water
Table 2: Definition of physical quantity symbols and units. The temperature of liquid solvent is 𝑇𝐵𝑆 :
𝑀 Mass (kg) 𝑑𝑇𝐵𝑆
𝑄 Heat (J) 𝑐𝑂 ⋅ 𝑀𝐵𝑆 = 𝑝𝐵 − 𝑔𝐵𝑂 ⋅ ℎ𝐵𝑆 . (4)
𝑑𝜏
𝑇 Temperature (K)
𝑞 Heat flow (w) The temperature of solvent vapor is 𝑇𝐵𝑂:
𝑚 Mass flow (kg/s)
𝑇𝐵𝑂 = 𝑇𝐵𝑆 . (5)
𝑔 Phase change rate (kg/s)
ℎ Phase change latent heat (J/kg) According to the test data provided by a manufacture
𝜏 Time (s) factory, the saturated vapor pressure of liquid solvent can be
𝐿 Length (m) computed by the following equation:
𝑟 Radius (m)
𝑆 Area (m2 ) 2153.6
lg𝑃𝐵𝑆 = 10.1 − . (6)
𝑐 Specific heat capacity (J/kg⋅K) 𝑇𝐵𝑆 − 34.3
𝑝 Power (W)
The pressure of solvent vapor is 𝑃𝐵𝑂 and it may be com-
𝑃 Pressure (Pa) puted by the perfect gas state equation under low pressure:
𝜌 Density (kg/m3 )
𝜂 Viscosity coefficient (Pa⋅s) 𝑀𝐵𝑂𝑅𝑇𝐵𝑂
𝑃𝐵𝑂 = . (7)
𝜕 Boiling heat transfer coefficient (w/m2 ⋅K⋅Pa) 𝑉𝐵 𝜇𝐵𝑂
𝛽 Condensation surface coefficient (w/m2 ⋅K⋅Pa)
𝜆 Convective heat transfer coefficient (w/m⋅K) 3.2.2. Unit Model of Transformer. The substances composi-
𝜇 Molar mass (kg/mol) tion concerned in the heat and mass transfer process in the
𝐶 Constant transformer 𝐷 consists of solvent, water, and air. There are
𝑘 Coefficient two phase states including liquid state and vapor state for sol-
vent and water. The solvent vapor flows into the transformer
Table 3: Definition of the first and second subscripts of each physi- tank through the logistics pipeline 2 and then condenses and
cal quantity. turns into liquid solvent. The heat released by phase change
ABCDEF Process unit label
will heat the transformer tank, iron cores, copper wires, and
The first subscript
insulation material. The mass conversion and heat exchange
123456789 Process logistics label
in this unit can be expressed by the following mathematical
I II III IV V Time phase
equations.
𝑆 Liquid solvent
𝑂 Solvent vapor (1) Air. The mass of air is 𝑀𝐷𝐴 :
𝑊 Liquid water
𝑉 Water vapor 𝑑𝑀𝐷𝐴
= 𝑚𝐴 − 𝑚3𝐴 , (8)
The second subscript 𝐴 Air 𝑑𝜏
𝑗 Insulation materials where 𝑚𝐴 is the amount of leakage air.
𝑖 Iron core The air temperature is 𝑇𝐷𝐴 and it is equal to the tempera-
𝑏 Transformer tank ture of solvent vapor 𝑇𝐷𝑂 in unit 𝐷:
𝑐 Copper wire
𝑇𝐷𝐴 = 𝑇𝐷𝑂. (9)
The mass of liquid solvent in 𝐵 is 𝑀𝐵𝑆 : The partial pressure of air in the transformer 𝐷 is 𝑃𝐷𝐴 :
𝑑𝑀𝐵𝑆 𝑀𝐷𝐴 𝑅𝑇𝐷𝐴
= 𝑚1𝑆 + 𝑚9𝑆 − 𝑔𝐵𝑂. (1) 𝑃𝐷𝐴 = . (10)
𝑑𝜏 𝑉𝐷𝜇𝐷𝐴
The mass of solvent vapor in 𝐵 is 𝑀𝐵𝑂:
(2) Water. The water in the transformer 𝐷 includes the liquid
𝑑𝑀𝐵𝑂 absorbed water existing in the insulation material and the
= 𝑔𝐵𝑂 − 𝑚2𝑂. (2)
𝑑𝜏 water vapor evaporating into the space. About the water vapor
see the following.
Phase change rate is 𝑔𝐵𝑂:
The mass of water vapor is 𝑀𝐷𝑉:
𝜇𝐵𝑆 𝑃 𝑃
𝑔𝐵𝑂 = 𝜕𝐵𝑆 ⋅ √ ⋅ ( 𝐵𝑆 − 𝐵𝑂 ) ⋅ 𝑆𝐵 𝑑𝑀𝐷𝑉
2𝜋𝑅 √𝑇𝐵𝑆 √𝑇𝐵𝑂 = 𝑔𝐷𝑉 − 𝑚3𝑉 . (11)
(3) 𝑑𝜏
𝑃𝐵𝑆 𝑃 The temperature of water vapor is 𝑇𝐷𝑉 and it is equal to
= 𝑘𝐵𝑆 ⋅ ( − 𝐵𝑂 ) ,
√𝑇𝐵𝑆 √𝑇𝐵𝑂 the temperature of solvent vapor 𝑇𝐷𝑂:
The partial pressure of water vapor is 𝑃𝐷𝑉: The condensation heat transfer coefficient of solvent
vapor in the tank is 𝛽𝐷𝑂:
𝑀𝐷𝑉 𝑅𝑇𝐷𝑉
𝑃𝐷𝑉 = . (13) 2 2
𝑉𝐷𝜇𝐷𝑉 𝛽𝐷𝑂 = 𝑘𝐷 (𝐶𝐷𝑆 − 𝑀𝐷𝑆 ). (22)
About the liquid water see the following. 𝛽𝐷𝑂 is related to the mass of condensed liquid solvent. The
The mass of liquid water is 𝑀𝐷𝑊: condensation heat transfer coefficient is decreasing with the
mass increasing of condensed liquid solvent. The maximum
𝑑𝑀𝐷𝑊 mass of condensed liquid solvent in the tank is defined as
= −𝑔𝐷𝑉. (14)
𝑑𝜏 𝐶𝐷𝑆 . When 𝑀𝐷𝑆 = 𝐶𝐷𝑆 , then 𝛽𝐷𝑂 = 0, which means that the
The temperature of liquid water is 𝑇𝐷𝑊 and it is equal to solvent vapor will no longer condense and the solvent vapor
the temperature of transformer body 𝑇𝐷: entering the tank is mainly taken out by the leakage pump
through the logistics pipeline 3. In this case, the intermediate
𝑇𝐷𝑊 = 𝑇𝐷. (15) pressure lower stage should be entered.
About the liquid solvent see the following.
The saturation pressure of liquid water is 𝑃𝐷𝑊 given by The mass of liquid solvent is 𝑀𝐷𝑆 :
the following equation [15, 16]:
𝑑𝑀𝐷𝑆
1657.46 = 𝑔𝐷𝑂 − 𝑚6𝑆 . (23)
lg𝑃𝐷𝑊 = 10.07406 − . (16) 𝑑𝜏
𝑇𝐷𝑊 − 45.98
The temperature of liquid solvent is 𝑇𝐷𝑆 :
The evaporation phase change rate of liquid water is 𝑔𝐷𝑉: 𝑇𝐷𝑆 = 𝑇𝐷. (24)
𝜇𝐷𝑊 𝑃 𝑃 The saturation pressure of liquid solvent is 𝑃𝐷𝑆 :
𝑔𝐷𝑉 = 𝜕𝐷𝑊 ⋅ √ ⋅ ( 𝐷𝑊 − 𝐷𝑉 ) ⋅ 𝑆𝐷
2𝜋𝑅 √𝑇𝐷 √𝑇𝐷𝑉
2153.6
(17) lg𝑃𝐷𝑆 = 10.1 − . (25)
𝑃 𝑃 𝑇𝐷𝑆 − 34.3
= 𝑘𝐷𝑊 ⋅ ( 𝐷𝑊 − 𝐷𝑉 ) ,
√𝑇𝐷 √𝑇𝐷𝑉 (4) Transformer Body. The transformer body includes iron
cores, copper wires, and insulation materials of transformer
where 𝑘𝐷𝑊 = 𝜕𝐷𝑊 ⋅ √𝜇𝐷𝑊/2𝜋𝑅 ⋅ 𝑆𝐷. core and transformer tank. The temperature of transformer
body is 𝑇𝐷:
(3) Solvent. The solvent includes the solvent vapor entering
through the logistics pipeline 2 and the liquid solvent pro- 𝑑𝑇𝐷
(𝑐𝐷𝑗 𝑀𝐷𝑗 + 𝑐𝐷𝑐 𝑀𝐷𝑐 + 𝑐𝐷𝑖 𝑀𝐷𝑖 + 𝑐𝐷𝑏 𝑀𝐷𝑏 )
duced by condensation phase change. Most of the liquid sol- 𝑑𝜏 (26)
vent flows out through the logistics pipeline 6, and the rest
retains on the wall surface of transformer tank and in the = 𝑔𝐷𝑂 ⋅ ℎ𝐷𝑆 − 𝑚6𝑆 ⋅ ℎ6𝑆 + 𝜆 𝐷𝑂 ⋅ 𝑆𝐷 ⋅ (𝑇𝐷𝑂 − 𝑇𝐷𝑆 ) .
insulation material.
About the solvent vapor see the following. 3.3. Process Logistics Model of the MVPD Equipment
The mass of solvent vapor is 𝑀𝐷𝑂:
(1) Logistics 2. There is only the solvent vapor in the logistics
𝑑𝑀𝐷𝑂 pipeline 2 and the flow from evaporator 𝐵 towards trans-
= 𝑚2𝑂 − 𝑔𝐷𝑂 − 𝑚3𝑂. (18) former 𝐷 is in accordance with the Poiseuille flow. The flux
𝑑𝜏
of logistics 2 depends on the pressure difference between the
The temperature of solvent vapor is 𝑇𝐷𝑂: solvent vapor pressure in the evaporator and the total pres-
sure of mixed gas in the transformer, as well as the structure
𝑇𝐷𝑂 = 𝑇𝐵𝑂 − Δ𝑇, (19)
of pipeline 2.
where Δ𝑇 is the temperature loss of pipeline 2, Δ𝑇 = 3 K. The volume flow of solvent vapor is 𝑄2𝑂:
The pressure of solvent vapor is 𝑃𝐷𝑂:
𝜋𝑟24 2 2
𝑄2𝑂 = ⋅ (𝑃𝐵𝑂 − 𝑃𝐷𝑇 ). (27)
𝑀𝐷𝑂𝑅𝑇𝐷𝑂 16𝜂2 𝐿 2
𝑃𝐷𝑂 = . (20)
𝑉𝐷𝜇𝐷𝑂 The mass flow of solvent vapor is 𝑚2𝑂:
The condense-phase-change rate of solvent vapor is 𝑔𝐷𝑂: 𝑚2𝑂 = 𝜌2𝑂 ⋅ 𝑄2𝑂
𝜇𝐷𝑂 𝑃 𝑃
𝑔𝐷𝑂 = 𝛽𝐷𝑂 ⋅ √ ⋅ ( 𝐷𝑂 − 𝐷𝑆 ) ⋅ 𝑆𝐷 𝑃𝐵𝑂𝜇𝐵𝑂 𝜋𝑟24 2 2
2𝜋𝑅 √𝑇𝐷𝑂 √𝑇𝐷𝑆 = ⋅ ⋅ (𝑃𝐵𝑂 − 𝑃𝐷𝑇 )
𝑅𝑇𝐵𝑂 16𝜂2 𝐿 2 (28)
(21)
𝑃𝐷𝑂 𝑃 𝑃𝐵𝑂
= 𝑘𝐷𝑂 ⋅ 𝛽𝐷𝑂 ⋅ ( − 𝐷𝑆 ) , = 𝑘𝐵𝑂 ⋅ 2
⋅ (𝑃𝐵𝑂 2
− 𝑃𝐷𝑇 ),
√𝑇𝐷𝑂 √𝑇𝐷𝑆 𝑇𝐵𝑂
where 𝑘𝐷𝑂 = √𝜇𝐷𝑂/2𝜋𝑅 ⋅ 𝑆𝐷. where 𝑘𝐵𝑂 = (𝜇𝐵𝑂/𝑅) ⋅ (𝜋𝑟24 /16𝜂2 𝐿 2 ).
Mathematical Problems in Engineering 7
(2) Logistics 3. The substance composition in logistics 3 (3) Other Logistics. The mass flow of condensed solvent liquid
includes solvent vapor, water vapor, and air. The flow of mix- in logistics 6 is directly related to the phase change rate of
ture from transformer 𝐷 towards condenser 𝐶 is in accor- solvent vapor in transformer 𝐷. Consider the following:
dance with Poiseuille flow. The flux of mixed gas depends
on the pressure difference between the total pressure in 𝑚6𝑆 = 𝑘6𝑆 ⋅ 𝑔𝐷𝑆 . (38)
condenser 𝐶 and the total pressure in the transformer, as well
as the structure of pipeline 3. The volume flow ratio of each The mass flow of condensed solvent liquid in logistics 5 is
gas is proportional to the partial pressure ratio of this gas in equal to the mass flow of solvent vapor in flow 3:
transformer 𝐷. The parameters 𝑚3𝐴, 𝑚3𝑉 , and 𝑚3𝑂 in (8), (11), 𝑚5𝑆 = 𝑚3𝑂. (39)
and (18) can be calculated through the following equation.
The total volume flow of mixed gas is 𝑄3𝑇 : The mass flow of liquid solvent in logistics 9 is equal to
the sum of logistics 5 and logistics 6:
𝜋𝑟34 2
𝑄3𝑇 = ⋅ (𝑃𝐷𝑇 − 𝑃𝐶2 ) , (29) 𝑚9𝑆 = 𝑚6𝑆 + 𝑚5𝑆 . (40)
16𝜂𝐷𝑇 𝐿 3
where the total pressure of mixed gas is equal to the sum of 4. Simulation and Discussion
each gas partial pressure; namely,
Based on the mathematical model given above, we can carry
𝑃𝐷𝑇 = 𝑃𝐷𝑂 + 𝑃𝐷𝐴 + 𝑃𝐷𝑉. (30) out the numerical simulation of MVPD process. Under the
limit of thesis space, the model established above and the
The equivalent viscosity coefficient of mixed gas depends
following simulation will focus on the heat and mass trans-
on the viscosity coefficient of each gas and the proportion of
formation and transfer process of various working medium
mixed gas; namely,
in evaporator 𝐵 and vacuum chamber 𝐷. Much attention is
𝜌𝐷𝐴 paid to the variations of the transformer body temperature,
𝜂𝐷𝑇 = 𝜂𝐷𝐴 ⋅ the pressure in the tank, and the water yield while the process
𝜌𝐷𝐴 + 𝜌𝐷𝑂 + 𝜌𝐷𝑉
parameters changes of other auxiliary process units and flows
𝜌𝐷𝑂 𝜌𝐷𝑉 are not concerned.
+ 𝜂𝐷𝑂 ⋅ + 𝜂𝐷𝑉 ⋅ .
𝜌𝐷𝐴 + 𝜌𝐷𝑂 + 𝜌𝐷𝑉 𝜌𝐷𝐴 + 𝜌𝐷𝑂 + 𝜌𝐷𝑉
(31) 4.1. Simulation Calculation Flow Chart. The simulation is in
accordance with time sequence and time is taken as the inde-
The volume flow of air in mixed gas is pendent variable; simulation process starts at a time point
𝑃𝐷𝐴 when all necessary technological parameters are given as the
𝑄3𝐴 = ⋅ 𝑄3𝑇 . (32) initial condition. The recurrence method is used to extrapo-
𝑃𝐷𝑇
late calculation value of every variable at the next moment in
The volume flow of water steam in mixed gas is turn (time step is 1 second). Every simulation cycle finishes a
whole MVPD process (about dozens of hours). In the process
𝑃𝐷𝑉 of simulation, the end criterion conditions of preparation
𝑄3𝑉 = ⋅ 𝑄3𝑇 . (33) stage, heating stage, intermediate pressure lower stage, and
𝑃𝐷𝑇
high vacuum stage are set. As simulation parameters reach
The volume flow of solvent vapor in mixed gas is the end criterion conditions, simulation program will shift to
the next stage automatically. The simulation calculation flow
𝑃𝐷𝑂 chart of MVPD process is shown in Figure 4.
𝑄3𝑂 = ⋅ 𝑄3𝑇 . (34)
𝑃𝐷𝑇
4.2. Calculation Example. A simulation example is given
If the volume flow is converted into the mass flow, the
below. This example comes from an actual field drying pro-
mass flow of air in mixed gas is
cess of a transformer. The voltage insulation class of this
𝑃𝐷𝐴 𝜇𝐷𝐴 large power transformer is 500 kV and the capacity of it is
𝑚3𝐴 = ⋅ 𝑄3𝐴 . (35) 750 MVA. The MVPD 8.2 equipment is adopted in actual pro-
𝑅𝑇𝐷
duction to perform dry processing.
The mass flow of water steam in mixed gas is According to the formulas of heat and mass transfer
model of MVPD system, simulation calculation procedure is
𝑃𝐷𝑉𝜇𝐷𝑉 programmed with MATLAB software. Based on actual pro-
𝑚3𝑉 = ⋅ 𝑄3𝑉 . (36)
𝑅𝑇𝐷 cess, the whole simulation time of heating stage and interme-
diate pressure lower stage is set as 66 hours, and time step is 1
The mass flow of solvent vapor in mixed gas is second. The structural and physical parameters of both trans-
former and MVPD equipment used in simulation calculation
𝑃𝐷𝑂𝜇𝐷𝑂
𝑚3𝑂 = ⋅ 𝑄3𝑂. (37) and the constant parameters in simulation process are listed
𝑅𝑇𝐷 in Table 4.
8 Mathematical Problems in Engineering
400
Start
380
Preparation stage 360
TD (K)
340
No
dTD 12000
=0
d𝜏
Manifold cycles 10000
Yes
8000
PDT (Pa)
Middle depressurization stage
6000
4000
dMDL No 2000
=0
d𝜏
0
0 10 20 30 40 50 60 70
Time (h)
Yes
High vacuum stage Figure 6: Simulation curve of pressure in the tank.
No
PD ≤ 10 Pa
Yes
In this simulation result, the simulation curve of total
water yield 𝑀𝐷𝐿 changing over time is shown in Figure 7.
Oiling stage During the actual drying process of 500 kV/750 MVA
transformer, the MVPD8.2 equipment can record the process
End parameters of different times by sensors. Summarizing the
data and program by VB, finally the actual process curve of
Figure 4: The flow chart of calculation program. MVPD can be got as shown in Figure 8. Figure 8 has shown,
respectively, the change law of key process parameters over
time such as temperature value, pressure value, and water
4.3. Simulation Results and Discussion. The simulation in this yield.
paper focuses on the variation of key process parameters By comparing the simulated curve and measured curve,
over time such as transformer body temperature, pressure in we can find that the simulation results are quite coincident
the tank, and water yield. The change curve of transformer with the measured values in the aspect of change trend
body temperature via time is shown in Figure 5. The change of key parameters. Both transformer body temperature and
curve of pressure in the tank via time is shown in Figure 6. pressure in the tank increase at heating stage and decrease at
Total water yield 𝑀𝐷𝐿 in the operation process via time is intermediate pressure lower stage. Besides, their range ability
the most intuitive index parameter reflecting transformer and change law are also very consistent.
drying effect. It should be equal to the difference between the Comparison results show that the mathematical model
initial moisture content 𝑀𝐷𝑊0 of insulation materials and the established and the calculation procedure programmed can
residual moisture content 𝑀𝐷𝑊. correctly reflect the heat and mass transfer law of MVPD
Namely, process and can predict quantitatively the change trend of key
physical parameters [17].
𝑀𝐷𝐿 = 𝑀𝐷𝑊0 − 𝑀𝐷𝑊. (41) The obvious difference between the simulation results and
the practical measured values is as follows. The simulation
The dehydration rate of insulation materials at a certain calculation shows that only 8 times heating stages and 7 times
time point is intermediate pressure lower stages are enough for the whole
drying process, while 11 times heating stages and 10 times
𝑑𝑀𝐷𝐿 intermediate pressure lower stages were used in the actual
= 𝑚𝐷𝑉 . (42)
𝑑𝜏 operation process. The reason for this difference is that
10 Mathematical Problems in Engineering
500 250
450
400 200
350
300
MDL (kg)
150
PB (kW)
250
200 100
150
100
50
50
0
0 10 20 30 40 50 60 70 0
0 10 20 30 40 50 60 70
Time (h) Time (h)
Figure 7: Simulation curve of total water yield.
Figure 9: Curve of evaporator heat power via time in the drying
process.
12-02-25
12-02-26
12-02-27
12-02-27
12-02-28
12-02-29
12-02-29
12-03-01
12-03-01
12-03-02
12-03-03
12-02-25
200
180
5. Conclusion
160
140 On the basis of analyzing in detail the working principle,
120 system composition, and technological process of the MVPD
100 equipment, a simplified process model reflecting the rela-
80 tionship between internal operation units and processes is
60
extracted in this paper. According to the principle of mass
40
20 and energy conservation, a heat mass transfer mathematical
0 model including 40 mathematical equations is established.
[Link]
[Link]
[Link]
[Link]
[Link]
[Link]
[Link]
[Link]
[Link]
[Link]
[Link]
[Link]
Acknowledgments
This research was supported by National Natural Sci-
ence Foundation of China (Grant nos. 31000665, 51176027,
31371873, and 31300408) and the Fundamental Research
Funds for the Central Universities of China (Grant no.
N130403001).
References
[1] Y. Li, The Study of Mobile Oil Vapor Phase Drying Equipment
on Oil-Immersed Transformer, North China Electric Power
University, 2009 (Chinese).
[2] B. G. Zhang, Q. W. Zhu, and F. Zhang, “Domestic oil vapor
phase drying equipment and application,” Transformer, no. 10,
1995 (Chinese).
[3] S. W. Ruan, “Oil vapor phase drying principle and process
improvement,” Process and Technology, vol. 9, 2013 (Chinese).
[4] W. Q. Li, Mechanism and Testing of Moisture Content Measure-
ment of Transformer Pressboard, Dalian University of Technol-
ogy, 2005, (Chinese).
[5] B. Z. Qiao and S. W. Zhang, “Mobile vacuum vapor phase drying
plant and its application,” Vacuum, vol. 50, no. 6, pp. 60–67, 2013
(Chinese).
[6] L. R. Liu, “Mobile oil vapor phase drying method and device,”
North China Electric Power, 2011 (Chinese).
[7] D. F. Garcı́a, B. Garcı́a, and J. C. Burgos, “Modeling power trans-
former field drying processes,” Drying Technology, vol. 29, no. 8,
pp. 896–909, 2011.
[8] S. D. Foss and L. Savio, “Mathematical and experimental anal-
ysis of the field drying of power transformer insulation,” IEEE
Transactions on Power Delivery, vol. 8, no. 4, pp. 1820–1828, 1993.
[9] W. W. Guidi and H. P. Fullerton, “Mathematical methods for
prediction of moisture take-up and removal in large power
transformers,” in Proceedings of the IEEE Winter Power Meeting,
pp. 242–244, 1974.
[10] S. W. Zhang and B. Z. Qiao, “Calculation of heat and mass trans-
fer in the process of transformer vacuum vapor phase drying,”
Drying Technology and Equipment, vol. 4, no. 3, 2006 (Chinese).
[11] C. H. Xu, S. W. Zhang, and K. Z. Guan, Vacuum Drying, Chem-
ical Industry Press, Beijing, China, 2004, (Chinese).
[12] B. Z. Qiao, Vacuum Vapor Phase Drying and Its Application,
Vacuum, 1988 (Chinese).
[13] S. H. Lin, “Prediction of the drying rate of transformer insula-
tion during the dry cycle,” Electric Power Systems Research, vol.
23, no. 3, pp. 227–231, 1992.
[14] H. Wildmoser, J. Scheiwiller, and E. J. Windhab, “Impact of
disperse microstructure on rheology and quality aspects of ice
cream,” LWT—Food Science and Technology, vol. 37, no. 8, pp.
881–891, 2004.
[15] M. J. Heathcote, J & P Transformer Book, pp. 597–763, Elsevier,
2007.
[16] J. H. Keenan and F. G. Keyes, Thermodynamic Properties of
Steam, John Wiley & Sons, New York, NY, USA, 1936.
[17] J. A. Almendros-Ibáñez, Transformer Field Drying Procedures:
A Theoretical Analysis, 1978.
Hindawi Publishing Corporation
Mathematical Problems in Engineering
Volume 2014, Article ID 794650, 15 pages
[Link]
Research Article
Research on the Impact of Wind Angles on the Residential
Building Energy Consumption
Copyright © 2014 Zou Huifen et al. This is an open access article distributed under the Creative Commons Attribution License,
which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Wind angles affect building’s natural ventilation and also energy consumption of the building. In winter, the wind direction in
the outdoor environment will affect heat loss of the building, while in summer the change of wind direction and speed in the
outdoor environment will affect the building’s ventilation and indoor air circulation. So, making a good deal with the issue of the
angle between local buildings and the dominant wind direction can effectively solve the winter and summer ventilation problems.
Thereby, it can enhance the comfort of residential person, improve indoor air quality, solve heat gain and heat loss problems in
winter and summer in the severely cold and cold regions, and reduce building energy consumption. The simulation software CFD
and energy simulation software are used in the paper. South direction of the building is the prototype of the simulation. The angle
between the direction of the building and the outdoor environment wind is changed sequentially. Energy consumption under
different wind angle conditions is compared with each other. Combined with natural ventilation under various wind angles, the
paper gives the best recommended solution of building direction in Shenyang.
− −
−
+ − + −
(a) (b)
−
+ −
+ −
+ −
(c) (d)
Figure 3: Pressure distributions under different roof heights. (a) Flat roof building (vertical section), (b) pitched roof construction with 30∘
slope, (c) pitched roof construction with 45∘ slope, and (d) building plans.
outdoors, the value of a pressure point of building envelope variation can be described by exponential function, and its
structure can be expressed as follows [10]: expression is shown as follows:
𝐻 𝛼
V𝑜 = V10 ( ) , (3)
𝜌𝑜 V𝑜2 10
𝑃=𝐾 , (2)
2
where V𝑜 is outdoor air velocity at the height of H, m/s; V10 is
outdoor wind speed at the standard height (10 m), m/s; 𝐻 is
where 𝑃 is the value of a pressure point on the building height of any point, m; 𝛼 is roughness coefficient of surface, in
envelope structure, Pa; 𝐾 is aerodynamic coefficient; V𝑜 is a built-up city center 𝛼 = 0.20, in offshore island and desert
outdoor air flow rate, m/s; 𝜌𝑜 is outdoor air flow rate, kg/m3 . 𝛼 = 0.12, and in rural and suburban 𝛼 = 0.16.
The aerodynamic coefficient K of different shapes of
buildings is different at different wind directions. When the 2.3. Natural Ventilation under the Action of Wind Pressure and
K value is positive, the wind pressure value of the point is Heat Pressure. In most cases, building natural ventilation is
positive, and when the K value is negative, the wind pressure not caused by a factor, but by the factors of wind pressure
value of the point is negative. Under the action of wind and heat pressure together. Their functions are varying and
pressure, the wind enters from the positive pressure side and their contributions to the whole building’s natural ventilation
gets out from the negative pressure side. Figure 3 represents are different. When the building is bearing heat and wind
the wind pressure distributions under different roof heights pressure at the same time, the internal and external pressure
[11]. In any case, the wind pressure is negative near the surface difference of each window hole is equal to the sum of
of flat roofs and irrelevant to the wind direction, and its the pressure difference of inside and outside window hole
variation is very small. When the pitched roof ’s slope is small, when the building bears the heat pressure and the wind
both the front and the back have negative pressure while, with pressure independently [12]. Effect of wind pressure and heat
the increase of slope, the front roof has positive pressure, and pressure on the natural ventilation in architecture is often
the back roof has negative pressure. inseparable, and the heat pressure of building is relatively
In the surface layer, outdoor wind speed of v has larger stable and easy to implement, but the action of the wind
variation with height; the average wind speed with the height is unstable and is often affected by the building’s shape, the
4 Mathematical Problems in Engineering
4200
Study
bedroom known as mass conservation equations, whose differential
expression is
13800
A
𝜕𝜌 𝜕
+ (𝜌𝑢𝑖 ) = 0,
2400
(4)
2400
13800
Lobby
𝜕𝑡 𝜕𝑥𝑖
1500
where 𝜌 is expressed as the density of the fluid and 𝜇𝑖 is
expressed as the fluid component along the i direction.
4500
𝜕 (𝜌𝑢𝑖 ) 𝜕 (𝜌𝑢𝑖 𝑢𝑗 )
+
1200
1200
Balcony 𝜕𝑡 𝜕𝑥𝑗
3000 3300 2100 1500 𝜕𝑝 𝜕 𝜕𝑢 𝜕𝑢𝑗
9 =− + [𝜂 ( 𝑖 + )] + 𝜌𝛽 (𝑇0 − 𝑇) 𝑔𝑖 .
𝜕𝑥𝑖 𝜕𝑥𝑗 𝜕𝑥𝑗 𝜕𝑥𝑖
(5)
Figure 4: Apartment A.
Energy equation is
(1) Mathematical Model Building. Mathematical model of (3) Select Turbulence Model. After the study of turbulence
the problem is the mathematical description of the study. by people for more than a century, now it has a variety of
According to thermal models, indoor gas meets the equation turbulence models, such as Reynolds stress model (RSM),
of the ideal gas state. Secondly, the pressure is usually treated algebraic stress model (ASM), and sticky coefficient model
as a constant flow of indoor air. Further, the common flow (EVM). Based on solving differential turbulent viscosity
of indoor air flowing is substantially vulgar, with flow often coefficient equations, EVM is divided into zero-equation
at 10∼20 m/s or less, and, therefore, the indoor air can be model, one-equation model, and two equation model. Their
treated as incompressible fluid [13]. Finally, the indoor air calculations have small amount and they are very suitable for
viscosity cannot be neglected; we must use the theory to study engineering applications where 𝑘-𝜀 model is most commonly
the dynamics of a viscous fluid. Further, the indoor air flow used [17].
is usually a turbulent flow, with the corresponding need to 𝑘-𝜀 model is a semiempirical formula model; it is assumed
simulate turbulence theory [14]. to be a fully developed turbulent flow, and viscosity derived
Mathematical Problems in Engineering 5
N
E
Normal
2.00
1.75
Entrance wind 1.50
1.25
Speed (m/s)
1.00
𝛼 Building surface 0.75
0.50
0.25
0.00
Figure 6: Wind angle.
𝜕 (𝜌𝜀) 𝜕 (𝜌𝑢𝑗 𝜀)
+
𝜕𝑡 𝜕𝑥𝑗
𝜕 𝜂 𝜕𝜀 𝜀
= ( − ) + + (𝑐1 𝐺𝑘 + 𝑐3 𝐺𝑏 − 𝑐2 𝜌𝜀) ,
𝜕𝑥𝑗 𝜎𝜀 𝜕𝑥𝑖 𝑘
(8) Figure 7: The wind speed distribution of the wind angle 𝛼 = 75∘ .
6 Mathematical Problems in Engineering
2
2.00 2.00
11.75 1.75
Speed (m/s)
11.50 1.50
Speed (m/s)
11.25 1.25
11.00 1.00
00.75 0.75
00.50 0.50
00.25 0.25
00.00 0.00
Figure 8: The wind speed distribution of the wind angle 𝛼 = 60∘ . Figure 9: The wind speed distribution of the wind angle 𝛼 = 45∘ .
where, 𝑐1 , 𝑐2 are the empirical constants and they are 1.44, (b) The preferred shape of the grid cell is a regular
1.92, respectively. 𝐾 is the turbulent dynamic coefficient. 𝜀 is hexahedron or tetrahedron.
the turbulent kinetic energy dissipation rate, the relationship
(c) A coarse mesh should be used in the place where
between them is
temperature gradient and velocity gradient are small.
𝑘3/2 However, mesh refinement should be used in the
𝜀 = 𝐶𝐷 . (9) place where temperature gradient and velocity gradi-
𝑙
ent are very big.
Laminar turbulent kinetic energy is generated by the
velocity gradient; the expression is 3.1.2. Physical Model. In this paper, the simulation object
is located in a residential district in Shenyang city, and it
𝜕𝑢𝑖 𝜕𝑢𝑗 𝜕𝑢𝑖 is seven layers of ordinary residence; its total construction
𝐺𝑘 = 𝜇𝑡 ( + ) . (10)
𝜕𝑥𝑗 𝜕𝑥𝑖 𝜕𝑥𝑗 area is 1669.18 m2 . Every floor has two households, and the
height of each household is 3 m. The function layout of every
𝐺𝑏 of buoyancy turbulent kinetic energy is generated; the room is the same with an apartment layout (see Figure 4). In
expression is the simulations, we headed for the building of the south as
the prototype of simulation, and the model simplified in the
𝜇𝑡 𝜕𝑇 premise of indoor thermal characteristics is not missing, as
𝐺𝑏 = 𝛽𝑔𝑖 ; (11)
𝜎𝑇 𝜕𝑥𝑖 shown in Figure 5.
for turbulent sticky factor, the expression is
3.2. The Setting of the Simulation Conditions. This paper
2 has studied the impact of the wind angle (see Figure 6) on
𝑐𝜇 𝜌𝑘
𝜇𝑡 = . (12) buildings’ energy consumption. Wind angle 𝛼 is the angle
𝜀 between the normal outdoor environment wind and the
(4) Meshing. CFD meshing software usually follows the normal direction of inlet plane (building surface). When the
following principles. flow is perpendicular to the surface of the building, 𝛼 is 0∘ ,
and when the flow is parallel to the surface of buildings, 𝛼 is
(a) A grid cell to another grid cell expansion coefficient 90∘ . As provisions, the clockwise direction is positive and the
ratio is generally maintained between 2 and 5; the counterclockwise direction is negative. The angles are shown
critical areas should be smaller. in Figure 6.
Mathematical Problems in Engineering 7
2.00 2.00
1.75 1.75
1.50 1.50
Speed (m/s)
Speed (m/s)
1.25 1.25
1.00 1.00
0.75 0.75
0.50 0.50
0.25 0.25
0.00 0.00
Figure 10: The wind speed distribution of the wind angle 𝛼 = 30∘ . Figure 12: The wind speed distribution of the wind angle 𝛼 = 0∘ .
Table 1: The residential construction and personnel structures under different modes.
184.571
568.068
161.500
497.060 Mean age of air (s)
Figure 14: The distribution of average air age of the wind angle 𝛼 =
Figure 13: The distribution of average air age of the wind angle 𝛼 = 60∘ .
75∘ .
people in the master bedroom, no more than 1 person in the Room Bedroom Living room Kitchen
second bedroom, and no more than 3 people in the living 7:00∼8:00
Lighting 22:00∼23:00 19:00∼22:00
room. The calorific value for one person in the room is 53 W. 18:00∼19:00
The largest calorific value of lighting, devices in living room, 7:00∼8:00
Devices All day is 0 19:00∼23:00
and devices in kitchen is 5 W/m2 , 9.3 W/m2 , and 48.2 W/m2 , 18:00∼19:00
respectively [21]. Changes of parameters are in Table 2. 7:00∼8:00
People 22:00∼7:00 19:00∼22:00
18:00∼19:00
(2) Air Condition Parameters Settings. Air conditioners and
heating devices would affect the indoor thermal environment
conditions. When room temperature is higher (higher than conditioners when the temperature is higher than the toler-
body’s tolerable temperature), people would turn on air able temperature rather than the desired temperature. It is
conditioners to adjust the room temperature and humidity found in the investigation that requirement of the thermal
in summer. In other words, people would turn on air comfort cannot be satisfied by natural ventilation when room
Mathematical Problems in Engineering 9
136.48 124.97
119.42 109.35
Figure 17: The distribution of average air age of the wind angle 𝛼 =
Figure 15: The distribution of average air age of the wind angle 𝛼 = 15∘ .
45∘ .
95.960
79.967
63.974 4. Analysis of Effects of Natural Ventilation
47.980
under Different Wind Angles
31.987
15.993 Laws of wind angles affecting natural ventilation are derived
0.0000 through the analysis of the effects of natural ventilation under
different wind angles. And they would provide guidance on
design of residential building orientation.
In the process of natural ventilation, the wind direction
of outdoor environment will change complicatedly. To make
better use of natural ventilation, when designing residential
buildings, the paper keeps an appropriate wind direction
angle of incidence between the local dominant wind direction
in summer and the residential building. The dominant wind
direction in summer in Shenyang is southwest. But the
orientation of conventional residential design is almost south
or north. South orientation is used as the standard in the
process of simulation. The variable factor of the simulation
is wind angle which would be 75∘ , 60∘ , 45∘ , 30∘ , 15∘ , and 0∘ .
The influence of wind speed at human height and air age
Figure 16: The distribution of average air age of the wind angle 𝛼 = distribution on human body is considered in the analysis of
30∘ . indoor ventilation’s effects. So, velocity field at 1.1 m height
10 Mathematical Problems in Engineering
600 (2) The Distribution of Age of Indoor Air under Different Wind
Average age of indoor air (s)
500 Angles. Figures 13, 14, 15, 16, 17, and 18 show the distribution
of the average age of indoor air at 1.1 m height (an average
400
sitting person’s nose height from the ground) in the cases of
300 different wind angles. Color in the picture transfers from the
200 cool tone to the warm tone. The color represents the average
100
age of indoor air. The cool tone means that the age value of
air is small. That means it needs a shorter time to complete
0
75∘ 60 45 30 15 0 the replacement under the condition of natural ventilation
Wind direction angle (∘ ) and fresh air is better; on the contrary, the warm tone means
that the age value of air is large, and it needs a longer time
Maximum value to complete the replacement under the condition of natural
Minimum value ventilation and the fresh air is poorer.
Average age
The distributions of indoor average air age of different
Figure 19: The contrast figures of average age of indoor air under wind angles are shown in Figures 13, 14, 15, 16, 17, and 18,
different working conditions. they show that (a) in the case of larger wind angle, the value
of indoor average air is generally higher. With the decrease
of the wind angle, the average age of indoor air in the same
area will significantly reduce and the degree of fresh air will
[22] (an average sitting person’s nose height from the ground)
increase, which improves the indoor air quality. (b) Because
and average air age field are chosen to be objects of study.
there is no outside window in the small indoor toilet, effective
(1) Wind Velocity Distribution under Different Wind Angles. air convection cannot form. Almost no wind flows into the
When outdoor air enters building through the opening in room and the airflow velocity is small. Thus, the average age
the surface, a velocity field will form in the interior. In of air is relatively bigger here than that in the other rooms.
the case of the same speed and direction of the outdoor The fresh air is bad, so we need to use auxiliary ventilation
wind and different outdoor wind angles, indoor velocity through the mechanical approach in this region. (c) Under
distribution differs. Figures 7, 8, 9, 10, 11, and 12 show airflow the predominant wind direction in summer, the ventilation
velocity distribution contours of 75∘ , 60∘ , 45∘ , 30∘ , 15∘ , and 0∘ blind spot will form in the master bedroom and in some
wind angles at 1.1 m height (an average sitting person’s nose areas of the residential study room. With greater wind angles
height from the ground). (except 0∘ ), the range of ventilation blind spot will be larger.
Mathematical Problems in Engineering 11
200
3
150
2
100
1
50
0 0
Jan. Feb. Mar. Apr. May Jun. Jul. Aug. Sep. Oct. Nov. Dec. Jan. Feb. Mar. Apr. May Jun. Jul. Aug. Sep. Oct. Nov. Dec.
Figure 20: The residential energy consumption distribution when the wind angle is 0∘ .
200
3
150
2
100
1
50
0 0
Jan. Feb. Mar. Apr. May Jun. Jul. Aug. Sep. Oct. Nov. Dec. Jan. Feb. Mar. Apr. May Jun. Jul. Aug. Sep. Oct. Nov. Dec.
Figure 21: The consumption distribution of the residential energy when the wind angle is 15∘ .
By the comparison analysis of average age of indoor air wind speed and the distributions of the air age in the case
under the different wind angles in Figure 19, it can be seen of different wind angles, it can be concluded that the indoor
that wind angle gets greater; the average age of indoor air will ventilation effect is better when the wind angle varies from 0∘
be bigger. When wind angle is 0∘ , the average age of indoor air to 45∘ .
will be the minimum and the fresh air will be the best. When
the range of wind angle varies from 15∘ to 45∘ , the change of
the average age of indoor air will not be much. It is basically 5. Residential Building Energy Simulation
stable within the range and the indoor air quality is relatively Analysis Performed under Different Wind
good.
Angles throughout the Year
In residential space, people mainly focus on the indoor
environment of the bedroom, study, and living room (people According to the predefined lifestyle and air conditioning
stay there for a long time). People (male or female) mainly operation mode, we use EQUEST software to
stay in these residential rooms, so whether indoor environ- perform residential building energy simulation
ment in these rooms is good or bad is directly related to analysis under different wind angles throughout the
the indoor personnel health. When we consider the indoor year. Figures 20, 21, 22, 23, 24, and 25 show the annual
12 Mathematical Problems in Engineering
200
3
150
2
100
1
50
0 0
Jan. Feb. Mar. Apr. May Jun. Jul. Aug. Sep. Oct. Nov. Dec. Jan. Feb. Mar. Apr. May Jun. Jul. Aug. Sep. Oct. Nov. Dec.
Figure 22: The consumption distribution of the residential energy when the wind angle is 30∘ .
200
3
150
2
100
1
50
0 0
Jan. Feb. Mar. Apr. May Jun. Jul. Aug. Sep. Oct. Nov. Dec. Jan. Feb. Mar. Apr. May Jun. Jul. Aug. Sep. Oct. Nov. Dec.
Figure 23: The consumption distribution of the residential energy when the wind angle is 45∘ .
energy consumption values of the residence under different of each month, they show that (1) consumption in winter
conditions. During the period of residence in the running, and summer is higher than that in transition seasons, and
energy is supplied by electricity and gas. In order to facilitate energy consumption in the hottest month of July and the
uniform comparison, the results are converted into a coldest month of January is higher than that in other months
unified dimension KJ; the specific conversion method is within the same seasons, respectively. (2) In individual energy
1 KWh = 3600 KJ, 1 Btu = 1.06 KJ. consumption, power consumption of devices and energy
Figure 20 shows the residential energy consumption consumption of residential interior lighting are almost the
distribution when the wind angle is 0∘ . Obviously, com- same each month. (3) In the process of annual energy con-
prehensive energy consumption of the residence includes sumption simulation, the energy consumption of residential
the energy consumption of lighting, household appliances, heating mainly occurs during the period from December
air conditioning, and heating. When the wind angle is 0∘ , to February and rarely in November and March. During
residential power consumption throughout the whole year is the period from November to January, it increases with the
27000 KWh (i.e., 97.2 GJ); the gas consumption for heating is decrease of outdoor temperature. During the period from
5.19 × 108 Btu (i.e., 571 GJ). Comparing energy consumption February to March, it reduces along with the rise of outdoor
Mathematical Problems in Engineering 13
200
3
150
2
100
1
50
0 0
Jan. Feb. Mar. Apr. May Jun. Jul. Aug. Sep. Oct. Nov. Dec. Jan. Feb. Mar. Apr. May Jun. Jul. Aug. Sep. Oct. Nov. Dec.
Figure 24: The consumption distribution of the residential energy when the wind angle is 60∘ .
temperature. The energy consumption of air conditioning of residential rooms reduces. In summer, the architectural
mainly occurs in July and August and rises with the decrease sunshine area is small when the building is facing south or
of outdoor temperature. near south. So, the solar radiant heat is relatively low. As a
The accumulation of electricity consumption and the gas result, the yearly comprehensive energy consumption of the
consumption under different angles will be drawn into the building is less when the building is facing south or near
chart months, as shown from Figures 26 to 28. south.
Figures 26 and 27 show the comparison of the residential
energy consumption of each month under different wind 6. Conclusion
angles. The changes in building orientation mainly affect
building heating and energy consumption of air conditioning (1) This paper analyzes the influence of wind direction
and have almost no effect on energy consumption in tran- angle for residential indoor natural ventilation on the
sition seasons. Because the building orientation affects both typical model of a residential district in Shenyang
the solar radiation’s heat gain and the infiltration’s heat loss, through CFD numerical simulation. According to
the heat gained from solar radiation in summer will increase the study of indoor air velocity and air age distri-
the cooling energy consumption of air conditioning, but in bution on different wind direction angles, when the
winter it will reduce the energy consumption of heating. wind direction angle ranges from 0∘ to 45∘ , cross-
Figure 28 shows the effect of the wind direction angle on ventilation will be formed. The wind velocity is large
comprehensive energy consumption of buildings throughout in the room, so the room does not have a wide area
the year. It can be seen that when the building rotates in the of ventilation dead corner basically and the average
positive direction, the comprehensive energy consumption of value of indoor air age is short, which means better
the building throughout the year presents a parabolic trend. freshness of the air. When the wind direction angle
When the wind direction angle is 15∘ , the comprehensive is larger than 45∘ , indoor wind velocity drops rapidly.
energy consumption of building is the minimum, 652.3 GJ. The ventilation dead corner increases gradually; this
As the wind direction angle increases gradually, the yearly is not beneficial for eliminating indoor polluted air
comprehensive energy consumption of the building increases by natural ventilation. Therefore, when the wind
gradually. This shows that the appropriate orientation of direction angle is from 0∘ to 45∘ , the effect of the
building should be south or 15∘ south by west. Being affected indoor natural ventilation is better.
by the changing rule of the solar altitude angle and the (2) This paper has studied the energy consumption of
azimuth angle, the building facing south or near south can the residence with different wind direction angles
make full use of solar energy, and the solar radiant heat is through EQUEST software. The conclusions are as
large, so the heating energy consumption can be reduced follows.
accordingly. While the dominant wind direction of Shenyang (a) In the simulation process of yearly energy
is north or northeast, the building facing south or near south consumption, the energy consumption of res-
can avoid the dominant wind direction in winter and reduce idential heating in Shenyang is mainly con-
the amount of cold air penetrating. As a result, the heat loss centrated from December to February, rarely
14 Mathematical Problems in Engineering
200
3
150
2
100
1
50
0 0
Jan. Feb. Mar. Apr. May Jun. Jul. Aug. Sep. Oct. Nov. Dec. Jan. Feb. Mar. Apr. May Jun. Jul. Aug. Sep. Oct. Nov. Dec.
Figure 25: The consumption distribution of the residential energy when the wind angle is 75∘ .
800
3
700
2
1 600
0 500
Jan. Feb. Mar. Apr. May Jun. Jul. Aug. Sep. Oct. Nov. Dec. 0 15 30 45 60 75
Wind direction angle (∘ )
𝛼 = 0∘ 𝛼 = 45∘
𝛼 = 15∘ 𝛼 = 60∘ Figure 28: The comprehensive energy consumption of residence
𝛼 = 30∘ 𝛼 = 75∘ each year under different wind angles.
𝛼 = 0∘ 𝛼 = 45∘
(c) When the wind direction angle of buildings
𝛼 = 15∘ 𝛼 = 60∘
𝛼 = 30∘ 𝛼 = 75∘
is from 0∘ to 45∘ , the comprehensive energy
consumption of buildings is relatively stable,
Figure 27: The comparison of the residential gas consumption of which is a relatively favorable building orien-
each month under different wind angles. tation. Among all the angles, when the wind
Mathematical Problems in Engineering 15
direction angle is 15∘ , the comprehensive energy [14] W. Yin, Potential estimation and energy efficient evaluation of
consumption of building is the lowest. So, in natural ventilation in building [Dissertation], Hunan University,
terms of energy saving, the orientation should 2009.
be the best solution in architectural design in [15] W. Tao, Numerical Heat Transfer, Xi’an Jiaotong University
Shenyang. Press, Xi’an, China, 2001.
[16] J. I. Woni, People, Climate, Building, China Building Industry
Press, Beijing, China, 1982.
Conflict of Interests [17] B. E. Launder and D. B. Spalding, “The numerical computation
The authors declare that there is no conflict of interests of turbulent flows,” Computer Methods in Applied Mechanics and
Engineering, vol. 3, no. 2, pp. 269–289, 1974.
regarding the publication of this paper.
[18] Y. Wu, Numerical Simulation and Experimental Study on Interior
Natural Ventilation of Residential Building, Tianjin University,
Acknowledgment 2009.
[19] Y. Wang and Y. Wu, “On the improved evaluation of the
The authors acknowledge the support of the China Ministry residential indoor air quality with natural ventilation in cold
of Housing and Urban-Rural Construction Research Project weather,” Journal of Safety and Environment, vol. 8, no. 5, pp.
“Study on the regulation strategies of heat and humidity 104–108, 2008.
environment in closed sheds in cold area,” 2013 (2013-K1-55). [20] J. Yiwen, The Research of Residential Thermal Performance
Evaluation Method, Tsinghua University, 2003.
References [21] ASHARE Standard 55-2002, Thermal Environmental Conditions
for Human Occupancy, American Society of Heating, Refriger-
[1] The annual development report on China building energy ating and Air-Conditioning Engineers, Inc., Atlanta, Ga, USA,
saving, 2009. 2002.
[2] G. Bo, Y. Nanyang, and W. Lei, “The strategy forms and [22] S. Zhu, C. Qi, and W. Li, “Experimental research on thermal
simulation analysis of natural ventilation,” Building Energy comfort of natural wind,” Journal of Beijing Forestry University,
Saving, vol. 264, no. 7, pp. 30–33, 2004. vol. 29, no. 4, pp. 54–58, 2007.
[3] J. P. Liu, Architectural Physics, China Building Industry Press,
Beijing, China, 2009.
[4] L. Yang, Climatic analysis and architectural design strategies for
bio-climatic design [Ph.D. thesis], Xi’an University of Architec-
ture and Technology, 2003.
[5] G. Hongliang, Study the Key Issues of Residential Areas Indoor
Natural Areas Indoor Natural Ventilation, Harbin Institute of
Technology, Shenzhen Graduate School, 2009.
[6] S. Duan, G. Zhang, P. Jianguo et al., “Development in research
of natural ventilation,” Heating Ventilating and Air Conditioning,
vol. 34, no. 3, pp. 22–28, 2004.
[7] L. Chuanzhi, F. Guohui, and X. Shuo, “Research into impact of
buildings of different height on wind-induced natural ventila-
tion,” Journal of Shenyang Jianzhu University (Natural Science),
vol. 23, no. 4, pp. 625–630, 2007.
[8] G. Guangcai, L. Hongxiang, and L. Yuguo, “The application and
research of natural ventilation,” Building Energy and Environ-
ment, vol. 4, pp. 4–7, 2003.
[9] K. Zhang, L. Wang, and N. Ren, “Analysis on the natural
ventilation design of residential buildings in Northern towns,”
Journal of Shenyang Jianzhu University (Social Science), vol. 13,
no. 3, pp. 274–277, 2011.
[10] Z. J. Cai and T. Y. Long, Fluid Mechanics Pumps and Fans, China
Building Industry Press, Beijing, China, 1997.
[11] Y. Lu, M. A. Zuiliang, and P. Zou, Heating Ventilation and Air
Conditioning, China Building Industry Press, Beijing, China,
1997.
[12] L. Wang, B. Gong, and N. Yu, “The research of the thermal
comfort in natural ventilation buildings,” Journal of Harbin
Institute of Technology, vol. 41, no. 12, pp. 254–258, 2009.
[13] W. A. Dalgliesh and D. Surry, “BLWT, CFD and HAM mod-
elling vs. the real world: bridging the gaps with full-scale
measurements,” Journal of Wind Engineering and Industrial
Aerodynamics, vol. 91, no. 12–15, pp. 1651–1669, 2003.
Hindawi Publishing Corporation
Mathematical Problems in Engineering
Volume 2014, Article ID 914725, 7 pages
[Link]
Research Article
Local Fractional Laplace Variational Iteration Method for
Nonhomogeneous Heat Equations Arising in Fractal Heat Flow
Shu Xu,1,2 Xiang Ling,1 Carlo Cattani,3 Gong-Nan Xie,4 Xiao-Jun Yang,5 and Yang Zhao6
1
School of Mechanical and Power Engineering, Nanjing University of Technology, Nanjing 210009, China
2
School of Mechanical Engineering, Huaihai Institute of Technology, Lianyungang 222005, China
3
Department of Mathematics, University of Salerno, Via Giovanni Paolo II, Fisciano, 84084 Salerno, Italy
4
School of Mechanical Engineering, Northwestern Polytechnical University, Xi’an, Shaanxi 710072, China
5
Department of Mathematics and Mechanics, China University of Mining and Technology, Xuzhou, Jiangsu 221008, China
6
Electronic and Information Technology Department, Jiangmen Polytechnic, Jiangmen 529090, China
Copyright © 2014 Shu Xu et al. This is an open access article distributed under the Creative Commons Attribution License, which
permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
The local fractional Laplace variational iteration method is used for solving the nonhomogeneous heat equations arising in the
fractal heat flow. The approximate solutions are nondifferentiable functions and their plots are also given to show the accuracy and
efficiency to implement the previous method.
arising in heat flow with local fractional derivative are inves- such that
tigated. Finally, conclusions are shown in Section 5.
∭ {𝑐𝛼 𝜌𝛼 𝑢𝑡(𝛼) (𝑥, 𝑦, 𝑧, 𝑡) − ∇𝛼
2. The Nonhomogeneous Heat
Equations Arising in Heat Flow with ⋅ [𝑘2𝛼 ∇𝛼 𝑢 (𝑥, 𝑦, 𝑧, 𝑡)] (8)
Local Fractional Derivatives
−𝑔 (𝑥, 𝑦, 𝑧, 𝑡) } 𝑑Ω(𝛾) = 0,
In this section we present the one-dimensional nonhomoge-
neous heat equations arising in heat flow with local fractional which leads to the nonhomogeneous local fractional heat
derivatives. equations [23]:
Let the local fractional volume integral of the function u
be defined as [19] 𝑐𝛼 𝜌𝛼 𝑢𝑡(𝛼) (𝑥, 𝑦, 𝑧, 𝑡) − ∇𝛼 ⋅ [𝑘2𝛼 ∇𝛼 𝑢 (𝑥, 𝑦, 𝑧, 𝑡)]
𝑁 (9)
∭u (𝑟𝑃 ) 𝑑Ω (𝛾)
=
(𝛾)
lim ∑ u (𝑟𝑃 ) ΔΩ𝑃 , (1) = 𝑔 (𝑥, 𝑦, 𝑧, 𝑡) .
𝑁→∞
𝑃=1
From (9) we obtain the nonhomogeneous heat equations in
(𝛾)
where the elements of the volume ΔΩ𝑃 → 0 as 𝑁 → the dimensionless case:
∞ and the fractal dimension of the volume 𝛾. The equality
𝑢 (𝑥, 𝑦, 𝑧, 𝑡) is the temperature at the point (𝑥, 𝑦, 𝑧) ∈ Ω, time 𝜙𝑡(𝛼) (𝑥, 𝑦, 𝑧, 𝑡) − ∇2𝛼 𝜙 (𝑥, 𝑦, 𝑧, 𝑡) = 𝜑 (𝑥, 𝑦, 𝑧, 𝑡) . (10)
𝑡 ∈ 𝑇, and the total amount of heat 𝐻 (𝑡) is described as
The two-dimensional case is [23]
𝐻 (𝑡) = ∭𝑐𝛼 𝜌𝛼 𝑢 (𝑥, 𝑦, 𝑧, 𝑡) 𝑑Ω(𝛾) , (2)
𝜙𝑡(𝛼) (𝑥, 𝑦, 𝑡) − ∇2𝛼 𝜙 (𝑥, 𝑦, 𝑡) = 𝜑 (𝑥, 𝑦, 𝑡) , (11)
where 𝑐𝛼 is the special heat of the fractal material and 𝜌𝛼 is the
density of the fractal material. and the one-dimensional case is [26]
The local fractional surface integral is defined as [19, 22]
𝜙𝑡(𝛼) (𝑥, 𝑡) − 𝜙𝑥(2𝛼) (𝑥, 𝑡) = 𝜑 (𝑥, 𝑡) . (12)
𝑁
(𝛽)
∬ u (𝑟𝑃 ) ⋅ 𝑑S(𝛽) = lim ∑ u (𝑟𝑃 ) ⋅ n𝑃 Δ𝑆𝑃 , (3)
𝑁→∞ 3. Local Fractional Laplace Variational
𝑃=1
where 𝑁 are elements of area with a unit normal local Iteration Method
(𝛽)
fractional vector n𝑃 , Δ𝑆𝑃 → 0 as 𝑁 → ∞ for 𝛾 = (3/2) 𝛽 = In this section, we give the idea of local fractional Laplace
3𝛼. variational method [29, 30] in order to investigate the one-
From (3) the local fractional Fourier law of the material dimensional nonhomogeneous heat equations arising in frac-
in fractal media [19, 23] was suggested as follows: tal heat flow.
𝑑𝛼 We present the following local fractional differential
𝐻 (𝑡) = ∯ 𝑘2𝛼 ∇𝛼 𝑢 (𝑥, 𝑦, 𝑧, 𝑡) ⋅ 𝑑S(𝛽) , (4) operator as follows:
𝑑𝑡𝛼 𝜕Ω(𝛽)
where 𝑑S(𝛽) is the fractal surface measure over Ω(𝛾) and 𝑘2𝛼 is 𝐿 𝛼 𝑢 − 𝑅𝛼 𝑢 = 0, (13)
the thermal conductivity of the fractal material.
In view of (4), the change in heat reads as follows [19, 23]: where the linear local fractional differential operator denotes
𝐿 𝛼 = 𝑑2𝛼 /𝑑𝑥2𝛼 and 𝑢 (𝑥) is a nondifferential function.
𝑑𝛼
𝐻 (𝑡) = ∭𝑐𝛼 𝜌𝛼 𝑢𝑡(𝛼) (𝑥, 𝑦, 𝑧, 𝑡) 𝑑Ω(𝛾) , (5) We can write the local fractional functional formula as
𝑑𝑡𝛼
𝑢𝑛+1 (𝑥) = 𝑢𝑛 (𝑥)
where 𝜕Ω(𝛽) is the boundary of Ω(𝛾) .
From (2) we suggest the following source term [23]: 𝜆(𝑥 − 𝑡)𝛼 (14)
+ 0 𝐼(𝛼)
𝑥 { [𝐿 𝑢 (𝑡) − 𝑅𝛼 𝑢𝑛 ]} .
Γ (1 + 𝛼) 𝛼 𝑛
𝐺 (𝑡) = ∭𝑔 (𝑥, 𝑦, 𝑧, 𝑡) 𝑑Ω(𝛾) . (6)
The local fractional Laplace transform is given as [29–32]
Making use of (4), (5), and (6), we have
̃ 𝛼 {𝑓 (𝑥)}
𝐿
∭𝑐𝛼 𝜌𝛼 𝑢𝑡(𝛼) (𝑥, 𝑦, 𝑧, 𝑡) 𝑑Ω (𝛾)
̃
= 𝑓𝑠𝐿,𝛼 (𝑠)
=∯ 𝑘2𝛼 ∇𝛼 𝑢 (𝑥, 𝑦, 𝑧, 𝑡) ⋅ 𝑑S(𝛽) (7) 1 ∞ (15)
𝜕Ω(𝛽) = ∫ 𝐸𝛼 (−𝑠𝛼 𝑥𝛼 ) 𝑓 (𝑥) (𝑑𝑥)𝛼 ,
Γ (1 + 𝛼) 0
+ ∭𝑔 (𝑥, 𝑦, 𝑧, 𝑡) 𝑑Ω(𝛾) 0 < 𝛼 ≤ 1,
Mathematical Problems in Engineering 3
and the inverse formula of local fractional Laplace transform such that
is suggested as [29–32]
𝛼
̃ −1 {𝑓𝐿,𝛼 (𝑠)}
𝑓 (𝑥) = 𝐿 𝛼 𝑠 ̃ 𝛼 { 𝜆(𝑥) } 𝑠2𝛼 = 0.
1+𝐿 (24)
Γ (1 + 𝛼)
𝛽+𝑖∞ (16)
1 ̃
= 𝛼 ∫ 𝐸 (𝑠𝛼 𝑥𝛼 ) 𝑓𝑠𝐿,𝛼 (𝑠) (𝑑𝑠)𝛼 ,
(2𝜋) 𝛽−𝑖∞ 𝛼 From (24) we get
𝛼
where 𝑓 (𝑥) is a local fractional continuous function, 𝑠 =
𝛽𝛼 + 𝑖𝛼 ∞𝛼 , Re (𝑠𝛼 ) = 𝛽𝛼 , and the local fractional integral of 𝛼
̃ 𝛼 { 𝜆(𝑥) } = − 1
𝐿 (25)
𝑓 (𝑥) of order 𝛼 in the interval [𝑎, 𝑏] is given as [23] Γ (1 + 𝛼) 𝑠2𝛼
𝑏
(𝛼) 1
𝑎 𝐼𝑏 𝑓 (𝑥) = ∫ 𝑓 (𝑡) (𝑑𝑡)𝛼
Γ (1 + 𝛼) 𝑎 such that local fractional iteration algorithm reads as
𝑗=𝑁−1
(17)
1 𝛼
= lim ∑ 𝑓 (𝑡𝑗 ) (Δ𝑡𝑗 ) , ̃ 𝛼 {𝑢𝑛+1 (𝑥)} = 𝐿
𝐿 ̃ 𝛼 {𝑢𝑛 (𝑥)}
Γ (1 + 𝛼) Δ𝑡 → 0 𝑗=0
1 ̃ (26)
with the partitions of the interval [𝑎, 𝑏] which is (𝑡𝑗 , 𝑡𝑗+1 ), − 𝐿 {(𝐿 𝛼 𝑢𝑛 (𝑥) − 𝑅𝛼 𝑢𝑛 (𝑥))} ,
𝑠2𝛼 𝛼
for Δ𝑡𝑗 = 𝑡𝑗+1 − 𝑡𝑗 , 𝑡0 = 𝑎, 𝑡𝑁 = 𝑏, and Δ𝑡 = max{Δ𝑡0 ,
Δ𝑡1 , Δ𝑡𝑗 , . . .}, 𝑗 = 0, . . . , 𝑁 − 1.
where the initial value is presented as follows:
From (15) the local fractional convolution of two func-
tions is defined as [29–32]
∞ ̃ 𝛼 {𝑢0 (𝑥)} = 𝑢 (0) .
𝐿 (27)
1
𝑓1 (𝑥) ∗ 𝑓2 (𝑥) = ∫ 𝑓1 (𝑡) 𝑓2 (𝑥 − 𝑡) (𝑑𝑡)𝛼 , (18)
Γ (1 + 𝛼) −∞
and we have Therefore, the local fractional series solution is given as
̃ 𝛼 {𝑓1 (𝑥)} 𝐿
𝐹𝛼 {𝑓1 (𝑥) ∗ 𝑓2 (𝑥)} = 𝐿 ̃ 𝛼 {𝑓2 (𝑥)} . (19)
̃ 𝛼 {𝑢} = lim 𝐿
𝐿 ̃ 𝛼 {𝑢𝑛 } . (28)
From (19) we obtain 𝑛→∞
̃ 𝛼 {𝑢𝑛+1 (𝑥)}
𝐿
From (28) we arrive at
̃ 𝛼 {𝑢𝑛 (𝑥)}
=𝐿 (20)
𝛼
̃ −1 {𝐿
𝑢 = lim 𝐿 ̃ 𝛼 𝑢𝑛 } . (29)
𝛼
̃ 𝛼 { 𝜆(𝑥) } 𝐿
𝑛→∞
+𝐿 ̃ 𝛼 {𝐿 𝛼 𝑢𝑛 (𝑥) − 𝑅𝛼 𝑢𝑛 (𝑥)} .
Γ (1 + 𝛼)
By the local fractional variation [23, 27, 29, 30], we obtain 4. The Nondifferentiable Solutions
̃ 𝛼 {𝑢𝑛+1 (𝑥)}}
𝛿𝛼 {𝐿 In this section, we discuss the one-dimensional nonhomoge-
neous heat equations arising in fractal heat flow.
̃ 𝛼 {𝑢𝑛 (𝑥)}}
= 𝛿𝛼 {𝐿 (21) Example 1. The nonhomogeneous local fractional heat equa-
𝛼
̃ 𝛼 { 𝜆(𝑥) } 𝐿
+ 𝛿𝛼 {𝐿 ̃ 𝛼 {𝐿 𝛼 𝑢𝑛 (𝑥) − 𝑅𝛼 𝑢𝑛 (𝑥)}} , tion with the nondifferentiable sink term is presented as
Γ (1 + 𝛼) follows:
which leads to
𝜕𝛼 𝑇 (𝑥, 𝑡) 𝜕2𝛼 𝑇 (𝑥, 𝑡)
𝛿𝛼 {𝐿̃ 𝛼 {𝑢𝑛+1 (𝑥)}} −
𝜕𝑡𝛼 𝜕𝑥2𝛼
̃ 𝛼 {𝑢𝑛 (𝑥)}}
= 𝛿𝛼 {𝐿 𝑥𝛼 𝐸𝛼 (−𝑡𝛼 )
(22) =− , 0 < 𝑥 < 1, 0 < 𝑡 ≤ 1, 0 < 𝛼 ≤ 1,
𝛼 Γ (1 + 𝛼)
̃𝛼 { 𝜆(𝑥) ̃ 𝛼 {𝐿 𝛼 𝑢𝑛 (𝑥)}}} = 0.
+𝐿 } {𝛿𝛼 {𝐿 (30)
Γ (1 + 𝛼)
From (22) we have subject to the initial-boundary value conditions
̃ 𝛼 {𝐿 𝛼 𝑢𝑛 (𝑥)}}
𝛿𝛼 {𝐿
𝜕𝛼 𝑇 (0, 𝑡)
̃ 𝛼 {𝑢𝑛 (𝑥)} − 𝑠𝛼 𝑢𝑛 (0) − 𝑢(𝛼) (0)}
= 𝛿𝛼 {𝑠2𝛼 𝐿 (23) = 𝐸𝛼 (−𝑡𝛼 ) ,
𝑛 𝜕𝑥𝛼 (31)
2𝛼 𝛼 ̃
= 𝑠 𝛿 𝐿 𝛼 {𝑢𝑛 (𝑥)} 𝑇 (0, 𝑡) = 0.
4 Mathematical Problems in Engineering
From (26) we obtain the local fractional iteration algorithm: Making use of (32) and (35), the third approximate term reads
as follows:
𝐿 ̃ 𝛼 {𝑇𝑛 (𝑥, 𝑡)} − 1 𝐿
̃ 𝛼 {𝑇𝑛+1 (𝑥, 𝑡)} = 𝐿 ̃ 𝛼
𝑠2𝛼 𝛼 ̃ 𝛼 {𝑇2 (𝑥, 𝑡)} − 1 𝜕 𝑇2 (𝑠, 𝑡)
̃ 𝛼 {𝑇3 (𝑥, 𝑡)} = 2𝐿
𝐿
𝑠2𝛼 𝜕𝑡𝛼
𝜕𝛼 𝑇𝑛 (𝑥, 𝑡) 𝜕2𝛼 𝑇𝑛 (𝑥, 𝑡)
⋅ {( −
𝜕𝑡𝛼 𝜕𝑥2𝛼 𝑇2(𝛼) (0, 𝑡) 𝐸𝛼 (−𝑡𝛼 )
− −
𝑠2𝛼 𝑠4𝛼
𝑥𝛼 𝐸𝛼 (−𝑡𝛼 ) (36)
+ )} 2𝐸 (−𝑡𝛼 ) 𝐸𝛼 (−𝑡𝛼 ) 𝐸𝛼 (−𝑡𝛼 )
Γ (1 + 𝛼) = 𝛼 2𝛼 − −
𝑠 𝑠4𝛼 𝑠2𝛼
̃ 𝛼 {𝑇𝑛 (𝑥, 𝑡)} − 1
=𝐿 𝐸𝛼 (−𝑡𝛼 ) 𝐸𝛼 (−𝑡𝛼 )
𝑠2𝛼 − = .
𝑠4𝛼 𝑠2𝛼
𝜕𝛼 𝑇𝑛 (𝑠, 𝑡) ̃ 𝛼 {𝑇𝑛 (𝑥, 𝑡)}
⋅{ − 𝑠2𝛼 𝐿 (32)
𝜕𝑡𝛼 From (32) and (36), the fourth approximate term can be
written as follows:
+ 𝑠𝛼 𝑇𝑛 (0, 𝑡) + 𝑇𝑛(𝛼) (0, 𝑡)
𝛼
𝐸 (−𝑡𝛼 ) 𝐿 ̃ 𝛼 {𝑇3 (𝑥, 𝑡)} − 1 𝜕 𝑇3 (𝑠, 𝑡)
̃ 𝛼 {𝑇4 (𝑥, 𝑡)} = 2𝐿
+ 𝛼 2𝛼 } 𝑠2𝛼 𝜕𝑡𝛼
𝑠
𝑇3(𝛼) (0, 𝑡) 𝐸𝛼 (−𝑡𝛼 )
𝛼 − −
̃ 𝛼 {𝑇𝑛 (𝑥, 𝑡)} − 1 𝜕 𝑇𝑛 (𝑠, 𝑡)
= 2𝐿 𝑠2𝛼 𝑠4𝛼
𝑠2𝛼 𝜕𝑡𝛼 (37)
2𝐸 (−𝑡𝛼 ) 𝐸𝛼 (−𝑡𝛼 ) 𝐸𝛼 (−𝑡𝛼 )
𝑇𝑛(𝛼) (0, 𝑡) 𝐸𝛼 (−𝑡𝛼 ) = 𝛼 2𝛼 − −
− − , 𝑠 𝑠4𝛼 𝑠2𝛼
𝑠2𝛼 𝑠4𝛼
𝐸𝛼 (−𝑡𝛼 ) 𝐸𝛼 (−𝑡𝛼 )
− = .
where the initial value is given as 𝑠4𝛼 𝑠2𝛼
𝐸𝛼 (−𝑡𝛼 ) Making the best of (32) and (36), we can write the fifth
̃ 𝛼 {𝑇0 (𝑥, 𝑡)} =
𝐿 . (33) approximate term as
𝑠2𝛼
𝛼
Using (32), we have the first approximation: ̃ 𝛼 {𝑇4 (𝑥, 𝑡)} − 1 𝜕 𝑇4 (𝑠, 𝑡)
̃ 𝛼 {𝑇5 (𝑥, 𝑡)} = 2𝐿
𝐿
𝑠2𝛼 𝜕𝑡𝛼
𝛼
𝐿 ̃ 𝛼 {𝑇0 (𝑥, 𝑡)} − 1 𝜕 𝑇0 (𝑠, 𝑡)
̃ 𝛼 {𝑇1 (𝑥, 𝑡)} = 2𝐿 𝑇4(𝛼) (0, 𝑡) 𝐸𝛼 (−𝑡𝛼 )
𝑠2𝛼 𝜕𝑡𝛼 − −
𝑠2𝛼 𝑠4𝛼
𝑇0(𝛼) (0, 𝑡) 𝐸𝛼 (−𝑡𝛼 ) (38)
− − 2𝐸 (−𝑡𝛼 ) 𝐸𝛼 (−𝑡𝛼 ) 𝐸𝛼 (−𝑡𝛼 )
𝑠2𝛼 𝑠4𝛼 = 𝛼 2𝛼 − −
(34) 𝑠 𝑠4𝛼 𝑠2𝛼
2𝐸 (−𝑡𝛼 ) 𝐸𝛼 (−𝑡𝛼 ) 𝐸𝛼 (−𝑡𝛼 )
= 𝛼 2𝛼 + − 𝐸𝛼 (−𝑡𝛼 ) 𝐸𝛼 (−𝑡𝛼 )
𝑠 𝑠4𝛼 𝑠2𝛼 − = .
𝑠4𝛼 𝑠2𝛼
𝐸𝛼 (−𝑡𝛼 ) 𝐸𝛼 (−𝑡𝛼 )
− = . Hence, we obtain the final term given as
𝑠4𝛼 𝑠2𝛼
𝛼
̃ 𝛼 {𝑇𝑛 (𝑥, 𝑡)} − 1 𝜕 𝑇𝑛 (𝑠, 𝑡)
= 2𝐿
𝑠2𝛼 𝜕𝑡𝛼
1.5
𝑇𝑛(𝛼) (0, 𝑡) cos𝛼 (𝑡𝛼 )
− + .
𝑠2𝛼 𝑠4𝛼
1
(43)
T (x, t)
Conflict of Interests
0 The authors declare that they have no conflict of interests in
this paper.
−0.1
−0.2 Acknowledgments
T (x, t)
−0.3
The work was supported by the Natural Science Foundation
−0.4 of Jiangsu Colleges and Universities (no. KK-12058) and
Jiangsu Province R&D Institute of Marine Resource (no.
−0.5
1 JSIUMR 201210).
1
0.5 0.8
0.6 References
t 0.4
0.2 x
0 0
[1] A. A. Kilbas, H. M. Srivastava, and J. J. Trujillo, Theory and
Figure 2: The nondifferentiable solution of the nonhomogeneous Applications of Fractional Differential Equations, vol. 204 of
local fractional heat equation with the nondifferentiable sink term North-Holland Mathematics Studies, Elsevier Science, Amster-
for 𝛼 = ln 2/ ln 3. dam, The Netherlands, 2006.
[2] B. J. West, M. Bologna, and P. Grigolini, Physics of Fractal Ope-
rators, Institute for Nonlinear Science, Springer, 2003.
In view of (43) and (47), the fifth approximate term is pre- [3] J. Klafter, S. C. Lim, and R. Metzler, Fractional Dynamics: Recent
sented as Advances, World Scientific, 2012.
𝛼 [4] D. Baleanu, J. A. T. Machado, and A. C. Luo, Fractional Dyna-
𝐿 ̃ 𝛼 {𝑇4 (𝑥, 𝑡)} − 1 𝜕 𝑇4 (𝑠, 𝑡)
̃ 𝛼 {𝑇5 (𝑥, 𝑡)} = 2𝐿 mics and Control, Springer, New York, NY, USA, 2012.
𝑠2𝛼 𝜕𝑡𝛼
[5] A. B. Alkhasov, R. P. Meilanov, and M. R. Shabanova, “Heat con-
𝑇4(𝛼) (0, 𝑡) cos𝛼 (𝑡𝛼 ) duction equation in fractional-order derivatives,” Journal of
− + Engineering Physics and Thermophysics, vol. 84, no. 2, pp. 332–
𝑠2𝛼 𝑠4𝛼 341, 2011.
(48)
2sin𝛼 (𝑡𝛼 ) cos𝛼 (𝑡𝛼 ) sin𝛼 (𝑡𝛼 ) [6] J. M. Angulo, M. D. Ruiz-Medina, V. V. Anh, and W. Grecksch,
= − − “Fractional diffusion and fractional heat equation,” Advances in
𝑠2𝛼 𝑠4𝛼 𝑠2𝛼
Applied Probability, vol. 32, no. 4, pp. 1077–1099, 2000.
𝐸𝛼 (−𝑡𝛼 ) sin𝛼 (𝑡𝛼 ) [7] Y. Z. Povstenko, “Fractional heat conduction equation and asso-
+ = .
𝑠4𝛼 𝑠2𝛼 ciated thermal stress,” Journal of Thermal Stresses, vol. 28, no. 1,
pp. 83–102, 2004.
Hence, we finally have
[8] H. M. Youssef, “Theory of fractional order generalized thermoe-
sin𝛼 (𝑡𝛼 ) lasticity,” Journal of Heat Transfer, vol. 132, no. 6, 7 pages, 2010.
̃ 𝛼 {𝑇𝑛 (𝑥, 𝑡)} =
𝐿 (49)
𝑠2𝛼 [9] M. A. Ezzat and A. S. El-Karamany, “Fractional order theory of
a perfect conducting thermoelastic medium,” Canadian Journal
so that the exact solution of nonhomogeneous local fractional of Physics, vol. 89, no. 3, pp. 311–318, 2011.
heat equation with nondifferentiable source term is
[10] M. A. Ezzat, “State space approach to thermoelectric fluid with
̃ −1 {𝐿
𝑇 (𝑥, 𝑡) = lim 𝐿 ̃ 𝛼 {𝑇𝑛 (𝑥, 𝑡)}} fractional order heat transfer,” Heat and Mass Transfer/Waerme-
𝛼
𝑛→∞ und Stoffuebertragung, vol. 48, no. 1, pp. 71–82, 2012.
(50)
𝑥𝛼 [11] H. H. Sherief, A. M. A. El-Sayed, and A. M. Abd El-Latief, “Frac-
= sin (𝑡𝛼 ) . tional order theory of thermoelasticity,” International Journal of
Γ (1 + 𝛼) 𝛼 Solids and Structures, vol. 47, no. 2, pp. 269–275, 2010.
For the fractal dimension 𝛼 = ln 2/ ln 3, the plot of the non- [12] L. Vázquez, J. J. Trujillo, and M. P. Velasco, “Fractional heat
differentiable solution of the nonhomogeneous local frac- equation and the second law of thermodynamics,” Fractional
tional heat equation with the nondifferentiable source term Calculus and Applied Analysis, vol. 14, no. 3, pp. 334–342, 2011.
is shown in Figure 2. [13] J. Hristov, “An inverse stefan problem relevant to boilover: heat
balance integral solutions and analysis,” Thermal Science, vol. 11,
no. 2, pp. 141–160, 2007.
5. Conclusions [14] J. Hristov, “A note on the integral approach to non-linear heat
conduction with Jeffrey’s fading memory,” Thermal Science, vol.
At the present work, the nonhomogeneous heat equations 17, no. 3, pp. 733–737, 2013.
arising in the fractal heat flow were investigated. The local [15] K. Davey and R. Prosser, “Analytical solutions for heat transfer
fractional Laplace variational iteration method was applied on fractal and pre-fractal domains,” Applied Mathematical Mod-
to obtain the nondifferentiable solutions for the nonhomo- elling, vol. 37, no. 1-2, pp. 554–569, 2013.
geneous local fractional heat equations with the nondiffer- [16] M. Ostoja-Starzewski, “Towards thermoelasticity of fractal
entiable source and sink terms. Finally, the graphs of the media,” Journal of Thermal Stresses, vol. 30, no. 9-10, pp. 889–
obtained solutions are also shown. 896, 2007.
Mathematical Problems in Engineering 7
Research Article
Numerical Simulation and Experimental Research on
Coal Ash Collecting and Grading System
Yuanhua Xie, Xianjin Li, Liwei Wang, Hong Yu, Bing Bai, Zhizhou Xu, and Tong Zhu
School of Mechanical Engineering and Automation, Northeastern University, 3-11 Wenhua Road, Heping District,
Shenyang 110004, China
Copyright © 2014 Yuanhua Xie et al. This is an open access article distributed under the Creative Commons Attribution License,
which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
The grading separation of coal ash can not only increase its economic value but also decrease its pollution to environment. Based
on the jet-attracting flow technology and the gas-solid two-phase flow theory, the force and motion of coal ash particles in airflow
were studied firstly. Focused on single coal ash particle, Matlab software was used to simulate the force conditions and separation
parameters of various diameter coal ash particles in airflow. Fluent software was used to simulate the nozzle fluidization domain
shape and to determine optimal jet flux. According to the theoretical results, a coal ash collecting and grading system was developed.
Using the separation efficiency as the evaluation index, the optimal experiment parameters of jet flux, attracting flux, and separation
time were obtained. At last, the calculated results and experimental results of coal ash particles median diameter from the first
grading separation exit under various attracting fluxes were compared. The reasons that could cause the errors were discussed. This
study has significant practical meaning and application value on coal ash collecting and grading separation.
1. Introduction bear its weight. The separated lighter material can be further
separated by cyclone separator [7]. Some problems in the
Coal ash is the solid waste of burning coal. Except for dust present airflow classification process are showed as follows
air pollution, coal ash can result in groundwater pollution [4]. The high fineness requirement to coal ash increases the
and soil pollution due to its compositions of trace amounts design difficulty of coal ash grading device. The performance
of heavy metals and tiny amounts of radioactive substances. and efficiency of separation device are severely affected by the
Coal ash can also cause geological disasters [1]. However, coal seal and abrasion. How can we combine the theory analysis
ash contains many kinds of useful substances, such as glass and experimental test? How can we solve the contradiction
microbead, high-aluminum powder, and carbon powder. between high work capacity and high grading fineness? How
Collecting and separating these substances effectively can can we solve low grading efficiency and low grading accuracy
greatly increase the economic value of coal ash, which can under high work capacity?
turn the waste into treasure and has more important meaning Based on the jet-attracting flow technology [8–11] and
on environment protection [2]. gas-solid two-phase flow theory [12–14], a new coal ash
At present, the coal ash grading separation methods collecting and grading device was developed in this study.
mainly include sieving, airflow classification, froth flotation, The new device could achieve the grading separation of coal
electrostatic separation, and magnetic separation [3, 4]. ash while collecting and transporting. The force and motion
For the separation of high fineness coal ash, the airflow of coal ash particles in airflow were studied in theory. Then
classification is often used [5, 6]. The basic principle of airflow the force and separation condition of various diameter coal
classification is to separate the lighter material out of the ash particles in airflow were simulated by Matlab software.
separation device upward or to take the lighter material to rel- The Fluent software was also used to simulate the shape of
atively distant place along the horizontal direction by airflow. nozzle fluidization domain and to determine the optimal
The heavier material will settle because the airflow cannot jet flux. Finally, the coal ash collecting and grading device
2 Mathematical Problems in Engineering
was used for the relevant experiments to study the relations was 1900∼2900 kg/m3 [14]; 𝑉𝑐 was the volume of coal ash
between the separation efficiency and jet flux, attracting flux particles (m3 ); 𝑑𝑐 was the equivalent volume diameter of coal
or separation time. The calculated and experimental values ash particles (m).
of coal ash particle median diameter in the first grading
separation exit under different attracting flux were compared.
2.1.2. Buoyancy. The buoyancy of coal ash particles in airflow
The reasons causing the errors were also discussed.
was shown as follows:
Single particle dynamics (SPD) is the simplest method in In this equation, 𝐹𝑓 was buoyancy (N); 𝜌𝑔 was the density of
studying the suspension two-phase flow. In SPD method, the airflow (kg/m3 ).
airflow field is regarded as known field, without considering
the particle fluctuation and the influence of particles on fluid 2.1.3. Pressure Gradient Force. For a single particle (or sus-
flow, but only considering the force and motion of a single pension system with very low concentration), the fluid flow
particle that is uncorrelated with other particles in flow field was not affected by the existing small particles. For the fluid
[15]. phase, the following equation existed:
4.296 × 10−4
𝐶𝐷 =
[𝑑𝑐 (V𝑔 − V𝑐 )] Figure 1: Velocity changes in turbulence.
(7)
6
+ + 0.4.
[1 + √(𝑑𝑐 (V𝑔 − V𝑐 )) / (1.79 × 10−5 )] After final transformation, it could be expressed as fol-
lows:
𝑑V𝑐 2
In this equation, V𝑐 was the particle velocity (m/s). =
𝑑𝑡 2𝜌𝑐 + 𝜌𝑔
10−5 10−7
Traction resistance (N)
10−8
10−9 10−12
0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0
Time (s) Time (s)
d = 250 𝜇m d = 250 𝜇m
d = 50 𝜇m d = 50 𝜇m
d = 10 𝜇m d = 10 𝜇m
(a) Traction resistance (b) False mass force
Figure 2: Changes of traction resistance and false mass force with time.
relationship between coal ash particle velocity and time Table 1: Gravity, buoyancy, and pressure gradient force of various
needed to use the numerical analysis computer software. diameter coal ash particles.
0.0 0.2
−0.1
−0.2 0.0
−0.3
Velocity (m/s)
−0.2
Velocity (m/s)
−0.4
−0.5 −0.4
−0.6
−0.7 −0.6
−0.8 −0.8
−0.9
−1.0 −1.0
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5
Time (s) Time (s)
d = 250 𝜇m d = 250 𝜇m
d = 50 𝜇m d = 50 𝜇m
d = 10 𝜇m d = 10 𝜇m
(a) V𝑔𝑝 = 0.06 m/s (b) V𝑔𝑝 = 0.3 m/s
0.4 0.7
0.6
0.2 0.5
0.4
Velocity (m/s)
d = 250 𝜇m d = 250 𝜇m
d = 50 𝜇m d = 50 𝜇m
d = 10 𝜇m d = 10 𝜇m
(c) V𝑔𝑝 = 1.38 m/s (d) V𝑔𝑝 = 1.98 m/s
1.2 2.5
1.0 2.0
0.8
Velocity (m/s)
Velocity (m/s)
1.5
0.6
0.4 1.0
0.2 0.5
0.0
0.0
−0.2
−0.4 −0.5
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5
Time (s) Time (s)
d = 250 𝜇m d = 250 𝜇m
d = 50 𝜇m d = 50 𝜇m
d = 10 𝜇m d = 10 𝜇m
(e) V𝑔𝑝 = 2.46 m/s (f) V𝑔𝑝 = 3.00 m/s
the speed of 0.2 m/s. However, 250 𝜇m particles were still in a velocity was 2.46 m/s, 250 𝜇m particles reached a suspending
quiescent state in theory (a small amount of 250 𝜇m particles steady state. And the velocity change range was still about
had already been drawn out by the airflow in the actual 0.2 m/s. However, the other two particles were in completely
situation). When the airflow velocity reached 1.98 m/s, 10 𝜇m unstable state. So, for the three different diameters of coal
particles began to present an erratic state. When the airflow ash particles (10 𝜇m, 50 𝜇m, and 250 𝜇m), the separation
6 Mathematical Problems in Engineering
1.0 3.0
2.5
0.5
2.0
1.5
0.0
1.0
−0.5 0.5
0.0
−1.0 −0.5
0.0 0.5 1.0 1.5 2.0 2.5 3.0 0.0 0.5 1.0 1.5 2.0 2.5 3.0
Attracting flow velocity (m/s) Attracting flow velocity (m/s)
d = 250 𝜇m d = 250 𝜇m
d = 50 𝜇m d = 50 𝜇m
d = 10 𝜇m d = 10 𝜇m
Figure 4: Relationships between coal ash particle average velocity Figure 5: Relationships between coal ash particle velocity pulsating
and airflow velocity. amplitude and airflow velocity.
airflow velocities (attracting velocity) in theory were 0.3 m/s, vertical direction. But, relatively speaking, it had a better
1.38 m/s, and 2.46 m/s. From Figure 4, it could be found that development in horizontal direction than in vertical direc-
the movement velocity of coal ash particles increased about tion. Though coal ash could be efficiently fractionated under
0.1∼0.5 m/s when the airflow velocity increased 0.5 m/s. The this flux (proved in Section 4.2), the throughput of coal ash
increasing amplitude of coal ash particle velocity increased fractionation was very low.
with the increase in airflow velocity. The increasing amplitude
As shown in Figure 6(b), when jet flux was 1.5 m3 /h, the
of big coal ash particles velocity was slightly larger than that
fluidization domain achieved the best shape. This type of
of the small particles.
fluidization domain developed well in both horizontal and
Figure 5 showed the relationships between the pulsating
vertical direction. The fluidized shape of coal ash was rather
amplitude of coal ash particle velocity and airflow velocity.
clear. On the premise of meeting higher separation efficiency,
The pulsating amplitude of coal ash particle velocity increased
this jet flux could achieve a higher throughput.
with the increase in airflow velocity. When the airflow
As shown in Figure 6(c), when jet flux was 3 m3 /h, the
velocity was less than 2.0 m/s, the pulsating amplitude of
fluidization domain still developed well in both horizontal
coal ash particle velocity was not obvious. When the airflow
and vertical direction. Though the jet flux was large, the
velocity was more than 2.0 m/s, the pulsating amplitude of
fractionation efficiency and the throughput were severely
three coal ash particles increased obviously. In this situation,
reduced. As can be seen in Figure 6(c), particle concentration
smaller particles presented larger pulsating amplitude and
of coal ash near the jet was particularly low under the strong
irregular local velocity. On the contrary, larger particles
jet airflow. Thereby, particles extracted by the attracting flow
showed smaller pulsating amplitude and regular local veloc-
were very little. At the same time, due to the large jet veloc-
ity. When the airflow velocity was more than 2.0 m/s, the
ity, the collision and interference among coal ash particles
average pulsating amplitudes were 2 m/s, 1 m/s, and 0.3 m/s
increased and the separation process was full of randomness.
for the particles of 10 𝜇m, 50 𝜇m, and 250 𝜇m.
As a result, fractionation efficiency and throughput were
reduced.
3.3. Effects of Jet Flux on Coal Ash Reclaiming. The coal The simulation results showed that the fluidization
ash collecting and grading system designed in this study domain of coal ash had a best fluidized state when jet flux
worked as the following principle. The jet airflow firstly was 1.5 m3 /h. So jet flux of 1.5 m3 /h was used as a reference in
stirred coal ash to form a fluidization domain. And then the the following experiments.
attracting pipe would extract coal ash. The coal ash would
be graded and separated in the device. As a result, the state
of fluidization domain played an important role in coal ash 4. Experimental Study of Coal
collecting and grading process. In the study, Fluent software Ash Collecting and Grading
under Eulerian-Eulerian model was used to simulate two-
phase flow. The fluidization domain shape of coal ash under The structure of coal ash collecting and grading system
different jet flux was simulated to determine the quality of designed in this study was shown in Figure 7. The system
fluidization performance and then determine the required jet was composed of air jet pipe, dust cover, suspension cham-
flux. Figure 6 shows the fluidization domain shape of coal ash bers, settling chambers, blinds, triangular storage rooms,
under different jet flux. The colorful coordinate presented the discharging pipes, attracting pipe, and other components.
volume fraction of coal ash. When worked, the airflow went into system from intake pipe
As shown in Figure 6(a), when jet flux was 0.05 m3 /h, 1, erupted from jet pipe 2, and then stirred up the coal ash
the fluidization domain developed in both horizontal and and formed a fluidization domain. Dust cover 12 limited
Mathematical Problems in Engineering 7
6.29e − 01 6.25e − 01
5.98e − 01 5.94e − 01
5.67e − 01 5.62e − 01
5.35e − 01 5.31e − 01
5.04e − 01 5.00e − 01
4.72e − 01 4.69e − 01
4.41e − 01 4.37e − 01
4.09e − 01 4.06e − 01
3.78e − 01 3.75e − 01
3.46e − 01 3.44e − 01
3.15e − 01 3.12e − 01
2.83e − 01 2.81e − 01
2.52e − 01 2.50e − 01
2.20e − 01 2.19e − 01
1.89e − 01 1.87e − 01
1.57e − 01 1.56e − 01
1.26e − 01 1.25e − 01
9.44e − 02 9.37e − 02
6.29e − 02 6.25e − 02
3.15e − 02 3.12e − 02
0.00e − 00 0.00e − 00
(a) Jet flux of 0.05 m3 /h (b) Jet flux of 1.5 m3 /h
6.23e − 01
5.92e − 01
5.61e − 01
5.30e − 01
4.99e − 01
4.67e − 01
4.36e − 01
4.05e − 01
3.74e − 01
3.43e − 01
3.12e − 01
2.80e − 01
2.49e − 01
2.18e − 01
1.87e − 01
1.56e − 01
1.25e − 01
9.35e − 02
6.23e − 02
3.12e − 02
0.00e − 00
(c) Jet flux of 3 m3 /h
Figure 6: Fluidization domain shape of coal ash under different jet flux.
the fluidized coal ash to stretch all around. The pressure in from discharge pipes (4, 7, and 10) at the particle size order
dust cover 12 would rise with the continuous entering of jet from big to small. The last remaining coal ash particles with
airflow. At the same time, the attracting pipe 21 showed an small size would be extracted through attracting pipe 21. Coal
attracting effect. The fluidized coal ash would move upwards ash particles extracted from discharge pipes (4, 7, and 10)
and enter into stair-like pipes. With the increasing diameter and attracting pipe 21 were separated into gas and solid by a
of the stair-like pipes, the coal ash flow velocity decreased. cyclone separator. Finally, this system achieved the collecting,
Meanwhile, particles with large density or big size carried grading, and transferring of coal ash simultaneously. Figure 8
by airflow could not move along the pipe any more. They is the schematic diagram of the whole separator system.
would suspend in suspension chambers (14, 16, 18, and 20) The simulation results in Section 3.2 showed that the
or settle along the pipe wall of settling chambers (15, 17, attracting airflow velocities to keep the coal ash particles
and 19). Sediment particles gathered in triangular storage of 10 𝜇m, 50 𝜇m, and 250 𝜇m at steady conditions were
0.3 m/s, 1.38 m/s, and 2.46 m/s, respectively. According to
rooms (5, 8, and 11) through blinds (3, 6, and 9). The
the continuity equation, the airflow in separation unit was
gathered coal ash particles were extracted slowly through
constant and the total airflow through each level was equal.
discharge pipes (4, 7, and 10) under the negative pressure. So the continuity equation was
This process accelerated the settling velocity of particles near
the settling chambers (15, 17, and 19) wall and the moving 𝐷 2
velocity of suspended particles to the pipe wall in suspension 𝑄1 = 𝑄2 = 𝑄3 = V𝜋( ). (13)
2
chambers (14, 16, 18, and 20). Three grading areas with
progressively increased diameter were designed in stair-like In this equation, 𝑄1 was the airflow flux in the first stage
pipe to separate coal ash. The particles would be drawn out (m3 /s); 𝑄2 was the airflow flux in the second stage (m3 /s);
8 Mathematical Problems in Engineering
5 6
8
9
3
10
2
4
11
1
12
Figure 8: Schematic diagram of coal ash separator system. 1: air compressor; 2: airflow control valve; 3, 7, and 11: air flow meters; 4: coal ash
separator; 5 and 9: cyclone separators; 6 and 10: filter units; 8 and 12: water ring vacuum pumps.
particles were affected by several forces and the low airflow jet flux and attracting flux, the coal ash separation efficiency
velocity could not attract large particle. Many small particles gradually increased with separation time. When separation
(≤45 𝜇m) were extracted out in the first grading stage and time reached 10 s, the separation efficiency was about 73% and
resulted in low separation efficiency. With the increase in kept steady. When separation time reached about 25 s, the
attracting flux, the updraft velocity increased. Many big separation efficiency presented another small increase, fol-
particles could be lift up and a clear increase in separation lowed by a sharp dropping to about 60%. This phenomenon
efficiency occurred. When attracting flux exceeded 70 m3 /h, could be explained as follows. With the increase in separation
the overlarge flux resulted in deterioration of particle dis- time, the coal ash concentration gradually increased, and
persibility and serious particle agglomeration. A lot of coal separation efficiency also gradually increased. Because of the
ash particles were pumped directly into the cyclone separa- fixing position of grading device, the fluidization domain of
tors and filter units without grading. Due to the instability coal ash gradually becomes a stable concave shape; therefore,
of particle agglomeration effects and the violent particle the separation efficiency was in a relatively stable value. With
collision, the separation effect of each grading stage dropped, the continuing of grading process, the amount of coal ash
resulting in the decrease of separation efficiency. In this study, reduced and less coal ash could be attracted into grading
the optimal attracting flux was 70 m3 /h. device, resulting in a small increase and a sharp decrease
in separation efficiency. In practical operation, the grading
4.4. Effect of Separation Time on Separation Efficiency. The device would do the horizontal or vertical movement to
attracting flux was set to 70 m3 /h and the jet flux was set achieve the stable separation efficiency and the continuous
to 1.5 m3 /h. The effects of separation time on separation coal ash collecting and grading effect. In this study, the
efficiency were examined, seen in Figure 13. Under a certain optimal separation time was 15 s.
10 Mathematical Problems in Engineering
100 10.0
90 9.0 Size (𝜇m) Cumu. (%)
80 8.0 1.179 0
70 7.0 2.389 0.54
Cumu. (%)
Diff. (%)
60 6.0 4.841 1.67
Diff. (%)
60 3.6 1.477 2.57
3.0 2.838 4.89
50
5.453 12.21
40 2.4
10.480 21.20
30 1.8 20.130 38.81
20 1.2 38.690 64.20
10 0.6 74.360 90.63
0 0.0 142.900 100.00
0.01 0.1 1 10 100 1000
Size (𝜇m)
(b) In the second stage
100 6.0
90 5.4 Size (𝜇m) Cumu. (%)
80 4.8 0.315 0
70 4.2 0.597 1.94
Cumu. (%)
Diff. (%)
60 3.6 1.132 6.44
50 3.0 2.147 12.26
40 2.4 4.070 23.18
30 1.8 7.718 39.42
1.2 14.630 63.15
20
27.740 82.26
10 0.6
52.600 96.41
0 0.0 99.750 100.00
0.01 0.1 1 10 100 1000
Size (𝜇m)
(c) In the third stage
Figure 10: Particle size distribution of coal ash (jet flux = 1.5 m3 /h) in different grading stages.
4.5. Comparison of Experimental and Theoretical Average of theoretical results. When the attracting velocity was less
Separated Particle Size. The median diameter of separated than 1.5 m/s, the experimental results were slightly larger
coal ash (average particle size) was an important indicator for than the theoretical results. When the attracting velocity
evaluating grading perfection in the coal ash grading process. was greater than 1.5 m/s, the experimental values were less
The jet flux was selected as 1.5 m3 /h and the separation time than the theoretical values, and the deviation became larger
was 15 s. The change of the average separated particle size with the attracting velocity increase. This phenomenon was
of the first separation exit was measured under different mainly caused by the different particle concentration in
attracting velocity (controlled by attracting flux) and was the grading device. When the attracting velocity was small,
compared with the theoretical value (in Figure 14). The mutual interference between the particles was not serious
reason to study the first separation exit was that the first and the particles were graded according to their own force.
separation exit was proved to be the smallest affected by tiny But a lot of small particles were in irregular movement and
particles less than 50 𝜇m and to have the most shallow particle were likely to be directly extracted out from the attracting
size distribution. pipe. Therefore, the experimental values were larger than the
As can be seen in Figure 14, under a certain jet flux, theoretical values. When the attracting velocity became large,
separated coal ash particle size gradually increased with particle concentration increased within the grading device.
the increase in attracting velocity. In the change process of Agglomeration and collision between particles enhanced.
attracting velocity from 0 to 3 m/s, the average separated Many small particles were mixed into the first discharged pipe
particle size gradually increased from 20 𝜇m to 260 𝜇m. The and then resulted in the relatively small average experimental
experimental results were consistent with the basic trend separated particle size.
Mathematical Problems in Engineering 11
80
70 400
300
60
200
50
100
40
0.0 0.5 1.0 1.5 2.0 2.5 3.0 0
0.0 0.5 1.0 1.5 2.0 2.5 3.0
Jet flux (m3 /h) Attracting velocity (m/s)
Figure 11: Separation efficiencies under different jet fluxes. Experimental values
Theoretical values
70
(2) The force of various sizes of coal ash particles was
60 obtained through Matlab numerical simulations. By
solving the coal ash particle motion mathematical
50 model using Matlab, the separation conditions of
different sizes of coal ash particles were simulated.
40
And the theoretical airflow velocities of different par-
0 20 40 60 80 100 120 140 ticle size classification required were also determined.
Attracting flux (m3 /h) For coal ash particles of 10 𝜇m, 50 𝜇m, and 250 𝜇m,
their theoretical classification airflow velocities were
Figure 12: Separation efficiencies under different attracting fluxes. 0.3 m/s, 1.38 m/s, and 2.46 m/s, respectively.
(3) The shape of nozzle fluidization domain under dif-
90
ferent jet flux was simulated using Fluent software.
80
Separation efficiency (%)
(1) In the base of jet-attracting flow technology and gas- Conflict of Interests
solid two-phase flow theory, the force and motion
of coal ash particles in gas flow were analyzed. The authors declare that there is no conflict of interests
The particle motion could be simplified to a one- regarding the publication of this paper.
dimensional vertical movement. The following force
as gravity, buoyancy, pressure gradient force, traction Acknowledgments
resistance, and false mass force had influence on
particle movement. A mathematical model of coal ash This study was jointly supported by National Natural Science
particle motion was established. Foundation of China (nos. 21107011, 51178098, and 31371873)
12 Mathematical Problems in Engineering
and the Fundamental Research Funds for the Central Univer- [16] K. F. Cen and J. R. Fan, “The analysis of the forces acting on coal
sities of China (nos. N100303006 and N130403001). particles and the trajectories in the gas flow,” Journal of Zhejiang
University, vol. 21, no. 6, pp. 1–11, 1987.
[17] D. Y. Liu, Two-Phase Hydrodynamics, Higher Education Press,
References Beijing, China, 1993.
[1] L. Sun, “Hazard and treatment of haze weather,” Environmental [18] M. Nakagawa, T. Matumi, H. Takeuchi, and N. Kokubo, “Mixing
Science and Management, vol. 37, no. 10, pp. 71–75, 2012. of the confined jet of mist flow,” JSME International Journal B:
Fluids and Thermal Engineering, vol. 39, no. 2, pp. 381–386, 1996.
[2] G. Q. Han, J. W. Li, B. Du, J. Huang, and R. J. Wang, “Present
status and future development of comprehensive utilization of [19] Y. S. Wang, Numerical Calculation and Experimental Study of
coal ash,” China Resource Comprehensive Utilization, vol. 24, no. Particle Motion in Jet Flow Sand Water Collectortor, Northeast-
7, pp. 12–14, 2006. ern University, Shenyang, China, 2010.
[3] Q. P. Li and G. X. Shao, “evelopment of air classification tech- [20] Y. Q. Xiong, M. Y. Zhang, Z. L. Yuan, and D. M. Xue, “Three-di-
nology,” Chemical Equipment Technology, vol. 23, no. 5, pp. 11– mensional numerical simulation on gas-solid two-phase flows
15, 2002. in gas-solid injector,” Proceedings of the Chinese Society of
Electrical Engineering, vol. 25, no. 20, pp. 77–82, 2005.
[4] J. C. Li, Study of Classification Principle Based on Air-Solid
Two-Phase Flow Theory and Its Application on SLK Coal Ash [21] K. Z. Qiu, X. D. Li, J. H. Yan, Y. Chi, M. J. Ni, and K. F. Cen,
Classifier, Southwest University of Science and Technology, “The experiment research of gas-solid interaction in horizontal
Mianyang, China, 2009. flow boundary layer,” Chinese Journal of Theoretical Applied
Mechanics, vol. 31, no. 4, pp. 456–465, 1999.
[5] M. Q. Li, C. B. Wu, and Q. N. Liu, “The application of frequency-
[22] J. R. Xu and Q. Luo, “The relative motion between the two
varied turbo-classifier in the classification of coal ash,” Coal Ash
phases of the solid-liquid hydrocyclone,” Chinese Journal of
Comprehensive Utilization, vol. 1, pp. 34–36, 2004.
Nonferrous Metals, vol. 8, no. 4, pp. 132–135, 1998.
[6] F. Shi, “The art of the state of coal ash classification technology in
[23] Y. Z. Meng, MATLAB 5.X Application Skills, Science Press,
china and probe into its developing prospects,” Coal Ash China,
Beijing, China, 1999.
vol. 3, pp. 41–44, 2004.
[24] Y. Cui, MATLAB 5.3 Projects Explanation, Aerospace Industry
[7] J. W. Leonard, Coal Preparation, American, Institute of Mining,
Press, Beijing, China, 1999.
Metallurgical and Petroleum Engineers, New York, NY, USA,
4th edition, 1979. [25] Y. Shi and X. C. Lin, “Three-dimensional numerical simulation
of gas-solid injector based on Fluent,” Coal Preparation Technol-
[8] Z. Tong, N. Tsutomu, H. Yonggang, I. XYuanhua, and H. Jin, ogy, vol. 4, pp. 19–23, 2011.
“Effect of particle properties on fluidization performance of
sedimentary layer based on the spraying jet type sand collector
the fundamental sutdy of a novel water bed pollutants collector,”
in Proceedings of the 4th IEEE Conference on Industrial Electron-
ics and Applications (ICIEA ’09), pp. 1093–1096, Xi’an, China,
May 2009.
[9] T. Zhu, T. Nozaki, Y. H. Xie et al., “Dynamic model, and nu-
merical simulation of particle motion in rotation flow field of
centrifuge,” in Proceedings of the 3rd International Symposium
on Advances in Computation and Intelligence (ISICA ’08), vol.
12, pp. 470–478, Wuhan, China, 2008.
[10] T. Zhu, T. Nozaki, Y. H. Xie, J. Han, and J. Jiang, “Numerical sim-
ulation and visualization experiment of solid particle motion
affected by parameters of flow in tapered drum rotating type
separator,” in Proceedings of the 7th International Conference on
System Simulation and Scientific Computing (ICSC ’08), pp. 867–
871, Beijing, China, October 2008.
[11] T. Zhu, J. Han, T. Nozaki, Y. H. Xie, and M. Y. You, “Charac-
teristics of collecting efficiency related to spraying jet flux and
suction flux in spraying jet type collector,” in Proceedings of the
6th International Conference on Material Handling (ICMH '08),
pp. 268–272, Shanghai, China, October 2008.
[12] D. Y. Fang, Two-Phase Flow Dynamics, National University of
Defense Technology Publishing, Changsha, China, 1988.
[13] L. Z. Zhang, Y. Z. Li, and G. J. Wang, “Study on particle flow
character in gas/solid two-phase wind tunnel,” Journal of System
Simulation, vol. 19, no. 14, pp. 3200–3202, 2007.
[14] K. F. Cen, Gas-Solid Separation Theory and Technology, Zhejiang
University Press, Hangzhou, China, 1999.
[15] Y. L. MA. Numerical simulation, Numerical Simulation of High
Concentrations of Gas-Solid Two Phase Flow, Zhejiang Uni-
versity, Hangzhou, China, 2001.
Hindawi Publishing Corporation
Mathematical Problems in Engineering
Volume 2014, Article ID 627528, 9 pages
[Link]
Research Article
Simulation of Microstructure during Laser Rapid Forming
Solidification Based on Cellular Automaton
Zhi-jian Wang,1 Shuai Luo,1 Hong-wu Song,2 Wei-dong Deng,1 and Wen-yi Li1
1
School of Mechatronics Engineering, Shenyang Aerospace University, Shenyang, Liaoning 110136, China
2
Institute of Metal Research, Chinese Academy of Sciences, Shenyang, Liaoning 110016, China
Copyright © 2014 Zhi-jian Wang et al. This is an open access article distributed under the Creative Commons Attribution License,
which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
The grain microstructure of molten pool during the solidification of TC4 titanium alloy in the single point laser cladding was
investigated based on the CAFE model which is the cellular automaton (CA) coupled with the finite element (FE) method. The
correct temperature field is the prerequisite for simulating the grain microstructure during the solidification of the molten pool.
The model solves the energy equation by the FE method to simulate the temperature distribution in the molten pool of the single
point laser cladding. Based on the temperature field, the solidification microstructure of the molten pool is also simulated with the
CAFE method. The results show that the maximum temperature in the molten pool increases with the laser power and the scanning
rate. The laser power has a larger influence on the temperature distribution of the molten pool than the scanning rate. During the
solidification of the molten pool, the heat at the bottom of the molten pool transfers faster than that at the top of the molten pool. The
grains rapidly grow into the molten pool, and then the columnar crystals are formed. This study has a very important significance
for improving the quality of the structure parts manufactured through the laser cladding forming.
1. Introduction the molten pool and the material features. With the devel-
opment of the computer science, the numerical simulation
Currently, the laser cladding technology has become the new technology is a good method to analyze the growth of grains
breakthrough of the complex structure parts for aviation in the solidification of molten pool. It contributes to revealing
manufacturing process. In the laser cladding technology, the evolution of microstructures and improving the quality of
the high input heat can make the surface of the substrate the workpiece in the laser cladding forming.
and the laser cladding material melt and then the molten The researches on the microstructure simulation have
pool is formed. It avoids the disadvantages of high machin- achieved many significant accomplishments since the end of
ing allowances in the traditional process, long computer the last century [5, 6]. Several kinds of calculation methods
numerically controlled (CNC) machining time, low material about nucleation and growth of the grains during the solid-
utilization, long production cycles and high manufacturing ification of materials were presented from these researches
cost [1–3]. Figure 1 is the principle of laser cladding forming. including deterministic modeling, stochastic modeling, and
In the process of laser cladding forming, the substrate is phase-field modeling [7]. Because of the determinacy of
heated by laser beam which results in melting the thin layer the deterministic modeling, it cannot simulate the random
of metal on the substrate surface. And then the molten pool growth of grains such as random nucleation and random
is formed with the powders offered by powder feeder [4]. crystal orientation. Phase-field modeling which avoids the
The solidification of the molten pool is a very complex difficulty of tracking the complex solid-liquid interface has
and unbalanced rapid solidification process. It is difficult disadvantages that slow computing speed, inefficiency, and
to observe the growth of grains during the solidification of small simulated area [8, 9].
molten pool experimentally because of the high temperature Cellular automaton (CA) which is a kinetic and stochastic
and the rapid solidification rate affected by the size of model scattered in time, space, and state has been widely
2 Mathematical Problems in Engineering
Workpiece
2.1. Heat Transfer Model. The temperature distribution has a
Control
significant impact on the grain microstructure in the molten
systerm Workbench pool. The heat transfer model [15] used in this paper is shows
as follows:
Figure 1: Principle of the laser cladding forming. 𝜕𝑇 𝜕 𝜕𝑇 𝜕 𝜕𝑇 𝜕 𝜕𝑇
𝜌𝑐 = (𝑘 ) + (𝑘 ) + (𝑘 ) , (1)
𝜕𝑡 𝜕𝑥 𝜕𝑥 𝜕𝑦 𝜕𝑦 𝜕𝑧 𝜕𝑧
where 𝜌 is the density, 𝑐 is the specific heat, 𝑇 is the
used in solidification and recrystallization of the material temperature, 𝑡 is the time, and 𝑥, 𝑦, and 𝑧 are the coordinates
science in recent decades [7]. Based on the solidification of the calculated area, respectively.
thermodynamic, grain nucleation, and growth kinetics, it can The enthalpy is defined in this paper to solve the release
simulate the growth of the grains in the solidification of the of the latent heat during the solidification of the molten pool:
molten pool and determine the locations of the nucleation 𝑇
and crystallization orientation. Cellular automaton has a 𝐻 (𝑇) = ∫ 𝑐𝑑𝑇 + 𝐿 (1 − 𝑓𝑠 ) , (2)
certain physical foundation. The simulated microstructure 0
is not dependent on the meshes of the model and the where 𝐻(𝑇) is the enthalpy, 𝐿 is the latent heat, and 𝑓𝑠 is the
computation speed is faster than other methods and the solid fraction.
computation area is larger.
Shang et al. [10] researched the various defects which 2.2. Nucleation. In terms of the nucleation of the solidi-
may be caused in forming parts on account of technical fication in the molten pool, the paper considers only the
parameters, equipment performance, material characteris- heterogeneous nucleation. For simulating the growth of the
tics, and other factors in the process of metal powder grain, a continuous nucleation model which was proposed
laser rapid forming. The results show that less or over by Rappaz and Gandin [16, 17] was used. The model assumes
accumulation resulting from powder feed delay is the main that the nucleation happens on a range of different nucle-
reason for affecting the dimensional accuracy of cladding ation locations. The continuous and nondiscrete distribution,
and the main reason which affects the appearance of surface 𝑑𝑛/𝑑(Δ𝑇), is used to describe the grain density increased
sticky powder is specific energy. Lin et al. [11] studied the with the increase of the undercooling. The total density of
solidification behavior and morphological evolution of 316L nucleation with a given Δ𝑇 is determined as follows:
stainless steel during laser rapid forming. It is found that
the sample completely consisted of columnar 𝛾 austenitic Δ𝑇
𝑑𝑛
dendrites which grow epitaxially from the base. Zhan et al. 𝑛 (Δ𝑇) = ∫ 𝑑 (Δ𝑇) . (3)
0 𝑑 (Δ𝑇)
[12] established a cellular automaton-finite difference model
and simulated the dendritic grains in the weld pool of Ni-Cr Gaussian distribution of nucleation sites is assumed to
binary alloy. The results reproduce the growth of secondary account for the heterogeneous nucleation in the molten pool:
and tertiary dendrite arms, the competitive growth and the
grain boundary segregation, and so forth. Seo et al. [13] 𝑑𝑛 𝑛max 1 Δ𝑇 − Δ𝑇max 2
= exp [− ( ) ], (4)
predicted the solidification of the grain structures in a Ni-base 𝑑 (Δ𝑇) √2𝜋Δ𝑇𝜎 2 Δ𝑇𝜎
superalloy based on a three-dimensional cellular automaton
model coupled with finite-element heat flow calculation. where 𝑛max is the maximum nucleation density, Δ𝑇𝜎 is the
An empirical relationship between the nucleation density standard deviation of the distribution, and Δ𝑇max is the
at the surface and the initial cooling time of the melt was mean nucleation undercooling. The locations of nucleation
proposed and applied to the input parameter in the model are chosen randomly among all the cellular automaton sites.
for describing the Gaussian distribution of nucleation sites Figure 2 shows the Gaussian distribution for nucleation sites.
at the mold surface. Yin and Felicelli [14] simulated the The parameters of the Gaussian distribution for nucle-
dendritic growth occurring in the molten pool during the ation sites have a deep effect on the growth of the grains. The
laser-engineered net shaping process. Based on the simu- area of the columnar crystal increases with the mean nucle-
lation results and experimental data, empirical expressions ation undercooling. The increase of the standard deviation
describing the relationship between the cooling rate and the undercooling results in the decrease of the grains which have
dendrite arm spacing were proposed. the minimum area but the uniformity of these grains deceases
In this paper, compared with the experimental micro- firstly and then increases. The maximum nucleation density
structures, a model combining the cellular automaton and contributes to the decrease of the grain size.
Mathematical Problems in Engineering 3
T i
dn/dT
Φi
ns [m−2 ]
2ΔTs,𝜎 ΔTs,max
Φj
n [m−3 ] j
Φk
2ΔT,𝜎 ΔT,max
Y
X
Figure 4: Geometry model and mesh of the single point laser cladding.
mesh on the base is 0.5 mm. And then the volume mesh where ℎ is the heat transfer coefficient, 𝑇𝑓 is the surrounding
can be generated automatically. The total number of nodes temperature, and 𝑇𝐵 is the surface temperature on the
is 32653 and the number of elements is 172478. boundary of the solving area.
40 4400
37 4340
34 4280
31 4220
Conductivity (W/(m·k))
Density (kg/m3 )
28 4160
25 4100
22 4040
19 3980
16 3920
13 3860
10 3800
0 200 400 600 800 1000 1200 1400 1600 1800 2000 0 200 400 600 800 1000 1200 1400 1600 1800 2000
Temperature (∘ C) Temperature (∘ C)
(a) (b)
2000 1.0
1700 0.9
1400 0.8
1100 0.7
Enthalpy (KJ/Kg)
800 0.6
Fraction-solid
500 0.5
200 0.4
−100 0.3
−400 0.2
−700 0.1
−1000 0
0 200 400 600 800 1000 1200 1400 1600 1800 2000 1650 1658 1666 1674 1682 1690 1698 1706 1714 1722 1730
Temperature (∘ C) Temperature (∘ C)
(c) (d)
Figure 5: Thermophysical perimeters of TC4 titanium alloy, (a) conductivity, (b) density, (c) enthalpy, and (d) fraction.
2500 2500
2000 2000
Temperature (∘ C)
Temperature (∘ C)
1500
1500
1000
1000
500
500
0
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 0
Time (s) 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0
1 4 The distance from the center of the molten pool
2 5 on Z direction (mm)
3 6 P = 1600 W P = 2000 W
P = 1800 W P = 2200 W
Figure 7: Temperature curves at different parts in the molten pool
during the single point laser cladding. Figure 9: Temperature distribution curves on the 𝑍 direction with
different laser powers.
2200
2100 2500
2000
1900 2000
Temperature (∘ C)
Temperature (∘ C)
1800
1700 1500
1600
1000
1500
1400
500
1300
1200 0
3.0 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 4.0
0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0
Time (s)
The distance from the center of the molten pool
1 4 on Z direction (mm)
t = 3s
2 5
t = 5s
3 6
t = 7s
Figure 8: Temperature curves at different positions in the molten
Figure 10: Temperature distribution curves on the 𝑍 direction with
pool during the solidification.
different laser scanning time.
Figure 11: Simulation results of the grain microstructure at different times, (a) 3.3 s, (b) 3.4 s, (c) 3.5 s, and (d) 3.6 s.
1.0
0.9
0.8
Z Z
0.7
X X
Y Y 0.6
(a) t = 3.0 s (b) t = 3.1 s
0.5
0.4
0.3
0.2
0.1
Z Z
X
Y X
Y
0.0
Figure 12: Fraction solid of the molten pool during the solidification at different times, (a) 3.0 s, (b) 3.1 s, (c) 3.3 s, and (d) 3.6 s.
the high cooling rate in the molten pool. Figure 13(b) is the Table 2: Depths of the molten pools on different laser process
experimental microstructure of the local molten pool. parameters (unit: mm).
In the laser cladding forming experiments, the input laser
Laser power
energy varies in different process parameters which have the Scanning time
great influence on the depths of the molten pools. Based on P = 1600 W P = 1800 W P = 2000 W P = 2200 W
the microstructure experiments; Table 2 shows the depths of t = 3s 0.305 0.358 0.579 0.595
the molten pools on different laser process parameters. The t = 5s 0.516 0.411 0.653 0.668
depths of the molten pools increase with the laser power and t = 7s 0.589 0.600 0.726 0.740
the scanning rate. The range of the depths is from 0.305 mm
to 0.740 mm.
On condition of the constant scanning time, the
6. Conclusions range of the maximum temperature is from 1894∘ C
to 2445∘ C with the increase of the laser power. In
(1) The paper established the CAFE model to simulate addition, the laser power has a larger influence on
the grain structure during the solidification of Ti-6Al- the temperature of the molten pool than the scanning
4V alloy in the molten pool of the single point laser rate.
cladding. The model solves the heat transfer equation (3) The growth of the grain has a close relationship with
in the molten pool by FE method to calculate the the temperature distribution of the molten pool in
temperature distribution. And then the growth of the the single point laser cladding. During the process
grains in the molten pool was simulated based on the of the grain growth, the faster growing grains are
temperature field. perpendicular to the temperature gradients in the
(2) The maximum temperature in the molten pool molten pool. As the solidification continues, the
increases with the laser power and the scanning rate. columnar crystals grow bigger and there is little sign
8 Mathematical Problems in Engineering
200 𝜇m 100 𝜇m
(a) (b)
Figure 13: (a) Experimental microstructure of the global molten pool; (b) experimental microstructure of the local molten pool.
of the isometric crystals. The mean radius of the method,” Journal of Inner Mongolia University of Science and
grains in the simulation is 0.138 mm. Technology, vol. 27, no. 1, 2003.
(4) This study has very important significance for [6] W. Qi, J. Zhang, B. Wang et al., “Three-dimensional simula-
tion of solidified microstructure of Fe-C alloy with cellular
improving the quality of the structure parts manufac-
automaton-finite element method,” The Chinese Journal of
tured from the laser cladding forming and the future Process Engineering, vol. 8, no. 1, pp. 64–67, 2008.
related research work. The CAFE model is a useful [7] S. Bo-wei, W. Lei, L. Xin, and W. Wei-dong, “Progresses in
tool for studying the solidification microstructure numerical simulation of solidification microstructure using
changing in the laser cladding forming. cellular automaton method,” Foundry, vol. 55, no. 5, pp. 439–
443, 2006.
Conflict of Interests [8] B.-C. Liu and T. Xing, Simulation and Quality Control in Casting
Engineering, China Machine Press, 2001.
The authors declare that there is no conflict of interests [9] Y. Yu and G. Yang, “Research development of microstructure
regarding the publication of this paper. numerical simulation in the solidification,” vol. 3, pp. 30–33,
2002.
[10] X. F. Shang, Z. J. Wang, D. X. Han, and Z. T. Xie, “Research
Acknowledgments on defects of titanium alloy laser rapid forming,” Advanced
Materials Research, vol. 189–193, pp. 3726–3730, 2011.
This research is supported by the National Natural Science
[11] X. Lin, H.-O. Yang, J. Chen, and W.-D. Huang, “Microstructure
Foundation of China (no. 51205261). The authors appreciate
evolution of 316L stainless steel during laser rapid forming,”
the Key Laboratory of Fundamental Science for National
Acta Metallurgica Sinica, vol. 42, no. 4, pp. 361–368, 2006.
Defense of Aeronautical Digital Manufacturing Process and
[12] X. Zhan, Z. Dong, Y. Wei, and Y. Wang, “Simulation of columnar
Shenyang Aerospace University. dendrite grain growth in weld pool of Ni-Cr binary alloy,” The
Chinese Journal of Nonferrous Metals, vol. 19, no. 8, pp. 1431–
References 1436, 2009.
[13] S.-M. Seo, I. Kim, C. Jo, and K. Ogi, “Grain structure prediction
[1] H.-M. Wang, S.-Q. Zhang, and H.-B. Shang, “The developments of Ni-base superalloy castings using the cellular automaton-
in laser rapid forming technology research of the large titanium finite element method,” Materials Science and Engineering A,
alloy structures,” Aviation Precision Manufacturing Technology, vol. 448-451, pp. 713–716, 2007.
vol. 44, no. 6, pp. 28–30, 2008. [14] H. Yin and S. D. Felicelli, “Dendrite growth simulation during
[2] J. Mazumder, J. Choi, K. Nagarathnam, J. Koch, and D. Hetzner, solidification in the LENS process,” Acta Materialia, vol. 58, no.
“The direct metal deposition of H13 tool steel for 3-D compo- 4, pp. 1455–1465, 2010.
nents,” JOM, vol. 49, no. 5, pp. 55–60, 1997. [15] L. Peng, S. Hai-Bo, L. Yang, and Z. Jia-Quan, “3D CAFE model
[3] F. Abe, K. Osakada, M. Shiomi, K. Uematsu, and M. Matsumoto, for simulating thesolidification microstructure of 430 stainless
“The manufacturing of hard tools from metallic powders by steel,” Journal of University of Science and Technology Beijing, vol.
selective laser melting,” Journal of Materials Processing Technol- 36, no. 3, pp. 315–322, 2014.
ogy, vol. 111, no. 1–3, pp. 210–213, 2001. [16] M. Rappaz and Ch.-A. Gandin, “Probabilistic modelling of
[4] Y. Huang, D. Zou, G. Liang, and J. Su, “Numerical simulation microstructure formation in solidification processes,” Acta Met-
on cladding track, fluid flow field and temperature field in laser allurgica et Materialia, vol. 41, no. 2, pp. 345–360, 1993.
cladding process with powder feeding,” Rare Metal Materials [17] Ch.-A. Gandin and M. Rappaz, “A coupled finite element-
and Engineering, vol. 32, no. 5, pp. 330–334, 2003. cellular automaton model for the prediction of dendritic grain
[5] W. Qi and J. Zhang, “Simulation of the microstructure evolution structures in solidification processes,” Acta Metallurgica et
of aluminum alloy by using cellular automaton-finite element Materialia, vol. 42, no. 7, pp. 2233–2246, 1994.
Mathematical Problems in Engineering 9
Research Article
Application of CFD, Taguchi Method, and
ANOVA Technique to Optimize Combustion and
Emissions in a Light Duty Diesel Engine
Copyright © 2014 Senlin Xiao et al. This is an open access article distributed under the Creative Commons Attribution License,
which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Some previous research results have shown that EGR (exhaust gas recirculation) rate, pilot fuel quantity, and main injection timing
closely associated with engine emissions and fuel consumption. In order to understand the combined effect of EGR rate, pilot fuel
quantity, and main injection timing on the NO𝑥 (oxides of nitrogen), soot, and ISFC (indicated specific fuel consumption), in this
study, CFD (computational fluid dynamics) simulation together with the Taguchi method and the ANOVA (analysis of variance)
technique was applied as an effective research tool. At first, simulation model on combustion and emissions of a light duty diesel
engine at original baseline condition was developed and the model was validated by test. At last, a confirmation experiment with
the best combination of factors and levels was implemented. The study results indicated that EGR is the most influencing factor
on NO𝑥 . In case of soot emission and ISFC, the greatest influence parameter is main injection timing. For all objectives, pilot fuel
quantity is an insignificant factor. Furthermore, the engine with optimized combination reduces by at least 70% for NO𝑥 , 20% in
soot formation, and 1% for ISFC, in contrast to original baseline engine.
1. Introduction in controlling the emissions and fuel economy [1–3, 10, 11].
However, the decreased experimental numbers still spend