Hattel, Jesper H.

Klingaa, Christopher G.

Rasmussen, Filip Salling

Sonne, Mads R.

The objective of the present report is to give a thorough description on how to implement a 3-Dimensional heat transfer solver applicable for the pultrusion process. The solver will use the Alternate Direction Implicit method modified by Douglass &amp; Gunn (3D-ADI-DG) [1]. A similar solver was introduced for the first time in this specific field of application by Baran et al. [2]. The current report does however serve as a thorough introduction for the reader to create theirown solver with the associated benefits, e.g. understanding and control of all aspects of the solver. The solver was developed, verified and validated following previous work by the authors [4] and used to investigate the effect of varying fiber volume fraction distribution on the cure behaviour during the pultrusion process [5]

Online Research Database In Technology

  General rights Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights.   Users may download and print one copy of any publication from the public portal for the purpose of private study or research.  You may not further distribute the material or use it for any profit-making activity or commercial gain  You may freely distribute the URL identifying the publication in the public portal  If you believe that this document breaches copyright please contact us providing details, and we will remove access to the work immediately and investigate your claim.         Downloaded from orbit.dtu.dk on: Mar 30, 2019Numerical Modelling of Heat Transfer using the 3D-ADI-DG Method - with Applicationfor Pultrusion.Technical reportRasmussen, Filip Salling; Klingaa, Christopher Gottlieb; Sonne, Mads Rostgaard; Hattel, Jesper HenriPublication date:2019Document VersionPublisher's PDF, also known as Version of recordLink back to DTU OrbitCitation (APA):Rasmussen, F. S., Klingaa, C. G., Sonne, M. R., & Hattel, J. H. (2019). Numerical Modelling of Heat Transferusing the 3D-ADI-DG Method - with Application for Pultrusion.: Technical report. Kgs. Lyngby: TechnicalUniversity of Denmark (DTU).Technical University of DenmarkDepartment of Mechanical EngineeringNumerical Modelling of Heat Transferusing the 3D-ADI-DG Method- with Application for Pultrusion.Technical ReportbyFilip Salling RasmussenCo-authors: Christopher G. Klingaa,Mads R. Sonne and Jesper H. HattelSection for Manufacturing Engineering1st Edition, March 2019Copyright c©, Filip Salling Rasmussen, 2019Department of Mechanical EngineeringTechnical University of DenmarkDK 2800 Kgs. Lyngby11 IntroductionThe objective of the present report is to give a thorough description on how to implement a3-Dimensional heat transfer solver applicable for the pultrusion process. The solver will use theAlternate Direction Implicit method modified by Douglass & Gunn (3D-ADI-DG) [1]. A similarsolver was introduced for the first time in this specific field of application by Baran et al. [2].The current report does however serve as a thorough introduction for the reader to create theirown solver with the associated benefits, e.g. understanding and control of all aspects of thesolver. The solver was developed, verified and validated following previous work by the authors[4] and used to investigate the effect of varying fiber volume fraction distribution on the curebehaviour during the pultrusion process [5].The Pultrusion ProcessPultrusion is a continuous process for manufacturing of Fiber Reinforced Polymers (FRP) withconstant cross-section. The fibers are either impregnated in a resin bath or by resin injectiondirectly in the die - the latter is illustrated in Fig. 1.1. The impregnated fibers are drawnthrough the heating/forming die where heaters ensure proper curing of the thermoset resin.Figure 1.1: 2D illustration of the resin injection pultrusion process [4].It will be assumed that the fibers are fully impregnated and no significant curing have happenedbefore the heating/forming die. Hence, the numerical solution is the discretized temperaturefield in the heating/forming die and the FRP composite. The reader is made aware that thepresent report will solely focus on heat transfer modelling and the reader is referred to [2] and[4] for a description on how to include the chemical modelling, i.e. curing.2 Heat Transfer Modelling - PultrusionThe governing equation for the convection-diffusion problem at hand is the energy equation∂ (ρcpT )∂t+∇ · (ρcpTu)−∇ · (λ∇T ) = S′′′T (2.1)where ρ is density in[kgm3], cp is specific heat capacity in[JkgK], T is temperature in [◦C], t istime in [s], u is velocity in[ms], λ is the thermal conductivity (tensor) in[WmK]and S′′′T is avolumetric source term in[Wm3]. The effective material properties are calculated using the ruleof mixtures and assumed constant. The following notation is equivalent to (2.1) and will be usedto ease the understanding of the time discretization in section 2.2.2Q˙+ qu − qλ = S′′′T (2.2)In (2.2): Q˙ is the change of heat content over time and qu, qλ denotes the convective anddiffusive heat flux, respectively. In the current work, energy conservation, as described in (2.1),will be modelled using a control volume based finite difference model (CV-FDM). The spatialdiscretization, the choice of time scheme, the implementation algorithm and an example ofsimulation results will be presented in the following four sections.2.1 Spatial DiscretizationThe principal of energy conservation is used to set up the discretized version of the energyequation. Fig. 2.1 illustrates the 3D notation used henceforth. The material properties ρ,cp and λ will be assumed constant. The discretization will account for a non-uniform butstructured mesh and will only consider convection in the pulling direction (x) and further assumethis velocity to be positive and constant. Hence, u = [ux, 0, 0] = const. and ux > 0. Thediffusive heat fluxes will be evaluated using a 3-point Central Difference Scheme (CDS), whilethe convective heat flux will be evaluated using a 2 point Upwind Difference Scheme (UDS).Hence, the discretized energy equation over a control volume ∆x∆y∆z becomesρcp∂T∂t ∆x∆y∆z+ρcpuxTP−TW∆xw∆x∆y∆z−λx∆y∆z∆xe(TE − TP )− λx∆y∆z∆xw (TP − TW ) +λy∆x∆z∆yb(TB − TP )− λy ∆x∆z∆yf (TP − TF ) +λz∆x∆y∆zn(TN − TP )− λz ∆x∆y∆zs (TP − TS)= S′′′T ∆x∆y∆z(2.3)The reader is encouraged to realise that (2.3) is equivalent to (2.1) - let (2.3) be divided withthe volume and subsequently let ∆x, ∆y and ∆z approach zero.The choice of time stepping scheme will be explained in detail in the following subsection. Hence,for now the time dependent term in (2.3), related to the change of heat content, will simply bediscretized as ∆x∆y∆z · ρcp · ∆TP∆t .For the purpose of numerical implementation, proper weight functions are introduced and (2.3)can be written asH0∆TP+Hu (TP − TW )−He (TE − TP )−Hw (TP − TW ) +Hb (TB − TP )−Hf (TP − TF ) +Hn (TN − TP )−Hs (TP − TS)= ST(2.4)where ST = S′′′T ∆x∆y∆z, and the weight functions are defined in the following way.3The capacity weight function:H0 =ρcp ·∆x∆y∆z∆t(2.5)The convective weight function:Hu = ρcpux∆y∆z (2.6)The conductivity weight functions in the x, y and z direction, respectively, using the notationillustrated in fig. 2.1 and 2.2:Hw = λx∆y∆z∆xwHe = λx∆y∆z∆xe(2.7)Hf = λy∆x∆z∆yfHb = λy∆x∆z∆yb(2.8)Hs = λz∆x∆y∆zsHn = λz∆x∆y∆zn(2.9)Figure 2.1: Schematic illustration of the 3D finite volume grid with labelling of the CV neighbours andCV faces.Figure 2.2: Schematic illustration of notation used for the structured non uniform mesh (x -direction).42.2 Time DiscretizationThe numerical discretization in time will be conducted using the Alternate Direction Implicitmethod modified by Douglas & Gunn (ADI-DG). The original 3D-ADI method implies dividingeach time step into three minor, consecutive time sub-steps equal in size [3]. In the modificationby Douglas & Gunn the three time sub-steps are not equal in size and should hence, to someextend, be interpreted as overlapping - this is illustrated in fig. 2.3. The advantage of usingthe modified 3D-ADI-DG method is good stability even for problems with a large convection todiffusion heat transfer ratio (i.e. large Peclet number, Pe > 2, see [2]) - which is typically thecase for pultrusion of FRPs [4].Figure 2.3: Illustration of the time sub stepping scheme for the ADI-DG algorithm.For simplicity the direction chosen to be implicit for each time sub-step are x, y and z inchronological order. Hence, for each of the three time sub-steps the energy equation (2.2) couldbe written as follows.Sub-step I:Q˙t→t+∆t3 + qt→t+ ∆t3u −(qt→t+ ∆t3x + qty + qtz)= StT (2.10)Sub-step II:Q˙t→t+2∆t3 + qt→t+ ∆t3u −(qt→t+ ∆t3x + qt→t+ 2∆t3y + qtz)= StT (2.11)Sub-step III:Q˙t→t+∆t + qt→t+∆t3u −(qt→t+ ∆t3x + qt→t+ 2∆t3y + qt→t+∆tz)= StT (2.12)It is important to notice the unequal size of the three sub-steps, which cause the capacity weightfunction (2.5) to change for each sub-step. Hence, introducing the following three sub-stepcapacity weight functions:H0x =ρcp ·∆x∆y∆z∆t/3(2.13)H0y =ρcp ·∆x∆y∆z2∆t/3(2.14)H0z =ρcp ·∆x∆y∆z∆t(2.15)The notation for the sub-step capacity weight functions is summarized in Table 1 together witha short notation for the sub-steps - this notation will be used henceforth.5old sub step I sub step II sub step IIITime t t+ ∆t3 t+2∆t3 t+ ∆tNotation t tx ty tzCapacity weight functions H0x H0y H0zTable 1: Notation used for the ADI-DG time scheme and the sub-step capacity weight functions.2.3 Implementation of 3D-ADI-DGThe ADI-DG method involves setting up a tri-diagonal matrix system for every row or column ofcontrol volumes (CV) in each of the three major directions. Hence, the following matrix systemis obtainedA · T = b (2.16)where A is the tri-diagonal coefficient matrix, T is a vector with the unknown temperaturesof the row/column of CVs and b is the the right hand side (RHS) vector. The coefficients inthe matrix A and the known values in the RHS vector b is obtained from (2.4) with the timestepping schemes illustrated in Fig. 2.3 and written in (2.10)-(2.12). The coefficients in the firstand last row of the system matrix A as well as the first and last element in b will be defined byeither a Dirichlet or a Neumann boundary condition as described in [3]. The final matrix system(2.16) is then readily solved using e.g. the backslash operator in MATLAB1.In the following let nx, ny, nz and i, j, k represent the number of CVs and the index in thethree spatial directions x, y, z, respectively.Sub-step I: Implicit formulation in the x -direction.This sub-step goes form t to t + ∆t3 ; see (2.10), Fig. 2.3 and Table 1. For each combination ofj = [2 : ny − 1] and k = [2 : nz − 1], running through all inner CVs in the x-direction, i = [2 :nx − 1], a coefficient matrix Ax of size nx × nx and a RHS vector bx with nx elements are setup by evaluating (2.4) using the time scheme (2.10). Hence,H0x(T txP − T tP)+ 12Hu{T txP − T txW + T tP − T tW}+ 12Hw{T txP − T txW + T tP − T tW}+ 12He{T txP − T txE + T tP − T tE}+Hf (TtP − T tF ) +Hb (T tP − T tB)+Hs (TtP − T tS) +Hn (T tP − T tN )= StT(2.17)where the lower-diagonal (LD), diagonal (D) and upper-diagonal (UD) elements of the coefficientmatrix Ax are directly given by the coefficients of TtxW , TtxP and TtxE , respectively. The RHS vectorbx are all the explicit parts of (2.17)2. Hence, the tri-diagonal coefficients in Ax and the RHSvector elements in bx are given by1https://se.mathworks.com/help/matlab/ref/mldivide.html2parts evaluated at t6ALDx = −12Hu − 12Hw (2.18)ADx = H0x +12Hu +12Hw +12He (2.19)AUDx = −12He (2.20)bx =H0xTtP− 12Hu {T tP − T tW }+ 12Hw {T tW − T tP }+ 12He {T tE − T tP }+Hf (TtF − T tP ) +Hb (T tB − T tP )+Hs (TtS − T tP ) +Hn (T tN − T tP )+StT(2.21)Sub-step II: Implicit formulation in the y-direction.This sub-step goes form t to t+ 2∆t3 ; see (2.11), Fig. 2.3 and Table 1. For each combination ofi = [2 : nx − 1] and k = [2 : nz − 1], running through all inner CVs in the y-direction, j = [2 :ny − 1], a coefficient matrix Ay of size ny × ny and a RHS vector by with ny elements are setup by evaluating (2.4) using the time scheme (2.11). Hence,H0y(TtyP − T tP)+ 12Hu{T txP − T txW + T tP − T tW}+ 12Hw{T txP − T txW + T tP − T tW}+ 12He{T txP − T txE + T tP − T tE}+ 12Hf{TtyP − T tyF + T tP − T tF}+ 12Hb{TtyP − T tyB + T tP − T tB}+Hs (TtP − T tS) +Hn (T tP − T tN )= StT(2.22)where the lower-diagonal (LD), diagonal (D) and upper-diagonal (UD) elements of the coefficientmatrix Ay are directly given by the coefficients of TtyF , TtyP and TtyB , respectively. The RHS vectorby are all the explicit parts of (2.22)3. Hence, the tri-diagonal coefficients in Ay and the RHSvector elements in by are given byALDy = −12Hf (2.23)ADy = H0y +12Hf +12Hb (2.24)AUDy = −12Hb (2.25)3parts evaluated at t and tx7by =H0yTtP− 12Hu{T txP − T txW + T tP − T tW}+ 12Hw{T txW − T txP + T tW − T tP}+ 12He{T txE − T txP + T tE − T tP}+ 12Hf {T tF − T tP }+ 12Hb {T tB − T tP }+Hs (TtS − T tP ) +Hn (T tN − T tP )+StT(2.26)Sub-step III: Implicit formulation in the z -direction.This sub-step goes form t to t + ∆t; see (2.12), Fig. 2.3 and Table 1. For each combination ofi = [2 : nx − 1] and j = [2 : ny − 1], running through all inner CVs in the z-direction, k = [2 :nz − 1], a coefficient matrix Az of size nz ×nz and a RHS vector bz with nz elements are set upby evaluating (2.4) using the time scheme (2.12). Hence,H0z(T tzP − T tP)+ 12Hu{T txP − T txW + T tP − T tW}+ 12Hw{T txP − T txW + T tP − T tW}+ 12He{T txP − T txE + T tP − T tE}+ 12Hf{TtyP − T tyF + T tP − T tF}+ 12Hb{TtyP − T tyB + T tP − T tB}+ 12Hs{T tzP − T tzS + T tP − T tS}+ 12Hn{T tzP − T tzN + T tP − T tN}= StT(2.27)where the lower-diagonal (LD), diagonal (D) and upper-diagonal (UD) elements of the coefficientmatrix Az are directly given by the coefficients of TtzS , TtzP and TtzN , respectively. The RHS vectorbz are all the explicit parts of (2.27)4. Hence, the tri-diagonal coefficients in Az and the RHSvector elements in bz are given byALDz = −12Hs (2.28)ADz = H0z +12Hs +12Hn (2.29)AUDz = −12Hn (2.30)bz =H0zTtP− 12Hu{T txP − T txW + T tP − T tW}+ 12Hw{T txW − T txP + T tW − T tP}+ 12He{T txE − T txP + T tE − T tP}+ 12Hf{TtyF − T tyP + T tF − T tP}+ 12Hb{TtyB − T tyP + T tB − T tP}+ 12Hs {T tS − T tP }+ 12Hn {T tN − T tP }+StT(2.31)4parts evaluated at t, tx and ty82.4 ResultsThe 3D-ADI-DG model has been used for a thick glass fiber-polyurethane composite beam witha square cross-section [5]. The steady state temperature and degree of cure predictions at chosencross-section positions are shown in Fig. 2.4. The reader should notice how the outer parts ofthe profile (corner and facecenter) is heated and cured predominantly in the first half of the diewhile curing of the inner parts (quarter diagonal and center) are taking place over the full dielength.Figure 2.4: Numerical steady state temperature and degree of cure predictions inside a glass fiber-polyurethane pultruded beam during manufacturing.3 ConclusionThe current report has presented a detailed guide on how to obtain a numerical solution to a3-dimensional heat transfer problem using the ADI-DG method (considering convection solelyin one of the major spatial directions). A 3-point central difference scheme and a 2-point upwinddifference scheme were used to discretized the diffusive and the convective heat fluxes, respec-tively. The proposed ADI-DG solver was applied for the first time for pultrusion by Baran et al.[2] and later implemented, validated and used by the authors in [4] and [5].9Bibliography[1] J. Douglas Jr., J.E. Gunn: A General Formulation of Alternating Direction Methods - PartI. Parabolic and Hyberbolic Problems. Numerische Mathematik, 6:428-453, 1964.[2] I. Baran, J.H. Hattel, Cem C. Tutum: Thermo-Chemical Modelling Strategies for the Pul-trusion Process. Appl. Compo. Mater., 20:1247-1263, 2013.[3] J.H. Hattel: ”Fundamental of Numerical Modelling of Casting Processes”, 1. edition 2005,Polyteknisk Forlag Kgs. Lyngby, ISBN 13: 978-87-502-0969-0, ISBN 10: 87-502-0969-8.[4] F.S. Rasmussen, M.R. Sonne, M. Larsen, J. Spangeberg, L.T. Lilleheden and J.H.Hattel:2D Numerical Modelling of thge Resin Injection Pultrusion Process Including ExperimentalResin Kinetics and Temperature Validation. Proceeding: ICCM21, 2017.[5] F.S. Rasmussen, C.G. Klingaa, M. Larsen, M.R. Sonne, J. Spangenberg and J.H. Hattel:Modelling the Effect of Non-Uniform Fibre Distribution on the Curing Behaviour in ResinInjection Pultrusion. Proceedings: ECCM18, 2018.10

Numerical Modelling of Heat Transfer using the 3D-ADI-DG Method - with Application for Pultrusion.:Technical report

http://orbit.dtu.dk/ws/files/171228339/F.S_Rasmussen_Numerical_Modelling_of_Heat_Transfer_using_the_3D_ADI_DG_Method_with_Application_for_Pultrusion.pdf

Numerical Modelling of Heat Transfer using the 3D-ADI-DG Method - with Application for Pultrusion.:Technical report

Abstract

Similar works

Full text

Available Versions

Online Research Database In Technology