Measurements and simulations of transient characteristics of heat pipes
- 格式:pdf
- 大小:485.10 KB
- 文档页数:7
Sediment Transport Modeling Review—Current andFuture DevelopmentsAthanasios N.͑Thanos ͒PapanicolaouAssociate Professor,Univ.of Iowa,IIHR-Hydroscience and Engineering,Dept.of Civil and Environmental Engineering,Iowa City,IA 52241.E-mail:apapanic@Mohamed ElhakeemResearch Associate Engineer,Univ.of Iowa,IIHR-Hydroscience and En-gineering,Dept.of Civil and Environmental Engineering,Iowa City,IA 52241.E-mail:melhakee@George KrallisERM Inc.,Exton,PA.E-mail:george.krallis@Shwet PrakashERM Inc.,Exton,PAJohn EdingerFaculty Research Associate,Bryn Mawr College,Dept.of Geology.E-mail:jeedgr@IntroductionThe use of computational models for solving sediment transport and fate problems is relatively recent compared with the use of physical models.Several considerations govern the choice be-tween physical and computational models;namely,the nature of the problem that needs to be solved,the available resources,and the overall cost associated with the problem solution.In some specific problems,a combination of physical and computational models can be used to obtain a better understanding of the pro-cesses under investigation ͑de Vries 1973͒.Using computational hydrodynamic/sediment transport models,in general,involves the numerical solution of one or more of the governing differential equations of continuity,momentum,and energy of fluid,along with the differential equation for sediment continuity.An advan-tage of computational models is that they can be adapted to dif-ferent physical domains more easily than physical models,which are typically constructed to represent site-specific conditions.An-other advantage of computational models is that they are not sub-ject to distortion effects of physical models when a solution can be obtained for the same flow conditions ͑identical Reynolds and Froude numbers,same length scale in the three directions,etc.͒as those present in the field.With the rapid developments in numerical methods for fluid mechanics,computational modeling has become an attractive tool for studying flow/sediment transport and associated pollutant fate processes in such different environments as rivers,lakes,and coastal areas.Representative processes in these environments in-clude bed aggradation and degradation,bank failure,local scour around structures,formation of river bends,fining,coarsening and armoring of streambeds,transport of point source and nonpointsource pollutant attached to sediments,such sediment exchange processes as settling,deposition,and self-weight consolidation;coastal sedimentation;and beach processes under tidal currents and wave action.Over the past three decades,a large number of computational hydrodynamic/sediment transport models have been developed ͑Fan 1988;Rodi 2006͒.Extensive reviews of different hydrodynamic/sediment transport models can be found in Nicollet ͑1988͒,Nakato ͑1989͒,Onishi ͑1994͒,Przedwojski et al.͑1995͒,Spasojevic and Holly ͑2000͒,and the ASCE Sedimentation Engi-neering Manual no.110͑2007͒.Broadly speaking,these models can be classified on the basis of the range of their applications ͑e.g.,suspended load versus bed-load;physical versus chemical transport ͒;and their formulation in the spatial and temporal con-tinua ͑e.g.,one-dimensional model ͑1D ͒;two-dimensional model ͑2D ͒;or three-dimensional model ͑3D ͒;and steady versus un-steady ͒.The choice of a certain model for solving a specific prob-lem depends on the nature and complexity of the problem itself,the chosen model capabilities to simulate the problem adequately,data availability for model calibration,data availability for model verification,and overall available time and budget for solving the problem.The objectives of this article are twofold.First,the article aims to trace the developmental stages of current representative ͑1D,2D,and 3D ͒models and describe their main applications,strengths,and limitations.The article is intended as a first guide to readers interested in immersing themselves in modeling and at the same time sets the stage for discussing current limitations and future needs.Second,the article provides insight about future trends and needs with respect to hydrodynamic/sediment transport models.In preparing this article,the authors may have uninten-tionally omitted some models,since including all the available models found in the literature is impossible.Finally,this article is mainly focused on multidimensional computational models ͑2D and 3D models ͒;however,a brief overview of the 1D models is also included for providing a rational comparison of the 1D model features with the main features of the 2D and 3D models.Description of ModelsThis section provides information about the model formulation,the spatial and temporal characteristics,the coupling/linkage of the hydrodynamic and sediment components,and the model’s predictive capabilities.Tables 1–3complement this description by providing useful information about the model capabilities to handle unsteady flows,bed load and suspended load,sediment exchange processes,type of sediment ͑cohesive versus cohesion-less ͒,and multifractional sediment rmation about model acronyms,language,availability,and distribution is also provided in Tables 1–3.Tables 4–6summarize examples of the different model applications.The reader can use these case stud-JOURNAL OF HYDRAULIC ENGINEERING ©ASCE /JANUARY 2008/1D o w n l o a d e d f r o m a s c e l i b r a r y .o r g b y Z h e j i a n g U n i v e r s i t y o n 01/26/13. C o p y r i g h t A S CE .F o r p e r s o n a l u s e o n l y ; a l l r i g h t s r e s e r v e d .ies as a reference guide for model setup,calibration,and verifi-cation.One-Dimensional ModelsSince the early 1980s,1D models have been used with some success in research and engineering practice.Most of the 1D models are formulated in a rectilinear coordinate system and solve the differential conservation equations of mass and momen-tum of flow ͑the St.Venant flow equations ͒along with the sedi-ment mass continuity equation ͑the Exner equation ͒by using finite-difference schemes.Some representative models that are developed on the basis of the previously mentioned equations include MOBED by Krishnappan ͑1981͒,IALLUVIAL by Karim and Kennedy ͑1982͒,CHARIMA by Holly et al.͑1990͒,SEDI-COUP by Holly and Rahuel ͑1990͒,3ST1D by Papanicolaou et al.͑2004͒.The HEC-6formulation by Thomas and Prashum ͑1977͒is also presented in a rectilinear coordinate and is dis-cretized by using finite-difference schemes;but it solves the differential conservation equation of energy instead of the mo-mentum equation.Among other 1D models that use different coordinate systems,equations or schemes of solution are FLUVIAL 11by Chang ͑1984͒,GSTARS by Molinas and Yang ͑1986͒,and OTIS by Runkel and Broshears ͑1991͒.Chang ͑1984͒used a curvilinear coordinate system to solve the governing equations of his model.Molinas and Yang ͑1986͒implemented the theory of minimum stream power to determine the optimum channel width and ge-ometry for a given set of hydraulic and sediment conditions.Runkel and Broshears ͑1991͒modified the 1D advection-diffusion equation with additional terms to account for lateral inflow,first-order decay,sorption of nonconservative solutes,and transient storage of these solutes.Most of the 1D models that are presented here can predict the basic parameters of a particular channel,including the bulk-velocity,water surface elevation,bed-elevation variation,and sediment transport load.All of them,except OTIS,can also pre-dict the total sediment load and grain size distribution of nonuni-form sediment.3ST1D by Papanicolaou et al.͑2004͒cannot differentiate the total sediment load into bed load and suspended load.HEC-6by Thomas and Prashum ͑1977͒and IALLUVIAL by Karim and Kennedy ͑1982͒are not applicable to unsteady flow conditions.Table 1contains the complete model reference and acronym explanation and summarizes the main features for each model.Some of these 1D models have additional specific features.HEC-6by Thomas and Prashum ͑1977͒,for example,decom-poses energy losses into form loss and skin friction loss.MOBED by Krishnappan ͑1981͒can predict the sediment characteristics of a streambed as a function of time and distance for different flow hydrographs.FLUVIAL 11by Chang ͑1984͒accounts for the presence of secondary currents in a curved channel by adjusting the magnitude of the streamwise velocity.The same model can predict changes in the channel bed profile,width,and lateral mi-gration in channel bends.CHARIMA by Holly et al.͑1990͒andHEC-6:Hydraulic Engineering Center;Thomas and Prashum ͑1977͒V .4.2͑2004͒Steady Yes Yes Yes No Entrainment and deposition PDPDF77MOBED:MObile BED;Krishnappan ͑1981͒—Unsteady Yes Yes Yes No Entrainment and deposition C C F90IALLUVIAL:Iowa ALLUVIAL;Karim and Kennedy ͑1982͒—Quasi-steady Yes Yes Yes No Entrainment and deposition C C FIV FLUVIAL 11;Chang ͑1984͒—Unsteady Yes Yes Yes No Entrainment and deposition C P FIV GSTARS:Generalized sediment transport models for alluvial River simulation͑Molinas and Yang,1986͒V .3͑2002͒Unsteady Yes Yes Yes No Entrainment and deposition PDPDF90/95CHARIMA:Acronym of the word CHARiage which means bedload in FrenchHolly et al.͑1990͒—Unsteady Yes Yes Yes Yes Entrainment and deposition C C F 77SEDICOUP:SEDIment COUPled;Holly and Rahuel ͑1990͒—Unsteady Yes Yes Yes No Entrainment and deposition C C F77OTIS:One-dimensional transport with inflow and storage;Runkel and Broshears ͑1991͒V .OTIS-P ͑1998͒Unsteady No Yes No No Advection-diffusion PDPDF 77EFDC1D:Environmental fluid dynamics code;Hamrick ͑2001͒—Unsteady Yes Yes Yes Yes Entrainment anddeposition PD PD F773STD1,steep stream sediment Transport 1D model;Papanicolaou et al.͑2004͒—Unsteady a Yes aYes Yes No Entrainmentand deposition C P F90Note:V ϭversion;C ϭcopyrighted;LD ϭlimited distribution;P ϭproprietary;PD ϭpublic domain;and F ϭFORTRAN.aTreated as a total load without separation.2/JOURNAL OF HYDRAULIC ENGINEERING ©ASCE /JANUARY 2008D o w n l o a d e d f r o m a s c e l i b r a r y .o r g b y Z h e j i a n g U n i v e r s i t y o n 01/26/13. C o p y r i g h t A S CE .F o r p e r s o n a l u s e o n l y ; a l l r i g h t s r e s e r v e d .OTIS by Runkel and Broshears ͑1991͒can treat transport and fate of conservative contaminants and heat.EFDC1D by Hamrick ͑2001͒can be applied to stream networks.3ST1D by Papanico-laou et al.͑2004͒is capable of capturing hydraulic jumps and simulating supercritical flows;therefore,it is applicable to un-steady flow conditions that occur over transcritical flow stream reaches,such as flows over step-pool sequences in mountain streams.Because of their low data,central processor unit ͑CPU ͒re-quirements,and simplicity of use,1D models remain useful pre-dictive tools even today,especially in consulting,for rivers and stream ecological applications where 2D or 3D models may not be needed and are computationally expensive.Table 4presents examples of different 1D model applications.Two-Dimensional ModelsSince the early 1990s,there has been a shift in computational research toward 2D models.Most of the 2D models are currently available to the hydraulic engineering community as interface-based software to allow easy data input and visualization of re-sults.This added capability has made these models user-friendly and popular.2D models are depth-averaged models that can pro-vide spatially varied information about water depth and bed elevation within rivers,lakes,and estuaries,as well as the mag-nitude of depth-averaged streamwise and transverse velocity com-ponents.Most 2-D models solve the depth-averaged continuity and Navier-Stokes equations along with the sediment mass bal-ance equation with the methods of finite difference,finite element,or finite volume.Table 2shows the complete model reference and explanation of the model acronym and summarizes the characteristics of selected 2D hydrodynamic/sediment trans-port models.The main specific features of each model are de-scribed below:SERATRA:A finite-element sediment-contaminant transport model developed by Onishi and Wise ͑1982͒.The model includes general advection-diffusion equations and incorporates sink/source terms.The model can predict overland ͑terrestrial ͒and in-stream pesticide migration and fate to assess the potential short-and long-term impacts on aquatic biota in receiving streams.SUTRENCH-2D:A finite-volume hydrodynamic and sediment transport model developed by van Rijn and Tan ͑1985͒for simu-SERATRA:SEdiment and RAdionuclide TRAnsport;Onishi and Wise ͑1982͒—Unsteadya Yes a Yes No Yes Advection-diffusion CC/LDFIVSUTRENCH-2D:SUspended sediment transport in TRENCHes;van Rijn and Tan ͑1985͒—Quasisteadya Yes a Yes No No Advection-diffusion C LD F90TABS-2;Thomas and McAnally ͑1985͒—Unsteadya Yes a Yes No Yes Entrainment and deposition C C F77MOBED2:MObile BED;Spasojevic and Holly ͑1990a ͒—Unsteady Yes Yes Yes No Entrainment and depositionC C F77ADCIRC:ADvanced CIRCulation;Luettich et al.͑1992͒—Unsteady a Yes aYes No Yes Advection-diffusionC/LD C/LD F90MIKE 21:Danish acronym of the word microcomputer;Danish Hydraulic Institute ͑1993͒—Unsteady a Yes aYes No Yes Entrainment and deposition CPF90UNIBEST-TC:UNIform BEach Sediment Transport—Transport Cross-shore;Bosboom et al.͑1997͒—Quasi-steadya Yes a Yes No No Entrainment and advection C LD F90USTARS:Unsteady Sediment Transport models for Alluvial Rivers Simulations;Lee et al.͑1997͒—Unsteady Yes Yes Yes No Entrainment and deposition P P F90FAST2D:Flow Analysis Simulation Tool;Minh Duc et al.͑1998͒—Unsteady Yes Yes No No Entrainment and deposition LD P F90FLUVIAL 12;Chang ͑1998͒—Unsteady Yes Yes Yes No Entrainment and deposition C P F77Delft 2D;Walstra et al.͑1998͒—Unsteady Yes Yes No Yes Advection-diffusion C LD F90CCHE2D:The National Center for Computational Hydroscience and Engineering;Jia and Wang ͑1999͒V .2.1͑2001͒Unsteady Yes Yes Yes No Advection-diffusion PD/CLDF77/F90Note:V ϭversion;C ϭcopyrighted;LD ϭlimited distribution;P ϭproprietary;PD ϭpublic domain;F ϭFORTRAN.aTreated as a total load without separation.JOURNAL OF HYDRAULIC ENGINEERING ©ASCE /JANUARY 2008/3D o w n l o a d e d f r o m a s c e l i b r a r y .o r g b y Z h e j i a n g U n i v e r s i t y o n 01/26/13. C o p y r i g h t A S CE .F o r p e r s o n a l u s e o n l y ; a l l r i g h t s r e s e r v e d .lating sediment transport and associated bed level change under conditions of combined quasi-steady currents and wind-induced waves over a sediment bed.The model solves the general advection-diffusion equations by incorporating a lag coefficient to account for the settling of sediments.TABS-2:A group of finite-element based hydrodynamic and sediment transport computer codes developed by the USACE Wa-terways Experimental Station ͑Thomas and McAnally 1985͒that currently operates by using the SMS v.9.0windows interface.These codes are applicable to rivers,reservoirs,and estuaries.The main components of TABS-2are the hydrodynamic component,RMA2;the sediment transport component,SED2D ͑formally STUDH ͒;and the water quality component,RMA4.MOBED2:A finite-difference hydrodynamic and sediment transport model used in a curvilinear coordinate system,devel-oped by Spasojevic and Holly ͑1990a ͒.The model can simulate water flow,sediment transport,and bed evolution in natural wa-terways such as reservoirs,estuaries,and coastal environments where depth averaging is appropriate.ADCIRC-2D:A finite-element hydrodynamic and sediment transport model developed by Luettich et al.͑1992͒in a rectilin-ear coordinate system for simulating large-scale domains ͑e.g.,the entire East Coast of the United States ͒by using 2D equationsfor the “external mode”but using the “internal mode”for obtain-ing detailed velocity and stress at localized areas.The internal mode is achieved by specifying the momentum dispersion and the bottom shear stress in terms of the vertical velocity profile.The wave-continuity formulation of the shallow-water equations is used to solve the time-dependent,free-surface circulation and transport processes.MIKE2:A finite-difference model in a rectilinear coordinate system developed by the Danish Hydraulic Institute ͑1993͒for simulating transport and fate of dissolved and suspended loads discharged or accidentally spilled in lakes,estuaries,coastal areas,or in the open sea.The system consists of four main model groups ͑modules ͒,namely,the hydrodynamic and wave models,the sediment process model,and the environmental hydrody-namic model groups.The hydrodynamic and wave models are relevant to the types of physical processes considered in flood-plain mapping.The sediment process models are used to simulate shoreline change and sand transport,whereas the environmental hydrodynamic models are used to examine water quality issues.UNIBEST-TC2:A finite-difference hydrodynamic and sedi-ment transport model in a rectilinear coordinate system to de-scribe the hydrodynamic processes of waves and currents in the cross-shore direction by assuming the presence of uniform meanECOMSED:Estuarine,Coastal,and Ocean Model—SEDiment transport;Blumberg and Mellor ͑1987͒V .1.3͑2002͒Unsteady aYesaYes No YesEntrainment and deposition PDPDF77RMA-10:Resource Management Associates;King ͑1988͒—Unsteady aYesaYes No Yes Entrainment and deposition C P F77GBTOXe:Green Bay TOXic enhancement;Bierman et al.͑1992͒—Unsteady No Yes No Yes Entrainment and deposition NA NA F77EFDC3D:Environmental Fluid Dynamics code;Hamrick ͑1992͒—Unsteady Yes Yes Yes Yes Entrainment and deposition PD P F77ROMS:Regional Ocean Modeling System;Song and Haidvogel ͑1994͒V .1.7.2͑2002͒Unsteady Yes Yes Yes No Entrainment and deposition LD LD F77CH3D-SED:Computational Hydraulics 3D-SEDiment;Spasojevic and Holly ͑1994͒—Unsteady Yes Yes Yes Yes Entrainment and deposition C C F90SSIIM:Sediment Simulation In Intakes with Multiblock options;Olsen ͑1994͒V .2.0͑2006͒Steady Yes Yes Yes No Advection-diffusion PD P C-Langua.MIKE 3:Danish acronym of the word Microcomputer;Jacobsen and Rasmussen ͑1997͒—Unsteady aYesaYes No Yes Entrainment and deposition C P F90FAST3D:Flow Analysis Simulation Tool;Landsberg et al.͑1998͒V .Beta-1.1͑1998͒Unsteady Yes Yes No No Entrainment and depositionLD P F90Delft 3D;Delft Hydraulics ͑1999͒V .3.25.00͑2005͒Unsteady Yes Yes No Yes Entrainment anddepositionC LD F77TELEMAC;Hervouet and Bates ͑2000͒—Unsteady a Yes a Yes No Yes Entrainment anddepositionC P F90Zeng et al.͑2005͒—UnsteadyYes Yes No No Entrainment anddepositionPPF90Note:V ϭversion;C ϭcopyrighted;LD ϭlimited distribution;P ϭproprietary;PD ϭpublic domain;F ϭFORTRAN.aTreated as a total load without separation.4/JOURNAL OF HYDRAULIC ENGINEERING ©ASCE /JANUARY 2008D o w n l o a d e d f r o m a s c e l i b r a r y .o r g b y Z h e j i a n g U n i v e r s i t y o n 01/26/13. C o p y r i g h t A S CE .F o r p e r s o n a l u s e o n l y ; a l l r i g h t s r e s e r v e d .longshore currents along the beach ͑Bosboom et al.1997͒.The bed load and suspended load transport processes are modeled by assuming local equilibrium conditions ͑no lag effects are consid-ered between flow and sediment ͒.USTARS:A modified form of ͑GSTARS ͒that is also based on the stream tube concept ͑Lee et al.1997͒.The hydrodynamic and sediment equations are solved with a finite-difference scheme in a rectilinear coordinate system.As in GSTARS,the theory of mini-mum stream power is used here to determine the optimum chan-nel width and geometry for a given set of hydraulic,geomorpho-logic,sediment,and man-made constraints.FAST2D:A finite-volume hydrodynamic and sediment model with boundary-fitted grids in a curvilinear coordinate system to simulate sediment transport and morphodynamic problems in al-luvial channels ͑Minh Duc et al.1998͒.The model accounts in-directly for secondary effects attributed to the complexity of the domain.FLUVIAL 12:A finite-difference hydrodynamic and sediment model in a curvilinear coordinate system developed by Chang ͑1998͒.The combined effects of flow hydraulics,sediment trans-port,and river channel changes can be simulated for a given flow period.The model is a mobile-bed model and simulates changes in the channel-bed profile,width,and sediment bed composition induced by the channel curvature.DELFT-2D:A finite-difference hydrodynamic and sediment transport model simulating waves and currents ͑Walstra et al.1998͒.The model couples the hydrodynamics with computed bot-tom morphological changes in a time-dependent way.The model can simulate bed-load and suspended load transport by using ei-ther a local equilibrium or a nonequilibrium ͑i.e.,the lag effects between flow and sediment ͒approach.The model can also show the effects of wave motion on transport magnitude and HE2D:A finite-element hydrodynamic and sediment model developed by Jia and Wang ͑1999͒.The model simulates the sus-pended sediment by solving the advection-diffusion equation and the bed-load transport by empirical functions ͑e.g.,Yalin 1972;van Rijn 1993͒.The model accounts for the secondary flow effect in curved channels.All the aforementioned models are applicable to unsteady flow conditions,except SUTRENCH-2D by van Rijn and Tan ͑1985͒and UNIBEST-TC2by Bosboom et al.͑1997͒.All models can predict the total sediment transport load;but only MOBED2,USTARS,FLUVIAL 12,and CCHE2D can handle multifrac-tional sediment transport and can decompose the total sediment load into bedload and suspended load.DELFT-2D and FAST2D can also separate the total sediment load into bedload and sus-HEC-6;Thomas and Prashum ͑1977͒Prediction of the flow and sediment transport along with the bed level change of the Saskatchewan River below Gardiner Dam,Canada ͑Krishnappan 1985͒Prediction of the bed profile for the eroded and redeposited delta sediment upstream from Glines Canyon Dam,Washington ͑U.S.Department of Interior,Bureau of Reclamation 1996͒MOBED;Krishnappan ͑1981͒Comparison of MOBED results with HEC-6results for the flow and sediment transport along with the bed-level change for the Saskatchewan River below Gardiner Dam,Canada ͑Krishnappan 1985͒Prediction of fine sediment transport under ice cover in the Hay River in Northwest Territories,Canada IALLUVIAL;Karim and Kennedy ͑1982͒Simulation of flow and sediment processes in the Missouri River,Nebraska ͑Karim and Kennedy 1982͒Simulation of flow and sediment processes downstream of the Gavins Point Dam on the Missouri River,Nebraska ͑Karim 1985͒FLUVIAL 11;Chang ͑1984͒Simulation of flow and sediment processes of the San Dieguito River,Southern California ͑Chang 1984͒Simulation of flow and sediment processes of the San Lorenzo River,Northern California ͑Chang 1985͒GSTARS;Molinas and Yang ͑1986͒Prediction of the scour depth and pattern at the Lock and Dam No.26replacement site on the Mississippi River,Illinois ͑Yang et al.1989͒Prediction of the variation in channel geometry for the unlined spillway downstream Lake Mescalero Reservoir,New Mexico ͑Yang and Simões 2000͒CHARIMA;Holly et al.͑1990͒Mobile-bed dynamics in the Missouri River from Ft.Randall to Gavins Point Dam,South Dakota ͑Corps of Engineers ͒Mobile-bed dynamics in the Missouri River from Gavins Point Dam to Rulo,Nebraska ͑National Science Foundation ͒SEDICOUP;Holly and Rahuel ͑1990͒Modeling of long-term effects of rehabilitation measures on bed-load transport at the Lower Salzach River,Germany ͑Otto 1999͒Long-term modeling of the morphology of the Danube River,Germany ͑Belleudy 1992͒OTIS;Runkel and Broshears ͑1991͒Simulation of field experiments conducted by Bencala and Walters ͑1983͒for the change in chloride concentration of the Uvas Creek,California.Estimation of the travel times and mixing characteristics of the Clackamas River,Oregon,using the slug of rhodamine data of Laenen and Risley ͑1997͒EFDC1D;Hamrick ͑2001͒Simulation of the flow and sediment transport processes in the Duwamish River and Elliott Bay,Washington ͑Schock et al.1998͒Development of a water quality model for the Christina River,Delaware ͑USEPA 2000͒3STD1;Papanicolaou et al.͑2004͒Prediction of the grain size distribution and bed morphology of the Cocorotico River,Venezuela ͑Papanicolaou et al.2004͒.Prediction of the grain size distribution and bed-load rate of the Alec River,Alaska ͑Papanicolaou et al.2006͒JOURNAL OF HYDRAULIC ENGINEERING ©ASCE /JANUARY 2008/5D o w n l o a d e d f r o m a s c e l i b r a r y .o r g b y Z h e j i a n g U n i v e r s i t y o n 01/26/13. C o p y r i g h t A S CE .F o r p e r s o n a l u s e o n l y ; a l l r i g h t s r e s e r v e d .pended load,but they are limited to uniform sediment sizes.Ex-amples of the different 2D model applications are summarized in Table 5.Three-Dimensional ModelsIn many hydraulic engineering applications,one has to resort to 3D models when 2D models are not suitable for describing cer-tain hydrodynamic/sediment transport processes.Flows in the vi-cinity of piers and near hydraulic structures are examples in which 3D flow structures are ubiquitous and in which 2D models do not adequately represent the physics.With the latest develop-ments in computing technology—such as computational speed,parallel computing,and data storage classification—3D hydrodynamic/sediment transport models have become much more attractive to use.Most 3D models solve the continuity and the Navier-Stokes equations,along with the sediment mass bal-ance equation through the methods of finite difference,finite ele-ment,or finite-volume.The Reynolds average Navier-Stokes ͑RANS ͒approach has been employed to solve the governing equations.The RANS models can be separated into hydrostatic and non-hydrostatic models.The hydrostatic models ͑e.g.,Gessler et al.1999͒are not able to accurately predict flow and transport phe-nomena in regions where the flow is strongly 3D and where large adverse pressure gradients or massive separation are present ͑e.g.,river bends containing hydraulic structures ͒.On the contrary,non-hydrostatic RANS models have been shown to adequately de-scribe intricate features of secondary flows in complex domains ͑e.g.,Wu et al.2000;Ruther and Olsen 2005͒.Table 3shows the complete model reference and explanation of the acronym and summarizes the characteristics of 12selected 3D hydrodynamic/sediment transport models.The main specific features of each model are described below:ECOMSED:A fully integrated 3D finite-difference hydrody-namic,wave,and sediment transport model in an orthogonal cur-SERATRA;Onishi and Wise ͑1982͒Investigation of the effects of sediment on the transport of radionuclides in Cattaraugus and Buttermilk Creeks,New York ͑Walters et al.1982͒Simulation of the hydrogeochemical behavior of radionuclides released to the Pripyat and Dnieper rivers from the Chernobyl Nuclear Power Plant in Ukraine ͑V oitsekhovitch et al.1994͒SUTRENCH-2D;van Rijn and Tan ͑1985͒Simulation of sand transport processes and associated bed-level changes along dredged pits and trenches at the lower Dutch coast,The Netherlands ͑Walstra et al.1998͒Modeling sediment transport and coastline development along the Iranian coast,Caspian Sea ͑Niyyati and Maraghei 2002͒TABS-2;Thomas and McAnally ͑1985͒Simulation of the flow and sediment transport processes in the Black Lake,Alaska ͑Papanicolaou et al.2006͒Evaluation of the hydraulic performance of different structures found in the Missouri River for creating new shallow water habitat ͑Papanicolaou and Elhakeem 2006͒MOBED2;Spasojevic and Holly ͑1990a ͒Simulation of mobile-bed dynamics in the Coralville Reservoir on the Iowa River,Iowa ͑Spasojevic and Holly 1990b ͒ADCIRC;Luettich et al.͑1992͒Simulation of the flow and sediment transport processes of the natural cap in the Matagorda Bay,Texas ͑Edge 2004͒Simulation of sand transport processes at Scheveningen Trial Trench,The Netherlands ͑Edge 2004͒MIKE 21;Danish Hydraulic Institute ͑1993͒Prediction of the spreading of dredged spoils in the Øresund Link,Denmark-Sweden Prediction of sediment transport rate at ebb flow in a tidal inlet,Grådyb,Denmark UNIBEST-TC;Bosboom et al.͑1997͒Coastal study for the impacts of constructing Kelantan Harbor,MalaysiaCoastal study for shoreline protection of Texel region,The NetherlandsUSTARS;Lee et al.͑1997͒Simulation of sand transport processes and associated bed-level changes of a reach in the Keelung River,Taiwan ͑Lee et al.1997͒Routing of flow and sediment of the Shiemen Reservoir,upstream Tan-Hsui River,Taiwan ͑Lee et al.1997͒FAST2D;Minh Duc et al.͑1998͒Simulation of sediment transport processes and associated bed level changes of a reach in the Bavarian Danube River,Germany ͑Minh Duc et al.1998͒Flood analysis and mitigation on the Orlice River,Poland ͑Beck et al.2003͒FLUVIAL 12;Chang ͑1998͒Simulation of flow and sediment processes of the San Dieguito River,Southern California ͑Chang 1994͒Simulation of flow and sediment processes of the Feather River,Northern California ͑Chang et al.1996͒Delft 2D;Walstra et al.͑1998͒Simulation of sand transport processes and associated bed-level changes along dredged pits and trenches at the lower Dutch coast,The Netherlands ͑Walstra et al.1998͒Simulation of the flow field and sediment transport processes of the Pannerdense Kop and IJssel Kop bifurcations in the Rhine River,The Netherlands ͑Sloff 2004͒CCHE2D;Jia and Wang ͑1999͒Investigation of the effects of the rock pile and the submerged dikes downstream of the Lock and Dam No.2of the Red River Waterway,LouisianaInvestigation of the effects of large woody debris structures on the fluvial processes in the Little Topashaw Creek,Mississippi ͑Wu et al.2005͒6/JOURNAL OF HYDRAULIC ENGINEERING ©ASCE /JANUARY 2008D o w n l o a d e d f r o m a s c e l i b r a r y .o r g b y Z h e j i a n g U n i v e r s i t y o n 01/26/13. C o p y r i g h t A S CE .F o r p e r s o n a l u s e o n l y ; a l l r i g h t s r e s e r v e d .。
Modeling of simultaneous heat and mass transfer during convective oven ring cake bakingMelike Sakin-Yilmazer a ,⇑,Figen Kaymak-Ertekin a ,Coskan Ilicali ba Ege University,Faculty of Engineering,Department of Food Engineering,35100Bornova,Izmir,TurkeybKyrgyzstan-Turkey Manas University,Faculty of Engineering,Department of Food Engineering,720044Bishkek,Kyrgyzstana r t i c l e i n f o Article history:Received 11November 2011Received in revised form 2February 2012Accepted 11February 2012Available online 20February 2012Keywords:Ring cakeBaking modeling SimulationFinite difference method Heat and mass transfera b s t r a c tThe convective oven ring cake baking process was investigated experimentally and numerically as a simultaneous heat and mass transfer process.The mathematical model described previously by the authors for cup cake baking was modified to simulate ring cake baking.The heat and mass transfer mech-anisms were defined by Fourier’s and Fick’s second laws,respectively.The implicit alternating direction finite difference technique was used for the numerical solution of the representative model.Prior to the utilization of the developed model in predicting the temperature and moisture profiles for ring cake bak-ing,the results of the numerical model were compared with analytical results involving only heat or mass transfer with constant thermo-physical properties.Excellent agreement was observed.The numerical temperature and moisture contents predicted by the model were compared with the experimental pro-files.They agreed generally reasonably well with the experimental temperature and moisture profiles.Ó2012Elsevier Ltd.All rights reserved.1.IntroductionBaking is a complex process since the physical,chemical and biochemical changes taking place are of complex nature and these changes are simultaneous and coupled.Quantitative understand-ing of processes like volume expansion,evaporation of water,for-mation of a porous structure,denaturation of protein,gelatinization of starch,crust formation and browning reaction are still limited (Sablani et al.,1998;Zhang and Datta,2006).Sim-plified models for the baking process can be developed by treating baking as simultaneous heat and moisture transfer.Even for such simplified models,the estimation of the thermo-physical and transport properties,and the transfer coefficients is a challenging task.Nevertheless,mathematical models supply valuable informa-tion on baking process which may be utilized to optimize the bak-ing conditions to reduce energy consumption.As a traditional process,optimization and process control to ob-tain the desired final product quality is still largely based on expe-rience and good craftsmanship rather than on predictive calculations (Feyissa et al.,2011).Experimental and theoretical studies for baking in general and cake baking in particular were gi-ven in details in Sakin (2005)and Sakin et al.(2007b)and will not be repeated here.Mathematical models representing the baking process reduce the need for the tedious trial and error baking experiments to achieve a high quality product and they serve as a quick and prac-tical tool for pre-design,optimization and process control (Sakin et al.,2007b ).Ring cake is a popular cake geometry all over the world.Selecting this geometry as the test geometry and developing accurate mathematical models will have practical implications to-wards optimizing the baking conditions for better quality and en-ergy saving.Since there is no published work in literature about ring cake baking and its modeling,a two dimensional simplified mathemat-ical model for simultaneous heat and moisture transfer will be developed.Volume rise during baking will be considered.The pre-dictions of the numerical model will be compared with experimen-tal temperature and moisture content profiles.2.Material and experimentalThe cake batter was prepared by thoroughly mixing of 49.4%(of total weight)ready to bake dry cake mix (Dr.Oetker,containing wheat flour,sugar,corn starch,and baking powder),24.7%pasteur-ized whole liquid egg,16.2%vegetable margarine and 9.7%water.Its preparation was given in details in Sakin et al.(2007b).The ini-tial moisture content of the batter was %0.50kg water/kg dry solid.Teflon coated Aluminum was chosen as the material for ring cake mould due to Aluminum’s high thermal conductivity value (206W/mK,Geankoplis,2003).The ring baking mould was an annular ring having an annular radius of 0.04m,radius of 0.11and a height of 0.05m.Baking experiments were carried out in an electrical baking oven (Teba High-01Inox)which was0260-8774/$-see front matter Ó2012Elsevier Ltd.All rights reserved.doi:10.1016/j.jfoodeng.2012.02.020⇑Corresponding author.Tel.:+902323113039;fax:+902323427592.E-mail addresses:melike.sakin@.tr (M.Sakin-Yilmazer),figen.ertekin@.tr (F.Kaymak-Ertekin),coskan.ilicali@ (C.Ilicali).0.39Â0.44Â0.35m(LÂWÂH)wide.The batter was baked at three different oven temperatures(165,185and205°C)under forced convection conditions.The oven was preheated before the baking experiment to obtain uniform oven baking conditions.The air was circulated by a fan in-stalled on the back side of the oven at a constant speed(0.56m/s, measured by Airflow anemometer,LCA6000,UK)and fresh air en-tered steadily into the oven cavity,through a hole(/:28mm)up in the oven to prevent moisture accumulation.Moisture and temperature profiles were recorded during baking experiments.The procedures for taking temperature and moisture content data were similar to that given in Sakin et al.(2007b). Moisture profiles were measured at the top,bottom,annular and side surface crusts as two parallel measurements.An infrared moisture analyzer(Ohaus,MB200,USA),with a precision of ±0.007g,was used for the moisture content determination of the samples.The results obtained by the instrument were previously validated by the standard air oven method for total solids and moisture in baked products,as1h at130°C(AOAC,1995).Temperature profiles were recorded at the top(r=0.063m), bottom(r=0.08m)and annular(z=0.015m)surfaces by using j-type thermocouples(wire size,/:1mm;accuracy:1–1.5°C)that were calibrated in water and ice baths against glass mercury ther-mometer,and by a multi-channel data-logger(Hanna Inst.,HI 98804,Portugal)having an accuracy of±0.5°C,excluding probe er-ror.All the temperature measurements were conducted as two parallels.Techniques used in thermocouple positioning were sim-ilar to Sakin et al.(2007b).The change in the height of the cake batter was measured by a digital caliper(Mitutoyo,Digimatic,Japan).3.Mathematical model descriptionSimultaneous heat and moisture transfer mechanisms in a ring cake during baking can be described by the two dimensional(in r and z directions)unsteady state heat and moisture transfer equa-tions with appropriate boundary conditions.The Fourier’s equation for unsteady state heat conduction for solids with constant thermo-physical properties in afinite cylinder geometry can be modified by the addition of a phase change term to account for the evaporation of the product moisture(Sakin, 2005;Sakin et al.,2003,2004,2005,2007b;Tong and Lund, 1993;Thorvaldsson and Janestad,1999;Shilton et al.,2002).@T¼a1Á@Tþ@2Tþ@2T!þkPÁ@Xð1Þwhere,T is the temperature(°C),t is the time(s),a is the thermal diffusivity(a=k/q c P,m2/s),r is the distance in the radial direction (m),z is the distance in the upward direction(m),k is latent heat of vaporization of water(J/kg),c p is the specific heat capacity(J/kgK),X is the moisture content(kg water/kg dry solid),k is the thermal con-ductivity(W/mK)and q is the density(kg/m3).Uniform initial temperature distribution throughout the cake batter was assumed:t¼0;T¼T ini at all r and zð2Þwhere,T ini is the initial batter temperature(°C).The top,bottom,side and annular surfaces boundary conditions for heat transfer were written as Eqs.(3)–(6),respectively as shown below.kÁ@T@zr¼r;z¼H¼hÁT aÀT r;HðÞð3ÞÀkÁ@T@zr¼r;z¼0¼hÁT aÀT r;0ðÞð4ÞkÁ@T@rr¼R;z¼z¼hÁT aÀT R;zðÞð5ÞNomenclatureA(i),B(i),C(i),AA(i),BB(i),CC(i)elements of coefficient matrixBi Biot number for heat transfer,Bi r¼h cÂR=k;Bi z¼h cÂðH=2Þ=kBi m Biot number for moisture transfer,Bi mr¼k cÂRk cÂR=D eff;Bi mz k¼k cÂðH=2ÞD effc p specific heat(J/kgK)D eff effective moisture diffusivity(m2/s)D(i),DD(i)elements of the right-hand side vectorH thickness of cake batter(m)h surface heat transfer coefficient(W/m2K)k thermal conductivity(W/mK)k c convective mass transfer coefficient(m/s)M number of intervals in z directionMP1M+1,the top surface layerN number of intervals in r directionNH number of intervals in the cake in r directionNH1NH+1,annular surfaceNP1N+1,the radial axisr distance in radial direction(m)R the radius(m)T temperature(°C)T r oven surface temperature(°C)t time(s)U b lower surface overall heat transfer coefficient(W/m2K) U l lateral surface overall heat transfer coefficient(W/m2K) X moisture content(kg water/kg dry solid)z distance in upward direction(m) Subscriptsa ambientannulus annular surfaceb batterbottom bottom surfacecup baking cupdc donecakeeff effectiveffinali the i th nodeini initialj the j th noder,z the r and z directionssides side surfacetop top surfaceGreek symbolsa thermal diffusivity(m2/s)k latent heat of vaporization of water(J/kg) D r,D z the spatial increments in r and z directions.D x thickness(m)D t the time step(s)q the product density(kg/m3)290M.Sakin-Yilmazer et al./Journal of Food Engineering111(2012)289–298ÀkÁ@Tr¼R a;z¼z¼hÁT aÀT Ra;zðÞwhere,h is a convective and radiative surface cient(W/m2K),T a is the baking oven thickness of product(m),R is the radius of the R a is the annular radius(m).Fick’s2nd law describes the unsteady statein solids.Infinite cylinder geometry,the change tent with time and position during baking (Crank,1975).@X¼D eff 1Á@Xþ@2Xþ@2X!Initially a uniform moisture distribution waspositionst¼0;X¼X ini at all r and zwhere,X ini is the initial moisture content of cakedry solid).The experimental data obtained in thisthe moisture contents in the top,bottom,sidedropped sharply at the start of the baking process.Although the top surface moisture content was the lowest,the moisture con-tents for the bottom,side and annular crusts were of comparable values.As the cake batter expands in the upward direction,sharp decrease in the moisture contents and sharp increase in tempera-tures were observed in the layers of the batter in contact with the cake mould.When crust formation begins on the surfaces in contact with the mould surfaces,it is believed that the crust sepa-rates slightly from the hot Aluminum surfaces enabling the loss of moisture also from these surfaces.Therefore,to account for the moisture loss in surface layers effective mass transfer coefficients were assigned to these surfaces.The top surface moisture transfer boundary condition was defined as below,by Eq.(9).ÀD effÁ@X@zr¼r;z¼H¼k c;topÁX r;HÀX1ðÞð9Þwhere,k c,top is an effective convective mass transfer coefficient(m/ s),X1is the equilibrium moisture content in the oven medium.The bottom,side and annulus boundary conditions were de-fined similarly as follows by Eqs.(10)–(12),respectively.D effÁ@X@zr¼r;z¼0¼k c;bottomÁX r;0ÀX1ðÞð10ÞÀD effÁ@X@rr¼R;z¼z¼k c;sidesÁX R;zÀX1ðÞð11ÞD effÁ@Xr¼R a;z¼z¼k c;annulusÁX Ra;zÀX1ðÞð12Þ4.Numerical solution4.1.Grid systemThe cake batter in a ring shaped baking mould was a ring of ra-dius R,annular radius R a and height H.A schematic diagram of the coordinate system used is given in Fig.1.Any rectangular cross-section from the annular surface represents the whole product.A grid system where terms i and j were used to index the nodes in r and z directions was formed(Sakin(2005)).D r and D z were spa-tial increments,and N and M were numbers of intervals in the r and z directions.M was taken as20and N was taken as44.The number of intervals within the cake in the radial direction NH was taken as 28.The volume increase and the subsequent decrease during bak-ing were taken into account as linear changes in height.In the numerical solution,the height of the baking product was refreshed at each time step by using the regression equations obtained from experimental height versus time data.4.2.Solution methodThe governing partial differential equations,defining heat and moisture transfer have been approximated byfinite differences. The implicit alternating direction method was used for the solu-tion.The tri-diagonal matrix systems were formed and solved by the Gaussian Elimination method for the temperature profile,T i,j and moisture profile,X i,j at all points from i=1to NH+1,j=1to M+1,with a time step(D t)of2s(Carnahan et al.,1969;Chandra and Singh,1995).The elements of the coefficients matrix and the right-hand side vector in tri-diagonal matrix system were defined as A(i,j),B(i,j), C(i,j)and D(i,j)for heat transfer;AA(i,j),BB(i,j),CC(i,j)and DD(i,j) for moisture transfer,respectively by Sakin(2005).These coeffi-cients were modified for ring cake baking.The simulation of the model has been done by a computer pro-gram written on Fortran95.Theflow schema of the program is similar to the one given in(Sakin et al.,2007b).4.3.Verification of numerical method against analytical solutions for simplified geometries and materialsThe mathematical model developed for ring cake baking and its numerical results for temperature and concentration profiles must be verified against well known analytical solutions before ring cake baking simulation.The solution of the given coupled partial differ-ential equations,defining the simultaneous heat and mass transfer case of baking,by analytical methods is essentially impossible (Franks,1972).Therefore,the developed model was compared with possible analytical solutions and independently for heat and mass transfer.For heat transfer verification Eq.(1)without the phase change term was used.To eliminate the effect of mass trans-fer,all the mass transfer coefficients were taken as zero.To trans-form the ring cake to an infinite slab,the height of the cake was taken to be much smaller than the radius.When this is the case,M.the ring cake becomes an infinite slab with thickness being equal to the height of the ring cake.Predictions of the numerical model were compared with analytical solutions (Cengel,2003)for a wide range of Biot numbers.Excellent agreement was observed.To the knowledge of the authors,no analytical solution exists for unstea-dy state heat conduction in an annulus,with convective boundary conditions.Therefore,for radial testing of the model,the height of the ring cake was taken to be much larger than the radius and the number of radial intervals in the ring cake,NH,was taken to be equal to the total number of intervals À1,N À1.Furthermore,the heat transfer coefficient for the annular surface was taken to be zero (adiabatic axis).Under these conditions,the resulting geome-try will closely simulate an infinite cylinder geometry.The predic-tions of the ring cake numerical model were compared with the analytical solutions for infinite cylinder.A very high level of agree-ment between analytical and numerical models was obtained.A fi-nal check for the heat transfer section was performed by comparing the predictions of the numerical model with analytical solutions for a finite cylinder.The method of superposition was used for the analytical solution of the two dimensional partial dif-ferential equations.Input data used in the verifications for a finite cylinder are given in Table parison of the predictions of the numerical model for the center temperature with analytical solu-tions for a finite cylinder at Bi r =Bi z =1and Bi r =Bi z =10are shown in Fig.2.Excellent agreements were observed.For mass transfer verifications of the model,different cases were considered:In all these cases,the ambient temperature was as-sumed to be equal to the initial temperature and the latent heat was taken to be zero to eliminate heat transfer effects.Initially,the radius of the cake was taken to be much larger than the height and only convective mass transfer on the top surface was consid-ered.Under these conditions,the ring cake becomes an infinite slabwith half thickness being equal to the height of cake.For the second case used in mass transfer verifications,the radius of the cake was taken to be much larger than the height and convective mass trans-fer was considered for the top and bottom surfaces.The same mass transfer coefficients were used for the both surfaces.For this case,the ring cake becomes an infinite slab with thickness equal to height of the cake.Numerical predictions for the surface,symmetry axis and average moisture content were compared with the analyt-ical solutions for a wide range of mass transfer Biot numbers (Bi m ).Very good agreements were observed for all cases.For radial testing of the model,the height of the ring cake was taken to be much lar-ger than the radius and the number of radial intervals in the ring cake,NH,was taken to be equal to the total number of intervals À1,N À1,similar to the heat transfer case.Furthermore,the mass transfer coefficient for the annular surface was taken to be zero.Un-der these conditions,the resulting geometry will closely simulate an infinite cylinder geometry.Again,very good results were ob-tained for these cases for center,surface and average moisture con-tents.Final verification of the model for mass transfer was carried out for an impermeable hollow short cylinder where the annular ra-dius was 1/40times the cylindrical radius.This geometry is practi-cally a short cylinder.The input parameters used for mass transfer verifications are also given in Table 1.Constant transport properties and constant product thickness were parison of the predictions of the numerical model with analytical solutions for Bi mr =Bi mz =1and Bi mr =Bi mz =10for center and average moisture contents are shown in Figs.3and 4,respectively.The average moisture content was calculated by the following equation;X ave ¼2R 2ÀR 2aZRR aXrdr ð13ÞAs can be observed from Figs.2–4,the numerical model devel-oped predicts accurately the analytical temperature and moistureprofiles.The slight differences in the average moisture contents in Fig.4were attributed to numerical integration error in Eq.(13).Therefore,it was concluded that the numerical model devel-oped for ring cake baking can be used to simulate experimental temperature and moisture content profiles.4.4.Input parameters to ring cake baking simulation and method of validation of simulation resultsThe ring cake baking process parameters used during the numerical solution were the product dimensions,thermo-physicalTable 1Input parameters used in the verification of the numerical model.Cake dimensions (m)H :0.0462,R =0.0231Initial temperature,(°C)T ini =25.7Ambient temperature,(°C)T a =185Heat transfer coefficients (W/m 2K)h :10.39,103.9Thermal conductivity (W/mK)k :0.24Density (kg/m 3)q =945Specific heat capacity (J/kg K)Cp =2650Initial moisture content (kg water/kg dry solids)X ini =0.56Ambient moisture content (kg water/kg dry solids)X 1:0Mass transfer coefficient (m/s)k c =1.3Â10À5Effective diffusivity (m 2/s)D ef =3Â10À7,3Â10À8292M.Sakin-Yilmazer et al./Journal of Food Engineering 111(2012)289–298and mass transport properties,baking oven temperatures,the sur-face heat and mass transfer coefficients,the initial temperatureand moisture content of the cake batter which are illustrated by Table 2.They were either experimentally determined,estimated or taken from the literature to simulate the experimental ring cake baking process.The mass transfer coefficients for oven tempera-ture of 205°C were taken to be 20%higher than the corresponding values at 165and 185°C since the baking time to approximately the same final moisture content was 50min at 205°C,compared to 60min at 165and 185°C.The product density and thermal con-ductivity was assumed to change linearly with cake height accord-ing to the following equations:q ¼q b þðq dc Àq b ÞðH ÀH ini ÞðH f ÀH ini Þð14Þk ¼k b þðk dc Àk b ÞðH ÀH ini ÞðH f ÀH ini Þð15ÞThe specific heat capacity was also evaluated separately and then the thermal diffusivities were calculated.The cake height was measured experimentally.The change in cake height with time was represented by regression straight lines.The %absolute deviation (p ,%)between the experimental and numerical results,for temperature and moisture profiles,were cal-culated by Eq.(16)to validate the simulation results (Boquet,Chi-rife,Iglesias,1978).Table 2Baking process parameters and input variables.Ring cake dimensions (m)H a :0.025,R a a :0.04,R a :0.11Baking temperatures (o C)T a a :165,185,205Initial moisture content (kg moisture/kg solids)X ini a :0.543–0.567Initial temperature (o C)T ini a :23.7–25.7Batter density (kg/m 3)q b a :947.8Donecake density (kg/m 3)q dc a :320Batter thermal conductivity(W/mK)k b b :0.24k b b :0.24Donecake thermalconductivity(W/mK)k dc b :0.121Specific heat capacity (J/kg K)C p b :2650Heat transfer coefficients (W/m 2K)h a :16.3Mass transfer coefficients (m/s)k c ,top a ,c :2.40Â10À5,k c ,bottom a :1.68Â10À5,k c ,sides a :1.56Â10À5,k c ,annulus a :1.56Â10À5k c ,205a =1.2k c,165Effective moisture diffusivity (m 2/s)D eff a ,d :1.7Â10À8a Experimentally measured or estimated values.b Baik et al.(1999).c Demirkol et al.(2006).dSakin et al.(2007a).M.Sakin-Yilmazer et al./Journal of Food Engineering 111(2012)289–298293p ð%Þ¼100N ÁX Ni ¼1X exp ÀXnum j jX expð16Þ5.Results and DiscussionTemperature and moisture profiles at certain positions of thering cake during baking were determined.The difference between the experimental and numerical moisture profiles were compared by p%values;where,a value below 10%indicates the good fit be-tween both results (Boquet,Chirife,Iglesias,1978).The experimentally observed and numerically predicted tem-perature profiles at the positions of the product top surfaceaccuracy of the predictions.For this reason,the experimental data for the center temperature and model predictions were not re-ported in this research.Althoughmeasurements were performed at three oven temperatures,165,185and 205°C,only the results at 165and 205°C were shown to make the figures less crowded.The top surface numerical temperature profile predictions were in good agreement with the experimental results,especially to-wards the completion of baking period (Fig.5).The p %values were calculated as 6.5%and 3.8%,for 165and 205°C,respectively.The specification of the heat transfer coefficient is crucial in heat trans-fer calculations especially at small Biot numbers.When we have a ring cake,the characteristic dimension to be used in the Biot num-ber is ambigous,since we have an annular geometry and the height of the cake also changes.However,if we take (R –R a )is the charac-294M.Sakin-Yilmazer et al./Journal of Food Engineering 111(2012)289–298h¼_mWkAðT aÀTÞtð17Þwhere_m W is the total moisture evaporated(kg)during a com-plete baking run,A is the total heat transfer area(m2),T is the sur-face temperature(o C)and t is the baking time(s);3600s for165°C and3000s for205°C.The surface temperature was assumed to be the mean of ambient temperature and100°C.A single heat trans-fer coefficient was used in all runs.Radiative contribution to heat transfer was not considered.Sakin et al.(2009)have shown that the oven wall temperature and the oven medium temperature dif-fer appreciably depending on oven medium temperature.For an oven medium temperature of160°C,the oven wall average tem-perature was181.6°C in the forced convection oven.At higher oven temperatures,the wall and medium temperature approaches to each other.In the present work,lack of experimental data on heat transfer coefficients will lead to inaccuracies in temperature predictions.The predicted bottom temperatures(Fig.6)were in good agree-ment with the experimental temperatures especially at an oven medium temperature of205°C.The p%values were calculated as 6.1%and3.9%,for165and205°C,respectively.The same heat transfer coefficient was used for all surfaces.However,the convec-tive and radiative heat transfer coefficients will be different for dif-ferent surfaces.Accurate determination of the heat transfer coefficients will be essential for more accurate predictions.Annular surface temperature predictions by the present model were com-pared with experimental values in Fig.7.Satisfactory agreement was observed for oven medium temperature of165°C,where the p%value was calculated as5.5%.However,experimental tempera-tures were greater than the predicted temperatures for205°C.The p%value was calculated as10.9%.This larger deviation was again attributed to the use of the same heat transfer coefficient for all surfaces.The estimation of thermo-physical properties of the ring cake and the transfer coefficients were very important in the numerical solution.They were mostly taken from the related liter-ature(Table2).The differences between the experimental temper-ature profiles and the numerical predictions can also be attributed to these inaccuracies in the thermo-physical properties.More com-prehensive studies on the thermo-physical property estimation of a baking cake batter and the estimation of heat transfer coefficients at different surfaces on a ring cake will result in more accurate numerical models.Similarly,the effective moisture diffusivity term and the mass transfer coefficients have a strong effect on the moisture profile prediction.The relative importance of the estimation of these two parameters will depend on the mass transfer Biot number.It may be shown from the data in Table2that mass transfer Biot numbers vary between11and66.For mass transfer,the internal resistance controls the diffusion of moisture.Therefore,accurate estimation of the effective diffusivity is essential for accurate mass transfer modeling.The available literature for this property for cake baking was limited mostly to the model equations generated for a thin layer of cake batter during drying/baking(Baik and Marcotte,2002;Sakin et al.,2005,2007a)and cup cake baking (Sakin et al.,2007b).As pointed out by Sakin et al.(2007b)the use of the proposed effective moisture diffusivity expressions by Baik and Marcotte(2002)and Sakin(2005)resulted in higher moisture profiles than experimental results.Although the effective diffusivity does not depend on the system geometry,for the same material,the geometry chosen in the estimation of the diffusivity affects the diffusion behaviour.Thus,the effective diffusivities cal-culated from experimental moisture data obtained from the drying or baking of a model system,will be indirectly affected by the geometry,in this case the ring cake.Rather than modifying these expressions,a constant effective diffusivity value was used for the sake of simplicity.The effective moisture diffusivity was taken as1.7Â10À8m2/s in line with the observations of Sakin(2005).When the moisture contents of different crusts of the ring cake were examined,it was observed that generally the top crust had the lowest moisture content and the bottom crust had the second lowest moisture.Moisture profiles for the side and annular crusts were higher than the top and bottom moisture contents and were very similar.Considering the moisture profiles,the top surface con-vective mass transfer coefficient was taken as 1.3Â10À5m/s (Demirkol et al.,2006).The effective convective mass transfer coef-ficients for the other surfaces were taken as shown in Table2.The experimentally observed and numerically predicted moisture pro-files at the top,bottom,side and annular surface positions for oven temperatures of165and205°C were shown in Figs.8–11.In almost all cases,thefinal numerical moisture contents were in close agreement with the experimental values.For the top crust moisture content(Fig.8)the numerical and experimental moisture profiles agreed quite well at165°C(p%value is9.1%).For oven temperature of205°C,the experimental profiles were lower than the numerical profiles;with a p%value larger than10%,indicating that the temperature dependence of the mass transfer coefficient and the effective diffusivity should be considered for this surface. For the bottom surface(Fig.9),the experimental and numerical moisture profiles agreed very well at205°C with a p%8.2.How-ever,at165°C,experimental profiles were higher than the numer-ical ones(the p%value equals to10.0).This can again beattributed M.Sakin-Yilmazer et al./Journal of Food Engineering111(2012)289–298295。
第52卷第12期表面技术2023年12月SURFACE TECHNOLOGY·169·实验参数对力曲线中Cassie-Wenzel状态转换信息的影响金卫凤a,李鑫a,李健b*,张家宇a,郑浩a(江苏大学 a.机械工程学院 b.材料科学与工程学院,江苏 镇江 212013)摘要:目的探明合适的实验参数范围,以获取稳定的润湿状态转换信息,考察实验系统误差、实验过程参数、表面润湿性能和润湿状态转换条件等对力曲线上润湿状态转换信息的影响规律。
方法根据力曲线法的理论,计算出采用不同实验参数时获取的含有润湿状态转换信息的力曲线,并采用文献的实验结果验证理论计算分析出的实验参数影响规律。
结果采用体积为0.1 mL的液滴进行挤压液滴实验时,1 µm的距离误差和0.8 mN的作用力误差即可保证力曲线上润湿状态转换信息较为明显。
对于一般精度的测试系统,能保证较为明显的润湿状态转换信息的液滴体积为0.050 mL以上,加载步长的优选值区间为10~25 µm,精度较高的系统可应用体积为0.010 mL的液滴。
加载表面的润湿性能对力曲线上的润湿状态转换信息影响较小,待测表面在Cassie状态的润湿性能影响也较小,但待测表面在Wenzel状态的接触角大可能使代表润湿状态转换信息的凸起宽度减小。
借助于润湿状态转换过程中液滴充填进超疏水表面微结构内导致的作用力下降信息,可增强力曲线上的润湿状态转换信息。
结论通过选用较小的距离误差和作用力误差,采用大体积液滴和可变加载步长,可以增强力曲线上的润湿状态转换信息,液滴充填微结构引起的力曲线上大幅度的凹陷也可以作为润湿状态转换的标志。
关键词:超疏水表面;状态转换;Cassie状态;Wenzel状态;力曲线法;实验参数中图分类号:TG176 文献标识码:A 文章编号:1001-3660(2023)12-0169-09DOI:10.16490/ki.issn.1001-3660.2023.12.015Effect of Experimental Parameters on Cassie-Wenzel StateTransition Information in Force Response CurveJIN Wei-feng a, LI Xin a, LI Jian b*, ZHANG Jia-yu a, ZHENG Hao a(a. School of Mechanical Engineering, b. School of Materials Science and Engineering,Jiangsu University, Jiangsu Zhenjiang 212013, China)ABSTRACT: Cassie-Wenzel wetting state transition is a typical phenomenon in the failure of superhydrophobic surfaces and has attracted much more attention recently. The force response curve method, which is carried out by squeezing a liquid droplet on the tested superhydrophobic surfaces, is a typical method to characterize the Cassie-Wenzel wetting state transition. In this收稿日期:2023-09-30;修订日期:2023-11-27Received:2023-09-30;Revised:2023-11-27基金项目:国家自然科学基金(51775248);江苏大学2023年大学生创新训练计划(202310299905X)Fund:The National Natural Science Foundation of China (51775248); Funded by the Innovation Training Program of Jiangsu University in 2023 for university students (202310299905X)引文格式:金卫凤, 李鑫, 李健, 等. 实验参数对力曲线中Cassie-Wenzel状态转换信息的影响[J]. 表面技术, 2023, 52(12): 169-177.JIN Wei-feng, LI Xin, LI Jian, et al. Effect of Experimental Parameters on Cassie-Wenzel State Transition Information in Force Response Curve[J]. Surface Technology, 2023, 52(12): 169-177.*通信作者(Corresponding author)·170·表面技术 2023年12月force response curve method, extracting the wetting state transition information from the force response curve is critical, which requires choosing suitable value of experimental parameters. In order to find out the suitable experimental parameter range for obtaining stable transition information of wetting state, the work aims to investigate the effect of experimental system error, experimental process parameters, surface wettability and wettability transition conditions on the wettability transition information in the force response curve. According to the theory of the force response curve method, a series of force response curves with the Cassie-Wenzel wetting state transition information were calculated for the squeezing droplet processes with different experimental parameters. The calculated results were verified by comparison with the force response curve from the former experiments. When the volume of droplet was 0.1 mL in the squeezing droplet experiment, a distance error of 1 μm and a force error of 0.8 mN could ensure the obvious wetting state transition information in the force response curve. More errors in both the distance and the force would make a considerable fluctuation on the force response curve, which might cover up the wetting state transition information. In the general situation of testing systems with a less precision, a droplet with volume above0.050 mL could be used in the experiments to obtain the wetting state transition information. A volume above 0.010 mL couldbe used in experiments on testing system with a higher precision. The optimum value range of loading step was 10~25 μm, which could guarantee a reasonable wetting state transition information over the fluctuation on the force response curve. The above calculated results were verified by comparison between the former experiments from two different groups. Moreover, it was observed that both the wettability of the loading surface and the wettability of the surface to be measured in Cassie state had few effect on the wetting state transition information on the force curve. However, increasing the contact angle of the surface to be measured in Wenzel state might decrease the width of the bulge that represented the wetting state transition information. The wetting state transition information on the force response curve can be enhanced by using a smaller distance error and a force error, a large droplet size and a variable loading step length. Another wetting state transition information, a much more deflection on the force response curve induced by the droplet filling into the microstructure of the tested superhydrophobic surface during the wetting state transition process, can be used in the squeezing droplet experiment. By optimizing the experimental parameters based on the guidance in this work, the mechanism of the Cassie-Wenzel wetting state transition may be further reliably explored.KEY WORDS: superhydrophobic surface; state transition; Cassie state; Wenzel State; force response curve method; experimental parameters固液-气液复合界面(即接触界面处于Cassie接触状态[1])是超疏水表面的本质特征,这种界面特征使超疏水表面具有减阻、防污等功能。
425 BibliographyH.A KAIKE(1974).Markovian representation of stochastic processes and its application to the analysis of autoregressive moving average processes.Annals Institute Statistical Mathematics,vol.26,pp.363-387. B.D.O.A NDERSON and J.B.M OORE(1979).Optimal rmation and System Sciences Series, Prentice Hall,Englewood Cliffs,NJ.T.W.A NDERSON(1971).The Statistical Analysis of Time Series.Series in Probability and Mathematical Statistics,Wiley,New York.R.A NDRE-O BRECHT(1988).A new statistical approach for the automatic segmentation of continuous speech signals.IEEE Trans.Acoustics,Speech,Signal Processing,vol.ASSP-36,no1,pp.29-40.R.A NDRE-O BRECHT(1990).Reconnaissance automatique de parole`a partir de segments acoustiques et de mod`e les de Markov cach´e s.Proc.Journ´e es Etude de la Parole,Montr´e al,May1990(in French).R.A NDRE-O BRECHT and H.Y.S U(1988).Three acoustic labellings for phoneme based continuous speech recognition.Proc.Speech’88,Edinburgh,UK,pp.943-950.U.A PPEL and A.VON B RANDT(1983).Adaptive sequential segmentation of piecewise stationary time rmation Sciences,vol.29,no1,pp.27-56.L.A.A ROIAN and H.L EVENE(1950).The effectiveness of quality control procedures.Jal American Statis-tical Association,vol.45,pp.520-529.K.J.A STR¨OM and B.W ITTENMARK(1984).Computer Controlled Systems:Theory and rma-tion and System Sciences Series,Prentice Hall,Englewood Cliffs,NJ.M.B AGSHAW and R.A.J OHNSON(1975a).The effect of serial correlation on the performance of CUSUM tests-Part II.Technometrics,vol.17,no1,pp.73-80.M.B AGSHAW and R.A.J OHNSON(1975b).The influence of reference values and estimated variance on the ARL of CUSUM tests.Jal Royal Statistical Society,vol.37(B),no3,pp.413-420.M.B AGSHAW and R.A.J OHNSON(1977).Sequential procedures for detecting parameter changes in a time-series model.Jal American Statistical Association,vol.72,no359,pp.593-597.R.K.B ANSAL and P.P APANTONI-K AZAKOS(1986).An algorithm for detecting a change in a stochastic process.IEEE rmation Theory,vol.IT-32,no2,pp.227-235.G.A.B ARNARD(1959).Control charts and stochastic processes.Jal Royal Statistical Society,vol.B.21, pp.239-271.A.E.B ASHARINOV andB.S.F LEISHMAN(1962).Methods of the statistical sequential analysis and their radiotechnical applications.Sovetskoe Radio,Moscow(in Russian).M.B ASSEVILLE(1978).D´e viations par rapport au maximum:formules d’arrˆe t et martingales associ´e es. Compte-rendus du S´e minaire de Probabilit´e s,Universit´e de Rennes I.M.B ASSEVILLE(1981).Edge detection using sequential methods for change in level-Part II:Sequential detection of change in mean.IEEE Trans.Acoustics,Speech,Signal Processing,vol.ASSP-29,no1,pp.32-50.426B IBLIOGRAPHY M.B ASSEVILLE(1982).A survey of statistical failure detection techniques.In Contribution`a la D´e tectionS´e quentielle de Ruptures de Mod`e les Statistiques,Th`e se d’Etat,Universit´e de Rennes I,France(in English). M.B ASSEVILLE(1986).The two-models approach for the on-line detection of changes in AR processes. In Detection of Abrupt Changes in Signals and Dynamical Systems(M.Basseville,A.Benveniste,eds.). Lecture Notes in Control and Information Sciences,LNCIS77,Springer,New York,pp.169-215.M.B ASSEVILLE(1988).Detecting changes in signals and systems-A survey.Automatica,vol.24,pp.309-326.M.B ASSEVILLE(1989).Distance measures for signal processing and pattern recognition.Signal Process-ing,vol.18,pp.349-369.M.B ASSEVILLE and A.B ENVENISTE(1983a).Design and comparative study of some sequential jump detection algorithms for digital signals.IEEE Trans.Acoustics,Speech,Signal Processing,vol.ASSP-31, no3,pp.521-535.M.B ASSEVILLE and A.B ENVENISTE(1983b).Sequential detection of abrupt changes in spectral charac-teristics of digital signals.IEEE rmation Theory,vol.IT-29,no5,pp.709-724.M.B ASSEVILLE and A.B ENVENISTE,eds.(1986).Detection of Abrupt Changes in Signals and Dynamical Systems.Lecture Notes in Control and Information Sciences,LNCIS77,Springer,New York.M.B ASSEVILLE and I.N IKIFOROV(1991).A unified framework for statistical change detection.Proc.30th IEEE Conference on Decision and Control,Brighton,UK.M.B ASSEVILLE,B.E SPIAU and J.G ASNIER(1981).Edge detection using sequential methods for change in level-Part I:A sequential edge detection algorithm.IEEE Trans.Acoustics,Speech,Signal Processing, vol.ASSP-29,no1,pp.24-31.M.B ASSEVILLE, A.B ENVENISTE and G.M OUSTAKIDES(1986).Detection and diagnosis of abrupt changes in modal characteristics of nonstationary digital signals.IEEE rmation Theory,vol.IT-32,no3,pp.412-417.M.B ASSEVILLE,A.B ENVENISTE,G.M OUSTAKIDES and A.R OUG´E E(1987a).Detection and diagnosis of changes in the eigenstructure of nonstationary multivariable systems.Automatica,vol.23,no3,pp.479-489. M.B ASSEVILLE,A.B ENVENISTE,G.M OUSTAKIDES and A.R OUG´E E(1987b).Optimal sensor location for detecting changes in dynamical behavior.IEEE Trans.Automatic Control,vol.AC-32,no12,pp.1067-1075.M.B ASSEVILLE,A.B ENVENISTE,B.G ACH-D EVAUCHELLE,M.G OURSAT,D.B ONNECASE,P.D OREY, M.P REVOSTO and M.O LAGNON(1993).Damage monitoring in vibration mechanics:issues in diagnos-tics and predictive maintenance.Mechanical Systems and Signal Processing,vol.7,no5,pp.401-423.R.V.B EARD(1971).Failure Accommodation in Linear Systems through Self-reorganization.Ph.D.Thesis, Dept.Aeronautics and Astronautics,MIT,Cambridge,MA.A.B ENVENISTE and J.J.F UCHS(1985).Single sample modal identification of a nonstationary stochastic process.IEEE Trans.Automatic Control,vol.AC-30,no1,pp.66-74.A.B ENVENISTE,M.B ASSEVILLE and G.M OUSTAKIDES(1987).The asymptotic local approach to change detection and model validation.IEEE Trans.Automatic Control,vol.AC-32,no7,pp.583-592.A.B ENVENISTE,M.M ETIVIER and P.P RIOURET(1990).Adaptive Algorithms and Stochastic Approxima-tions.Series on Applications of Mathematics,(A.V.Balakrishnan,I.Karatzas,M.Yor,eds.).Springer,New York.A.B ENVENISTE,M.B ASSEVILLE,L.E L G HAOUI,R.N IKOUKHAH and A.S.W ILLSKY(1992).An optimum robust approach to statistical failure detection and identification.IFAC World Conference,Sydney, July1993.B IBLIOGRAPHY427 R.H.B ERK(1973).Some asymptotic aspects of sequential analysis.Annals Statistics,vol.1,no6,pp.1126-1138.R.H.B ERK(1975).Locally most powerful sequential test.Annals Statistics,vol.3,no2,pp.373-381.P.B ILLINGSLEY(1968).Convergence of Probability Measures.Wiley,New York.A.F.B ISSELL(1969).Cusum techniques for quality control.Applied Statistics,vol.18,pp.1-30.M.E.B IVAIKOV(1991).Control of the sample size for recursive estimation of parameters subject to abrupt changes.Automation and Remote Control,no9,pp.96-103.R.E.B LAHUT(1987).Principles and Practice of Information Theory.Addison-Wesley,Reading,MA.I.F.B LAKE and W.C.L INDSEY(1973).Level-crossing problems for random processes.IEEE r-mation Theory,vol.IT-19,no3,pp.295-315.G.B ODENSTEIN and H.M.P RAETORIUS(1977).Feature extraction from the encephalogram by adaptive segmentation.Proc.IEEE,vol.65,pp.642-652.T.B OHLIN(1977).Analysis of EEG signals with changing spectra using a short word Kalman estimator. Mathematical Biosciences,vol.35,pp.221-259.W.B¨OHM and P.H ACKL(1990).Improved bounds for the average run length of control charts based on finite weighted sums.Annals Statistics,vol.18,no4,pp.1895-1899.T.B OJDECKI and J.H OSZA(1984).On a generalized disorder problem.Stochastic Processes and their Applications,vol.18,pp.349-359.L.I.B ORODKIN and V.V.M OTTL’(1976).Algorithm forfinding the jump times of random process equation parameters.Automation and Remote Control,vol.37,no6,Part1,pp.23-32.A.A.B OROVKOV(1984).Theory of Mathematical Statistics-Estimation and Hypotheses Testing,Naouka, Moscow(in Russian).Translated in French under the title Statistique Math´e matique-Estimation et Tests d’Hypoth`e ses,Mir,Paris,1987.G.E.P.B OX and G.M.J ENKINS(1970).Time Series Analysis,Forecasting and Control.Series in Time Series Analysis,Holden-Day,San Francisco.A.VON B RANDT(1983).Detecting and estimating parameters jumps using ladder algorithms and likelihood ratio test.Proc.ICASSP,Boston,MA,pp.1017-1020.A.VON B RANDT(1984).Modellierung von Signalen mit Sprunghaft Ver¨a nderlichem Leistungsspektrum durch Adaptive Segmentierung.Doctor-Engineer Dissertation,M¨u nchen,RFA(in German).S.B RAUN,ed.(1986).Mechanical Signature Analysis-Theory and Applications.Academic Press,London. L.B REIMAN(1968).Probability.Series in Statistics,Addison-Wesley,Reading,MA.G.S.B RITOV and L.A.M IRONOVSKI(1972).Diagnostics of linear systems of automatic regulation.Tekh. Kibernetics,vol.1,pp.76-83.B.E.B RODSKIY and B.S.D ARKHOVSKIY(1992).Nonparametric Methods in Change-point Problems. Kluwer Academic,Boston.L.D.B ROEMELING(1982).Jal Econometrics,vol.19,Special issue on structural change in Econometrics. L.D.B ROEMELING and H.T SURUMI(1987).Econometrics and Structural Change.Dekker,New York. D.B ROOK and D.A.E VANS(1972).An approach to the probability distribution of Cusum run length. Biometrika,vol.59,pp.539-550.J.B RUNET,D.J AUME,M.L ABARR`E RE,A.R AULT and M.V ERG´E(1990).D´e tection et Diagnostic de Pannes.Trait´e des Nouvelles Technologies,S´e rie Diagnostic et Maintenance,Herm`e s,Paris(in French).428B IBLIOGRAPHY S.P.B RUZZONE and M.K AVEH(1984).Information tradeoffs in using the sample autocorrelation function in ARMA parameter estimation.IEEE Trans.Acoustics,Speech,Signal Processing,vol.ASSP-32,no4, pp.701-715.A.K.C AGLAYAN(1980).Necessary and sufficient conditions for detectability of jumps in linear systems. IEEE Trans.Automatic Control,vol.AC-25,no4,pp.833-834.A.K.C AGLAYAN and R.E.L ANCRAFT(1983).Reinitialization issues in fault tolerant systems.Proc.Amer-ican Control Conf.,pp.952-955.A.K.C AGLAYAN,S.M.A LLEN and K.W EHMULLER(1988).Evaluation of a second generation reconfigu-ration strategy for aircraftflight control systems subjected to actuator failure/surface damage.Proc.National Aerospace and Electronic Conference,Dayton,OH.P.E.C AINES(1988).Linear Stochastic Systems.Series in Probability and Mathematical Statistics,Wiley, New York.M.J.C HEN and J.P.N ORTON(1987).Estimation techniques for tracking rapid parameter changes.Intern. Jal Control,vol.45,no4,pp.1387-1398.W.K.C HIU(1974).The economic design of cusum charts for controlling normal mean.Applied Statistics, vol.23,no3,pp.420-433.E.Y.C HOW(1980).A Failure Detection System Design Methodology.Ph.D.Thesis,M.I.T.,L.I.D.S.,Cam-bridge,MA.E.Y.C HOW and A.S.W ILLSKY(1984).Analytical redundancy and the design of robust failure detection systems.IEEE Trans.Automatic Control,vol.AC-29,no3,pp.689-691.Y.S.C HOW,H.R OBBINS and D.S IEGMUND(1971).Great Expectations:The Theory of Optimal Stop-ping.Houghton-Mifflin,Boston.R.N.C LARK,D.C.F OSTH and V.M.W ALTON(1975).Detection of instrument malfunctions in control systems.IEEE Trans.Aerospace Electronic Systems,vol.AES-11,pp.465-473.A.C OHEN(1987).Biomedical Signal Processing-vol.1:Time and Frequency Domain Analysis;vol.2: Compression and Automatic Recognition.CRC Press,Boca Raton,FL.J.C ORGE and F.P UECH(1986).Analyse du rythme cardiaque foetal par des m´e thodes de d´e tection de ruptures.Proc.7th INRIA Int.Conf.Analysis and optimization of Systems.Antibes,FR(in French).D.R.C OX and D.V.H INKLEY(1986).Theoretical Statistics.Chapman and Hall,New York.D.R.C OX and H.D.M ILLER(1965).The Theory of Stochastic Processes.Wiley,New York.S.V.C ROWDER(1987).A simple method for studying run-length distributions of exponentially weighted moving average charts.Technometrics,vol.29,no4,pp.401-407.H.C S¨ORG¨O and L.H ORV´ATH(1988).Nonparametric methods for change point problems.In Handbook of Statistics(P.R.Krishnaiah,C.R.Rao,eds.),vol.7,Elsevier,New York,pp.403-425.R.B.D AVIES(1973).Asymptotic inference in stationary gaussian time series.Advances Applied Probability, vol.5,no3,pp.469-497.J.C.D ECKERT,M.N.D ESAI,J.J.D EYST and A.S.W ILLSKY(1977).F-8DFBW sensor failure identification using analytical redundancy.IEEE Trans.Automatic Control,vol.AC-22,no5,pp.795-803.M.H.D E G ROOT(1970).Optimal Statistical Decisions.Series in Probability and Statistics,McGraw-Hill, New York.J.D ESHAYES and D.P ICARD(1979).Tests de ruptures dans un mod`e pte-Rendus de l’Acad´e mie des Sciences,vol.288,Ser.A,pp.563-566(in French).B IBLIOGRAPHY429 J.D ESHAYES and D.P ICARD(1983).Ruptures de Mod`e les en Statistique.Th`e ses d’Etat,Universit´e deParis-Sud,Orsay,France(in French).J.D ESHAYES and D.P ICARD(1986).Off-line statistical analysis of change-point models using non para-metric and likelihood methods.In Detection of Abrupt Changes in Signals and Dynamical Systems(M. Basseville,A.Benveniste,eds.).Lecture Notes in Control and Information Sciences,LNCIS77,Springer, New York,pp.103-168.B.D EVAUCHELLE-G ACH(1991).Diagnostic M´e canique des Fatigues sur les Structures Soumises`a des Vibrations en Ambiance de Travail.Th`e se de l’Universit´e Paris IX Dauphine(in French).B.D EVAUCHELLE-G ACH,M.B ASSEVILLE and A.B ENVENISTE(1991).Diagnosing mechanical changes in vibrating systems.Proc.SAFEPROCESS’91,Baden-Baden,FRG,pp.85-89.R.D I F RANCESCO(1990).Real-time speech segmentation using pitch and convexity jump models:applica-tion to variable rate speech coding.IEEE Trans.Acoustics,Speech,Signal Processing,vol.ASSP-38,no5, pp.741-748.X.D ING and P.M.F RANK(1990).Fault detection via factorization approach.Systems and Control Letters, vol.14,pp.431-436.J.L.D OOB(1953).Stochastic Processes.Wiley,New York.V.D RAGALIN(1988).Asymptotic solutions in detecting a change in distribution under an unknown param-eter.Statistical Problems of Control,Issue83,Vilnius,pp.45-52.B.D UBUISSON(1990).Diagnostic et Reconnaissance des Formes.Trait´e des Nouvelles Technologies,S´e rie Diagnostic et Maintenance,Herm`e s,Paris(in French).A.J.D UNCAN(1986).Quality Control and Industrial Statistics,5th edition.Richard D.Irwin,Inc.,Home-wood,IL.J.D URBIN(1971).Boundary-crossing probabilities for the Brownian motion and Poisson processes and techniques for computing the power of the Kolmogorov-Smirnov test.Jal Applied Probability,vol.8,pp.431-453.J.D URBIN(1985).Thefirst passage density of the crossing of a continuous Gaussian process to a general boundary.Jal Applied Probability,vol.22,no1,pp.99-122.A.E MAMI-N AEINI,M.M.A KHTER and S.M.R OCK(1988).Effect of model uncertainty on failure detec-tion:the threshold selector.IEEE Trans.Automatic Control,vol.AC-33,no12,pp.1106-1115.J.D.E SARY,F.P ROSCHAN and D.W.W ALKUP(1967).Association of random variables with applications. Annals Mathematical Statistics,vol.38,pp.1466-1474.W.D.E WAN and K.W.K EMP(1960).Sampling inspection of continuous processes with no autocorrelation between successive results.Biometrika,vol.47,pp.263-280.G.F AVIER and A.S MOLDERS(1984).Adaptive smoother-predictors for tracking maneuvering targets.Proc. 23rd Conf.Decision and Control,Las Vegas,NV,pp.831-836.W.F ELLER(1966).An Introduction to Probability Theory and Its Applications,vol.2.Series in Probability and Mathematical Statistics,Wiley,New York.R.A.F ISHER(1925).Theory of statistical estimation.Proc.Cambridge Philosophical Society,vol.22, pp.700-725.M.F ISHMAN(1988).Optimization of the algorithm for the detection of a disorder,based on the statistic of exponential smoothing.In Statistical Problems of Control,Issue83,Vilnius,pp.146-151.R.F LETCHER(1980).Practical Methods of Optimization,2volumes.Wiley,New York.P.M.F RANK(1990).Fault diagnosis in dynamic systems using analytical and knowledge based redundancy -A survey and new results.Automatica,vol.26,pp.459-474.430B IBLIOGRAPHY P.M.F RANK(1991).Enhancement of robustness in observer-based fault detection.Proc.SAFEPRO-CESS’91,Baden-Baden,FRG,pp.275-287.P.M.F RANK and J.W¨UNNENBERG(1989).Robust fault diagnosis using unknown input observer schemes. In Fault Diagnosis in Dynamic Systems-Theory and Application(R.Patton,P.Frank,R.Clark,eds.). International Series in Systems and Control Engineering,Prentice Hall International,London,UK,pp.47-98.K.F UKUNAGA(1990).Introduction to Statistical Pattern Recognition,2d ed.Academic Press,New York. S.I.G ASS(1958).Linear Programming:Methods and Applications.McGraw Hill,New York.W.G E and C.Z.F ANG(1989).Extended robust observation approach for failure isolation.Int.Jal Control, vol.49,no5,pp.1537-1553.W.G ERSCH(1986).Two applications of parametric time series modeling methods.In Mechanical Signature Analysis-Theory and Applications(S.Braun,ed.),chap.10.Academic Press,London.J.J.G ERTLER(1988).Survey of model-based failure detection and isolation in complex plants.IEEE Control Systems Magazine,vol.8,no6,pp.3-11.J.J.G ERTLER(1991).Analytical redundancy methods in fault detection and isolation.Proc.SAFEPRO-CESS’91,Baden-Baden,FRG,pp.9-22.B.K.G HOSH(1970).Sequential Tests of Statistical Hypotheses.Addison-Wesley,Cambridge,MA.I.N.G IBRA(1975).Recent developments in control charts techniques.Jal Quality Technology,vol.7, pp.183-192.J.P.G ILMORE and R.A.M C K ERN(1972).A redundant strapdown inertial reference unit(SIRU).Jal Space-craft,vol.9,pp.39-47.M.A.G IRSHICK and H.R UBIN(1952).A Bayes approach to a quality control model.Annals Mathematical Statistics,vol.23,pp.114-125.A.L.G OEL and S.M.W U(1971).Determination of the ARL and a contour nomogram for CUSUM charts to control normal mean.Technometrics,vol.13,no2,pp.221-230.P.L.G OLDSMITH and H.W HITFIELD(1961).Average run lengths in cumulative chart quality control schemes.Technometrics,vol.3,pp.11-20.G.C.G OODWIN and K.S.S IN(1984).Adaptive Filtering,Prediction and rmation and System Sciences Series,Prentice Hall,Englewood Cliffs,NJ.R.M.G RAY and L.D.D AVISSON(1986).Random Processes:a Mathematical Approach for Engineers. Information and System Sciences Series,Prentice Hall,Englewood Cliffs,NJ.C.G UEGUEN and L.L.S CHARF(1980).Exact maximum likelihood identification for ARMA models:a signal processing perspective.Proc.1st EUSIPCO,Lausanne.D.E.G USTAFSON, A.S.W ILLSKY,J.Y.W ANG,M.C.L ANCASTER and J.H.T RIEBWASSER(1978). ECG/VCG rhythm diagnosis using statistical signal analysis.Part I:Identification of persistent rhythms. Part II:Identification of transient rhythms.IEEE Trans.Biomedical Engineering,vol.BME-25,pp.344-353 and353-361.F.G USTAFSSON(1991).Optimal segmentation of linear regression parameters.Proc.IFAC/IFORS Symp. Identification and System Parameter Estimation,Budapest,pp.225-229.T.H¨AGGLUND(1983).New Estimation Techniques for Adaptive Control.Ph.D.Thesis,Lund Institute of Technology,Lund,Sweden.T.H¨AGGLUND(1984).Adaptive control of systems subject to large parameter changes.Proc.IFAC9th World Congress,Budapest.B IBLIOGRAPHY431 P.H ALL and C.C.H EYDE(1980).Martingale Limit Theory and its Application.Probability and Mathemat-ical Statistics,a Series of Monographs and Textbooks,Academic Press,New York.W.J.H ALL,R.A.W IJSMAN and J.K.G HOSH(1965).The relationship between sufficiency and invariance with applications in sequential analysis.Ann.Math.Statist.,vol.36,pp.576-614.E.J.H ANNAN and M.D EISTLER(1988).The Statistical Theory of Linear Systems.Series in Probability and Mathematical Statistics,Wiley,New York.J.D.H EALY(1987).A note on multivariate CuSum procedures.Technometrics,vol.29,pp.402-412.D.M.H IMMELBLAU(1970).Process Analysis by Statistical Methods.Wiley,New York.D.M.H IMMELBLAU(1978).Fault Detection and Diagnosis in Chemical and Petrochemical Processes. Chemical Engineering Monographs,vol.8,Elsevier,Amsterdam.W.G.S.H INES(1976a).A simple monitor of a system with sudden parameter changes.IEEE r-mation Theory,vol.IT-22,no2,pp.210-216.W.G.S.H INES(1976b).Improving a simple monitor of a system with sudden parameter changes.IEEE rmation Theory,vol.IT-22,no4,pp.496-499.D.V.H INKLEY(1969).Inference about the intersection in two-phase regression.Biometrika,vol.56,no3, pp.495-504.D.V.H INKLEY(1970).Inference about the change point in a sequence of random variables.Biometrika, vol.57,no1,pp.1-17.D.V.H INKLEY(1971).Inference about the change point from cumulative sum-tests.Biometrika,vol.58, no3,pp.509-523.D.V.H INKLEY(1971).Inference in two-phase regression.Jal American Statistical Association,vol.66, no336,pp.736-743.J.R.H UDDLE(1983).Inertial navigation system error-model considerations in Kalmanfiltering applica-tions.In Control and Dynamic Systems(C.T.Leondes,ed.),Academic Press,New York,pp.293-339.J.S.H UNTER(1986).The exponentially weighted moving average.Jal Quality Technology,vol.18,pp.203-210.I.A.I BRAGIMOV and R.Z.K HASMINSKII(1981).Statistical Estimation-Asymptotic Theory.Applications of Mathematics Series,vol.16.Springer,New York.R.I SERMANN(1984).Process fault detection based on modeling and estimation methods-A survey.Auto-matica,vol.20,pp.387-404.N.I SHII,A.I WATA and N.S UZUMURA(1979).Segmentation of nonstationary time series.Int.Jal Systems Sciences,vol.10,pp.883-894.J.E.J ACKSON and R.A.B RADLEY(1961).Sequential and tests.Annals Mathematical Statistics, vol.32,pp.1063-1077.B.J AMES,K.L.J AMES and D.S IEGMUND(1988).Conditional boundary crossing probabilities with appli-cations to change-point problems.Annals Probability,vol.16,pp.825-839.M.K.J EERAGE(1990).Reliability analysis of fault-tolerant IMU architectures with redundant inertial sen-sors.IEEE Trans.Aerospace and Electronic Systems,vol.AES-5,no.7,pp.23-27.N.L.J OHNSON(1961).A simple theoretical approach to cumulative sum control charts.Jal American Sta-tistical Association,vol.56,pp.835-840.N.L.J OHNSON and F.C.L EONE(1962).Cumulative sum control charts:mathematical principles applied to their construction and use.Parts I,II,III.Industrial Quality Control,vol.18,pp.15-21;vol.19,pp.29-36; vol.20,pp.22-28.432B IBLIOGRAPHY R.A.J OHNSON and M.B AGSHAW(1974).The effect of serial correlation on the performance of CUSUM tests-Part I.Technometrics,vol.16,no.1,pp.103-112.H.L.J ONES(1973).Failure Detection in Linear Systems.Ph.D.Thesis,Dept.Aeronautics and Astronautics, MIT,Cambridge,MA.R.H.J ONES,D.H.C ROWELL and L.E.K APUNIAI(1970).Change detection model for serially correlated multivariate data.Biometrics,vol.26,no2,pp.269-280.M.J URGUTIS(1984).Comparison of the statistical properties of the estimates of the change times in an autoregressive process.In Statistical Problems of Control,Issue65,Vilnius,pp.234-243(in Russian).T.K AILATH(1980).Linear rmation and System Sciences Series,Prentice Hall,Englewood Cliffs,NJ.L.V.K ANTOROVICH and V.I.K RILOV(1958).Approximate Methods of Higher Analysis.Interscience,New York.S.K ARLIN and H.M.T AYLOR(1975).A First Course in Stochastic Processes,2d ed.Academic Press,New York.S.K ARLIN and H.M.T AYLOR(1981).A Second Course in Stochastic Processes.Academic Press,New York.D.K AZAKOS and P.P APANTONI-K AZAKOS(1980).Spectral distance measures between gaussian pro-cesses.IEEE Trans.Automatic Control,vol.AC-25,no5,pp.950-959.K.W.K EMP(1958).Formula for calculating the operating characteristic and average sample number of some sequential tests.Jal Royal Statistical Society,vol.B-20,no2,pp.379-386.K.W.K EMP(1961).The average run length of the cumulative sum chart when a V-mask is used.Jal Royal Statistical Society,vol.B-23,pp.149-153.K.W.K EMP(1967a).Formal expressions which can be used for the determination of operating character-istics and average sample number of a simple sequential test.Jal Royal Statistical Society,vol.B-29,no2, pp.248-262.K.W.K EMP(1967b).A simple procedure for determining upper and lower limits for the average sample run length of a cumulative sum scheme.Jal Royal Statistical Society,vol.B-29,no2,pp.263-265.D.P.K ENNEDY(1976).Some martingales related to cumulative sum tests and single server queues.Stochas-tic Processes and Appl.,vol.4,pp.261-269.T.H.K ERR(1980).Statistical analysis of two-ellipsoid overlap test for real time failure detection.IEEE Trans.Automatic Control,vol.AC-25,no4,pp.762-772.T.H.K ERR(1982).False alarm and correct detection probabilities over a time interval for restricted classes of failure detection algorithms.IEEE rmation Theory,vol.IT-24,pp.619-631.T.H.K ERR(1987).Decentralizedfiltering and redundancy management for multisensor navigation.IEEE Trans.Aerospace and Electronic systems,vol.AES-23,pp.83-119.Minor corrections on p.412and p.599 (May and July issues,respectively).R.A.K HAN(1978).Wald’s approximations to the average run length in cusum procedures.Jal Statistical Planning and Inference,vol.2,no1,pp.63-77.R.A.K HAN(1979).Somefirst passage problems related to cusum procedures.Stochastic Processes and Applications,vol.9,no2,pp.207-215.R.A.K HAN(1981).A note on Page’s two-sided cumulative sum procedures.Biometrika,vol.68,no3, pp.717-719.B IBLIOGRAPHY433 V.K IREICHIKOV,V.M ANGUSHEV and I.N IKIFOROV(1990).Investigation and application of CUSUM algorithms to monitoring of sensors.In Statistical Problems of Control,Issue89,Vilnius,pp.124-130(in Russian).G.K ITAGAWA and W.G ERSCH(1985).A smoothness prior time-varying AR coefficient modeling of non-stationary covariance time series.IEEE Trans.Automatic Control,vol.AC-30,no1,pp.48-56.N.K LIGIENE(1980).Probabilities of deviations of the change point estimate in statistical models.In Sta-tistical Problems of Control,Issue83,Vilnius,pp.80-86(in Russian).N.K LIGIENE and L.T ELKSNYS(1983).Methods of detecting instants of change of random process prop-erties.Automation and Remote Control,vol.44,no10,Part II,pp.1241-1283.J.K ORN,S.W.G ULLY and A.S.W ILLSKY(1982).Application of the generalized likelihood ratio algorithm to maneuver detection and estimation.Proc.American Control Conf.,Arlington,V A,pp.792-798.P.R.K RISHNAIAH and B.Q.M IAO(1988).Review about estimation of change points.In Handbook of Statistics(P.R.Krishnaiah,C.R.Rao,eds.),vol.7,Elsevier,New York,pp.375-402.P.K UDVA,N.V ISWANADHAM and A.R AMAKRISHNAN(1980).Observers for linear systems with unknown inputs.IEEE Trans.Automatic Control,vol.AC-25,no1,pp.113-115.S.K ULLBACK(1959).Information Theory and Statistics.Wiley,New York(also Dover,New York,1968). K.K UMAMARU,S.S AGARA and T.S¨ODERSTR¨OM(1989).Some statistical methods for fault diagnosis for dynamical systems.In Fault Diagnosis in Dynamic Systems-Theory and Application(R.Patton,P.Frank,R. Clark,eds.).International Series in Systems and Control Engineering,Prentice Hall International,London, UK,pp.439-476.A.K USHNIR,I.N IKIFOROV and I.S AVIN(1983).Statistical adaptive algorithms for automatic detection of seismic signals-Part I:One-dimensional case.In Earthquake Prediction and the Study of the Earth Structure,Naouka,Moscow(Computational Seismology,vol.15),pp.154-159(in Russian).L.L ADELLI(1990).Diffusion approximation for a pseudo-likelihood test process with application to de-tection of change in stochastic system.Stochastics and Stochastics Reports,vol.32,pp.1-25.T.L.L A¨I(1974).Control charts based on weighted sums.Annals Statistics,vol.2,no1,pp.134-147.T.L.L A¨I(1981).Asymptotic optimality of invariant sequential probability ratio tests.Annals Statistics, vol.9,no2,pp.318-333.D.G.L AINIOTIS(1971).Joint detection,estimation,and system identifirmation and Control, vol.19,pp.75-92.M.R.L EADBETTER,G.L INDGREN and H.R OOTZEN(1983).Extremes and Related Properties of Random Sequences and Processes.Series in Statistics,Springer,New York.L.L E C AM(1960).Locally asymptotically normal families of distributions.Univ.California Publications in Statistics,vol.3,pp.37-98.L.L E C AM(1986).Asymptotic Methods in Statistical Decision Theory.Series in Statistics,Springer,New York.E.L.L EHMANN(1986).Testing Statistical Hypotheses,2d ed.Wiley,New York.J.P.L EHOCZKY(1977).Formulas for stopped diffusion processes with stopping times based on the maxi-mum.Annals Probability,vol.5,no4,pp.601-607.H.R.L ERCHE(1980).Boundary Crossing of Brownian Motion.Lecture Notes in Statistics,vol.40,Springer, New York.L.L JUNG(1987).System Identification-Theory for the rmation and System Sciences Series, Prentice Hall,Englewood Cliffs,NJ.。
The Q400EM is a high-performance, research-grade thermomechanical analyzer (TMA), with unmatched flexibility in operating modes, test probes, fixtures, and available signals. For standard TMA applications, the Q400 delivers the same performance and reliability. It is ideal f or research, teaching, and quality control applications, with perf ormance equivalent to competitive research models.Temperature Range (max)-150 to 1,000°C-150 to 1,000°CTemperature Precision + /- 1°C+ /- 1°CFurnace Cool Down Time <10 min from 600°C to 50°C <10 min from 600°C to 50°C (air cooling)Maximum Sample Size - solid 26 mm (L) x 10 mm (D)26 mm (L) x 10 mm (D)Maximum Sample Size - film/fiber 26 mm (L) x 0.5 mm (T)26 mm (L) x 0.5 mm (T)x 4.7 mm (W)x 4.7 mm (W)Measurement Precision +/- 0.1 %+/- 0.1 %Sensitivity15 nm15 nmDynamic Baseline Drift <1 µm (-100 to 500°C)<1 µm (-100 to 500°C)Force Range 0.001 to 1 N 0.001 to 1 N Force Resolution 0.001 N 0.001 N Frequency 0.01 to 2 Hz Not Available Mass Flow Control Optional Optional AtmosphereInert, Oxidizing, Inert, Oxidizing, (static or controlled flow)or Reactive Gasesor Reactive GasesNote: The Q400 can be field upgraded to the Q400EM.Operational Modes Standard Included Included Stress/Strain Included Not Available CreepIncluded Not Available Stress Relaxation Included Not Available Dynamic TMA (DTMA)IncludedNot AvailableModulated TMA ™(MTMA ™) Included Not Available1The Q400 features a rugged and reliable furnace. Its customized electronics provide excellent heating rate control and rapid response over a wide temperature range. Furnace raising and lowering is soft-ware controlled. The design ensures long life and performance consistency. The excellent heating rate control provides for superior baseline stability and improved sensitivity, while the rapid response permits Modulated TMA™operation. Furnace movementprovides operational convenience, and easy access to the sample chamber.2Located in the furnace core, the easily accessed chamber provides complete temperature and atmosphere control for sample analysis. Purge gas regulation is provided by an optional digital mass flow controller. These include enhanced data quality, ease-of-use, and productivity. The open design simplifies installation of available probes (see Modes of Deformation), sample mounting, and thermocouple placement. Data precision is enhanced by mass flow control of the purge gas.2 13Force Motor A non-contact motor provides aprecisely controlled, friction-free, calibrated force to the sample via the measurement probe or fixture. The force is programmable from 0.001 to 1 N, and can be increased to 2 N by addition of weights to a special tray. A precision sine wave generator provides a set of ten individual frequencies for use in dynamic experiments. Benefits:The motor smoothly generates the accurate and precise static, ramped, or oscillatory dynamic force necessary for quality measurementsin all modes of operation. The choice of frequencies allows optimization of dynamic TMA (DTMA) experiments in compression, 3-point bending, or tension modes of deformation.4Linear Variable Differential Transducer The heart of the Q400TMA sample measurement system is the precision, moveable-core, linear variable differential transducer (LVDT). Benefits:It generates an accurate output signal that is directly proportional to a sample dimension change. Its precise and reliable response over a wide temperature range (–150 to 1,000°C) makes for reproducible TMA results. Its location below the furnace protects it from unwanted temperature effects and ensures stable baseline performance.3421Expansion measurements determine a material’scoefficient of thermal expansion (CTE), glass transi-tion temperature (Tg), and compression modulus.A flat-tipped standard expansion probe (Figure 1)is placed on the sample (a small static force may be applied), and the sample is subjected to a temperature program. Probe movement records sample expansion or contraction. This mode is used with most solid samples. The larger surface area of the macro-expansion probe (Figure 2)better facilitates analysis of soft or irregular samples, powders, and filmsP ENETRATIONPenetration measurements use an extended tip probeto focus the drive force on a small area of the samplesurface (Figure 3). This provides precise measurementof Tg, softening, and melting behavior. It is valuablefor characterizing coatings without their removal froma substrate. The probe operates like the expansionprobe, but under a larger applied force. The hemi-spherical probe (Figure 4)is an alternate penetrationprobe for softening point measurements in solids.C OMPRESSIONIn this mode, the sample is subjected to either a static, linear ramp, or dynamic oscillatory force, while under a defined temperature program, and atmosphere. Sample displacement (strain) is recorded by either expansion / penetration experiments to measure intrinsic material properties, or dynamic tests to determine viscoelastic parameters (DTMA), to detect thermal events, and to separate overlapping transitions (MTMA™).Figure 2Figure 43-P OINT B ENDINGIn this bending deformation (also known as flexure), the sample is supported at both ends on a two-point, quartz anvil atop the stage (Figure 7). A fixed static force is applied vertically to the sample at its center, via a wedge-shaped, quartz probe. Material properties are determined from the force and the measured probe deflection. This mode is considered to represent “pure” deformation, since clamp-ing effects are eliminated. It is primarily used to determine bending properties of stiff materials (e.g., composites), and for distortion temperature measurements.Dynamic (DTMA) measurements are also available with the Q400EM, where aspecial low-friction metallic anvil replaces the quartz version.S PECIALTY P ROBE / F IXTURE K ITSAdditional sample measurement probes and fixtures are available for use with both the Q400 and Q400EM in specialty TMA applications. These include:Dilatometer Probe Kit –for use in volume expansion coefficient measurementsParallel Plate Rheometer –for the measurement of low shear viscosity of materials (10 to 107Pa.s range)under a fixed static force.The expansion, macro-expansion, and penetration probes are supplied with the Q400. These probes, plus the flexure probe, and the low-friction bending fixture, are included with the Q400EM module. Data analysis programs relevant to each of the measurements described are provided in our Thermal Advantage ™for Q Series ™software.T ENSIONTension studies of the stress/strain properties of films and fibers are performed using a Film/Fiber probe assembly (Figure 5). An alignment fixture (Figure 6)permits secure, and reproducible, sample positioning in the clamps. The clamped sample is placed in tension between the fixed and moveable sections of the probe assembly. Application of a fixed force is used to generate stress/strainand modulus information. Additionalmeasurements include Tg, softening temperatures, cure, and cross-link density. Dynamic tests (e.g. DTMA,MTMA™) in tension can be performed to determine viscoelastic parameters (e.g., E |, E ||, tan δ), and to separate overlapping transitions.Figure 6Figure 7TMA measures material deformation under controlled conditions of force, atmosphere, time, and temperature. Force can be applied in compression, flexure, or tension modes using probes previously described. TMA measures intrinsic material properties (e.g., expansion coefficient, glass transition temperature,Young’s modulus), plus processing / product performance parameters (e.g., softening points). These measure-ments have wide applicability, and can be performed by the Q400/Q400EM.TMA can also measure polymer viscoelastic properties using transient (e.g., creep, stress relaxation)or dynamic tests. These require the Q400EM module. In creep, a known stress is applied to the sample, and its deformation is monitored. After a period, the stress is removed, and strain recovery is recorded. In stress relaxation, a fixed strain is applied, and stress decay is monitored.In Dynamic TMA (DTMA), a known sinusoidal stress and linear temperature ramp are applied to the sample, and the resulting sinusoidal strain, and sine wave phase difference (δ), are measured . From this data, storage modulus (E |), loss modulus (E ||), and tan δ(E ||/E |) are calculated as functions of temperature, time, or stress.In Modulated TMA ™(MTMA ™), the sample experiences the combined effects of a linear ramp, and a sinusoidal temperature of fixed amplitude and period . The net signals, after Fourier transformation of the raw data, are total displacement and change in thermal expansion coefficient. Both can be resolved into their reversing and non-reversing component signals.The reversing signals contain events attributable to dimension changes, and are useful in detecting related events (e.g., Tg). The non-reversing signals contain events that relate to time dependent kinetic processes (e.g., stress relaxation).The Q400 and 400EM operating modes permit multiple material property measurements.The Q400 features the Standard mode, while the Q400EM additionally offers Stress/Strain,Creep, Stress Relaxation, Dynamic TMA, and Modulated ™TMA modes.Temperature (Time)Force StrainT(F o r c e )Force (Time)TFS TANDARD M ODE (Q400/Q400EM)Force is constant, and displacement is monitored undera linear temperature ramp. Provides intrinsic property measurements.Strain is constant, and the force required to maintain it ismonitored under a temperature ramp. Permits assessment of shrinkage forces in films/fibers.Force is ramped, and strain measured at constant temperature togenerate force/displacement plots, and modulus information.S TRESS /S TRAIN M ODE (Q400EM)Stress or strain is ramped, and the resulting strain or stress is measured at constant temperature. Both provide stress / strain plots and related modulus information.Strain (Stress)TStrain Stress(S t r a i n )D YNAMIC TMA M ODE (Q400EM):A sinusoidal force (stress) is applied during a temperature ramp. Analysis of the resulting strain and phase data provides viscoelastic property parameters (e.g., E |, E ||tan δ).Timet 2t 1S t r a i n / S t r e s sTemperature (time)STTemperatureTM o d u l a t e d L e n g t hM o d u l a t e d T e m p e r a t u r eC REEP /S TRESS R ELAXATION M ODES (Q400EM)In Creep, stress is held constant, and strain is monitored. In Stress Relaxation,strain is held constant, and stress decay is monitored. Both are transient tests used to assess material deformation and recovery properties.M ODULATED TMA M ODE (Q400EM):Temperature is programmed linearly, and simultaneously modulated at constant stress to generate signals relating to total displacement, CTE, and their reversing and non-reversing components. These permit detection of thermal transitions,and separation of overlapping events (e.g., Tg and stress relaxation).804020090-40-80-12050402030-1001080At a Point 127.3˚C α=25.8µm/m˚CPoint-to-Point Method α=27.6µm/m˚C Average Method α=26.8µm/m˚C230.0˚C45.0˚CAluminumExpansion Probe Size: 7.62mm Prog.: 5˚C/min Atm.: N26040140120100240220200180160260Temperature (˚C)60-20-4071.24˚C -17.48µmSize: 0.492 x 5.41 x 5.08 mm Force: 78.48 mNDeflection: -17.48 µm70605040302080Temperature (˚C)FIGURE 12D i m e n s i o n C h a n g e (µm )I NTRINSIC AND P RODUCT P ROPERTY M EASUREMENTSshows expansion and penetration probe measurements of Tg, and softening point of a synthetic rubber using a temperature ramp at constant force. The large CTE changes in the expan-sion plot indicate the transition temperatures. In penetration, they may be detected by the sharp movement of the loaded probe into the changing material structure.A CCURATE C OEFFICIENT OF T HERMAL E XPANSION (CTE) M EASUREMENTSFigure 11demonstrates the use of the expansion probe to accurately measure small CTE changes in an aluminum sample over a 200˚C temperature range. Advantage ™software permits analysis of the curve slope using an “at point”, “straight line” or “best fit” method to compute the CTE (α) at a selected temperature, or over a range.M ATERIAL P ERFORMANCE ANDS ELECTIONis an example of a 3-point bending mode (flexure probe) experiment on a polyvinyl chloride (PVC) sample, using the ASTM International T est Method E2092 to determine the distortion temper-ature. This test specifies the temperature at which a sample of defined dimensions produces a certain deflection under a given force. It has long been used for predicting material performance.-140-120-100-80-60-40-200102.54˚C-93.22 µm257.71˚C-108.0 µm50100150200250300-50350Temperature (˚C)FIGURE 13D i m e n s i o n C h a n g e (µm )202005202520202015201075250.20.10.30.0-251020304050Time (min)FIGURE 14D i m e n s i o n C h a n g e (µm )T e m p e r a t u r e (˚C )F o r c e (N )20300.0000.0150.0100.005510Slope = Modulus1520Strain (%)FIGURE 15S t r e s s (M P a )0.020M ULTILAYER F ILM A NALYSISFigure 13shows a compression mode analysis, using a penetration probe, of a double layer PE / PET film sample, supported on a metal substrate. The sample temperature was linearly ramped from ambient to 275 ˚C at 5 ˚C/min. The plot shows probe penetra-tions of the PE layer (93.22 µm) at 102 ˚C, and the PET layer (14.78 µm) at 257 ˚C respectively.F ILM P ROPERTY T ESTINGillustrates a classic isostrain experiment, in the tension mode, on a food wrapping film. The film was strained to 20% at room temperature for 5 minutes, cooled to -50 ˚C and held for 5 more minutes, then heated at 5 ˚C/min to 40 ˚C. The plot shows the force variation required to maintain a set strain in the film. The test simulates its use from the freezer to the microwave.F ILM T ENSILE T ESTINGFigure 15displays a strain ramp experiment, at a constant temperature, on a proprietary film in tension. The plot shows an extensive region where stress and strain are linearly related, and over which a tensile modulus can be directly determined. The results show the ability of the Q400EM to function as a mini tensile tester for films and fibers.01230.050.100.150.200.250.300.350.40Force (N)Yield RegionElastic Region40.40.20.60.020406080100120140160180200Temperature (˚C)As ReceivedCold Drawn0.00.20.40.60.81.0-1123456789101112Time (min)FIGURE 18CreepRecoveryS t r a i n (%)1.2F IBER S TRESS /S TRAIN M EASUREMENTSStress/strain measurements are widely used to assess,and compare, materials. shows the different regions of stress/strain behavior in a polyamide fiber (25 µm) in tension, when subjected to a force ramp at a constant temperature. The fiber undergoes instantaneous deformation, retardation,linear stress/strain response, and yield elongation.Other parameters (e.g., yield stress; Young’s modulus) can be determined.T HERMAL S TRESS A NALYSISOFF IBERSdisplays a tension mode experiment,using a temperature ramp at a constant strain (1%), to perform a stress analysis on a polyolefin fiber, as received, and after cold drawing. The plot shows the forces needed to maintain the set strain as a func-tion of temperature. The data has been correlated with key fiber industry, processing parameters, such as shrink force, draw temperature, draw ratio,elongation at break, and knot strength.C REEP A NALYSISCreep tests help in materials selection for end-uses where stress changes are anticipated. Figure 18illustrates an ambient temperature creep study on a polyethylene film in tension. It reveals the instantaneous deformation, retardation, and linear regions of strain response to the set stress, plus its recovery with time on stress removal. The data can also be plotted as compliance, and recoverable compliance, versus time.1301351401450.010.110.00110Time (min)FIGURE 19R e l a x a t i o n M o d u l u s (M P a )150050010001500200025000.100.080.060.040.02200050100150406080100120140160Temperature (˚C)FIGURE 20S t o r a g e M o d u l u s (M P a )T a n D e l t aL o s s M o d u l u s (M P a )3000-20202004000206080100120131.68˚C140160180200Temperature (˚C)FIGURE 21D i m e n s i o n C h a n g e (µm )N o n -R e v D i m e n s i o n C h a n g e (µm )R e v D i m e n s i o n C h a n g e (µm )40S TRESS R ELAXATION A NALYSISshows a stress relaxation test in tension on the same polyolefin film used for the creep study. A known strain is applied to the film, and maintained,while its change in stress is monitored. The plot shows a typical decay in the stress relaxation modulus. Such tests also help engineers design materials for end uses where changes in deformation can be expected.V ISCOELASTIC P ROPERTYD ETERMINATION – D YNAMIC TMAillustrates a dynamic test, in which a semi-crystalline polyethylene terephthate (PET) film in tension is subjected to a fixed sinusoidal stress during a linear temperature ramp. The resulting strain and phase data are used to calculate the material’s viscoelastic properties (E |, E ||, and tan δ). The plotted data shows dramatic modulus changes as the film is heated through its glass transition temperature.S EPARATING O VERLAPPINGT RANSITIONS - M ODULATED ™ TMAFigure 21shows a MTMA™ study to determine the Tg of a printed circuit board (PCB). The signals plotted are the total dimension change, plus its reversing, and non-reversing components. The total signal is identical to that from standard TMA, but does not uniquely define the Tg. The component signals, however, clearly separate the actual Tg from the stress relaxation event induced by non-optimum processing of the PCB.•conduct experiments and simultaneously analyzes data•operates up to 8 modules simultaneously•Wizards – guides and prompts in setting up experiments•provides a real-time display of the progress of the experiment •Autoqueuing– permits pre-programmed set-up of planned experiments •Autoanalysis– permits pre-programmed data analysis of planned experiments•– provides extensive, context sensitive, assistance•– terminates a test upon attaining a specified value (e.g., CTE)UNIVERSAL A NALYSIS ATA A NALYSIS•analyzes data from all TA Instruments modules•provides easy one plot analysis of large and small events•–analyzes data “as it arrives”••within UA 2000 using Microsoft Word™& Excel™templates •– for quick retrieval of previously analyzed data filesI NNOVATIVE E NGINEERINGTA Instruments is the recognized leader for supplying innovative technology,investing twice the industry average in research and development. Our new Q Series™ Thermal Analysis modules are the industry standard. The Q400TMA provides innovative technology suitable for research as well as QC laboratories. The Q400EM includes Dynamic TMA and also Modulated TMA ™, a technique unavailable from other manufacturers.T ECHNICAL S UPPORTCustomers prefer TA Instruments because of our reputation for after-sales support. Our worldwide technical support staff is the largest and most experienced in the industry. They are accessible daily by telephone, email, or via our website. Multiple training opportunities are available including on-site training, seminars in our application labs around the world, and convenient web-based courses.ALESANDERVICEWe pride ourselves in the technical competence and professionalism of our sales force, whose only business is thermal analysis and rheology. TA Instruments is recognized worldwide for its prompt, courteous, and knowledgeable service staff. Their specialized knowledge and experience are major reasons why current customers increasingly endorse our company and products to their worldwide colleagues.Q UALITY P RODUCTSAll thermal analyzers and rheometers are manufactured to ISO 9002 procedures in our New Castle, DE (USA) or our Leatherhead, UK facilities. Innovative flow manufacturing procedures and a motivated, highly skilled, work force ensure high quality products with industry leading delivery times.130******** 33130489460 3227060080 441372360363 31765087270 49602396470 390227421283 81354798418 34936009300 61395530813 46859469200。
CFD COMPUTATIONS OF EMISSIONS FOR LDI-2 COMBUSTORSWITH SIMPLEX AND AIRBLAST INJECTORSKumud AjmaniCFD Nexus, LLCCleveland, Ohio, USAHukam Mongia Purdue University West Lafayette, Indiana, USA Phil LeeWoodward FST, IncZeeland, MI, USAABSTRACT An effort was undertaken to perform CFD analysis of fluid flow in Lean-Direct Injection (LDI) combustors with axial swirl-venturi elements for the next-generation LDI-2 design. The National Combustion Code (NCC) was used to perform non-reacting and reacting flow computations on several LDI-2 injector configurations in a thirteen-element injector array for different operating conditions. All computations were performed with a consistent approach of mesh-optimization, spray-modeling, ignition and kinetics-modeling with the NCC. Computational predictions of emissions (EINOx, EICO and UHC) were compared with the two sets of experimental data representing respectively low-pressure and high-pressure engine cycle conditions. INTRODUCTION Emissions targets set for NOx in NASA‘s Environmentally Responsible Aviation program has revived interest in LDI injection concepts and the attendant low emissions levels achieved with these injectors in previous technology demonstration efforts. The results of the previous generation (LDI-1) efforts were summarized by [Lee 2007] and [Tacina 2008], and some highlights of current generation (LDI-2) efforts have been reported by [Lee 2013] and [Suder 2013]. Experimental measurements of axial-swirler LDI-2 configurations based on a combination of thirteen simplex-type pressure-atomizing injectors and airblast injectors as designed by Woodward FST (WFST) were performed at WFST (for low P3) and NASA Glenn Research Center (for high P3). The LDI-2 design used the lessons learnt from experimental studies of nine-element LDI-1 arrays with axial-swirlers, as performed and reported by [Tacina 2005]. In more recent work, [Tacina 2014a] and [Hicks 2014] have reported optical measurements of unsteady effects and comparison of cold and reacting flows, respectively, of LDI-1 multipoint configurations with bladed axial-swirlers and simplex injectors. The National Combustion Code (NCC) has been developed at NASA Glenn Research Center (GRC) through extensive validation with available LDI data, making it very attractive as a tool to help guide technology development of the 2nd generation Lean Direct Injection, LDI-2. Previous work of [Ajmani2013a] and [Mongia, 2008] have shown the importance of establishing CFD best practices for correlating emissions data from state of the art combustions systems. This paper provides an overview of the efforts undertaken with the NCC to compute heat release, NOx and CO emissions (using RANS reacting flow CFD) for LDI-2 designs proposed by Woodward FST, Inc. The CFD results are compared with experimental measurements of NASA GRC [Tacina 2014], and Woodward FST, Inc, at various engine power conditions.COMPUTATIONAL APPROACH WITH THE NCCThe National Combustion Code (NCC) was used to perform simulations of a Woodward, FST, Inc., LDI-2 configuration. The NCC is a state-of-the-art computational tool that is capable of solving the time-dependent, Navier-Stokes equations with chemical reactions. The NCC is being developed primarily at the NASA Glenn in order to support combustion simulations for a wide range of applications, and has been extensively validated and tested for low-speed chemically reacting flows.D o w n l o a d e d b y C R A N F IE L D U N I V E R S I T Y o n M a y 11, 2016 | h t t p ://a r c .a i a a .o r g | D O I : 10.2514/6.2014-3529 50th AIAA/ASME/SAE/ASEE Joint Propulsion ConferenceJuly 28-30, 2014, Cleveland, OHAIAA 2014-3529Propulsion and Energy ForumThe NCC uses second-order accurate central-differences for the convective and diffusion flux discretization, and a Jameson operator (a blend of 2nd and 4th-order dissipation terms) for numerical stability. The second and fourth order dissipation parameters are typically set to 10-4 and 0.05, respectively ([Swanson 1997]). The value of k 2, the constant that scales the second order dissipation gradient switch, is typically set to 0.25. In order to enhance convergence acceleration in pseudo-time, implicit residual smoothing is used to smooth the computed residuals in NCC RANS. Turbulence closure is obtained by using a two-equation, cubic k-ε model with variable Cµ ([Shih 1998]) and dynamic wall functions with pressure gradient effects ([Shih 2000]).Liquid Phase and Spray ModelingThe specification of the fuel injector exit condition plays a major role in the fidelity of the NCC simulations. Spray injection of fuel particles needs to be specified slightly downstream (1.0mm) of the injector-exit plane of the simplex injectors. In earlier work with NCC RANS, it has been reported that spray initial conditions have a significant impact on reacting-flow results ([Iannetti 2008]). A droplet initial temperature of 300K, and fuel/air ratio dependent initial SMD and injection velocity, as supplied by the injector manufacturer Woodward FST, were specified as inputs for the spray solver. These inputs were based on measurements made by Woodward FST for different fuel-flow rates through the simplex and airblast injectors at simulated atmospheric flow conditions.The best-practice spray inputs were used for the NCC for 60o hollow cone spray with the annular angle width of 100 (viz. droplets spread within 55o to 65o region), 10 droplet groups, discretized into 32 spatial streams along the 3600 circumference. The typical spray integration time-steps were 2e-7s (local time-step, dtml ) and 4e-5s (global time-step, dtgl ), which translates to 200 local time-steps for each global time-step for the spray solver. At each spray time-step, the spatial streams were permitted a stochastic variation of the stream location within the 10o cone thickness. CFD trials with larger number of streams (64, 96) and larger droplet groups (12, 16) with a single element LDI configuration, showed no significant impact on the NCC reacting-flow predictions [Ajmani 2013b].The liquid spray (Jet-A fuel) was modeled by tracking spray particles in a Lagrangian framework, where each particle represents a group of actual spray droplets [Raju 2012]. The governing equations for the liquid phase are based on a Lagrangian formulation where the spray particle position and velocity are described by a set of ordinary differential equations. The Lagrangian solution process used for this study employs a best-practice unsteady spray model such that droplet groups are only integrated for a fraction of their lifetime (but restarted at this point for the next iteration), rather than to a completely steady-state solution. An inflow droplet size distribution is prescribed by the correlation equation: dn n =4.21x 106d d 32!"#$%&3.5e−16.98d d 32()*+,-0.4ddd 32Here n is the total number of droplets, d 32 is the Sauter mean diameter (SMD), and dn is the number of droplets in the size range between d and d + dd. A user-specified number of ‘droplet groups’ is used to represent the drop size distribution among a finite number of droplet classes. As the typical SMD for the conditions computed here was less than 10µm, the specified droplets undergo evaporation without any primary or secondary breakup modeling.Chemical Kinetics and Ignition ModelingA finite-rate chemistry model was used to compute the species source-terms for Jet-A/air chemistry. Reacting flow computations were performed with the chemical-kinetics model described in Table 1. The chemistry model incorporates 14 species and 18 chemical reaction steps. Jet-A fuel is modeled as a surrogate mixture of decane (73%), benzene (18%) and hexane (9%). The kinetics mechanism was validated by matching adiabatic flame temperature, flame-speed and ignition-delay with experimental shock-tube data in the equivalence ratio range of 0.5 to 1.0. The CHEMKIN code was used to find an appropriate activation energy for the fuel breakup step which would produce a close fit between computed D o w n l o a d e d b y C R A N F I E L D U N I V E R S I T Y o n M a y 11, 2016 | h t t p ://a r c .a i a a .o r g | D O I : 10.2514/6.2014-3529ignition delay and experimental ignition delay for a particular equivalence ratio, initial pressure and, initial temperature [Ajmani 2010].The kinetics model uses A (pre-exponential factor), n (temperature exponent) and E (activation energy, cal/mol) to compute the Arrhenius rate coefficient, k = A (T/T 0)n e (-E/RT), for a given temperature, T (K). (R = universal gas constant, T 0 (K) is a reference temperature). Note that reaction steps 1-3 are irreversible, and reaction steps 4-18 are formulated as reversible reactions. The kinetics for NOx prediction includes an extended Zeldovich mechanism (four steps for NO) and an additional four steps for N 2O species. The inclusion of N 2O is expected to improve the NOx predictions in the small local regions where fuel-rich burning is occurring in the flow.Reaction A n E1 C11H21 + O2 => 11CH + 10H + O2 1.00E+12 0.00 3.10E+04GLO / C11H21 0.8/GLO / O2 0.8/2 CH + O2 => CO + OH 2.00E+15 0.00 3.00E+033 CH + O => CO + H 3.00E+12 1.00 0.00E+004 H2 + O2 <=> H2O + O 3.98E+11 1.00 4.80E+045 H2 + O <=> H + OH 3.00E+14 0.00 6.00E+036 H + O2 <=> O + OH 4.00E+14 0.00 1.80E+047 H2O + O2 <=> 2O + H2O 3.17E+12 2.00 1.12E+058 CO + OH <=> CO2 + H 5.51E+07 1.27 -7.58E+029 CO + H2O <=> CO2 + H2 5.50E+04 1.28 -1.00E+0310 CO + H2 + O2 <=> CO2 + H2O 1.60E+14 1.60 1.80E+0411 N + NO <=> N2 + O 3.00E+12 0.30 0.00E+0012 N + O2 <=> NO + O 6.40E+09 1.00 3.17E+0313 N + OH <=> NO + H 6.30E+11 0.50 0.00E+0014 N + N + M <=> N2 + M 2.80E+17 -0.75 0.00E+0015 H + N2O <=> N2 + OH 3.50E+14 0.00 7.55E+0216 N2 + O2 + O <=> N2O + O2 1.00E+15 0.00 3.02E+0217 N2O + O <=> 2NO 1.50E+15 0.00 3.90E+0418 N2O + M <=> N2 + O + M 1.16E+15 0.00 3.32E+0418 N2O + M <=> N2 + O + M 1.16E+15 0.00 3.32E+04Table 1: Kinetics mechanism for Jet-A fuel surrogateA computationally affordable kinetics mechanism (of fewer than 20 species) for liquid spray simulations of Jet-A fuel which can provide NOx and CO predictions, and lean blow-out, at high P 3, high T 3, and low equivalence ratio (< 0.5) conditions for LDI combustor design, remains an open challenge for the chemical kinetics community. Current work at NASA Glenn is focusing on optimizing the current kinetics mechanism for emissions predictions by comparing EINOx and EICO values of the optimized reduced mechanism with emissions values of ‘detailed mechanisms’ and experimental data. More details of the optimization efforts of this kinetics mechanism for emissions predictions for LDI type configurations is provided in [Ajmani 2014].Ignition modeling was performed by introducing artificial ignition source terms in a region 2-3mm downstream of each venturi-exit plane. The ignition sources were turned off in the NCC when every computational cell in the ignition zone reached a ‘cutoff’ temperature of 1600K, or if 1000 iterations with the ignition sources were reached. No “re-ignition” of the mixture was allowed, once the NCC solver had turned off the ignition sources. This practice ensured consistency of CFD computations between the different geometries and different fuel/air ratios computed in this study.D o w n l o a d e d b y C R A N F IE L D U N I V E R S I T Y o n M a y 11, 2016 | h t t p ://a r c .a i a a .o r g | D O I : 10.2514/6.2014-3529NO Computations and EINOx ComparisonsReacting flow computations with the finite-rate chemistry model for Jet-A/Air described in Table 1 were performed until the mass-weighted averaged EINOx value across the exit plane (located 150mm downstream of the combustor dump plane) was within a range of ±5% of previous values, over 1000 steady-state iterations. A computational averaging of the values of NOx across the entire exit plane was deemed an acceptable best-practice for comparison with experimental data of [Tacina 2014b ].Chemistry-turbulence interaction is an important physical phenomenon to couple the effect of turbulent fluctuations with chemical reactions. The developers of the NCC have implemented and validated several different modeling approaches for chemistry-turbulence interaction. Some results from these efforts have been reported in [Liu 2013]. A computational solution using a Time Filtered Navier Stokes (TFNS) approach with an Eulerian-PDF turbulence-chemistry interaction and detailed chemical kinetics for a single-element LDI configuration was also undertaken and reported in [Ajmani 2013c].LDI-2 GEOMETRY, MESHING AND CFD SETUPThe National Combustor Code (NCC) was used to perform simulations of a multiple elements of the thirteen-element WFST LDI-2 configuration with a combination of airblast and simplex injection elements . A candidate LDI-2 arrangement with no recessed elements (also referred to as ‘baseline’) was chosen for these computations. The Pilot and all three Main stages have their venturi-exit plane in-line with combustor dump-plane. This arrangement, identified as Baseline Config 10, had the following setup for the elements:Pilot: Simplex (pressure atomizer) tip, 55oCCW Main1: 4 elements, Simplex tips, 45o CCWMain2, Main3: 8 Elements, Airblast tips, OAS/IAS 45o CW/45 o CWOAS/IAS – Outer/Inner Air SwirlersCW – Clockwise orientation, CCW – Counter Clockwise orientationThe LDI-2 geometry for the Baseline Config 10 as supplied by WFST (see figure 1) was imported into the CUBIT mesh-generation software, to create a fully tetrahedral mesh with 17M elements. Each blade passage and venturi was meshed as an individual block, and these blocks were then “imprinted” or merged with connecting volumes at their respective common surfaces. This ensured consistency of meshing across similar geometric elements, and also allowed for ‘drop-in’ replacement of single (or multiple) swirlersand/or injectors without needing to regenerate the complete mesh for all thirteen elements.Figure 1. Geometry and mesh for LDI-2 Baseline configuration (Config 10) D o w n l o a d e d b y C R A N F I E L D U N I V E R S I T Y o n M a y 11, 2016 | h t t p ://a r c .a i a a .o r g | D O I : 10.2514/6.2014-3529NCC RANS Computations were performed for three different conditions (as tested at NASA GRC) designated Case A, Case B and Case C, as shown in figure 2. These three conditions may be nominally viewed as low, medium, and high power cycle conditions. Some highlights of the computational analysis with the NCC were:• Fully tetrahedral mesh consists of approximately 17M+ elements. Best-practices for meshing of LDI-1 injector elements reported in [Ajmani, 2013 (ASM)] were leveraged for LDI-2 meshing.• Boundary Conditions– Fixed mass flow rate and static-temperature at inflow– Fixed static-pressure (P 3-dP) at outflow– Adiabatic, no-slip conditions at all wallsPower T 3 avg P 3 avg FAR, P ilot FAR, M ain1 FAR, Main2 FAR, M ain3 FAR, T otal (F) (psi) Case A Low 425 83 0.0654 0.0403 0.0000 0.0000 0.0172 Case B Medium 1001 264 0.0260 0.0261 0.0261 0.0263 0.0251 Case C High 1085 234 0.0257 0.0351 0.0349 0.0352 0.0329dP Air f low AC d P 4 EINOx EIHC EICO T 4 (Expt)(psi) (lb/sec) (in 2) (psi) (g /kg F uel) (g /kg F uel) (g /kg F uel) (F)Case A 2.8 1.0604 1.91 79 2.8 11.2 94.3 1637Case B 10.4 3.0352 2.06 253 3.8 0.0 0.0 2512Case C 7.3 2.3122 2.04 226 11.9 0.0 1.0 2989Figure 2. Three power conditions for Woodward FST LDI-2 design, tested at NASA GRCThe first stage of each simulation consists of non-reacting RANS computations till a ‘mass-imbalance (outflow-inflow) convergence of 0.1% over 500 consecutive NCC RANS iterations is achieved. Typically, 100,000 RANS iterations at a CFL of 1.95 are run, to obtain a converged, steady-state, non-reacting flowfield. The CFL is then lowered to 0.95, and the Lagrangian spray is then initiated for 100 iterations, followed by introduction of artificial ignition sources, to ignite the mixture. An 18 step, 14 species finite-rate chemical kinetics (as discussed earlier in Table 1) is used for the reacting flow simulations. Thereacting RANS computations are run until the value of EINOx at the computational exit plane “converges”. The EINOx is considered converged when its variation over 1000 successive iterations is within a 5% band at the exit plane. For comparison with experimental (five-hole probe) data, computed EINOx values are based on the mass-weighted average of EINOx across the entire exit plane.The specification of liquid fuel flow for each element of the Pilot, Main1, Main2 and Main3 for all cases is listed below:– Each Simplex (Pilot, Main1) injector is modeled with SMD=8.8µm, Vinj=38.6m/s, 60o hollow cone, 8 droplet groups, 32 streams– Each Airblast injector is modeled as 8 or 16 ‘discrete’ circumferential injection holes with SMD=7.5µm, Vinj=5m/s, 10o solid cone, 8 droplet groups, 8 streams. The modeling with 16 ‘discrete’ holes provides better results than 8 holes, as the circumferential film is better resolved with greater number of holes.NCC predictions for effective-area (AC d ) at CASE A conditionsNon-reacting flow: ΔP = 17446Pa = 3.29% of P 3; Computed AC d = 1.98in 2Reacting flow: ΔP = 18908Pa = 3.63%; Computed AC d = 1.90in 2 (Experiment = 1.91in 2) NCC predictions for AC d at Case B conditionsNon-reacting flow: ΔP = 63013Pa = 3.46%; Computed AC d = 2.08in 2Reacting flow: ΔP = 76428Pa = 4.2%; Computed AC d = 1.87in 2 (Experiment = 2.06in 2) NCC predictions for AC d at CASE C conditionsNon-reacting flow: ΔP = 56997Pa = 3.53%; Computed AC d = 1.88in 2Reacting flow: ΔP = 55917Pa = 3.47%; Computed AC d = 1.90in 2 (Experiment = 2.04in 2) D o w n l o a d e d b y C R A N F I E L D U N I V E R S I T Y o n M a y 11, 2016 | h t t p ://a r c .a i a a .o r g | D O I : 10.2514/6.2014-3529Reacting Flow: Axial Velocity Predictions with NCC for Case A (Low Power)Figures 3 and 4 show axial -velocity contour comparisons for NCC predictions for a typical CASE A computation. Cross-sections in four y-planes for the pilot and three Main stages, and six cross-sections in the (axial) x-planes near the dump plane, are plotted in figures 3 and 4, respectively. Figure 3 shows strong primary recirculation zones behind all the pressure atomizer elements, and no recirculation behind any of the air-blast elements. Figure 4 shows that these recirculation regions behind the four pressure atomizers of Main 1 stage are present at the 20mm plane downstream of the dump plane. The experimental measurement plane (114mm) is also shown in figure 4, with a central core of positive axial flow. A five hole-probe was used in the experiment to collect and report emissions data at the 114mm location [Tacina 2014].Figure 3. Axial velocity (m/s) contours for the four planes in line with the Pilot and Main element arrays for Case A.Figure 4. Axial velocity (m/s) contours for six axial planes for Case A; dump plane is at 0mm. D o w n l o a d e d b y C R A N F I E L D U N I V E R S I T Y o n M a y 11, 2016 | h t t p ://a r c .a i a a .o r g | D O I : 10.2514/6.2014-3529Temperature Predictions with the NCC for Case A (Low Power)Figures 5 and 6 show temperature contour comparisons for NCC predictions for Case A. Note that only the Pilot (see Y=0 plane of fig 5) and Main 1 stage (pressure atomizers at Y=0.014 and 0.034) were fueled for this case. The elongated flame behind the pilot (at FAR=0.065) and the shorter flames behind the Main 1 elements (FAR=0.045) are clearly seen in figure 5. The flames for all the PA elements sit in the diverging section of their respective venturi cups. Figure 6 shows that the individual flames from the Main 1 elements start merging with the Pilot flame at 20mm downstream of the dump plane. The 114mm location shows a very distinct hot region along the centerline, and the temperature gradually tapers towards the walls. This trend matches the five-point probe measurements by the traversing probe in the NASA CE5 test cell. The NCC computed mass-averaged value of temperature at the 114mm location is 1082K, as compared to the CEA computed (equilibrium) T 4 of 1165K. The experiment does not measure temperature but uses the CEA computed T 4 (based on FAR, T 3 and P 3) as the ‘experimental’ reference value for T 4.Figure 5. Temperature (K) contours for the four planes in line with the Pilot and Main element arrays for Case A.Figure 6. Temperature (K) contours at six axial planes for Case A; dump plane is at 0mm. D o w n l o a d e d b y C R A N F I E L D U N I V E R S I T Y o n M a y 11, 2016 | h t t p ://a r c .a i a a .o r g | D O I : 10.2514/6.2014-3529NO mass-fraction Predictions with the NCC for Case A (Low Power)Figures 7 and 8 shows NO mass-fraction contour comparisons for NCC predictions for Case A. The NO production from the pilot (Y=0 plane) dominates the overall NO for the system, as seen in figure 7 (Y=0) and the high NO values along the centerline in figure 8 (+10mm, +20mm, +114mm). The 114mm location shows a very high NO concentration at the centerline, which gradually tapers towards the walls. This trend matches the five-point probe measurements by the traversing probe in the experiment. The NCC predicted mass-averaged values of EINOx, EICO and EIHC at the exit are 1.87, 1.7 and 22.8, as compared to the experimental values of 2.81, 94.3 and 11.2, respectively. The EICO prediction is very poor, while the EIHCand EINOx predictions are reasonable.Figure 7. NO mass-fraction contours for the four planes in line with the Pilot and Main element arrays for Case A.Figure 8. NO mass-fraction contours at six axial planes for Case A; dump plane is at 0mm. D o w n l o a d e d b y C R A N F I E L D U N I V E R S I T Y o n M a y 11, 2016 | h t t p ://a r c .a i a a .o r g | D O I : 10.2514/6.2014-3529Reacting Flow: Axial Velocity Predictions with NCC for Case B (Medium Power)Figures 9 and 10 show axial velocity contour comparisons for NCC predictions for a typical Case B computation. Cross-sections in four y-planes for the pilot and three Main stages, and six cross-sections in the (axial) x-planes near the dump plane, are plotted in figures 9 and 10, respectively. Figure 9 shows attached, primary recirculation zones behind all the PA elements (pilot, Main 1), and weak, corner recirculation zones in the venturis of the AB elements (Main 2, Main 3). Figure 10 shows that the recirculation (dark blue) regions behind the Pilot and the four PA elements of Main 1 remain strong at 10mm and 20mm downstream of the dump plane. The experimental measurement plane at 114mm shown in figure 10 shows a fairly well mixed out distribution of axial velocity at this location.Figure 9. Axial velocity (m/s) contours for the four planes in line with the Pilot and Main element arrays for Case B.Figure 10. Axial velocity (m/s) contours for six axial planes for Case B; dump plane is at 0mm. D o w n l o a d e d b y C R A N F I E L D U N I V E R S I T Y o n M a y 11, 2016 | h t t p ://a r c .a i a a .o r g | D O I : 10.2514/6.2014-3529Temperature Predictions with the NCC for Case B (Medium Power)Figures 11 and 12 show temperature contour comparisons for NCC RANS predictions for the Case B computation. The attached flames behind the pressure-atomizing pilot (FAR=0.026) and Main 1 elements (FAR=0.026) are distinct from the weaker, detached flames behind Main 2 and Main 3 airblast elements, as shown in figure 11. Figure 12 shows that the flames from the four Main 1 elements dominate the heat release at the 10mm and 20mm locations. The 114mm location shows some non-mixedness in thetemperature distribution at this location. The computed mass-averaged value of T 4 at the exit plane of the computational domain is 1698K, which is 3% higher than the equilibrium temperature of 1651K computed by the CEA code.Figure 11. Temperature (K) contours for the four planes in line with the Pilot and Main element arrays for Case B.Figure 12. Temperature (K) contours at six axial planes for Case B; dump plane is at 0mm. D o w n l o a d e d b y C R A N F I E L D U N I V E R S I T Y o n M a y 11, 2016 | h t t p ://a r c .a i a a .o r g | D O I : 10.2514/6.2014-3529NO mass-fraction Predictions with the NCC for Case B (Medium Power)Figures 13 and 14 shows NO mass-fraction contour comparisons for NCC RANS predictions for the CaseB conditions. The NO production from pilot element dominates the overall NO for the system, as seen in figure 13 (Y=0.0m, top-left) and figure 14 (see central areas of -10mm, +10mm, +20mm planes). The +20mm location in figure 14 shows that the NCC predicts two high NO regions adjacent to the central pilot. These regions correspond to the Main1 stage, which is fueled by pressure atomizing injectors. The Main2 and Main3 stages, fueled by airblast injectors, produce insignificant amount of NO, even though they are operating at the same FAR as the pilot and Main1 stage. The NCC prediction at 114mm(experimental data location), shows an unmixed profile for NO mass-fraction. The NCC predicted mass-averaged value of EINOx is 5.4, which is 40% higher than the experimentally reported EINOx value of 3.8.Figure 13. NO mass fraction contours for four planes in line with Pilot and Main element arrays (Case B).Figure 14. NO mass fraction contours at six axial planes for Case B; dump plane is at 0mm. D o w n l o a d e d b y C R A N F I E L D U N I V E R S I T Y o n M a y 11, 2016 | h t t p ://a r c .a i a a .o r g | D O I : 10.2514/6.2014-3529Reacting Flow: Axial Velocity Predictions with NCC for CASE C (High Power)Figures 15 and 16 show axial -velocity contour comparisons for NCC predictions for a simulated CASE C condition. Cross-sections in four y-planes for the pilot and three Main stages, and six cross-sections in the (axial) x-planes near the dump plane, are plotted in figures 15 and 16, respectively. Figure 15 shows attached, strong primary recirculation zones behind all the PA elements (pilot, Main 1), and no primary recirculation zones behind the AB elements (Main 2, Main 3). Some ‘corner’ recirculation zones are observed in the diverging sections of the AB venturis. Figure 16 shows that the recirculation (dark blue) regions behind the four PA elements of Main 1 remain strong at +20mm downstream, with a weaker recirculation zone behind the Pilot. The 114mm plane is also shown in figure 16, with a weak central core of positive axial flow surrounded by higher axial velocity near the four Main-3 element locations.Figure 15. Axial velocity (m/s) contours for the four planes in line with the Pilot and Main element arrays for Case C.Figure 16. Axial velocity (m/s) contours at six axial planes for Case C; dump plane is at 0mm. D o w n l o a d e d b y C R A N F I E L D U N I V E R S I T Y o n M a y 11, 2016 | h t t p ://a r c .a i a a .o r g | D O I : 10.2514/6.2014-3529Temperature Predictions with the NCC for CASE C (High Power)Figures 17 and 18 show temperature contour comparisons for NCC predictions for Case C. The attached flames behind the pilot (FAR=0.026) and the much stronger, attached flames behind the Main 1 elements (FAR=0.035) are distinct from the weaker, detached flames behind Main 2 and Main 3 (FAR=0.035) airblast elements (see figure 17). The four Main 1 elements dominate the heat release at the 10mm and 20mm locations downstream of the dump plane, with much smaller heat release from Main2 and Main3 (figure 18). The 114mm location shows four very distinct cold spots near the corners, and considerable non-mixedness in the temperature distribution. The computed mass-averaged value of temperature at the exit is 1922K, which compares excellently with the CEA computed T 4 of 1916K.Figure 17. Temperature (K) contours for the four planes in line with the Pilot and Main element arrays forCase C.Figure 18. Temperature (K) contours at six axial planes for Case C; dump plane is at 0mm. D o w n l o a d e d b y C R A N F I E L D U N I V E R S I T Y o n M a y 11, 2016 | h t t p ://a r c .a i a a .o r g | D O I : 10.2514/6.2014-3529。
Journal of Measurement Science and InstrumentationJournal of Measurement Science and Instrumentation(CN14-1357/TH,ISSN1674-8042)is published by North University of China.It is a comprehensive academic journal,aiming to present scientific research papers in the fields of measurement science and instrumentation,including general principles,technologies and instrumentation of measurement and control applied to such academic and industrial fields as mechanics,electric and electronic engineering,magnetics,optics,chemistry,biology,and so on.We accept papers on the following topics:1)General theory and principles of measurement science and instrumentation,including measurement methods and practical developments such as precision measurements,new measurement principles, advanced technology of measurement,new generation instruments and systems.2)New techniques and systems in measurement:nanometer test techniques,remote test and calibration, automated test&diagnostics systems,calibration&self-calibration,virtual measurement systems,non-invasive measurement systems,distributed measurement systems.3)Measurement science and techniques in biological applications,medical and life sciences,material research.4)Measurement methods and apparatus for astronomy and astrophysics.5)Sensors,sensor systems and their applications(sensing principles and mechanisms,optoelectronic sensors,mechanical sensors,thermal sensors,magnetic sensors,optical sensors,fibre optic sensors,smart sensors,chemical and biological sensors,wireless network sensors,sensors for extreme environments,etc.). 6)Signal processing techniques(A/D&D/A converters,data acquisition,data output,signal transmission, data preprocessing&post-processing,image processing&pattern recognition,signal detection& classification,inverse problems&signal reconstruction,etc.).7)Optical and laser based techniques(optical metrology;optical methods for process control,optical microelectromechanical systems,optical techniques in micro-mechanics;microscopy and adaptive optics;laser material processing,laser beam delivery and diagnostics,laser remote sensing and environmental monitoring, laser safety,laser applications in medicine and biology,etc.).8)Imaging techniques(microscopy,tomography,holography,THz technology,etc.).9)Spectroscopy(optical,acoustic,dielectric principles,mass spectroscopy,fluorescence,x-ray,engineering applications of spectroscopy,etc.).10)Novel instrument systems and technology(virtual instrument system,design and manufacture of instruments,detectors,instruments and methods for non-destructive tests,instruments and techniques for dosimetry,monitoring and radiation damage,etc.).We accept manuscripts that meet the requirements stated in the instructions below:1)Contribution article must be the original work and has not been published before or under consideration for another publication,or a translation version of any paper published in other language.2)The manuscript should contain items arranged as follows:title,author,affiliation,abstract(including objective,method,result and conclusion),keywords(between5and8keywords),text body, acknowledgements,references・3)You can send your paper as an attachment by e-mail to************.cn or nu可************.4)It is necessary for Chinese authors to fill the Inquiry Form on the Protection of Secrets(download from the website:/tgfs.htm)about your manuscript and stamp the official seal of security department on it,and then send it to the Editorial Office when submitting the paper.Copyright:As the manuscript has been accepted,all authors should transfer the copyright of the article, including that of the printing and online version,to the publisher.This will ensure the widest possible dissemination.。
Measurements and simulations of transient characteristicsof heat pipesJarosław Legierski a ,Bogusław Wie ßceka,*,Gilbert de MeybaTechnical University of Ło ´dz ´,Institute of Electronics,ul.Wo ´lczanska 223,90-924Ło ´dz ´,PolandbUniversity of Ghent,Department of Electronics and Information Systems,Sint Pieternieuwstraat 41,9000Ghent,BelgiumReceived 1April 2004;received in revised form 13May 2005Available online 5August 2005AbstractThe rejection of heat generated by components and circuits is a very important aspect in design of electronics sys-tems.Heat pipes are very effective,low cost elements,which can be used in cooling systems.This paper presents the modelling and measurements of the heat and mass transfer in heat pipes.The physical model includes the effects of liquid evaporation and condensation inside the heat pipe.The internal vapour flow was fully simulated using compu-tational fluid dynamics.The theory has been compared with experimental measurements using thermographic camera and contact thermometers.The main purpose of this study is to determine the effective heat pipe thermal conductivity in transient state during start up the pipe operation and temperature increase.Ó2005Elsevier Ltd.All rights reserved.1.IntroductionThe investigations of heat pipes and their applica-tions into thermal management are known for years,but lately they have become more attractive to transport and dissipate heat in electronics [2,3,6,7].It is mainly due to the high effectiveness,small size and compact con-struction.Most of the previous research was made for static thermal conditions,and it resulted in getting the effective thermal resistance [6,7,9–11].Heat pipe is a passive heat transfer device with a high effective thermal conductivity.The device is built up as hermetic container that is partially filled with a fluid (typically water for electronic cooling design).The inner surfaces of the heat pipe have a capillary wicking mate-rial.The liquid in the pipe saturates the capillary wick.When one heats up the end of the pipe,the liquid at this place evaporates.Because evaporated medium in the heated section has a higher pressure than in the other part of the pipe,the vapour moves to the colder section,where it is cooled down and condensed.The condensed liquid is then transported back to the hot area thought the capillary wicks.This process causes a transport of heat and mass along the pipe [7,12].In such a heat sink one can define three sections:–evaporator,–condenser,–adiabatic region in the middle of the pipe.Usually,the heat pipe is modelled as a very high ther-mal conductivityelement (k =50000W/m K)[5].The purpose of the measurements and simulations presented in this paper is to determine the effective thermal0026-2714/$-see front matter Ó2005Elsevier Ltd.All rights reserved.doi:10.1016/j.microrel.2005.06.003*Corresponding author.E-mail address:wiecek@p.lodz.pl (B.Wie ßcek).Microelectronics Reliability 46(2006)109–115/locate/microrelconductivity of the heat pipe during the transient process of temperature rise in a typical cooling system (Fig.1).2.Measurements of transient characteristicsThe experimental layout for the measurement of heat pipe effective thermal conductivity is presented in Fig.2.The heat pipe was placed between two containers.One of them (container 2)was filled with water at ambient temperature.Container 1was initially empty.The ther-mal input function was done by filling in container 1with hot water (about 90°C).The measurement proce-dure is started when the hot water reached the heat pipe edges.The results of the temperature measurements were saved on the hard disc in the computer as a film from the thermal camera.The investigated heat pipe was insulated by low conductive material in order to keep adiabatic condition in the middle region of the heat sink and to have the considerable amount of energytransported to the cold container.At the beginning and the end of heat pipe,the insulation has been re-moved in order to measure temperature T B and T E by the thermographic camera.Additionally,the tempera-ture of the hot and cold water was measured by two con-tact thermometers,as shown in Fig.2.During the measurement,the temperatures of the fol-lowing characteristic thermal points at the surface of the heat pipe and in containers were recorded:T H mean temperature of a thermal source in hot water container 1T B temperature at the beginning of the heat pipe (evaporator)T E temperature at the end of the pipe (condenser)T Cmean temperature of water in cold container 2Temperature results obtained from the measurement (experiment)and simulation were described with the index exp and sim ,respectively.Fig.1.Heat pipe construction and principle ofoperation.Fig.2.Transient thermal characteristic measurement setup.110J.Legierski et al./Microelectronics Reliability 46(2006)109–115A water heat pipe,of 200mm length,4mm diameter and 25mm evaporator and condenser lengths was used in the investigation [8].Two containers,hot and cold ones were insulated by foamed polystyrene of 10mm thickness.Hot and cold reservoirs were filled with about 2l and 250ml of water,respectively.The transient ther-mal characteristics obtained from the measurement are shown in Figs.3–5display the temperature fields as observed with the thermal camera.As shown in Fig.3,during the first 20s,the heat pipe is in a transient state.This period can be explained as a time constant of the heat pipe itself.During starting time the evaporation and condensing rates are rising and will reach the steady state conditions in which the heat pipe operates optimally.One can notice a certain anomaly in Fig.3.After some time (25s)from starting point,the temperatures T B,exp (t )and T E,exp (t )are slightly reduced.It seems that the temperature profile gets the maximum value at time t =25s.It is due to a slight cooling of the hot water container T H,exp (t ),what is clearly visible in Fig.3(upper curve).Water temperature in cold container was changing from T C,exp (0)=22.7°C to T C,exp (130s)=23.2°C dur-ing 130s of the experiment.In order to calculate the average power transported through the heat pipe into the cold container filled with 250ml of water,the tem-perature measurement was performed using digital ther-mometer equipped with thermocouples.Assuming the accuracy of the temperature readout better then 0.1°C,one can evaluate the relative accuracy of temper-ature and average power measurement at about 10%.Power transported through the pipe during the experi-ment was evaluated as P ¼m Ác S ÁT C ;exp ð130ÞÀT C ;exp ð0Þ%4Wð1Þwhere m is the water mass in cold container,c S is the specific heat of water.The temperature T B,exp and T E,exp was measured with a thermographic camera,while T H,exp was recorded using digital electronic thermometer.In steady state,temperature drop T B,exp ÀT E,exp along the pipe is about 0.7°C,which is an experimental proof of the high thermal conductivity of such a device.3.Simulations of transient processes in the heat pipe The simulation of heat and mass transfer inside the pipe was performed using CFD (Computer Fluid Dynamics)software—Fluent 6.0Òwith pre-processor Gambit 2.0Ò.3D model of the heat pipe consists of three elements:evaporator,condenser and central adia-batic ing C-based programmingtools—UserFig. 3.Thermal characteristics of the heat pipe (measure-ments).Fig.4.Temperature field for time t =8s.Fig.5.Temperature field for time t =39s.J.Legierski et al./Microelectronics Reliability 46(2006)109–115111Define Functions (UDF),one can easily introduce evap-oration and condensation into the modelling.The sur-faces where evaporation takes place were modelled using build-in mass-flow-inlet elements in Fluent Òenvi-ronment.In the condenser,an outflow element was implemented.The evaporation and condensation heat fluxes (negative and positive)were placed on two small areas in the evaporator and condenser section to simu-late absorption and dissipation of energy for both ends of the heat pipe as shown in Fig.6.The thermal excita-tion function was implemented as temperature (T S )in UDF file,equal to the temperature of the evaporator section obtained from the measurement T B,exp .Cold water container was simulated as a box filled with water at T C temperature.The rectangular mesh with 10666elements (cells)was used in the simulations for all parts of the cooling system.Wick material and its influence of heat transfer in a pipe were neglected.In the preliminary investigations presented in this paper,one assumed that the evapora-tion,vapour transport and condensation are phenomena of the highest significance in the operation of the heat pipe,and the effective thermal conductivity can be eval-uated without considering the liquid return in a wick.The further more detailed research can include a wick in the simulations.The mass transfer in the pipe was described by well-known the Navier–Stokes equations,while heat transfer was solved using Fourier–Kirchoffexpressions.Hertz–Knudsen model was applied to determine mass flux during the evaporation [1,4]u EVAP ¼e ÁffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiM 2Áp ÁR r p 0ffiffiffiffiffiT 0p Àp p ffiffiffiffiffiffiT p p!ð2Þwhere u EVAP is the mass flux during the evaporation[kg/m 2s];M is the molar mass of liquid [kg/mol];R =8.314is the universal gas constant [J/mol K];p 0is the saturation vapour pressure at the interface [Pa];p p is the saturation vapour pressure in the end region of a pipe [Pa];e is the evaporation coefficient (06e 61);T 0is the temperature of the vapour at the interface[K];T p is the temperature of the vapour in the end range of a pipe [K].The evaporation coefficient e describes the rate of evaporation.A smooth surface will have a low e value,while a rough surface will have a higher value because more ‘‘kernels’’for initiating evaporation are available.Because air is not presented in the simulated heat pipe,the saturation vapour pressure over the interface is equal to the static pressure at this place.The evaporating liquid absorbs a thermal heat flux U EVAP from the ambient:U EVAP ¼u EVAP Áh fgð3Þwhere h fg is the latent heat,describing vaporization [J/kg].One should assume that the entire evaporating mass flux is condensing in condenser section,what leads to the equation u COND ¼u EVAPð4ÞThe heat flux U COND corresponding to phase change in the condenser can be expressed as U COND ¼u COND Áh fgð5Þwhere u COND is the condensation mass flux [kg/m 2s].Simulations have been performed for different values of evaporation coefficient e .Although theadiabaticFig.6.Model of the heatpipe.Fig.7.Thermal characteristic for e =1.58·10À3and a =0.25.112J.Legierski et al./Microelectronics Reliability 46(2006)109–115section of the heat pipe is insulated thermally,the heat losses may not be neglected as compared to the power of 4W transported into the cold container.Therefore,a heat transfer coefficient a along the adiabatic section has been introduced in the simulation.The coefficient a stands both for the conduction in the insulation and the natural convection from the outer surface to the ambient.The simulation results were compared with the measurement.In Fig.7,the thermal characteristics obtained from the simulation are shown.Figs.8–10graphically present the results of the simulations (tem-perature and vapour velocity distributions).The images present only one half of the simulated area because of the cylindrical symmetry of the modelled heat pipe.paring simulation results with measurement The evaporation coefficient values e used in the sim-ulations were between 10À4and 10À2.Figs.11–13pres-ent the simulation results for increased evaporation rate.The time constant values obtained from the simu-lation agree very well with the experimental data,as shown in Figs.12and 13.The temperature valuesalsoFig.8.Temperature in condenser and cold water container for e =1.58·10À3and t =90s.Fig.9.Velocity distribution in evaporator section [m/s]:e =1.58·10À3and t =20s.Fig.10.Velocity distribution in condenser section [m/s]:e =1.58·10À3and t =90s.paring results of simulation with measurement for e =10À4and a =0.25(P =0.79W).J.Legierski et al./Microelectronics Reliability 46(2006)109–115113agree very well.The only discrepancy is for the time de-lay for T E,sim if the evaporation coefficient is e =10À4(Fig.11)or e =1.58·10À3(Fig.12).For e =10À2(Fig.13)this delay has disappeared but for longer times,the temperature drop T B,sim ÀT E,sim =0.09°C,which is much less than the experimentally found T B,exp ÀT E,exp =0.7°C.In this case the power transported by the heat pipe is six times higher (24W)than in the exper-iment (4W).This makes e =10À2an unacceptable va-lue.For e =1.58·10À3(Fig.12)the simulated temperature drop T B,sim ÀT E,sim =0.73°C is in perfect agreement with the experimental observations,as well as heat transported by the pipe (4.046W)corresponds precisely with the measurements.Hence one must con-clude that e =1.58·10À3should be the correct evapora-tion coefficient.The delay between T B,sim and T E,sim still visible in Fig.12but it is relatively small and acceptable.Resuming the best results obtained by comparing simulation with measurement were for evaporation coef-ficient e =1.58·10À3.For e =10À4heat pipe is ‘‘too slow’’and evaporator heating time to temperature T =55°C is 50s.This result is opposite to experimental data.For evaporation coefficient value e =10À2temper-ature gradient along the pipe is too small what results in disagreement of the power transported through the pipe.5.Effective thermal conductivity of the pipeThe effective thermal conductivity of the pipe k was determined by using the simulation result as a function of time.By definition the thermal conductivity can by expressed as k ¼À/L eff D Tð6Þwhere /is the heat flux along the pipe [W/m 2];D T is the temperature gradient;L effis the effective heat pipe length [m].Effective heat pipe length can be approximated by the equation [5]L eff ¼L C 2þL A þL E2ð7Þwhere L C ,L E ,L A are condenser,evaporator and adia-batic section lengths,respectively.Both heat flux /and the total power P T transferred to the cold container as well as temperature gradient D T were calculated by Fluent software.The heat flux and power were evaluated in the condenser section on the contact surface between the pipe and water,as shown in Fig.6.One can notice that the total power transferred to the cold reservoir includes both the power generated by the vapour condensation and the heat flux on the contact surface between the hot vapour to the cold water.Results are presented in Table 1.One may argue whether a time dependent thermal conductivity is the right way to represent the phenome-non observed here.A network involving capacitances and resistors is more common to deal with transientTable 1Effective thermal conductivity k (e =1.58·10À3)t (s)1510152030405060708090P T (W)0.39 1.41 2.27 4.25 5.54 5.15 4.34 4.20 4.02 4.38 4.23 4.10D T (°C) 4.023 2.3353 1.7582 1.2396 1.01850.60280.79860.78080.67840.63230.77060.7607U (W/m 2)7823281064517684608110129102387863908355479939871758414581508k (W/m K)340.32106.24496.511944.518922.429724.318931.118726.820621.124127.119109.018751.1paring results of simulation with measurement for e =1.58·10À3and a =0.25(P =4.05W).paring results of simulation with measurement for e =10À2and a =0.25(P =24.61W).114J.Legierski et al./Microelectronics Reliability 46(2006)109–115phenomena.However,in the beginning of heat pipe operation,the velocity of the vapour particles is quite small resulting in lower thermal conductivity.One must also bear in mind that the thermal capacity of the va-pour in the adiabatic section is rather small,especially at low pressure.Hence,it makes sense to represent the heat pipe by a time dependent conductor or alternatively by an equivalent time dependent thermal resistance net-work.One can remark that the maximum measured value found here is around25000W/m K which is less than the value recommended in Ref.[5].This may be due to the fact that the heat pipe was placed horizontally in our experiments.It is well known that heat pipes operate more efficiently in the vertical position,in which the evaporating section is at the bottom.It is also clear from Fig.14that the heat pipe,or at least its adiabatic section becomes a very good conductor only after20s.6.ConclusionAs a result of the measurements and the simulations of heat pipe,one was able to determine the evaporation coefficient.The value is about1.58·10À3for a pipe pre-sented in this paper.The effective thermal conductivity of the heat pipe was found to be time dependent with a maximum value in the range:15000–30000W/m K, which is close to the value obtained from technical paper e.g.[8].The effective thermal conductivity reaches its steady state value after%20–30s.This result was ob-tained for the water heat pipe of200mm length and 4mm diameter,but general model and overall approach of measurement and simulation can be extended to any other heat pipe.References[1]Eames IW,Marr NJ,Sabir H.The evaporation coefficientof water:a review.Int J Heat Mass Transfer1997;40(12): 2963–73.[2]Terpstra M,Van Veen JG.Heat pipes construction andapplications—study of patents.London,NY:Elsevier Applied Science;1987.[3]Garner SD.Heat pipes for electronics cooling applications.Electronics Cooling,September1996.[4]Wajman T,Wiecek B,Felczak M.Vaporation fromcapillaries in cooling electronics devices—evaporation coefficient evaluation.Mixed Design of Integrated Cir-cuits and Systems(MIXDES),Lodz/Poland,June26–28, 2003.[5]Thyrum G.Critical aspects of modelling heat pipe assistedheat sinks.Thermacore Inc.,.[6]DeHoffR,Grubb K.Heat pipe application guidelines.Thermacore Inc.,.[7]Toth J,DeHoffR,Grubb K.Heat pipes:the silent wayto manage desktop thermal problems.Thermacore Inc., I-THERM Conference Seattle,1998.[8]Miniature Heat Pipes HP-NB Product Data Guide,Thermacore Inc.,4/8/1998.[9]Legierski J,Wiecek B.Steady state analysis of coolingelectronic circuits using heat pipes.IEEE Trans Comp Pack Technol2001;24(4).[10]Wiecek B.Heat Pipe Cooling Fins for Integrated Circuits,Workshop Microtechnology and Thermal Problems in Electronics,MicrothermÕ98,Zakopane,21–27September 1998.[11]Legierski J,Wiecek B.Modelling and measurements ofpower circuits cooled by heat pipes.Proc.Quantitative Infrared Thermography,Reims,July18–21,2000. [12]Legierski J,Wiecek B.Transient thermal characteristics ofheat pipes.Measurements and Simulation.Proc.Quanti-tative Infrared Thermography,Dubrovnik,September24–27,2002.J.Legierski et al./Microelectronics Reliability46(2006)109–115115。