Thermal Modelling,Simulation and Evaluation of a high power battery cell for automotive applications
ROHM BD9G500EFJ-LA DC DC 转换器解决方案用户指南说明书
User’s GuideROHM Solution SimulatorDC/DC Converter BD9G500EFJ-LA Thermal SimulationThis document contains electrical simulations of the DC/DC converter BD9G500EFJ-LA and introduces and describes the use of a simulation environment that allows simultaneous thermal simulation of devices including Schottky Barrier Diodes (SBD: RB088BM100TL). By changing the parameters of the components, it is possible to simulate a wide range of conditions.1 Simulation circuitFigure 1. Simulation circuit (BD9G500EFJ-LA)In Figure 1, the area within the green line shows the thermal simulation circuit and the rest of the figure shows the electrical simulation circuit.This circuit is an application circuit based on a 1-channel buck DC/DC converter with a current output of up to 5 A using the BD9G500EFJ-LA.The thermal simulation circuit feeds the device losses and SBD losses calculated in the electrical simulation into the thermal simulation model, and calculates the IC and SBD temperatures.2 Simulation methodSimulation settings such as simulation time and convergence options can be set from “Simulation Settings” shown in Figure 2, and the initial simulation settings are shown in Table 1.If you are having problems with the convergence of the simulation, you can change the advanced options to fix the problem. The simulation temperature and various parameters of the electrical circuit are defined in “Manual Options”.Figure 2. Simulation Settings and executionTable 1. Initial values for Simulation SettingsParametersInitial valuesRemarksSimulation Type Time-Domain Do not change the simulation type End time7 msecs Advanced Options More SpeedManual Options .PARAM …See Table 2 for details3 Simulation conditions3.1 Definition of parametersThe parameters for the components shown in blue in Figure 3 are defined in the manual options as they need to be set in the simulation conditions. Table 2 shows the initial values for each parameter. These values are written in a text box in the “Manual Options” section of the simulation settings, as shown in Figure 4.Figure 3. Definition of component parametersSimulation SettingsSimulateTable 2. Simulation conditionsParameters VariablenamesInitial values Unit DescriptionTemperature Ta25°C Ambient temperatureVoltage V_VIN48V Input voltage Set in the range of 7 to 76 VVoltage V_VOUT5V Output voltage Set in the range of 1 V to (0.97 × V_VIN)Current I_IOUT1A Output current 5 A (MAX)Inductance L_PRM33µH Smoothing inductorFigure 4. Definition of parameters3.2 Setting of component constantsFor the method of setting switching frequency, output LC filter constant, output voltage, etc., refer to “Selection of Components Externally” in the data sheet or the calculation sheet.BD9G500EFJ-LA Data sheetCalculation-Sheet For The Circuit Theoretical Formula – BD9G500EFJ-LAWrite parameters3.3 Thermal circuitThe “BD9G500EFJ_LA” symbol in Figure 5 is the thermal simulation model of the BD9G500EFJ-LA. The nodes shown inred in Figure 5 can be used to check the temperature of the junction, the mold surface and the FIN surface. Detailedinformation for each node is shown in Table 3.You can check the temperature bytouching the red node with a probeFigure 5. BD9G500EFJ-LA thermal simulation modelTable 3. Description of the nodes in Figure 5Node name DescriptionBD9G500EFJ_Tj Monitors the junction temperature of BD9G500EFJ-LASBD_Tj Monitors the junction temperature of RB088BM100BD9G500EFJ_Tt Monitors the top center temperature of BD9G500EFJ-LASBD_Tt Monitors the top center temperature of RB088BM100SBD_Tfin Monitors the FIN center temperature of RB088BM1003.4 Selecting a thermal simulation modelThere are a number of thermal simulation models to choose from and their components are shown in Table 4. Figure 6 shows how to select one. First, right-click on the BD9G500EFJ-LA component and select “Properties”. In the “Property Editor”, set the value of the “SpiceLib Part” to the value you selected from Table 4 to change the thermal simulation model.Figure 6. How to select a thermal simulation modelTable 4. List of available componentsComponent name SpiceLib Part valueDescriptionBD9G500EFJ-LA2s Thermal simulation model for a two-layer board 2s2pThermal simulation model for a four-layer boardFor more information on the board, see “Reference: About the BD9G500EFJ-LA thermal simulation model” on page 7.Changing the value of the SpiceLib Part allows you to select a different thermal model4 Links to related documents4.1 ProductsBD9G500EFJ-LARB088BM1004.2 User’s GuideSingle Buck Switching Regulator BD9G500EFJ-LA EVK User’s GuideReference: About the BD9G500EFJ-LA thermal simulation modelAn image of the 3D model used to create the thermal simulation model is shown in Figure A. Structural information is also shown in Table A.Figure A. BD9G500EFJ-LA 3D imageTable A. Structural informationStructural parts DescriptionBoard outline dimensions114.3mm × 76.2mm, t=1.6mmBoard material FR-4Layout pattern Refer to “Single Buck Switching Regulator BD9G500EFJ-LA EVK User’s Guide”2-layer board Layer structure Top Layer : 70µm ( 2oz ) Bottom Layer : 70µm ( 2oz )4-layer board Layer structure Top Layer : 70µm ( 2oz )Middle1 & Middle2 Layer : 35µm ( 1oz ) Bottom Layer : 70µm ( 2oz )BD9G500EFJ-LA (HTSOP-J8)RB088BM100 (TO-252)BoardNoticeROHM Customer Support System/contact/Thank you for your accessing to ROHM product informations.More detail product informations and catalogs are available, please contact us.N o t e sThe information contained herein is subject to change without notice.Before you use our Products, please contact our sales representative and verify the latest specifica-tions :Although ROHM is continuously working to improve product reliability and quality, semicon-ductors can break down and malfunction due to various factors.Therefore, in order to prevent personal injury or fire arising from failure, please take safety measures such as complying with the derating characteristics, implementing redundant and fire prevention designs, and utilizing backups and fail-safe procedures. ROHM shall have no responsibility for any damages arising out of the use of our Poducts beyond the rating specified by ROHM.Examples of application circuits, circuit constants and any other information contained herein areprovided only to illustrate the standard usage and operations of the Products. The peripheral conditions must be taken into account when designing circuits for mass production.The technical information specified herein is intended only to show the typical functions of andexamples of application circuits for the Products. ROHM does not grant you, explicitly or implicitly, any license to use or exercise intellectual property or other rights held by ROHM or any other parties. ROHM shall have no responsibility whatsoever for any dispute arising out of the use of such technical information.The Products specified in this document are not designed to be radiation tolerant.For use of our Products in applications requiring a high degree of reliability (as exemplifiedbelow), please contact and consult with a ROHM representative : transportation equipment (i.e. cars, ships, trains), primary communication equipment, traffic lights, fire/crime prevention, safety equipment, medical systems, servers, solar cells, and power transmission systems.Do not use our Products in applications requiring extremely high reliability, such as aerospaceequipment, nuclear power control systems, and submarine repeaters.ROHM shall have no responsibility for any damages or injury arising from non-compliance withthe recommended usage conditions and specifications contained herein.ROHM has used reasonable care to ensur e the accuracy of the information contained in thisdocument. However, ROHM does not warrants that such information is error-free, and ROHM shall have no responsibility for any damages arising from any inaccuracy or misprint of such information.Please use the Products in accordance with any applicable environmental laws and regulations,such as the RoHS Directive. For more details, including RoHS compatibility, please contact a ROHM sales office. ROHM shall have no responsibility for any damages or losses resulting non-compliance with any applicable laws or regulations.W hen providing our Products and technologies contained in this document to other countries,you must abide by the procedures and provisions stipulated in all applicable export laws and regulations, including without limitation the US Export Administration Regulations and the Foreign Exchange and Foreign Trade Act.This document, in part or in whole, may not be reprinted or reproduced without prior consent ofROHM.1) 2)3)4)5)6)7)8)9)10)11)12)13)。
PerkinElmer STA 8000 同步热分析仪简介说明书
IntroductionWhile determining the composition of alloys has traditionally been the domain of differential scanning calorimeters ordifferential thermal analyzers,1 the STA 8000 has shown itself to be capable of this type of demanding analysis as well. The requirements for melt analysis are accurate temperature and melt energy measurement and the ability to exclude oxygen – or if necessary, nitrogen – during the analysis. This note shows examples of two high temperature melting systems, including iron-nickel alloys which require oxygen exclusion.ExperimentalThe STA 8000 (Figure 1) provides the analysis of sample sizes typically in the 10 to 200 milligram range with the samples heated by a small furnace (~20 cc volume) capable of operating from 15 to 1600 degrees Celsius.2 The STA sensor is a double pan differential temperature sensor which is calibrated to generate data for heat flow to the sample with an accuracy of 5% or better. The weight sensor is located remotely below the furnace where it is isolated from sample decomposition products by inert purge through a narrow channel. This provides microgram-level weight change detectability of sample loss or oxidative gain. Unless otherwise noted the purge gas used for this note was nitrogen at a flow rate of 100 cc/min. The use of argon instead of nitrogen is supported in the Pyris ™ software which controls the flow rate of gas through the furnace chamber and provides for gas switching and flow rate changes.a p p l i c a t i o n n o t eAuthorsBruce Cassel Kevin P. Menard PerkinElmer, Inc. Shelton, CT USA Professor Charles Earnest Berry CollegeDepartment of Chemistry Mount Berry, GA USAUse of the STA 8000 Simultaneous Thermal Analyzer for MeltAnalysis of AlloysFigure 1. STA 8000 Simultaneous Thermal Analyzer.2temperature after correcting for thermal lag. This is done by using the slope of the leading edge of the DTA thermal curve of the pure primary component (Figure 4). The peak of the observed melting endotherm is the point at which the last of the solid material melts; however, since the latent heat of melting causes the temperature of the sample to lag behind the sensor, the peak maximum should be corrected for this effect. This thermal lag can be evaluated from the melting peak of a high purity sample, since when it is melting, the temperature of the sample remains constant at the melting point until all the sample has melted. Since the furnace temperature is increasing linearly, the leading edge of the peak should be linear. Notice that the peak shape for pure gold is linear and faithfully follows the theoretical shape expected from melting in a quantitative calorimeter. The inverse slope of that leading edge is just the ratio of apparent temperature change per unit heat flow change. Thus, the solidus of the alloy is seen to be the peak temperature minus the product of its peak height and the inverse slope (as calculated using Pyris software) of the leading edge of the melt of the pure material. The results obtained were in excellent agreement with literature.5The samples chosen for analysis included elementalmaterials of high purity and industrial samples of uncertain composition. The two component phase diagrams for these systems (which indicate the solidus and liquidus temperatures versus composition) are well documented in the literature, so that a determination of the melting onset from the thermal curve is sufficient to determine the composition. Conversely, the analyzer can be used to generate the phase diagram for accurately prepared, multiple component alloy systems.ResultsGold-Copper Alloys. With applicability to jewelry production, dental fillings, aerospace, and especially electronics, the gold-copper alloy is used as a reliable braising system to bond metallic surfaces or as an inert coating. The materials whose STA data is shown in Figure 2 were produced by weighing out reagent copper granules and 99.99% purity gold wire into the STA pan, then heating the mixture and holding isothermally above the gold melting point to allow alloy formation. The alloy was then cooled, then analyzed in the STA, all these steps being part of one analytical method using Pyris Software. As can be seen from Figure 2, small increases in the copper content result in significant changes in the melting point. Thus the STA, which has a temperature accuracy specification of better than a degree, is able to determine the composition of a material in a known alloy family. The two-component phase diagram for this system taken from the literature 4 is shown in Figure 3 with the placement of the data points (marked by vertical slashes, not to be taken as error bars).The extrapolated onset temperature from the leading edge of the DTA thermal curve can be taken as the determination of the solidus, the onset of melting when heating the solid solution. In compositions where there is a melting range, i.e., where the liquidus is measurably higher than thesolidus, the liquidus can be taken as the alloy melting peakFigure 2. Melting endotherms for gold and three gold-copper alloys.Figure 3. Gold-copper phase diagram showing compositions of endotherms in Figure 2.Figure 4. Determination of solidus and liquidus (1.7% Cu: 98.3% Au sample). The slope calculation from the gold melt (the higher temperature peak) allows lag correction for the peak parameter of the alloy.3Invar and Other Iron-Nickel Alloys. The Invar alloys of Iron and nickel have uniquely low coefficients of expansion which make them useful in fine watches, sensitive instrumentation, aerospace and especially in electronicapplications where small differences in expansion may cause failure. A sample of Invar 36 (36% nickel) was obtained from the PerkinElmer manufacturing group where it is used in construction in analytical instruments. The sample was analyzed in the STA 8000 at 20, 10 and 5 ˚C/min using a sample purge of 100 cc/min nitrogen. The resulting thermal curve can be seen in Figure 5. Using a strong magnet over the STA furnace allowed the Curie temperature (the temperature above which all attraction to a magnet disappears), to be recorded (Figure 6). At the end of the analysis the weight of the Invar sample was determined to not have changed by more than two tenths of a percent, an indication that any oxidation from back diffusion from air into the exit port was minimal at this purge rate. Figure 7 shows the DTA thermal curve for Invar when heated at 5 ˚C/min. This rate allows a better resolution of the solidus and liquidus, which were calculated from the peak maxima and corrected using the leading edge slope of melting reagent iron.SummaryThe STA 8000 is designed to give accurate melting characteristics over a temperature range from ambient to 1600 ˚C. The quality of the calorimetry is sufficient to allow quantitative corrections for thermal lag as is normally done in a differential scanning calorimeter. Using a purge rate of 100 cc/min of house nitrogen there was sufficient oxygen exclusion to limit oxidation of iron and iron alloys to less than 0.2% of original weight after extended analysis to 1600 ˚C. The use of argon purge and faster, or slower, purge rates is also supported. This should qualify the STA 8000 for quantitative characterization of alloy systems up to 1600 ˚C. Examples of using a DTA approach to characterize ternary alloy systems can also be found in the literature.6,7Figure 5. Invar36 at three heating rates, and cooling.Figure 6. Invar36 alloy at 10 ˚C/min showing the melt and the step up in apparent weight at the Curie temperature.Figure 7. Invar 36 at 5 ˚C/min showing values of the liquidus and solidus.For a complete listing of our global offices, visit /ContactUsCopyright ©2012, PerkinElmer, Inc. All rights reserved. PerkinElmer ® is a registered trademark of PerkinElmer, Inc. All other trademarks are the property of their respective owners.010824_01PerkinElmer, Inc. 940 Winter StreetWaltham, MA 02451 USA P: (800) 762-4000 or (+1) 203-925-4602References1. W.J. Boettinger, U.R. Kattner, K.-W. Moon and J.H.Perepezko, “DTA and Heat-Flux DSC Measurements of Alloy Melting and Freezing,” NIST ® Special Pub. 950-15 (2006).2. STA Specification sheet #010619_01.3. STA Brochure # 010452_01.4. H Okamoto, et. al., Bulletin of Phase Diagrams, Vol. 8, No. 5, 1987.5. M.H. Sloboda, Gold Bulletin 1971 Vol. 4, No. 1, Chamber of Mines of So. Africa. Pub.6. S.K. Lin, C.F. Yang, S.H. Wu, S.W. Chen “Liquidusprojection and solidification of the Sm-In-Cu ternary alloys”, Journel of Electronic Materials, 2008, Springer.7. S. Chen, H. Hsu, C. Lin, “Liquidus projection of the ternary Ag-Sn-Ni System” J. Material Research, 2004 Cambridge Univ Press.。
NPT型IGBT电热仿真模型参数提取方法综述_徐铭伟
电力自动化设备Electric Power Automation Equipment Vol．33No．1 Jan.2013第33卷第1期2013年1月0引言近年来，绝缘栅双极型晶体管IGBT（Insulated Gate Bipolar Transistor）因其不断改善的电压、电流承受能力和工作频率、功率损耗等性能指标而被广泛应用到机车牵引、开关电源、新能源发电等电能变换和处理领域中［1］，因此IGBT的可靠性受到国内外科研工作者的广泛关注。
研究表明，与IGBT器件结温（T j）相关的热循环过程和器件封装材料热膨胀系数不一致是致其故障的主要诱因［2-3］，IGBT的电热仿真模型可以估计结温的变化情况，从而可用于IGBT可靠性的评估。
国内外对IGBT的电热仿真模型开展了大量研究工作［4-6］，其中基于半导体物理并考虑自热效应（Self-heating）的IGBT A.R.Hefner器件模型［6］和反映其封装传热过程的Cauer网络［7-9］联合组成的IGBT 电热模型准确度较高，并已在Saber、Pspice等电路仿真软件中得到应用［10-11］，但是，仿真软件有限的器件模型库无法满足仿真需要，同时出于技术保密的缘故，半导体制造商并不会提供建立电热模型需要的模型参数，因此如何建立一种有效并准确的参数提取方法就显得十分必要。
IGBT电热仿真模型参数同半导体物理、器件以及封装结构直接相关，无法直接测量，只能通过一定的技术方法和手段获取。
一个有效的参数提取过程是获得有效的电热模型的前提条件；此外，实现模型参数的准确提取对于分析IGBT的性能、优化驱动电路的设计、指导其应用以及选型都具有重要意义。
在参数提取之后，有效性验证也至关重要，可以让使用者合理选择器件的工作范围。
由于非穿通（NPT）型IGBT目前在工业领域中已获得了广泛而成熟的应用［12］，本文将以其作为参数提取的研究对象。
本文从NPT型IGBT电热仿真模型的工作原理出发，首先将模型参数分为电参数和热参数两大类。
SOLIDWORKS Flow Simulation 产品说明书
OBJECTIVESOLIDWORKS® Flow Simulation is a powerful Computational Fluid Dynamics (CFD) solution fully embedded within SOLIDWORKS. It enables designers and engineers to quickly and easily simulate the effect of fluid flow, heat transfer and fluid forces that are critical to the success of their designs.OVERVIEWSOLIDWORKS Flow Simulation enables designers to simulate liquid and gas flow in real-world conditions, run “what if” scenarios and efficiently analyze the effects of fluid flow, heat transfer and related forces on or through components. Design variations can quickly be compared to make better decisions, resulting in products with superior performance. SOL IDWORKS Flow Simulation offers two flow modules that encompass industry specific tools, practices and simulation methodologies—a Heating, Ventilation and Air Conditioning (HVAC) module and an Electronic Cooling module. These modules are add-ons to a SOLIDWORKS Flow Simulation license. BENEFITS• Evaluates product performance while changing multiple variables at a rapid pace.• Reduces time-to-market by quickly determining optimal design solutions and reducing physical prototypes.• Enables better cost control through reduced rework and higher quality.• Delivers more accurate proposals.CAPABILITIESSOLIDWORKS Flow SimulationSOLIDWORKS Flow Simulation is a general-purpose fluid flow and heat transfer simulation tool integrated with SOLIDWORKS 3D CAD. Capable of simulating both low-speed and supersonic flows, this powerful 3D design simulation tool enables true concurrent engineering and brings the critical impact of fluid flow analysis and heat transfer into the hands of every designer. In addition to SOL IDWORKS Flow Simulation, designers can simulate the effects of fans and rotating components on the fluid flow and well as component heating and cooling. HVAC ModuleThis module offers dedicated simulation tools for HVAC designers and engineers who need to simulate advanced radiation phenomena. It enables engineers to tackle the tough challenges of designing efficient cooling systems, lighting systems or contaminant dispersion systems. Electronic Cooling ModuleThis module includes dedicated simulation tools for thermal management studies. It is ideal for companies facing thermal challenges with their products and companies that require very accurate thermal analysis of their PCB and enclosure designs.SOLIDWORKS Flow Simulation can be used to:• Dimension air conditioning and heating ducts with confidence, taking into account materials, isolation and thermal comfort.• Investigate and visualize airflow to optimize systems and air distribution.• Test products in an environment that is as realistic as possible.• Produce Predicted Mean Vote (PMV) and Predicted Percent Dissatisfied (PPD) HVAC results for supplying schools and government institutes.• Design better incubators by keeping specific comfort levels for the infant and simulating where support equipment should be placed.• Design better air conditioning installation kits for medical customers.• Simulate electronic cooling for LED lighting.• Validate and optimize designs using a multi-parametric Department of Energy (DOE) method.SOLIDWORKS FLOW SIMULATIONOur 3D EXPERIENCE® platform powers our brand applications, serving 12 industries, and provides a rich portfolio of industry solution experiences.Dassault Syst èmes, t he 3D EXPERIENCE® Company, provides business and people wit h virt ual universes t o imagine sust ainable innovat ions. It s world-leading solutions transform the way products are designed, produced, and supported. Dassault Systèmes’ collaborative solutions foster social innovation, expanding possibilities for the virtual world to improve the real world. The group brings value to over 220,000 customers of all sizes in all industries in more than 140 countries. For more information, visit .Europe/Middle East/Africa Dassault Systèmes10, rue Marcel Dassault CS 4050178946 Vélizy-Villacoublay Cedex France AmericasDassault Systèmes 175 Wyman StreetWaltham, Massachusetts 02451-1223USA Asia-PacificDassault Systèmes K.K.ThinkPark Tower2-1-1 Osaki, Shinagawa-ku,Tokyo 141-6020Japan©2018 D a s s a u l t S y s t èm e s . A l l r i g h t s r e s e r v e d . 3D E X P E R I E N C E ®, t h e C o m p a s s i c o n , t h e 3D S l o g o , C A T I A , S O L I D W O R K S , E N O V I A , D E L M I A , S I M U L I A , G E O V I A , E X A L E A D , 3D V I A , B I O V I A , N E T V I B E S , I F W E a n d 3D E X C I T E a r e c o m m e r c i a l t r a d e m a r k s o r r e g i s t e r e d t r a d e m a r k s o f D a s s a u l t S y s t èm e s , a F r e n c h “s o c i ét é e u r o p ée n n e ” (V e r s a i l l e s C o m m e r c i a l R e g i s t e r # B 322 306 440), o r i t s s u b s i d i a r i e s i n t h e U n i t e d S t a t e s a n d /o r o t h e r c o u n t r i e s . A l l o t h e r t r a d e m a r k s a r e o w n e d b y t h e i r r e s p e c t i v e o w n e r s . U s e o f a n y D a s s a u l t S y s t èm e s o r i t s s u b s i d i a r i e s t r a d e m a r k s i s s u b j e c t t o t h e i r e x p r e s s w r i t t e n a p p r o v a l .• Free, forced and mixed convection• Fluid flows with boundary layers, including wall roughness effects• Laminar and turbulent fluid flows • Laminar only flow• Multi-species fluids and multi-component solids• Fluid flows in models with moving/rotating surfaces and/or parts• Heat conduction in fluid, solid and porous media with/without conjugate heat transfer and/or contact heat resistance between solids• Heat conduction in solids only • Gravitational effectsAdvanced Capabilities• Noise Prediction (Steady State and Transient)• Free Surface• Radiation Heat Transfer Between Solids • Heat sources due to Peltier effect• Radiant flux on surfaces of semi-transparent bodies• Joule heating due to direct electric current in electrically conducting solids• Various types of thermal conductivity in solid medium • Cavitation in incompressible water flows• Equilibrium volume condensation of water from steam and its influence on fluid flow and heat transfer• Relative humidity in gases and mixtures of gases • Two-phase (fluid + particles) flows • Periodic boundary conditions.• Tracer Study• Comfort Parameters • Heat Pipes • Thermal Joints• Two-resistor Components • PCBs•Thermoelectric Coolers• Test the heat exchange on AC and DC power converters.• Simulate internal temperature control to reduce overheating issues.• Better position fans and optimize air flux inside a design.• Predict noise generated by your designed system.Some capabilities above need the HVAC or Electronic Cooling Module.SOLIDWORK Design Support• Fully embedded in SOLIDWORKS 3D CAD• Support SOLIDWORKS configurations and materials • Help Documentation • Knowledge base• Engineering database• eDrawings ® of SOLIDWORKS Simulation results General Fluid Flow Analysis• 2D flow • 3D flow • Symmetry• Sector Periodicity • Internal fluid flows • External fluid flowsAnalysis Types• Steady state and transient fluid flows • Liquids • Gases• Non-Newtonian liquids • Mixed flows• Compressible gas and incompressible fluid flows •Subsonic, transonic and supersonic gas flowsMesher• Global Mesh Automatic and Manual settings • Local mesh refinementGeneral Capabilities• Fluid flows and heat transfer in porous media • Flows of non-Newtonian liquids • Flows of compressible liquids •Real gases。
thermal duct model derived by fin-theory approach and applied on an unglazed solar collector
A steady state thermal duct model derived by fin-theory approachand applied on an unglazed solar collectorB.Stojanovic´*,D.Hallberg,J.Akander Building Materials Technology,KTH Research School,Centre for Built Environment,University of Ga ¨vle,SE-80176Ga ¨vle,SwedenReceived 12November 2008;received in revised form 28January 2010;accepted 15June 2010Available online 24August 2010Communicated by:Associate Editor Brian NortonAbstractThis paper presents the thermal modelling of an unglazed solar collector (USC)flat panel,with the aim of producing a detailed yet swift thermal steady-state model.The model is analytical,one-dimensional (1D)and derived by a fin-theory approach.It represents the thermal performance of an arbitrary duct with applied boundary conditions equal to those of a flat panel collector.The derived model is meant to be used for efficient optimisation and design of USC flat panels (or similar applications),as well as detailed thermal analysis of temperature fields and heat transfer distributions/variations at steady-state conditions;without requiring a large amount of computa-tional power and time.Detailed surface temperatures are necessary features for durability studies of the surface coating,hence the effect of coating degradation on USC and system performance.The model accuracy and proficiency has been benchmarked against a detailed three-dimensional Finite Difference Model (3D FDM)and two simpler 1D analytical models.Results from the benchmarking test show that the fin-theory model has excellent capabilities of calculating energy performances and fluid temperature profiles,as well as detailed material temperature fields and heat transfer distributions/variations (at steady-state conditions),while still being suitable for component analysis in junction to system simulations as the model is analytical.The accuracy of the model is high in comparison to the 3D FDM (the prime benchmark),as long as the fin-theory assumption prevails (no ‘or negligible’temperature gradient in the fin perpendicularly to the fin length).Comparison with the other models also shows that when the USC duct material has a high thermal conductivity,the cross-sectional material temperature adopts an isothermal state (for the assessed USC duct geometry),which makes the 1D isothermal model valid.When the USC duct material has a low thermal conductivity,the heat transfer course of events adopts a 1D heat flow that reassembles the conditions of the 1D simple model (for the assessed USC duct geometry);1D heat flow through the top and bottom fins/sheets as the duct wall reassembles a state of adiabatic condition.Ó2010Elsevier Ltd.All rights reserved.Keywords:Unglazed solar collector;Roof integrated;Duct;Modelling;Fin-theory;Benchmarking1.IntroductionHeat transfer simulations and calculations on absorption of heat in flat-plate solar collectors,have generally been based on simplified one-dimensional (1D)heat transfer modelling (e.g.see:Duffie and Beckman,2006;Fischer et al.,2004).Traditionally,a fluid element flowing through the collector (or arbitrary duct)is assumed to absorb heatfrom its ambient in a 1D manner,along its flow path(Duffie and Beckman,2006).This approach renders a model that requires small computational resources and time,which is beneficial in system simulations when sea-sonal or annual performances are analysed.The efficiency and accuracy of the model is adequate for analysing energy performance,calculating bulk or outlet fluid temperatures,or ordinary absorber plate and fluid temperature distribu-tions (Duffie and Beckman,2006).If detailed analysis of temperature fields and heat transfer distributions/variations at steady-state or dynamic conditions are required,more0038-092X/$-see front matter Ó2010Elsevier Ltd.All rights reserved.doi:10.1016/j.solener.2010.06.016*Corresponding author.Tel.:+4626648137;fax:+4626648181.E-mail address:bojan.stojanovic@hig.se (B.Stojanovic´)./locate/solenerSolar Energy 84(2010)1838–1851sophisticated and complex models have to be used.For instance,this is needed when solar collectors have a duct/ tubefluid volume circumference that is significantly larger ()2)than the heat-absorbing surface.In this case regarding the duct/tube as a point in the absorber plate(e.g.Duffie and Beckman,2006;Hilmer et al.,1999)is not appropriate. Or in other cases,the thermal bridging between the absorb-ing surface and the backside,via the duct walls,cannot be neglected(Ammari,2003;Cristofari et al.,2002),nor the heat exchange between a duct wall and the collectorfluid (Ammari,2003;Hachemi,1999;Yeh et al.,2002).Detailed thermal modelling(of e.g.solar collectors)usually result in a multi-dimensional numerical model(e.g.Hassan and Bel-iveau,2007).The model can vary in complexity,depending on the desired level of detail.However,a detailed numerical model requires a substantial amount of computational power and time,which primarily makes it suitable and use-ful for specific and detailed case studies;not suitable as a component in system simulations.Furthermore,multi-dimensional numerical models have afixed geometrical lay-out of nodes/elements/volumes(depending on mathemati-cal discretisation principle),which makes the task of changing and comparing different solar collector designs tedious.In some cases there is a need to have a model which can calculate/simulate detailed thermal distributions and varia-tions,without requiring a large amount of pre-processing or computational power and time;thereby attaining a model suitable for component analysis while being useful in system simulations.The advantages are that solar collec-tor analysis and optimisation can be performed in junction to system operation at different geographical locations, during long-term simulation scenarios.As an example,this procedure is useful in durability assessments of solar collector absorber surfaces(e.g.see Carlsson et al.,1994).Annual simulations of solarNomenclatureA area(m2)bfin thickness(m)C equation constant(unit)C[T f(z)]equation constant as a function offluid temper-ature(unit)c heat capacity(J/K)c p specific heat capacity(J/kg K)d diameter(m)G global solar irradiance on a surface(W/m2)h heat transfer coefficient(W/m2K)K parameter for simplifying thefin temperature equation constants(–)Lfin length(m)q[x i,T f(z)]heatflux in an infinite smallfin element as a function of position in afin andfluid tempera-ture(W/m2)Q heat transfer rate(W)Q[x i,T f(z)]heat transfer rate in an infinite smallfin ele-ment as a function of(W)position in afin andfluid temperatureQ[T f(z)]heat transfer rate as a function offluid temper-ature(W/m)R thermal resistance(m2K/W)t time(s)T temperature(K)T(z)temperature as a function of USC duct length position(K)T(in)inlet temperature(K)T[x i,T f(z)]temperature as a function of position in afin andfluid temperature(K)V ambient air speed(m/s)U total heat loss coefficient(W/m2K)x coordinate parameter(–)y coordinate parameter(–)z coordinate parameter(–)_m massflow(kg/s)Nu Nusselt number(–)Greek symbolsa solar thermal absorptance(–)b parameter for simplifying thefin temperatureequation constants(unit)c[T f(z)]parameter as a function offluid temperature for simplifying thefin temperature equation con-stants(unit)D difference(–)e long wave(IR)emittance(–)j thermal conductance(W/K)k thermal conductivity(W/m K)q density(kg/m3)r Stefan–Boltzmann constant(W/m2K4) Subscriptsa ambientb backsideffluidh hydraulici nodal numberj nodal numberk nodal numberm materialn number of time-stepr radiations surfacefinfinsky skyUSC unglazed solar collectorB.Stojanovic´et al./Solar Energy84(2010)1838–18511839collectors operating in system solutions are performed atdifferent geographical locations,in order to assess the absorber surface microclimatic exposure to relevant degra-dation agents,such as temperatureet al.,1994;Van der Linden et al.,In general,durability issues of to a number of problems:lation or freezing (Wennerholm,issues are the exposure to various temperature,wind,rain/humidity,lutants)as a result of outdoor These contribute to a degradation tion of the collector material performance due to change in Carlsson et al.,1994;Stojanovic´assess this degradation,the tude of degradation agents (e.g.has to be detailed.Degradation nying degradation mechanisms and be strongly affected by for actively heated and cooled only be physicochemical activation values are passed.If only being degradation agents and assessment on degradation will in –time of wetness (TOW)is an degradation experiments exposure testing,e.g.see Carlsson that degradation agents and and that dose levels under a certain ute to a negligible rate of threshold values are exceeded.These changes in optical posed throughout the system,ciency (Hollands et al.,1992),economical and environmental a result of these types of changes,assessment is needed,so that the the performance lowering assessed.Work on assessing the mance due to solar collector tion of microclimate in/on solar of the most important steps (Carlsson et al.,1994;Hollands 2009).1.1.Background:a roof-integrated An unglazed solar collector the EU project Endothermic cient Housing in the EU unched in strated the use of a solar-assisted 2008;Stojanovic´,2009; e.g.see 2007)to provide the annual heating,cooling and hot water houses in different regions of the USC is a dual-purpose component that integrates with the building to form its roof (Virk,2008).The collector propi-tiously blends into its surroundings (different shapes and col-1840 B.Stojanovic´et al./Solar Energy 84(2010)1838–1851prevailing local climatic conditions.These modificationsconsisted in a variation of system:build-up,operation and control strategy.Regarding the USCs,no changes in the design of the component was made for the different systems,except for:coating colour,length of the extruded USC and number of USC panels used in the assembled system and rate of flow of the heat transfer fluid.The USC consists of an extruded aluminium profile which comes in two shapes,flat and bold rolled (see Fig.1).After extrusion,the ends are sealed and in-/outlet pipes are attached by welding.Thereafter,the panels are painted with a polyester powder coating.Each panel has a fixed:width (0.22m),number of ducts and material thickness (see also Table 2),while the length can vary up to a maximum of 6m.Feet and folds enable fixation onto the roof and inter-locking between adjacent panels.The following panel is interlocked to the previous and thereafter fixated.This forms an assembling chain that integrates the USCs into the roof and forms a homogenous surface.The heat transfer liquid flows through the panel in one homogenous direction,where each duct makes up an individual parallel flow streak (see Fig.2).The USC circuit has a flow regime that can vary from serial to parallel,depending on how the panels are con-nected.The fluid is distributed to and from the panels by a manifold.The fluid flow in the USC roof circuit is obtained via one or several circulation pumps.1.2.ObjectivesThis paper presents thermal modelling of the Endohous-ing USC flat panel that uses heat transfer liquids.The aimis to produce a detailed yet swift thermal steady-state model.The model is analytical,1D and derived by fin-the-ory approach.It represents the thermal performance of an arbitrary duct with applied boundary conditions equal to the flat panel USC;hence the model is also applicable for other thermal duct calculations.As the USC predomi-nantly consists of angular ducts,the traditional thermal modelling approach of the duct/tube being a point in the absorber is not appropriate (if detail analysis of tempera-ture fields and heat transfer distributions/variations is required);in this case the ducts constitute the absorber,hence the USC flat panel model development.The derived model is to be used for efficient optimisation and design of the USC flat panels (or similar applications),as well as detailed thermal analysis at steady-state conditions;with-out requiring a large amount of computational power and time.Detailed surface temperatures are necessary fea-tures of the fin-theory based USC thermal duct model in order to be useful for durability studies of the surface coat-ing,hence the effect of coating degradation on USC andsystem performance (see Section 1and reference Stojanovic´et al.,2008).The model accuracy and proficiency is bench-marked against a detailed 3D numerical model as well as two 1D analytical models.The paper objective is also to present results,discussion and conclusions of the different benchmarking models comparison,as these are commonly used for this type of application;in line with the work pre-sented by Schnieders (1997).1.3.Model assumptionsThis section represents a summary of assumptions applied in the derivation of the various USC thermal mod-els.Table 1displays which assumption is applicable to each model.1.The fluid temperature in each duct is equal,as the flow in the USC panel is parallel (see Section 1.1).The USC roof installation should also be seen as having a parallel circuit,as is the case in theSwedish Endosite system (Virk,2008;Stojanovic´,2009).2.The thermal effects of the feet and folds on the USC flat panel are neglected (see Figs.1and 3).Panel feet affect the USC by thermally bridging the duct assem-bly (especially the adjacent ducts)to the substrate by conduction.A roof integrated USC installation con-sists of a number of panels that are interlocked with each other,the last duct in one panel and the firstTable 2The geometrical dimensions of the modelled USC half duct;representing the actual dimensions of the flat panel USC developed and used in the Endohousing project.Note that the wall fin/sheet thickness is not eligible for the 1D simple steady-state C models Half duct inner with (m)Half duct inner height (m)Top and bottomfin/sheet thickness (m)Wall fin/sheet thickness (m)USC duct length (m)For all models0.01250.0130.0020.0015Table 1A summary of assumptions applied in the derivation of the various USC models.Assumption Fin-theory steady state 3D steady state FDM 1D isothermal steady state 1D simple steady state 1x x x x 2x x x x 3x x x x 4x x x x 5x x x x 6x x x 7x x x x 8x xx x 9xx x 10x11xNote :A detailed thermal analysis at the specific USC flat panel boundary conditions of the panel feet and the end duct wall is beyond the scope of this paper or the presented thermal duct models.B.Stojanovic ´et al./Solar Energy 84(2010)1838–18511841in a neighbouring panel will have a different build-up than a centrally placed duct.This is especially appar-ent at the duct wall,which will have a different heat transfer than the wall of a centrally placed duct.3.Boundary effects at the perimeter of the roof are neglected due to marginal impact on overall USC roof performance.4.The duct fluid temperature is isothermal (no temper-ature gradient,fully mixed)in each cross-section,except at the thermal boundary layer which gives rise to the convective heat transfer coefficient (Holman,2002).5.The average Nusselt number for(cross-section aspect ratio 1/2)wall heat flux and fully applied for the duct convection (1978).6.A linearisation of the equation energy (heat)transfer rate by tion (Duffie and Beckman,2006the discussion in Section 2.1.7.No heat exchange by long wave within the duct since the heat (see Section 1.1).8.Material,fluid and heat transfer exception of long wave stants,independent of 9.The model assumes that there is pendicular to the USC duct 10.The USC duct consists of a material cross-section.11.The USC duct consists of a temperature difference over and and the lower USC duct wall sheet/fin is adiabatic.The sheets/fins have each a their respective fin,which is length.2.Fin-theory based USC thermal model derivation A 1D steady state thermal model based on fin-theory isderived in order to obtain a tool that swiftly calculates and assesses the USC’s performance:material temperature and heat transfer distribution/variation and fluid temperature.This section presents the model derivation procedure.The basic concept of the 1D steady state fin-theory is that heat is transported from the fin-base throughout the fin-material in length direction;there is no ‘or negligible’temperature gradient in the fin perpendicularly to the fin length.While conduction prevails,heat is either emitted or absorbed by the fin depending on the interaction with the ambient (Holman,2002).The USC flat panel cross-section represented in Fig.3shows that a thermal symmetry can be applied on a USC flat panel duct,dividing it into a half duct.The USC half duct becomes a system of three fin sheets.Fig.4presents the thermal system with applied boundary conditions and defined energy flows used in the fin-theory based model.These boundary conditions represent the thermal course of event in a central duct on flat panel USC.As previously discussed,by adopting that statements 1,2and 3repre-sented in Section 1.3(see also Table 1),the heat transfer occurring in a centrally placed duct can be seen as repre-senting an entire USC panel or roof (see also Fig.3).Addi-tionally,please also observe the notification at the end of Section 1.3.Analytical thermal modelling of ducts in similar applica-tions has also been presented by Ammari (2003),Hachemi (1999)and Yeh et al.(2002).Ammari (2003)presents a Fig.3.A cross-section of the Endohousing flat panel USC.The figure shows the applied thermal symmetry lines on a centrally placed duct that makes a half duct,which is seen as representing the USC thermal behaviour.The models neglect the thermal effects of the two lower feet and the folds at each edge.Energy 84(2010)1838–1851air-heaters.In their models,the duct walls are seen as fac-tors that increase the thermal convection from duct sur-faces to the airflow.No thermal bridging by thermal conduction via the duct walls is applied.2.1.Mathematical derivationThe following example presents the derivation of the temperature distribution expression of Fin1(see Fig.5), which is the topfin/sheet of the USC half duct.The heat balance of over an infinite small element in a cross-section of Fin1is according to Figs.4and5as:Q1½x1;T fðzÞ þQ3½x1;T fðzÞ þQ4þQ5½x1;T fðzÞ¼Q6½x1;T fðzÞ þQ2½x1þdx1;T fðzÞ ð1ÞThe rate of heat conduction to the element at x1isQ1½x1;T fðzÞ ¼Àk1ÁA1ÁdT1½x1;T fðzÞdx1ð2Þand heat conduction from the element at position x1+dx1isQ2½x1þdx1;T fðzÞ ¼Àk1ÁA1ÁdT1½x1;T fðzÞdx1þddx1Àk1ÁA1ÁdT1½x1;T fðzÞdx1dx1ð3ÞHeat transfer rate by ambient(outdoor)convection:Q3½x1;T fðzÞ ¼h aÁA1:2ÁðT aÀT1½x1;T fðzÞ Þð4ÞThe outdoor wind convection coefficient h a used through-out this paper is according to Watmuffet al.(1977)as:h a¼2:8þ3:0ÁVð5ÞEnergy(heat)transfer rate by solar radiation:Q4¼aÁA1:2ÁG USCð6ÞEnergy(heat)transfer rate by long wave(IR)radiation (Duffie and Beckman,2006):ture T1[x1,T f(z)]is a function of bothfin length/position(x1) and USC ductfluid temperature along the duct length/posi-tion(z),a representative mean value of T1[x1,T f(z)]derived in an iterate manner is used in order to calculate h r.Heat transfer rate byfluid convection:Q6½x1;T fðzÞ ¼h fÁA1:2ÁðT1½x1;T fðzÞ ÀT fðzÞÞð8ÞThe ductfluid convection coefficient h f used throughout this paper is according to Holman(2002)as:h f¼NuÁk fd hð9Þwhere the average Nusselt number for a rectangular duct(cross-section aspect ratio1/2)with constant axial wall heatflux and fully developed laminarflow is Nu= 4.123(Shah and London,1978).The combination of the above-presented equations in accordance to Eq.(1)will give:d2T1½x1;T fðzÞdx1Àðh fþh aþh rÞk1Áb1|fflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflffl}b1ÁT1½x1;T fðzÞþh fÁT fðzÞþh aÁT aþaÁG USCþh rÁT skyk1Áb1|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}c1½T fðzÞ¼0ð10ÞThe partial differential equation of the Fin1tempera-ture distribution becomes:d2T1½x1;T fðzÞdx1Àb1ÁT1½x1;T fðzÞ þc1½T fðzÞ ¼0ð11ÞThe solution to the partial differential equation in Eq.(11)is:T1½x1;T fðzÞ ¼C11½T fðzÞ ÁeÀffiffiffiffib1pÁx1þC12½T fðzÞ Áeffiffiffiffib1pÁx1þc1½T fðzÞb1ð12ÞThe procedure for deriving thefin-temperature distribu-tion for Fins2and3is in direct analogy with the above presented example.Fins2and3are applied with the boundary conditions and heat balances that are displayed in Fig.4and presented in Appendix A.In order to be able to calculate the cross-sectional tem-perature distribution in Fins1,2and3(Eqs.(12),(A3), and(A5)),the unknown constants in thefin-temperature equations have to be solved.By using boundary conditions as illustrated in Fig4and combining the heat transfer rates exiting from onefin and entering another,provides the pos-sibility to solve the unknown constants.Fin1dT1½x1¼0;T fðzÞ1¼0;Àk1ÁA1ÁdT1½x1¼L1;T fðzÞdx1¼Àk2ÁA2ÁdT2½x2¼0;T fðzÞdx2B.Stojanovic´et al./Solar Energy84(2010)1838–18511843Fin2T2½x2¼0;T fðzÞ ¼T1½x1¼L1;T fðzÞ T2½x2¼L2;T fðzÞ ¼T3½x3¼L3;T fðzÞFin3dT3½x3¼0;T fðzÞdx3¼0;Àk3ÁA3ÁdT3½x3¼L3;T fðzÞdx3¼k2ÁA2ÁdT2½x2¼L2;T fðzÞdx2The solved unknown constants of thefin-temperature distribution equation in a USC duct cross-section are defined and presented in Appendix A.By integrating the heatflux(containing thefin-tempera-ture distribution)over the entirefin length attains the heat transfer rate from afin cross-section.Eq.(13)presents the heat transfer rate from Fin1to the USCfluid(see also Fig.4).An example of the derivation of the heat transfer rate equations to and from the USC half duct is presented in Appendix A.Q6½T fðzÞ ¼Àh fÁÀC11½T fðzÞ Áb1þC12½T fðzÞ Áb1þC11½T fðzÞ ÁeðÀffiffiffiffib1pÁL1ÞÁb1ÀC12½T fðzÞ Áeðffiffiffiffib1pÁL1ÞÁb1Àc1½T fðzÞ ÁL1Áffiffiffiffiffib1pþT fðzÞÁL1Ábð3=2Þ1 0@1Ab1ð13ÞHaving derived heat transfer rates from thefin cross-sec-tions to the USCfluid,provides the possibility to calculate thefluid temperature along the USC duct,hence thefin-temperature distribution.Eq.(14)presents the steady state fluid temperature along the USC duct.The equation assumes that there is no heat transfer perpendicular to the USC duct cross-section.T fðzÞ¼T fðinÞþC2 C1Áe zÁC1_mÁc p fÀC2C1ð14ÞBy attainment of Eq.(14),thefin-temperature distribu-tion at a USC duct cross-section,as well as along the USC duct,can be calculated at steady-state conditions.The der-ivation of Eq.(14)and definition of equation constants are presented in Appendix A.3.Benchmarking modelsThe prime objective of the benchmarking models is to be used as a comparison for thefin-theory derived USC ther-mal duct model,in order to investigate calculated/simulated results deviations,but also as a comparison amongst themselves.3.1.1D isothermal steady-state modelThis model is based on the assumption that the thermal conductivity of the duct material is infinitely large,thus making the USC duct material cross-section isothermal. Other assumptions applied are listed in Table1;please also see the discussion in the beginning of Section2.The model is basically‘equivalent’to the Hottel–Whiller–Bliss model as presented by Duffie and Beckman(2006).As the USC duct(in this case)constitutes the absorber plate,there is no absorber plate surface that is without direct contact with the ductfluid.By applying this case to the Hottel–Whiller–Bliss model(as presented by Duffie and Beckman (2006))and neglecting the bond thermal resistance between the absorber plate and tube/duct,brings that the duct material(fin)temperature is isothermal;as the model regards the duct as a point in the absorber plate cross-sec-tion.Fig.6illustrates the model set-up.The heat balance of afluid element,having the massflow_m in z-direction and temperature T f(z),gives that the increase influid tempera-ture dT f(z)due to heat supply corresponds to:h fÁðL1þL2þL3ÞÁdzÁðT s½T fðzÞ ÀT fðzÞÞ|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}Q6½T fðzÞ þQ9½T fðzÞ þQ12½T fðzÞ¼h aÁL1ÁdzÁðT aÀT s½T fðzÞ Þ|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}Q3½T fðzÞþaÁG USCÁL1Ádz|fflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflffl}Q4þh rÁL1ÁdzÁðT skyÀT s½T fðzÞ Þ|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}Q5½T fðzÞþU bÁL3ÁdzÁðT bÀT s½T fðzÞ Þ|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}Q13½T fðzÞð15ÞThe equation above can be reformulated with constants C1,C2and C3as expressed below.C1¼aÁL1ÁG USGþh aÁL1ÁT aþh rÁL1ÁT skyþU bÁL3ÁT bh aÁL1þh rÁL1þU bÁL3þh fÁðL1þL2þL3Þð16ÞC2¼Àðh aÁL1þh rÁL1þU bÁL3Þðh aÁL1þh rÁL1þU bÁL3þh fÁðL1þL2þL3ÞÞð17ÞC3¼h fÁðL1þL2þL3Þ_mÁc pfð18ÞLocal cross-sectional material/surface isotherm can accordingly be found such that:T s½T fðzÞ ¼aÁG USCÁL1þh aÁL1ÁT aþh rÁL1ÁT skyþU bÁL3ÁT bþh fÁðL1þL2þL3ÞÁT fðzÞf123a1r1b3ð19ÞConsidering that the inlet temperature T f(in)will increase to be the outlet temperature T f(z)at the duct length z,the heat balance of thefluid element can be inte-grated such that:Z T fðzÞT fðinÞdT fðzÞC1þC2ÁT fðzÞ¼C3ÁZ zdzð20ÞThe solution for Eq.(20)isfinallyT fðzÞ¼T fðinÞþC1C2Áe C2ÁC3ÁzÀC1C2ð21Þ1844 B.Stojanovic´et al./Solar Energy84(2010)1838–1851。
FEA simulation
Technical ReportFinite element simulation of the temperature and stress fields in single layers built without-support in selective lasermeltingAhmed Hussein,Liang Hao ⇑,Chunze Yan,Richard EversonCollege of Engineering,Mathematics and Physical Sciences,University of Exeter,Exeter,EX44QF Devon,United Kingdoma r t i c l e i n f o Article history:Received 5February 2013Accepted 23May 2013Available online 7June 2013a b s t r a c tOverhanging and floating layers which are introduced during the build in selective laser melting (SLM)process are usually associated with high temperature gradients and thermal stresses.As there is no underlying solid material,less heat is dissipated to the powder bed and the melted layer is free to deform resulting undesired effects such as shrinkage and crack.This study uses three-dimensional finite element simulation to investigate the temperature and stress fields in single 316L stainless steel layers built on the powder bed without support in SLM.A non-linear transient model based on sequentially coupled thermo-mechanical field analysis code was developed in ANSYS parametric design language (APDL).It is found that the predicted length of the melt pool increases at higher scan speed while both width and depth of the melt pool decreases.The cyclic melting and cooling rates in the scanned tracks result high VonMises stresses in the consolidated tracks of the layer.Ó2013Elsevier Ltd.All rights reserved.1.IntroductionOne of the limitations of SLM is the thermal distortion of the part during forming,which may lead to undesired shrinkage and cracks.The elevated temperature gradients developed yield a non-homogeneous permanent strains and residual stresses distrib-uted within consolidated layers [1].As a result,the dimensional size,shape,and mechanical properties of the fabricated parts are affected considerably.The processing parameters such as laser power,scanning speed,laser spot size,and scanning strategy all play a crucial role on the development of temperature gradients and residual stresses in the part [2].These parameters are usually optimised through experimental means for specific machines and materials.However,detailed investigation of all different parame-ters and materials for SLM through experiments can be time con-suming and costly.Therefore,numerical methods are mainly used as a tool to study the role of these parameters on temperature distribution,residual stress and other thermal mechanisms.The finite element method is the most commonly used numer-ical method for predicting temperature and stress fields in SLM.Childs et al.[3]investigated the influence of process parameters on the mass of melted single layers in SLM and found that melted mass increased with increasing scanning speed.The effect of layer thickness on part deformation in SLM parts was investigated by Zaeh and Branner [4].A thinner layer thickness resulted in a higher deformation due to the effect of higher temperature variation concentrated within thin layers.Matsumoto et al.[5]developed a finite element analysis (FEA)method for single layer parts on the loose powder in SLM.A 2D non-linear heat transfer with volume internal heat source problem is numerically solved based on the coupling of Matlab and ANSYS FEM models [6].The surface tem-perature distribution during SLM of 90W–7Ni–3F materials was predicted by Zhang et al.[7].Dai and Shaw [8]use the finite ele-ment analysis to investigate the effect of laser scanning strategy on residual thermal stresses and distortion.Nickel et al.[9]developed FEA model to simulate the laser deposition process and found that the deposition pattern has a significant effect on the part stresses and deflections.Yin et al.[10],Ibraheem et al.[11],Shuai et al.[12]report their research results of temperature fields in single metallic layer SLM processes by using element birth and death technique.In [13],a more comprehensive understanding of the SLM thermal field has been achieved by creating a 3D model and considering the interval time for new powder recoating.Ma and Bin [14]proposed a 3D FEA model with fixed tempera-ture heat source for calculating the evolution of temperature and thermal stresses within a single metallic layer formed on the pow-der bed using two different scanning patterns in SLS.It was found that the distortion and transient stresses of a layer processed by a moving laser beam decreased with fractal scanning pattern.From this review,it is evident that using 2D analysis with generalised plane strain conditions seems to be convenient with less computer processing requirements,but 3D analysis remains absolutely nec-essary to fully understand the problem.An in-depth understanding of the materials and laser interaction in overhanging regions and its associated thermal and stress mechanisms is essential step for the design of efficient and reliable support structure and improving surface quality of overhanging geometries.Corresponding author.Tel.:+44(0)1392723665.E-mail addresses:Ayh203@ (A.Hussein),l.hao@ (L.Hao).In this study,a3D non-linear transientfinite element model based on sequentially coupled thermo-mechanicalfield analysis was developed in ANSYS programme to predict the temperature distribution,thermal stresses and melt pool dimensions of laser scanned single layers built on the powder bed.Simulation of the moving heat source and changing boundary conditions are accom-plished by a user written subroutine implemented in ANSYS para-metric design language(APDL).Since the heat energy is transported well below the surface of the powder bed in SLM, the laser energy density was applied as a volumetric heat source rather than a surface heatflux.Temperature dependent physical properties of316L stainless steel powder material are taken into account and latent heat of fusion is considered.2.3Dfinite element modelANSYS provides the capability of performing indirect sequen-tially coupled thermo-mechanical analysis for both heat and stress analysis.In this study,a multiple physics environment was used in a single database as this allows a quick switching between physics environments for subsequent stress analysis.A non-linear tran-sient thermal analysis was performedfirst to obtain the global temperature history generated during the laser melting.A tran-sient stress analysis is then developed with an automatic exchange of the element type from thermal-to-structural,and applying thethe heat source from laser over few elements at the time in the la-ser path.The time period for which the laser beam is retained on each step,is given as,T step¼d xVð1Þwhere d x is the length of elements under the laser spot in mm and V is the scanning speed of the laser beam in mm/s.Table1summarises the parameters used in thefinite element simulation.2.1.Thermal modellingThe thermal equilibrium equation satisfies the following classi-cal3D heat conduction equation given by[15],q c d Td t¼dd xkd Td xþdd ykd Td yþdd zkd Td zþQð2Þwhere q is the material density(kg/m3);c is the specific heat capac-ity(J/kg K);T is the temperature;t is the interaction time;k is ther-mal conductivity(W/mK);and Q=(x,y,z,t)is the volumetric heat generation(W/m3).The effective thermal conductivity is a function of porosity of the powder[16].The porosity of the powder/can be calculated as, /¼qbulkÀq powderqð3ÞA.Hussein et al./Materials and Design52(2013)638–647639H¼Zq cðTÞd Tð5ÞFig.2shows the calculated temperature dependant thermal conductivity and enthalpy of the material up to the melting tem-perature.The sudden change of the thermal conductivity from low value calculated from Eq.(4)to higher value in the curve is due to transition of material from powder state to liquid state,in which case the corresponding bulk material properties is used[18].2.1.1.Boundary conditionsThe initial condition of uniform temperature distribution throughout the powder bed prior to laser melting at time t=0 can be applied as,Tðx;y;z;0Þ¼T0ðx;y;zÞð6ÞT0is the ambient temperature taken as293K(20°C).Laser axis direction at z=0,Àk@T@zz¼0¼QÀh T0ÀT surfðÞð7Þwhere h is the heat transfer coefficient at the powder surface which is taken as(10W/mK);and T surf is the temperature of the powder bed surface.Since the layers are built on powder bed with large thickness, the heat transfer at the bottom of loose powder can be assumednegligible. @Twhere r is the radial distance from the beam centre;I0is the inten-sity of the beam at r=0;and x is the radius of the beam at which I¼I0eÀ2.This can be written as,IðrÞ¼2APp x2expÀ2r2x2ð10ÞA is the absorptivity of the powder material which can be calcu-lated if the reflectivity of the material k is known.A¼1Àk(a reflectivity of iron=0.7was considered for316L stainless steel).2.3.Mechanical analysisThe same FE mesh used in thermal analysis is employed here, except for the element type and the boundary conditions.To calcu-late the distribution of stress,the elastic FEA simulation is used. Stress is related to strain by[19],f r g¼½D f e e gð11Þwhere{r}is the stress vector,[D]is the elasticity matrix and,f e e g¼f e gÂf e th gð12Þwhere{e}the total is strain vector and{e th}is the thermal strain vector.Eq.(11)may be written as,f e g¼½D Â1f r gÂf e th gð13ÞFor isotropic material,the above stress–strain relationship can be written in Cartesian co-ordinates as follows[19,20],e xx¼1E½r xxÂmðr yyþr zzÞ þa e D Te yy¼1E½r yyÂmðr xxþr zzÞ þa e D Te zz¼1E½r zzÂmðr xxþr yyÞ þa e D Te xy¼1þmr xy;e xz¼1þme xz;e yz¼1þme yzð14ÞTable1Finite element simulation parameters.Parameter ValueLaser power,P100WScanning speed,V100,200,300mm/sTrack length,L10mmNumber of tracks scanned,N5tracksPower bed thickness,T1mmHatch spacing,H s75l mLaser spot size,D150l mThermal element type,3D SOLID70Structural element type,3D SOLID185640 A.Hussein et al./Materials and Design52(2013)638–647where T ref is the reference temperature at t=0.a is a function of temperature and can be written as,e th¼Z TT refaðTÞd Tð16ÞThe effective stress can be given as,r eff ¼ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffir21þr22þr23qþ2mðr2r2þr1r3þr2r2Þð17Þwhere r1,r2,and r3are the three principal stresses.The equivalent stress or VonMises can be computed as,r eq m ¼ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi1½ðr2Âr2Þ2þðr2Âr3Þ2þðr3Âr1Þ2rð18ÞThe equivalent stress is related to the equivalent strain by this relationship,r0¼E e0ð19ÞFor stress analysis,in addition to density q,the following ther-mo-structural material properties depending on temperature are required are thermal expansion coefficient a,elastic modulus E, Poisson’s ratio m,and yield strength r y,as listed in Table2.3.Results and discussion3.1.Temperature distributionIn SLM,the temperature distribution in the powder bed and consolidated layers changes rapidly with time and space.Fig.3a shows the temperature at the beginning of the laser scanning,from which,the very high temperature gradients in the vicinity of the la-ser spot on the powder bed can be clearly seen due to an applied Gaussian heat source.The temperature of the powder particles is elevated rapidly under the action of absorbed energy,causing a molten pool when the temperature exceeds the melting tempera-ture and heat affected zones in the surrounding loose powder.Note that the energy intensity of the source might also be high enough to cause the material to evaporate[16].The highest predicted tem-perature corresponding to the molten zone of the powder material is2600K for[P=100W,V=100mm/s]and exceeds the melting temperature of316L stainless steel(1672K).However,this maxi-mum temperature at the start of track1is reduced at the end of first track scan to2392K and so at the end offifth track to 2225K as shown in Fig.3b and c respectively.The drop of maxi-mum temperature can be attributed to the increased conductivity of the previously solidified regions of the track compared to the low thermal conductivity available initially in the powder bed. The thermalfield changes as the laser source moves along the track and the melt pool moves along with the laser source.It is further observed that the temperature gradient in the front side of the moving laser beam is much steeper than that in the rear side. The melt pool shape resembles as in comet tail profile(see Fig.3c).This trend of skewed temperature distribution towards the rear of the laser was also reported in other temperature simu-lations[13,22].This can be attributed to the fact that the rapidly cooling molten material has greater conductive properties than the untreated powder in front of the laser.The temperature distribution in the layer is very much affected by the energy density which is a function of laser power,spot size, scanning speed,hatch spacing and scanning strategy.Additionally, the temperature gradient in the layer is similarly influenced by conductivity of the material underneath the deposited layer being a loose powder bed,previously re-melted layer,support structure or a solid substrate.Fig.4a and b shows temporal variation of tem-perature of single track deposited on powder bed and solid sub-strate respectively.For the powder bed,the melted layer thickness depends on the energy input from the laser,while for substrate a thin layer of75l m is deposited on solid material and then melted with the laser beam.A fast cooling rate can be ob-served when the track is built on a solid substrate.Because the so-lid bulk material has higher conductivity compared to powder bed allowing more heat sink through conduction in a shorter period of time and subsequent solidification to room temperature.Fig.5a further shows that the temperature of the scannedfirst track for powder bed and solid substrate decreases as the distance from the laser position is increased.The region directly under the laser beam is in the molten state,whereas the remaining portion of the track length experiences a rapid fall in temperature and undergoes solidification.The solidification begins when the tem-perature drops below the melting temperature of the material and approaches to room temperature.The temperature gradient is higher at the start of solidification and approaches minimum va-lue at the end of solidification.The rapid cooling of tracks built on solid substrate result a more stable molten zone as well as heat af-fected zone at the vicinity of laser spot.Note that multiple reflec-tions effect between the powder particles leads to higher optical penetration depths compared to bulk materials.The various heating and cooling cycles taking place in various tracks which form the layer during laser scan are shown in Fig.5b and c.The peak temperatures represent the start of each track forming the layer and are taken from neighbouring nodes in the Y-axis separated by a hatch spacing of75l m between the tracks.For track1the recorded temperature is taken from a node located at(X=1.225,Y=1.225,Z=0),track2at(X=1.225, Y=1.3,Z=0),track3at(X=1.225,Y=1.375,Z=0),track4at (X=1.225,Y=1.45,Z=0),and track5at(X=1.225,Y=1.525, Z=0).Due to overlap between subsequent tracks which depends on the hatch spacing distance,some sections are exposed to multiple scans and melted twice to bond adjacent tracks and form a fully dense layer.The temperature in half portion of the1st track is raised above its melting temperature while scanning the2nd track.How-ever,as more tracks are deposited in3rd and other subsequent tracks,the temperature falls to below the melting temperature of the material.The lower secondary peaks correspond to the laser scanning the same position on successive tracks.This is also exhib-ited in multi-layer build when upper layers are deposited and sug-gests that the addition of layers and subsequent scanning are of significance to the temperaturefield in the model[13].The rate of overlap between track and layers influence the surface roughness, porosity,and mechanical properties of SLM fabricated parts[23].The cyclic melting/heating and cooling continues in all tracks forming the layer.The heating and subsequent cooling to the ambi-ent(chamber)temperature occur within a few tenths of a millisec-ond of each other,thus suggesting that the irradiated spots areTable2Temperature-dependant mechanical properties of316L stainless steel[21].Temperature(K)29847367387397310731573 Thermal expansion,a(1°CÀ1Â10À6)14.315.616.917.718.519.119.5 Elastic modulus,E(GPa)198.518717215714110610A.Hussein et al./Materials and Design52(2013)638–647641subjected to rapid thermal cycles.These rapid cycles are associated with commensurate thermal stress changes[8].Long track length of the laser beam highly influences the thermal gradients which occur within the layer,while shorter tracks tend to produce more homogeneous temperature gradients but with increased melting depth.As the vector length increases,the higher delay period be-tween successive irradiation leads to a decrease in the amount of energy stored on the surface.3.2.Melt pool dimensions surement was achieved through temperature distribution plots re-corded at the instant when the laser had travelled along its path of scanning.The length represents the distance of the molten mate-rial along the X-axis and parallel to the scan direction while width is taken as the molten region along the Y-axis and perpendicular the scanning direction of the laser.The melted depth is measured from the powder surface to the molten depth inside the powder bed along the Z-axis.Fig.6b and c shows the width and depth pro-files of the melt pool for scanning speeds of100mm/s,200mm/s, and300mm/s with afixed laser power of100W.With liquidus line imposed onto the plots,which represent regions above thetemperature distribution during layer melting(a)at the beginning offirst track scan,(b)at the end of thefirst track scan at time=0.091 scan.and subsequent tracks where peaks were found in areas at the bot-tom surface of the scanned layer.This can be contributed to the higher cooling rates and solidification at the bottom side of the layer through conduction in the powder compared to the upper surface resulting higher residual stresses at the bottom surface.It is expected that the stress in the molten zone is tensile and it will transform to compressive stress as the distance from molten zone increases.According to the equilibrium of force and momentum of the part,the irradiated zone will become surrounded by a zone of compressive stress.Other studies are on agreement that tensile stresses are developed at vicinity of the laser spot while compres-sive stresses are more towards the edges[5].To further understand the stress variation and history in X,Y, and Z directions,two different nodes(node1located at the start of the1st track,and node2located at the start of the5th track of the layer were investigated(see Fig.8).As indicated in Fig.8a, the X-directional stress of node1is cyclic compressive-tensile-compressive stresses;tensile stresses are developed during the la-ser melting and turned into compressive during cooling cycle (solidification).However,as the heating source moves away,the parts cools down and the node remains compressive throughout. The X-directional stress of node2is tensile most of the time and occasionally compressive.Similar trend were also found in Y and Z directional stresses of the nodes as presented in Fig.8b and c. The highest stress was found in X-component along the direction of laser scan of node1and reached(612MPa),whereas the Y and Z components of the stress reach(172MPa)and(162MPa) respectively.The high stress value in node1arises from the fact that it expe-riences an early heating and cooling cycles and thus causes higher temperature gradients and related thermal stresses.The tempera-ture gradients are higher at the beginning of the layer during laser processing due to the lower temperature of the surrounding pow-der,and therefore have large effect on the distribution and inten-sity of residual stresses.Thefirst melted layer shrinks during cooling cycle,but as more layers are deposited,the already solidi-fied layers constrains and prevents further shrinking in the top layer.Since this mechanism occurs for each layer at each step of the SLM process,residual stresses may develop inside the manu-factured component.For this reason,SLM parts are stress relieved through heat treated process prior to otherfinishing operations such as the removal of the part from the solid substrate or support-ing structures.Otherwise,the residual stresses will cause part deflection due to stress relaxation.For layers freely built on the powder bed,there is no underlying solid material to constrain the thermal stresses observed in the FEA model.The material expansion(i.e.tensile stresses)during heating cycle and material contractions during cooling(pressive stresses)cause shrinkages and cracks to the layer.After the laser beam leaves that area,the irradiated zone will cool and tends to shrink.The shrinkage is partially inhibited as a consequence of the plastic deformation developed during heating,yielding a resid-ual tensile stress condition at the irradiated zone.The high tensile stresses generated in the x-direction(along the laser scanning direction)may also lead to transverse cracking of the layer.Crack-ing can be avoided by preheating or by using shorter scanning tracks which reduce the cooling rate.The mechanism that prevents cracking by preheating increases the ductility of the material and enhances the possibility of stress relief by plastic deformation. 4.ConclusionThree dimensional transientfinite element model was devel-process.The model results can be summarised in the following points.Highest temperature gradient was recorded at the start offirst track scan and drops subsequently for all scanning speeds.High cooling rates were predicted when the layer is scanned over a solid substrate compared to when layer is scanned on loose powder bed.The predicted length of the melt pool increases at higher scan speeds while both width and depth of the melt pool decreases. High VonMises stresses was predicted within the layer caused by the stepwise increase in the temperature with each succes-sive overlapping laser tracks and alternating compressive and tensile residual stresses within each track.Cracks usually found in overhang layers not only initiate in areas of compressive stress,but if tensile stress is present which occur if compressive yield stress for the material is exceeded.Simulation results reveals the magnitude of the temperature and stress in different locations of the layer.It provide us insight into regions requiring more support structures so that the design aspects of reliable structures can be implemented to withstand residual stress forces and dissipate heat during part building.AcknowledgementThe authors acknowledge funding by the UK Technology Strat-egy Board(TSB)under research project grant(BA036D). References[1]Kruth J,Froyen L,Van Vaerenbergh J,Mercelis P,Rombouts M,Lauwers B.Selective laser melting of iron-based powder.J Mater Process Technol 2004;149(1):616–22.[2]Mercelis P,Kruth J.Residual stresses in selective laser sintering and selectivelaser melting.Rapid Prototyp J2006;12(5):254–65.[3]Childs T,Hauser C,Badrossamay M.Selective laser sintering(melting)ofstainless and tool steel powders:experiments and modelling.Proc Inst Mech Eng,Part B:J Eng Manuf2005;219(4):339–57.[4]Zaeh M,Branner G.Investigations on residual stresses and deformations inselective laser melting.Prod Eng2010;4(1):35–45.[5]Matsumoto M,Shiomi M,Osakada K,Abe F.Finite element analysis of singlelayer forming on metallic powder bed in rapid prototyping by selective laser processing.Int J Mach Tools Manuf2004;42(1):61–7.[6]Patil RB,Yadava V.Finite element analysis of temperature distribution in singlemetallic powder layer during metal laser sintering.Int J Mach Tools Manuf 2007;47(7–8):1069–80.[7]Zhang DQ,Cai QZ,et al.Select laser melting of W–Ni–Fe powders:simulationand experimental study.Int J Adv Manuf Technol2010;51(5–8):649–58. [8]Dai K,Shaw L.Distortion minimization of laser-processed components throughcontrol of laser scanning patterns.Rapid Prototyp J2002;8(5):270–6.[9]Nickel A,Barnett D,Prinz F.Thermal stresses and deposition patterns inlayered manufacturing.Mater Sci Eng2001;A317(1):59–64.[10]Yin J,Zhu H,et al.Simulation of temperature distribution in single metallicpowder layer for laser put Mater Sci2012;53(1):333–9.[11]Ibraheem Ameer K,Derby Brian,et al.Thermal and residual stress modelling ofthe selective laser sintering process.Mater Res Soc2003;758:47–52.[12]Shuai C,Feng P,et al.Simulation of dynamic temperaturefield during selectivelaser sintering of ceramic powder.Mathe Computer Modell Dynam Syst 2012;1(11).[13]Roberts IA,Wang CJ,et al.A three-dimensionalfinite element analysis of thetemperaturefield during laser melting of metal powders in additive layer manufacturing.Int J Mach Tools Manuf2009;49(12):916–23.[14]Ma L,Bin H.Temperature and stress analysis and simulation in fractalscanning-based laser sintering.Int J Adv Manuf Technol2007;34(9):898–903.[15]Carslaw H,Jaeger J.Conduction of heat in solids,2nd ed.vol.1;1959.[16]Dong L,Makradi A,Ahzi S,Remond Y.Three-dimensional transientfiniteelement analysis of the selective laser sintering process.J Mater Process Technol2009;29:700–6.[17]Thummler F,Oberacker R.An introduction to powder metallurgy.Cambridge(London):The University Press;1993.[18]Mills,Kenneth C.Recommended values of thermophysical properties forselected commercial alloy.Woodhead Publishing;2002.646 A.Hussein et al./Materials and Design52(2013)638–647[20]Yilbas BS,Arif AFM.Material response to thermal loading due to short pulselaser heating.Int J Heat Mass Transfer2001;44:3787–98.[21]Deng D,Kiyoshima S.Numerical simulation of residual stresses induced bylaser beam welding in a SUS316stainless steel pipe with considering initial residual stress influences.Nucl Eng Des2010;240:688–96.[22]Shen X-F et al.Finite element simulation of thermal distribution in directmetal laser multi-track sintering.J Sichuan Univ(Eng Sci Ed) 2005;37(1):47–51.[23]Guan K,Wang Z,Gao M,Li X,Zeng X.Effects of processing parameters ontensile properties of selective laser melted304stainless steel.Mater Des 2013;50:581–6.[24]Gu D,Shen Y.Balling phenomena in direct laser sintering of stainless steelpowder:metallurgical mechanisms and control methods.Mater Des 2009;30(8):2903–10.[25]Simchi A,Pohl H.Direct laser sintering of iron–graphite powder mixture.Mater Sci Eng A2004;383:191–200.[26]Brandl E,Heckenberger U,Holzinger V,Buchbinder D.Additive manufacturedAlSi10Mg samples using Selective Laser Melting(SLM):microstructure,high cycle fatigue,and fracture behavior.Mater Des2012;34:159–69.[27]Agarwala M,Bourell D,Beaman J,Marcus H,Barlow J.Direct selective lasersintering of metals.Rapid Prototyp J1995;1:26–36.[28]Gu DD,Shen YF,Xiao J.Influence of processing parameters on particulatedispersion in direct laser sintered WC–Cop/Cu MMCs.Int J Refract Met Hard Mater2008;26:411–22.[29]Boccalini M,Goldenstein H.Solidification of high speed steels.Int Mater Rev2001;46:92–115.[30]Simchi A,Asgharzadeh H.Densification and microstructural evaluation duringlaser sintering of M2high speed steel powder.Mater Sci Technol 2004;20:1462–8.A.Hussein et al./Materials and Design52(2013)638–647647。
55simnmo钢热变形行为实验研究及数值模拟
55SiMnMo钢热变形行为实验研究及数值模拟作者：于恩林， 赵玉倩， 闫涛， 许小林作者单位：于恩林,闫涛,许小林(燕山大学机械工程学院)， 赵玉倩(燕山大学机械工程学院 东北大学秦皇岛分校)相似文献(2条)1.期刊论文杜林秀.姚圣杰.熊明鲜.王国栋.DU Lin-xiu.YAO Sheng-jie.XIONG Ming-xian.WANG Guo-dong低碳结构钢的奥氏体晶粒超细化-东北大学学报（自然科学版）2007,28(11)以两种不同成分的低碳微合金结构钢为研究对象,结合热模拟实验与实验室热轧实验,研究原始组织、化学成分及部分加热条件对低碳钢加热过程奥氏体晶粒超细化的影响规律.结果表明,以热轧态铁素体/珠光体经温轧并且冷变形的组织为原始组织最有利于获得超细晶奥氏体(1 μm);此外适量添加合金元素Nb,Ti,V,适当提高加热速度均有利于细化奥氏体,而当加热速度大于100 ℃/s时,对奥氏体的超细化效果不明显;另外,加热前预变形可以显著细化奥氏体晶粒,且提高其尺寸均匀性.2.学位论文张永军现代汽车用钢生产技术开发2009汽车工业是一个国家的重要支柱产业，其发展水平是衡量一个国家工业发展水平的重要标志。
在所有的钢铁产品中，汽车用钢是技术含量最高、附加值最大的品种之一。
在工业发达国家，汽车用钢占钢材总产量的17％左右。
因此，对汽车用钢的研究具有十分重要的意义。
钢铁材料是汽车制造中的基本材料，尤其是结构钢或合金结构钢，不仅用量大，而且主要用于关键零件，如发动机构件、齿轮、各类轴件、紧固件等，是保证汽车运行性能的核心部件的制造材料。
<br> 近年来，邢台钢铁有限责任公司实施了“做精、做专、做强”发展战略，通过大力推进技术创新，全面优化产品结构和工艺技术结构，着力提升装备技术水平，具备了开发生产汽车零部件用特殊钢的技术条件。
为扩大产品钢在汽车零部件中的市场份额，公司开展了汽车零部件用钢的开发生产。
CHEMICALEORTHEPASTDOESITHAVEFUTURE化学EOR过去它有未来
KYPAM Polymers
Comb like with short branched chain to maintain effectiveness in high salinity brinesWide MW range for reservoirs with different permeabilitiesSuccessfully applied in some reservoirs in China
ASP Pilot Test Well Pattern, 2Z-B9-3 Well Group
Karamay Oil Field (SPE 64726)
ASP Slug Design and Injection Sequence
(SPE 64726)
PAM
PAM
Modelling Core Flood - Karamay ASP Project
Harry L. Chang
Oil Production by Polymer Flooding
Daqing Oilfield, China
A Typical PF Field Performance
Typical Pressure and Polymer Production
Oil Production by Polymer Flooding
Harry L. Chang
EOR Chemicals
Polymers and related chemicalsSurfactantsCo-surfactantsCo-solvents
Harry L. Chang
Improvements on Polymers
纬编针织物导热性能的有限元仿真
纬编针织物导热性能的有限元仿真作者：杨恩惠初曦邱华来源：《丝绸》2020年第01期摘要：為了提出一种用于优化、评估、预测针织物导热性能的方法，文章探索了纬编针织物热量传递的有限元仿真模拟，尤其是夏季汗衫面料，因为其导热性能与人体舒适性密切相关。
首先使用织物厚度仪和超景深显微镜获取织物结构参数，通过插入三次B样条曲线建立纬平针织物的三维几何模型。
然后将单元线圈模型导入有限元分析软件STARCCM+中，通过网格划分、边界设置、迭代计算步骤进行数值计算，并和实验值进行了比较分析，误差在4%以内。
表明此方法具有一定的实用性，可用于预测纬编针织物的导热性能。
关键词：纬编针织物;导热性能;三维模型;有限元仿真;网格划分中图分类号： TS186.2文献标志码： A文章编号： 10017003（2020）01003106引用页码： 011106DOI： 10.3969/j.issn.10017003.2020.01.006Finite element simulation of thermal conductivity of the weft knitted fabricYANG Enhui， CHU Xi， QIU HuaAbstract： In order to propose a method for optimizing， evaluating and predicting the thermal conductivity of knitted fabrics， the finite element simulation of heat transfer of weft knitted fabrics was explored in this paper. Especially for summer Tshirt fabric， its thermal conductivity is closely related to the comfort of human body. Firstly， fabric structure parameters were obtained by means of fabric thickness meter and ultradepthoffield microscope. And threedimensional geometric model of weftflat knitted fabric was established by inserting cubic Bspline curve. Then， the element loop model was imported into the finite element analysis software STARCCM+. The numerical calculation was carried out by mesh division， boundary setting， and iterative calculation steps. Compared with the experimental data， the error was less than 4%. It is shown that this method is practical and can be used to predict the thermal conductivity of weft knitted fabrics.Key words： weft knitted fabric; thermal conductivity; threedimensional model; finite element simulation; mesh division服装织物的热舒适性一般指人体在不同环境条件下或者不同活动条件下，外界环境与内在皮肤产生热量的交换，达到平衡时使人感到舒适的特性。
(GOOD)Mathematical modeling and thermal performance analysis
Mathematical modeling and thermal performance analysisof unglazed transpired solar collectorsM.Augustus Leon,S.Kumar*Energy Field of Study,Asian Institute of Technology,P.O.Box 4,Klong Luang,Pathumthani 12120,ThailandReceived 21June 2005;received in revised form 7March 2006;accepted 14June 2006Available online 1September 2006Communicated by:Associate Editor Charles KutscherAbstractUnglazed transpired collectors or UTC (also known as perforated collectors)are a relatively new development in solar collector tech-nology,introduced in the early nineties for ventilation air heating.These collectors are used in several large buildings in Canada,USA and Europe,effecting considerable savings in energy and heating costs.Transpired collectors are a potential replacement for glazed flat plate collectors.This paper presents the details of a mathematical model for UTC using heat transfer expressions for the collector com-ponents,and empirical relations for estimating the various heat transfer coefficients.It predicts the thermal performance of unglazed transpired solar collectors over a wide range of design and operating conditions.Results of the model were analysed to predict the effects of key parameters on the performance of a UTC for a delivery air temperature of 45–55°C for drying applications.The parametric stud-ies were carried out by varying the porosity,airflow rate,solar radiation,and solar absorptivity/thermal emissivity,and finding their influence on collector efficiency,heat exchange effectiveness,air temperature rise and useful heat delivered.Results indicate promising thermal performance of UTC in this temperature band,offering itself as an attractive alternate to glazed solar collectors for drying of food products.The results of the model have been used to develop nomograms,which can be a valuable tool for a collector designer in optimising the design and thermal performance of UTC.It also enables the prediction of the absolute thermal performance of a UTC under a given set of conditions.Ó2006Elsevier Ltd.All rights reserved.1.IntroductionDrying of fruits and vegetables require hot air in the temperature range of 45–60°C.With abundant solar radi-ation in the tropical climates of Asia,an unglazed tran-spired solar collector (UTC)system could readily provide hot air at this temperature range for almost 300days of the rge roof-mounted installations could offer a cost-effective alternative to the expensive glazed collectors to supply significant quantities of hot air for drying appli-cations.A number of ventilation air heating systems based on UTC collectors have been installed in Canada,USA and Europe.As of 2003,more than 80UTC systems cov-ering a collector area of more than 35,000m 2have been installed around the world,mostly for preheating ventila-tion air.Several UTC-based drying systems have also been installed recently,for drying fruits,nuts,coffee beans,wood chips and barks,wool,and chicken manure (SOLARWALL,2001/2006).A list of UTC installations used for drying applications is given in Table 1,to highlight the range of size,location and application.The unglazed transpired solar collector has a dark,per-forated vertical/inclined sheet metal absorber is fixed to another parallel surface or wall,with a gap of 10–15cm between them,with all sides closed and sealed.Ambient0038-092X/$-see front matter Ó2006Elsevier Ltd.All rights reserved.doi:10.1016/j.solener.2006.06.017*Corresponding author.Tel.:+6625245410;fax:+6625245439.E-mail address:kumar@ait.ac.th (S.Kumar)./locate/solenerSolar Energy 81(2007)62–75air,pulled through the perforations using a blower, absorbs the heat available at the absorber,and delivers hot air at the blower outlet.These collectors reportedly offer the lowest cost and highest efficiency(60–75%)for air heating(IEA,1999;Christensen,1998).The geometry of the UTC absorber and the operating conditions influence the thermal performance of unglazed transpired solar collectors.Providing transpiration in a UTC serves to enhance the convective heat transfer coeffi-cient between the absorber and the air stream.The suction captures the boundary layer,and significantly reduces wind heat losses;The design allows for elimination of the glazing with its associated cost and optical reflectance.Studies on UTC to investigate the heat and mass trans-fer,efficiency,airflow distribution and pressure drop have been carried out since1991.These include theoretical and empirical models to predict its thermal performance over a wide range of operating conditions.The basic heat loss theory for unglazed transpired collector is explained in detail by Kutscher et al.(1993)and Hollands(1998). Cao et al.(1993)and Golneshan and Hollands(1998), based on experimental results,reported correlations for heat exchange effectiveness for a plate with an array of slits as perforations,with windflow transverse to the slits. An empirical model for thin plates with circular holes on a triangular layout has been presented by Kutscher(1994), with the wind parallel to the plate.Kutscher(1994)and Van Decker et al.(1996)measured heat exchange effective-ness on thin and thick plates with circular holes on a square or triangular layout over a range of parameters,NomenclatureA s collector area(m2)C p.air specific heat capacity of air(J/kg K)C p.bp specific heat capacity of back plate material(J/kg K)C p.col specific heat capacity of absorber material(J/kg K)D perforation diameter(m)d plen plenum depth(m)F cs collectorÀsky view factorF cg collectorÀground view factorF cl cloud factorH absorber height(m)h cloud base height(km)h conv convective heat transfer coefficient(W/m2K) I T solar radiation incident on the collector(W/m2) K air thermal conductivity of air(W/m K)m air.out massflow rate of air through the collector(kg/s) m bp mass of back plate(kg)m col mass of absorber plate(kg)n sky fractional area of sky covered by cloudsNu Nusselt numberP pitch of perforations(m)P atm atmospheric pressure at the collector location (mbar)Q conv.air$bp convection heat transfer from air to back plate(W)Q conv.bp$amb convection heat transfer from back plate to surrounding air(W)Q conv.col$air convection heat transfer from absorber to air(W)Q rad.col$bp radiation heat transfer from absorber to back plate(W)Q rad.col-sur radiation heat transfer from absorber to ambient(W)Q rad.bp-sur radiation heat transfer from back plate to surrounding(W)Re Reynolds number T amb ambient air temperature(K)T air.out exit air temperature(K)T bp temperature of back plate(K)T dp dew point temperature(K)T col average absorber plate temperature(K)T sky sky temperature(K)T sur temperature of the surrounding(K)t hour of the day(1–24)v app approach velocity(m/s)v plen plenum air velocity(m/s)v wind free stream wind velocity over the absorber and back plate(m/s)W absorber width(m)Greek symbolsa col solar absorptance of the collector surfaceb absorber porosityD P total pressure drop across the collector(Pa)D P abs pressure drop across the absorber plate(Pa)D P fric frictional pressure drop in the plenum(Pa)D P buoy buoyancy pressure drop in the plenum(Pa)D P acc acceleration pressure drop in the plenum(Pa) e bp thermal emissivity of back plate outer surface e col thermal emissivity of the absorber surface facingthe sky/grounde col.in thermal emissivity of the absorber surface facingthe back platee csky clear sky emissivitye hc hemispherical cloud emissivitye HX heat exchange effectivenesse sky sky emissivity at collector locationm kinematic viscosity(m2/s)q air density of air(kg/m3)q col density of absorber material(kg/m3)q bp density of back plate material(kg/m3)l air viscosity of ambient air(kg/m s)g col collector efficiencyM.A.Leon,S.Kumar/Solar Energy81(2007)62–7563and presented a correlation for heat exchange effectiveness byfitting the measured data.They also correlated their results with equations involving dimensionless groups such as Reynolds number,ratio of approach velocity to wind velocity,and hole pitch to diameter ratio.Arulanandam et al.(1999)obtained a correlation for heat exchange effec-tiveness using a CFD model for a plate with circular holes on a square layout,under no-wind conditions.Van Decker et al.(2001)investigated thick and thin plates with circular holes on a square or triangular layout over a range of approach velocities and wind speeds,and pre-sented a relationship for heat exchange effectiveness. Unglazed transpired solar collectors have found most of their applications in ventilation air heating in the high-lat-itude countries of Northern America and Europe.The studies conducted so far on UTC have mainly concen-trated on its performance in such operating conditions, and very few studies have been carried out in tropical weather conditions where the solar input and ambient air temperatures are generally higher,while its use for ven-tilation air heating is limited.Also,the application of UTC for drying,which means operating the collector for higher delivery air temperatures,has not been studied in detail.Considering the potential for unglazed transpired solar collectors for drying applications,a mathematical model has been developed for predicting the thermal performance of UTC over a range of design and operating conditions suitable for drying of food products.This article presents the model developed,and an analysis of the performance of the unglazed transpired solar collector.The model pre-dicts the air temperature rise for specific operating condi-tions,and estimates the UTC design parameters for specific hot air requirements.2.Development of a mathematical model of unglazed transpired solar collector2.1.Collector configurationThe configuration of the UTC under analysis is illus-trated in Fig.1.The collector,mounted vertically,has a perforated absorber and a back plate.An air gap separates the absorber sheet from the back plate(known as the ple-num).The gaps at the four sides along the absorber edge are closed to form a box.A blowerfixed at the top of the box provides the required suction during collector operation.Table1UTC-based air heating system installations for drying applicationsLocation/year of installation Collector size(m2)Delivery air temperature/flowrateProductKorina Farms,California,USA/200346527°C/n.a.a Pecan nutsXin Zheng Feng Li Food Co.,/China,200356n.a.Jujube fruit Coopeldos R.L.,Costa Rica/2003860n.a.Coffee beanCarriere&Sons,California/200330043°C/142m3/h m2WalnutSunsweet Dryers,California,USA/2002110n.a./182m3/h m2PruneKreher’s Poultry Farms Clarence,New York,USA/200250n.a./146m3/h m2Chicken manureHoyt Farm,Casleton,New York,USA/2002130n.a.WoolCafe Duran,Panama/2002900n.a.Coffee beanBiowarme Klein St.Paul Power Plant,Austria/200110060°C/72m3/h m2Shredded wood and bark Gelee Chicken Farm,Cookshire,Quebec,Canada/2001158n.a.Chicken manure Malabar,Indonesia/199960035°C/125m3/h m2TeaSynthite Industrial Chemical Ltd.,Coimbatore,India/1998130070°C Marigoldflowers ASEAN-Canada Project on Solar Energy,Malaysia/199437060°C/n.a.CocoaCaribbean/199990,170,140(three collectors)n.a./133,44,and55m3/h m2TobaccoSource:SOLARWALL(2001/2006);GPEKS(2006);IEA(1999);Hollick(1999).a Not available.64M.A.Leon,S.Kumar/Solar Energy81(2007)62–75The heat transfer in the unglazed transpired solar collec-tor is studied by considering the overall energy balance between the components of the transpired collector.Rate equations have been used to estimate the convective and radiative heat transfer rates between the components.The model incorporates several empirical relations to cal-culate the various heat transfer coefficients used in the rate equations.2.2.AssumptionsThe assumptions used in developing the model are listed below:(i)The absorber and back plate temperatures are assumed to be uniform throughout their respective surfaces (isothermal).While metal absorbers are mostly isothermal from hole-to-hole,non-metal absorbers show some non-isothermality.However,studies by Gawlik (1995),Gawlik and Kutscher (2002)and Gawlik et al.,2005)have shown that the non-isothermality does not have a major impact on collector performance.(ii)Airflow through the perforations is assumed homoge-neous.In reality,due to buoyant flow (driven by the absorbed solar energy,especially in collectors with large height)and friction (due to forced flow)in the plenum,airflow profile could be slightly non-homo-geneous,depending on whether the buoyant flow or forced flow dominates (Gunnewiek et al.,1996;Dymond and Kutscher,1997).(iii)Convection losses from the absorber plate to theambient air are considered negligible.In large area collectors,providing a minimum pressure drop of 25Pa has been shown to reduce absorber convection losses to insignificant levels (Kutscher,1994).(iv)Flow reversal through the absorber is assumed tobe negligible.Flow reversal is a phenomenon driven by buoyancy and wind,when airflow at the top of the absorber is out of the collector rather than into it.This can however be minimized by providing a minimum pressure drop of 25Pa across the absor-ber plate (Kutscher,1994,1997;Kutscher et al.,2003).(v)The absorber is considered to be diffuse,(i.e.,onewithout directional characteristics)and gray (i.e.,one independent of wavelength)for all absorbed and emitted radiation.(vi)Losses along the plenum edge are generally not signif-icant in large area collectors (Summers David,1996),and hence are neglected.(vii)An absorber configuration with circular perforationsin triangular pitch is assumed.Studies by Kutscher (1994)and Van Decker et al.,1996)have assumed similar absorber configuration.These assumptions are consistent with studies conducted earlier.2.3.Energy balance equationsFor the present analysis,the energy balance equations for the three collector components—absorber plate,air and the back plate—have been written by considering the solar input,the heat and mass flow rates across the collec-tor,and thermal losses.The heat transfer modes and heat transfer exchanges in the collector are presented in Fig.2.The energy input to the system is from the solar radiation received on the absorber surface.The net losses from the system are due to convection and radiation losses from the absorber sur-face and backplate.M.A.Leon,S.Kumar /Solar Energy 81(2007)62–75652.3.1.Absorber platem colÃC p:colÃðd T p=d tÞ¼ða colÃI TÃA sÞÀðQ conv:col$airþQ rad:col$bpþQ rad:col$surÞð1ÞThe term Q conv.col$air refers to the heat gain from the ab-sorber plate to air,which includes heat transfer from the absorber front surface,hole and the back surface to the plenum air.Q rad.col$bp gives the radiative heat transfer from the absorber to the back plate.Q rad.col$sur is the term representing the radiative heat loss from the absorber sur-face to the surrounding.2.3.2.Plenum airðm air:outÃd tÞÃC p:airÃðd T air:out=d tÞ¼Q conv:col$airÀQ conv:air$bpð2ÞQ conv.air$bp in Eq.(2)represents the heat exchange between the back plate and airflowing through the collector.2.3.3.Back platem bpÃC p:bpÃðd T bp=d tÞ¼Q conv:air$bpþQ rad:col$bpÀQ conv:bp$ambÀQ rad:bp$surð3ÞWhile estimating the radiation loss from the back plate to the surrounding,the temperature of the surrounding T sur is calculated from the sky temperature and ground surface temperature.2.4.Rate equations2.4.1.Convection heat transferThe heat transfer rates due to convection were estimated using the convection heat transfer coefficient,the heat transfer surface area,and the temperature difference between the surface and the surroundingfluid.The follow-ing correlations were used to estimate the average convec-tion heat transfer coefficients.An average value has been used for the plenum velocity since the velocity varies from zero at the bottom to maximum at the top(Summers David,1996).2.4.1.1.Absorber plate to plenum rmation on con-vection heat transfer with suction from a hot perforated flat plate is rather limited(Kutscher et al.,1993).The empirical correlation reported by Kutscher(1994)is used in the present model to estimate the Nusselt number.Nu1¼2:75ýðP=DÞÀ1:21ÃRe0:431þ0:011ÃbÃRe1Ãðv wind=v appÞ0:48 ð4ÞwhereRe1¼ðq airÃv holeÃDÞ=l airð5ÞThe convection heat transfer coefficient between the absor-ber plate and plenum air is then given byh conv:col$air¼ðNu1ÃK airÞ=Dð6ÞNow,Q conv.col$air=h conv.col$air*(T colÀT air.plen).2.4.1.2.Plenum air to back plate.The Nusselt number for convection heat transfer between the plenum air and back plate,Nu2¼0:664ÃðRe2Þ0:5ÃðPr2Þ0:333ð7ÞwhereRe2¼ðq airÃv plenÃHÞ=l airð8ÞPr2¼ðC p:airÃl airÞ=K airð9ÞConvective heat transfer coefficient between the air and the back plate is estimated fromh conv:air$bp¼ðNu2ÃK airÞ=d plenð10ÞThe convection heat transfer between the plenum air and back plateQconv:air$bp¼h conv:air$bpÃðT air:plenÀT bpÞð11Þ2.4.1.3.Absorber plate to ambient.Convection heat losses from the absorber surface of an UTC to the ambient air are negligible during its normal operation as the convective boundary layer is continuously sucked off(Kutscher,1992, 1994;Arulanandam et al.,1999).Providing a minimum pressure drop of25Pa and a minimum approach velocity (estimated as theflow rate per unit area of collector at which outside air is drawn into the collector)of0.02m/s has been shown to reduce absorber convection losses to insignificant levels in large area collectors(Kutscher, 1994;Summers David,1996).As the present analysis con-siders25Pa and0.02m/s as the lower limits for absorber pressure drop and approach velocity respectively,the con-vective heat loss component has been neglected.2.4.1.4.Back plate to ambient.The Nusselt number in the above expression is calculated using the following correla-tion(Cengel and Turner,2001):Nu3¼0:664ÃðRe3Þ0:5ÃðPr3Þ0:333ð12ÞRe3¼ðq airÃv windÃWÞ=l airð13ÞConvection heat transfer coefficient:h conv:bp$amb¼ðNu3ÃK airÞ=L cð14ÞNow;Qconv:bp$amb¼h conv:bp$ambÃðT bpÀT ambÞð15Þ2.4.2.Radiation heat transferThe heat transfer rates due to radiation were estimated using the Stefan–Boltzmann law,involving the total radiat-ing area,absolute temperature of the radiating body,and Stefan–Boltzmann constant.66M.A.Leon,S.Kumar/Solar Energy81(2007)62–752.4.2.1.Absorber plate to ambient.Radiation heat loss from the absorber surface occurs both to the sky and to the ground,with the view factor depending on the absorber tilt.Noting that the absorber has been considered as gray and diffuse,the radiation heat loss from the absorber isQrad:col$sur ¼e colÃr sbÃA sÃðT4colÀF cs T4skyÀF cg T4gndÞð16ÞSeveral correlations are available for the estimation of sky temperature T sky(Bliss,1961;Swinbank,1963;Berger et al.,1984;Martin and Berdahl,1984).In the present study,a sub-routine based on the correlation developed by Martin and Berdhal has been used.T sky¼e0:25skyÃT ambð17Þwhere sky emissivitye sky¼e cskyþð1Àe cskyÞÃn skyÃe hcÃF clð18Þand the clear sky emissivitye csky¼0:711þ0:56T dp100þ0:73T dp1002þ0:013cosp t12h iþ0:00012ðP atmÀ1000Þð19Þwhere t is the hour of the day.The cloud factor F cl is com-puted by the sub-routine internally from the solar radiation data.The ground surface temperature T gnd depends on several factors such as the surface exposure to insolation,wind speed,soil moisture,vegetation cover and vegetation height,and its determination is complex(Signorelli and Kohl,2004).Considering that the ground surface tempera-ture variation coincides with the local air temperature var-iation(Bi et al.,2004),and that the the daily mean surface temperature is commonly near the mean air temperature (GSFC,2006),T gnd has been taken as the ambient air tem-perature T amb.Radiative heat transfer rates from absorber plate to back plate,and from back plate to the surroundings were estimated using the temperature and emissivity values of the collector components.2.4.2.2.Absorber plate to back plate.The radiation heat transfer between the absorber plate and back plate was esti-mated using the following relation:Qrad:col$bp ¼r sbÃA sÃðT4colÀT4bpÞ=ð1=e col:inþ1=e bpÀ1Þð20ÞThe emissivity of the inner surface of the absorber is taken to be different from that of the outer surface since the outer surface generally has a coating of high absorptivity.The in-ner surface of the absorber is normally bare,free of any coating.2.4.2.3.Back plate to surrounding.The radiation heat trans-fer from the back plate to the surrounding was estimated using the following relation:Qrad:bp$sur¼e bpÃr sbÃA sÃðT4bpÀT4surÞð21ÞThe emissivity of both the surfaces of the back plate are as-sumed to be equal.2.5.Collector efficiencyThe efficiency of the unglazed transpired solar collector is the ratio of the useful heat delivered by the solar collector to the total solar energy input on the collector surface.The useful heat delivered by the collector can be estimated from the temperature andflow rate of the outlet air.Therefore,gcol¼m air:out C p:airðT air:outÀT ampÞI T A sð22Þ2.6.Heat exchange effectivenessHeat exchange effectiveness of an UTC depends on the overall heat transfer coefficient for the air passing through the absorber,and is defined as the ratio of the actual tem-perature rise of air as it passes through the absorber plate to the maximum possible temperature rise(Kutscher, 1994).Mathematically,e HX¼T air:plenÀT ambT colÀT ambð23ÞHeat exchange effectiveness e HX between the absorber and plenum air can be estimated using the relationship based on logarithmic mean temperature difference(LMTD)ap-plied for heat exchangers(Kutscher,1992,1994).e HX¼1Àexp½Àh conv:col$airÃA s=ðm air:outÃC p:airÞ ð24ÞThese two expressions provide a relation between the exit air temperature,absorber temperature and ambient temperature.2.7.Pressure dropTotal pressure drop across the collector(D P)is a sum of pressure drop across the absorber plate(D P abs),and pres-sure drops in the plenum,which include frictional pressure drop(D P fric),buoyancy pressure drop(D P buoy)and accel-eration pressure drop(D P acc)(Kutscher,1994;Dymond and Kutscher,1997).i:e:;D P¼D P absþD P fricÀD P buoyþD P accð25ÞBuoyancy force tends to push the plenum air up,acting in the direction opposite to that of frictional force.Pressure drop across the absorber plate has to be at least 25Pa to ensure a uniformflow and temperature distribu-tion over the collector(Kutscher,1997).If temperature dis-tribution is not uniform,hot spots could develop on the collector surface,which will increase the radiation loss to the surroundings(Kokko and Marshal,1992).M.A.Leon,S.Kumar/Solar Energy81(2007)62–75673.Data and solution procedureThe energy balance and rate equations described in Section2were solved tofind the exit air temperature, heat exchange effectiveness,collector efficiency and heat delivered for a given set of input values,and for a given time period.The intermediate values would include the estimation of the heat transfer coefficients and heatflux between the UTC components.For subsequent times, the initial temperatures are specified equal to the temper-ature values of the respective previous sections.The iterative process is continued until the end time is reached.As initial condition,the absorber plate,back plate and plenum air are all assumed to be at ambient temperature.As noted by earlier studies(Kutscher et al.,1993;Kut-scher,1994;Van Decker et al.,1996,2001),the key param-eters of an UTC system include perforation diameter, pitch,airflow rate/approach velocity,collector absorptiv-ity/emissivity,solar radiation,and wind velocity.A para-metric analysis was conducted for a range of perforation diameter–pitch combinations,radiation and airflow rates. Triangular pitch was assumed for the absorber perfora-tions.The input parameters considered for the analysis are:(a)porosity(involving perforation diameter and pitch),(b)approach velocity,(c)solar radiation,and(d) solar absorptivity/thermal emissivity.The range of parameters selected for the present analysis was as below.(i)Approach velocity:For drying of most food products,a drying air temperature of50–60°C is required,which can be attained only at low approach veloci-ties.Tests at NREL(National Renewable Energy Laboratory,USA)showed that below the approach velocity of0.02m/s,performance is worse than theo-retical prediction perhaps because of natural convec-tion effects,or non-homogeneous suction.Theseflow rates result in high absorber surface temperatures and associated low collector efficiency.Considering this minimum requirement,in the present study,an approach velocity range of0.02–0.03m/s was selected,which is equivalent to an airflow rate of 72–108m3/h m2of collector area.However,in some cases,the range has been increased to0.0125–0.0375m/s.(ii)Solar radiation:A solar radiation range of400–900W/m2was selected for the analysis,which reason-ably corresponds to the daily average solar radiation in tropical locations.(iii)Absorber pressure drop:A minimum absorber pres-sure drop of25Pa is required to ensure uniformflow distribution through the collector.However,a slightly lower pressure drop may be allowable if the approach velocity is maintained above0.02m/s (Summers David,1996).The upper limit to absorber pressure drop is guided by the power consumption bythe suction blower,which will become excessive at pressure drops higher than80Pa.A pressure drop range of25–80Pa was therefore selected for the cur-rent study.(iv)Collector plenum:The frictional pressure drop across the collector was estimated using the model over a plenum(distance between the absorber and back plate)range of50–150mm.For the selected collector parameters,a plenum of120mm was found to be suf-ficient to keep the frictional pressure drop to a mini-mum.A higher value may be required for larger airflow rates as in the case of ventilation air heating applications.(v)Pitch and perforation diameter:Considering past studies on the pitch of perforations(Kutscher,1992, 1994;Van Decker and Hollands,1999;Arulanandam et al.,1999),a pitch range of12–24mm was selected for the present study,which corresponds to a perfora-tion diameter of0.80–1.55mm.Lower values of pitch dictate absorbers with very small perforation diame-ters,and could pose manufacturing difficulties.A per-foration diameter of0.8mm has been used in commercial rger values of pitch tend to lower the collector efficiency and heat exchange effectiveness.Table2summarises the input parameters and the range of their values used in the present study.The output parameters estimated were(a)collector effi-ciency,(b)heat exchange effectiveness,(c)air temperature rise(or delivery air temperature),and(d)useful heat deliv-ered.The effect of varying the input parameters on these were also studied.For the collector with mild steel absorber coated with flat black paint,within the selected range of solar radiation and airflow rate,the UTC simulation was performed for a number of combinations of pitch and perforation diameter. The results are discussed in the following sections.Table2Input parameters and their values used in the studyS.No.Input parameter Range1Approach velocity0.02–0.03m/s 2Solar radiation400–900W/m2 3Ambient temperature30°C4Wind velocity 1.2m/s5Pressure drop across the absorber25–80Pa6Plenum depth120mm7Pitch(triangular)12–24mm8Perforation diameter0.80–1.55mm 9Absorber material Mild steel10Design parameters used for reference collectorSolar absorptance0.95Thermal emittance0.85Pitch20mmPerforation diameter 1.25mm68M.A.Leon,S.Kumar/Solar Energy81(2007)62–75。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
Thermal Modelling,Simulation and Evaluation of a High Power Battery Cell for AutomotiveApplicationsDragan SimicAIT Austrian Institute of Technology,Giefinggasse2,1210Vienna,Austriafax:+43–50550–6347,e-mail:dragan.simic@ait.ac.atAbstract—This paper deals with the thermal and elec-trical modelling and simulation of an automotive high power battery cell.The electrical modelling is imple-mented for determining the electrical losses.The elec-trical losses are converted to heat losses and therefore a heatflow in the cell.The thermal modelling is based on a1D implementation by means of discrete volume elements.The resulting simulation model includes the thermal and electrical behavior of a high power cell.All simulation models of the investigated cell are imple-mented and developed in the simulation environment Dymola,using Modelica programming language,[1]. Thermal and electrical measurements have been car-ried out to validate the simulation results of the thermal and electrical modelling.Three tests setups have been developed to measure electrical and thermal signals for determination of the parameters for the simulation model.The measurement and simulation results will be presented in this paper.I.IntroductionA T the moment electrical vehicles and electrical pow-ertrain components are allocated in niche markets. However,they will considerably gain market share when electrical component development and production get cheaper.Furthermore the oil prize will increase constantly since global oil resources are said to go down and end soon.Full electric vehicles that could meet consumers’expectations are still far from realization since suitable electric energy storages that facilitate acceptable cruising ranges have not been developed yet.Considering these facts,a thermal simulation of a high power cell was developed by Austrian Institute of Technol-ogy(AIT).Some simulations using the investigated high power cell will be executed and simulation results will be provided.Thefinal thermal and electrical results of the high power cell will be compared with measurements. The investigation of new vehicle components at AIT is based on the interconnection of three phases;simulation, methods and components and testing and validation.The development process starts with modelling and simulation of the thermal,electrical and mechanical sub-components as well as the entire vehicle.Subsequently,methods and components have to be developed andfinally testing and validation on components and system level has to be carried out.Each thermal,electrical and mechanical component can be optimized,validated and realized very effectively using this process.II.Simulation EnvironmentsLibraries developed for simulations of automotive appli-cations at the AIT usually are based on the Modelica simulation language.This libraries consist of elementary components for developing vehicle drive systems.Two simulation libraries,the SmartElectricDrives(SED)li-brary and the SmartPowerTrains(SPT)library can be used to model the components of an electrical vehicle or vehicle systems.SmartElectricDrives LibraryThe SED library is a tool for modeling,simulation and tuning of electric drives in electromechanical systems,ac-cording to[2],[3],[4].This library is particularly designed for simulations of automotive applications with a Modelica development platform.The SED library contains elementary components for developing drive systems and the powerful’ready-to-use’models.These’ready-to-use’library models include the characteristics of fully controlled modern electric drives in only one component.Furthermore,the SED library supports modeling and simulation on different levels of abstraction.The users can choose between transient and quasi stationary models.In the latter case,electrical transient effects are neglected.SmartPowerTrain LibraryThe SPT library includes models and blocks for modeling and simulation of mechanical components of a vehicle, in this work were used only basic components from this library,according to[5],[6].The library allows longitu-dinal dynamic simulations of vehicles.All driving effects such as aerodynamic and rolling resistances are included in models of this library.Figure1.Simulation model of the high power cell in Modelica Additionally,this library contains components for mod-eling and simulation of thermal effects of a vehicle.A virtual driver with operating strategy for controlling the acceleration pedal position,the brake pedal position and the clutch pedal position is implemented in this library, too.The communication between these two developed libraries as well as the independent control of each model allows the highestflexibility in the design of new vehicle concepts. All’ready-to-use’models from both libraries include a bus connector extended from ModelicaStandard library for providing external communication.This type of model-ing allows a high harmonizing level between developed libraries.III.Thermal and Electrical Model of the HighPower CellThe modelling was conducted in two steps.Thefirst step is the electrical modelling of the high power cell.The electrical model of the high power cell comprises algebraic and ordinary differential equation.The electrical model of the high power cell is presented infigure1.The cell is modelled with two resistances(Rs and Rd)and one capacitor(Cd)model,they are connected with the positive and negative electrical pins(pin_p and pin_n) of the battery.The determination and calculation of the state of charge(SOC)is based on the integration of the battery current.The open circuit voltage(OCV)is provided as characteristic curve(OCVTable).The second step is the thermal modelling of the cell.The model is developed by means of discrete volume elements, [7].In this thermal model the coefficient of heat transfer for each discrete volume is calculated(thermal infigure 1),according to[8],[6].The heatflow inside the high power cell was regarded in all three directions.The model is parametrized using geometrical and thermal measurement data,such as thermal conductivity,specific Figure2.Thermal model of the high power cell in Modelicathermal capacity,density and thickness of the layer mate-rials(anode,cathode,separator,etc),etc.Heat losses in the thermal model are generated using equation(1).˙Qflow=Rd·i2Rd+Rs·i2Rs(1)Figure2shows the distribution of the discrete volume of investigated cubical-shaped cell in all three directions,x, y and z.The discrete volume number can be defined in each direction.For investigation of the thermal behaviour of the battery package which includs cylindrical cells,first a detailed thermal model of one cell should be implemented.This base model includes the same equations and models as described in equation(1).The heat losses distribution of the anode and cathode in this cell is implemented homogeneously.This heatflow in x and y direction,cross directions,has the same size and conductivity coefficents.Therefore the model of the cell is reduced from a3D to a2D problem. This reduced model of the cubical-shaped cell is used in this paper for evaluation of the simulation model.The discrete volume(discrete volume infigure2)consists of six thermal connectrors(heat ports).In each direction exist two connectors to connect the discrete volume with the next.Each of the connector includs two variables,temperature and heatflow rate.This connectors allow the heatflow in all three,x,y and z directions.The connvection to air(heat exchange with environment)is ensured based on connection the face discrete volumes with environ-ment components using outside connector(heatPort in figure1)IV.Experimental Setup and MeasurementThe thermal measurement was performed in three stages. Thefirst steps includes the determination of the specificFigure3.The testing scheme for determination of the specific heat capacity of the high power cellheat capacity.A cubical-box is build and bandaged with an isolation material(isolation),figure3.This box is filled with water at a known temperature(in experimental sepup20◦C).The water volume was exactlly defined. The battery is cooled in a climate chamber to get a high temperature difference between water and battery. When the battery is immersed in the water a steady state temperature is approached.The stationaryfinal temperature is then measured.These three temperatures (start temperature of the cell,T1c,start temperature of the water,T1w,and the stationaryfinal temperature, T end)can used to calculate the specific heat capacity,c pc, of the high power cell,equation(2)and(3).c pc=m w·c pw·(T1w−T2w)m c·(T1c−T2c)(2) T2w=T2c=T end(3)The calculated specific heat capacity is used for parame-terization of the thermal model of the high power cell in Modelica.This specific heat capacity was also compared with the data sheet of the cell to minimize the parame-terization error.The second step deals with the thermal measurement of the high power cell.The cell is inserted into a ther-mal chamber with a constant air temperature of about 41◦C.A thermal sensor is mounted on the cell face (temperature sensor),figure4.During the measurement the thermal and electrical values like cell temperature,en-vironment temperature,current and voltage are recorded. With the third step the thermal conductances of the high power cell are determined.Three high power cells are stacked in a package.One temperature sensor is mounted on the cell face of thefirst cell and one temperature sensor on the cell face of the third cell to measure the temperature difference between this twofaces.Figure4.The testing scheme for measuring the heating of the high power cellThefirst cell was heated with a heat source.The heatflow of the heat source was constant.The third cell was cooled with a metal plate.The temperature of the metal plate was constant.Based on the measured temperatures and heatflow rate of the heat source was determined the ther-mal conductance of this three high power cells.Therefore the thermal conductance was re-calculated and prepared to characteristic thermal conductance of a material,in this case for high power cell in z direction(figure2).V.Simulation Results and ValidationFor the validation of the high power cell model the mea-sured and simulated values have been compared.Figure5 shows the congruence between measured(i measured)and simulated(i simulated)current profile.The current profile for charging and discharging of the high power cell is used as a typical test profile.In thisfigure5are the both current profile the same.The current profile were depicted to verify that the current profile of the cell-testing bench and simulated current profile are the same.The comparison of the measured and simulated voltages is presented infigure6.It can be seen that the measured and simulated voltage correspond very well in the overall regarded time frame.This small deviation verifies that the quasi stationary electrical behavior of the high power cell was implemented in a satisfying way.The transient electrical behavior of the high power cell was not modelled in this work.The transient modelling is not necessary to implement if the electrical model is only used for simulation the heat losses.The thermal time constants of a high power cell are relatively higher as the electrical time constants.The diference between the electrical quasi stationary and electrical transient modelling,in this case, is very low according to incluence of the simulated thermal results.Afterwards,the simulation results of the thermal model were compared with measurement results.Figure7showsFigure 5.Measured and simulated current profile of the high powercell Figure 6.Measured and simulated voltage of the high power cellthe temperature profiles of the measured and simulated high power cell,where Tambience is the temperature of the environment,Tmeasured is the face temperature of the cell and Tsimulated is the face temperature of the simulated cell.The simulated temperature behavior is al-most identical with the measured temperature.This good congruence between measured and simulated temperature of the cell proves that the thermal behavior of the high power cell was also modelled in a satisfying way.VI.ConclusionIn this contribution an electrical and a thermal model of a high power cell were presented.The results of the distinct approaches were shown.Electrical and thermal measure-ments have been carried out to validate the simulation of the thermal modelling and were presented.The thermal model of the high power cell was implemented in Modelica and simulated using the Dymola simulation environment.The described cell simulation allows the determination of the temperature behavior in a large number of applica-tions.Due to accurate simulations a significant accelera-tion of the development process of electricalcomponentsFigure 7.Measured and simulated temperature of the high power cell and measured temperature of the environmenthas been achieved and effort and costs have been drasti-cally reduced.References[1]Dymola,Dynamic Modeling Laboratory,User‘s Manual ,Dy-nasim AB,2004,.[2]T.Baeuml,H.Giuliani,D.Simic,and F.Pirker,“An advancedsimulation tool based on physical modelling of electric drives in automotive applications”,IEEE Vehicle Power and Propulsion Conference,VPPC,Arlington,USA ,2007.[3]J.V.Gragger,D.Simic,C.Kral,H.Giuliani,V.Conte,andF.Pirker,“A simulation tool for electric auxiliary drives in hevs -the SmartElectricDrives library”,World Automotive Congress,FISITA ’06,Yokohama,Japan ,2006.[4]J.V.Gragger, D.Simic, C.Kral,and F.Pirker,“A fullvehicle simulation of an HEV starter-generator concept with the smartelectricdrives library”,SAE World Congress,2007,Detroit,USA ,April 2007.[5]H.Giuliani,D.Simic,J.V.Gragger,C.Kral,and F.Pirker,“Optimization of a four wheel drive hybrid vehicle by means of the SmartElectricDrives and the SmartPowerTrains library”,International Battery,Hybrid and Fuel Cell Electric Vehicle Symposium &Exposition,EVS22,Yokohama,Japan ,pp.1983–1993,2006.[6]Dragan Simic,Hannes Lacher,Christian Kral,and Franz Pirker,“Evaluation of the smartcooling(sc)library for the simulation of the thermal management of an internal combustion engine”,SAE World Congress,Detroit,USA ,2007.[7]H.Oberguggenberger and D.Simic,“Thermal modeling ofan automotive nickel metall hydrid battery in modelica using dymola”,International Modelica Conference,6th,Bielefeld,Germany ,,no.485,pp.77–81,2008.[8] D.Simic,C.Kral,and F.Pirker,“Simulation of the coolingcircuit with an electrically operated water pump”,IEEE Vehicle Power and Propulsion Conference,VPPC,Chicago,USA ,2005.Definition,Acronyms and Abbreviations AIT Austrian Institute of Technology SED SmartElectricDrives SPT SmartPowerTrains SOC state of chargeOCVopen circuit voltage。