A Thermodynamic Framework for Additive Manufacturing, using Amorphous Polymers, Capable of Predicting Residual Stress, Warpage and Shrinkage
aa r X i v : . [ phy s i c s . c l a ss - ph ] F e b A Thermodynamic Framework for AdditiveManufacturing, using Amorphous Polymers, Capable ofPredicting Residual Stress, Warpage and Shrinkage
P Sreejith a , K Kannan a , K R Rajagopal b a Department of Mechanical Engineering, Indian Institute of Technology Madras,Chennai 600036, India b Department of Mechanical Engineering, Texas A&M University, College Station, TX77804, US bstract A thermodynamic framework has been developed for a class of amorphouspolymers used in fused deposition modeling (FDM), in order to predict theresidual stresses and the accompanying distortion of the geometry of theprinted part (warping). When a polymeric melt is cooled, the inhomogeneousdistribution of temperature causes spatially varying volumetric shrinkage re-sulting in the generation of residual stresses. Shrinkage is incorporated intothe framework by introducing an isotropic volumetric expansion/contractionin the kinematics of the body. We show that the parameter for shrinkage alsoappears in the systematically derived rate-type constitutive relation for thestress. The solidification of the melt around the glass transition temperatureis emulated by drastically increasing the viscosity of the melt.In order to illustrate the usefulness and efficacy of the derived constitutiverelation, we consider four ribbons of polymeric melt stacked on each othersuch as those extruded using a flat nozzle: each layer laid instantaneouslyand allowed to cool for one second before another layer is laid on it. Eachlayer cools, shrinks and warps until a new layer is laid, at which time theheat from the newly laid layer flows and heats up the bottom layers. Theresidual stresses of the existing and newly laid layers readjust to satisfy equi-librium. Such mechanical and thermal interactions amongst layers result ina complex distribution of residual stresses. The plane strain approximationpredicts nearly equibiaxial tensile stress conditions in the core region of thesolidified part, implying that a preexisting crack in that region is likely topropagate and cause failure of the part during service. The free-end of theinterface between the first and the second layer is subjected to the largest2agnitude of combined shear and tension in the plane with a propensity fordelamination.
Keywords:
Additive manufacturing, Residual stress, Shrinkage, Warpage,Flat nozzle, Helmholtz free energy, Rate of entropy production, Meannormal stress, Magnitude of stress, Nature of stress state
1. Introduction
Manufacturing processes such as fused deposition modelling (FDM), ex-trusion molding, injection molding and blow molding induce residual stressesin the thermoplastic components. Determining these stresses is crucial, sinceit alters the mechanical response of the final solidified component. Followinginstances demonstrate the same: (1) Plastic pipes manufactured by extru-sion, develop compressive residual stresses on the outer surfaces and tensileresidual stresses on the inner surfaces. Creep tests carried out to studycrack propagation on v-notched specimens cutout from these plastic pipes,confirmed that compressive stresses hindered crack growth whereas tensilestresses aided it (see Chaoui et. al. [1]). (2) Similarly, quenching of ther-moplastics induce compressive residual stresses in the material. Three pointbending tests and Izod impact tests, conducted to determine the fatigue lifeand impact strength of the quenched specimens, respectively, imply that,the strength increases with the magnitude of the residual stresses (see Horn-berger et. al. [2] and So et. al. [3]).In this paper, we attempt to develop a theory to capture the thermo-dynamic process associated with FDM during the layer by layer fabrication3f a polymeric component. We aim to achieve the following: (1) Develop athermodynamic framework that can represent the FDM process. (2) Use thedeveloped model to: (a) Determine the residual stresses induced in the finalsolidified component due to rapid cooling and (b) capture the dimensionalchanges (shrinking and warping) of the part, which manifests as a conse-quence of the residual stresses.In recent years, FDM has been proving to be useful in the fabricationof patient specific transplants, tissue scaffolds and pharmaceutical products(see Melchels et. al. [4]; Melocchi et. al. [5]). The use of this technologyin the medical industry has been growing fast due to the flexibility it of-fers in fashioning products, such as tools for surgical planning and medicalimplants. For instance, it has replaced solvent casting, fiber-bonding, mem-brane lamination, melt molding, and gas forming, which were previouslyused to fabricate medical grade polymethylmethacrylate (PMMA) implantsfor craniofacial reconstruction and regeneration surgeries (see Espanil et.al.[6]; Nyberg et. al. [7]). One of the most important advantages of FDM isthe possibility to architect the microstructure of the part. This flexibilitywas demonstrated by Muller et. al. [8], wherein they fabricated an implantfor a frontal-parietal defect with branching vascular channels to promote thegrowth of blood vessels.Stress-shielding effect caused by metal implants have led researchers toturn to polymers such as polyetheretherketone (PEEK) for orthopedic im-plants (see Han et.al. [9]). PEEK displays mechanical behaviour similar tothat of the cortical bone (see Vaezi and Yang [10] ), due to the presenceof ether and ketone groups along the backbone of it’s molecular structure.4here have been surgeries such as cervical laminoplasty (see Syuhada et. al.[11] ), total knee replacement surgery(see Swathi and Devadath [12] ; Borgeset. al. [13] ), patient specific upper limb prostheses attachment (see Kateet. al. [14] ) and patient specific tracheal implants (see Freitag et. al. [15])recorded in the literature, which use PEEK, polylactic acid (PLA), acryloni-trile butadiene styrene (ABS), polyurethane (PU) and other medical gradepolymer implants that have been fabricated using FDM. Recently, additivemanufacturing has received a huge impetus in the healthcare sector due tothe ongoing Covid-19 pandemic, and the technology is being used to fabri-cate products such as ventilator valves, face shields, etc.During fabrication of such components, thermoplastics are heated abovethe melting point ( θ m ) and then rapidly cooled below the glass transition tem-perature ( θ g ), which causes drastic temperature gradients. Cooling reducesconfigurational entropy of the polymer molecules, and the sharp tempera-ture gradients induce differential volumetric shrinkage, consequently leadingto the formation of residual stresses. As previously stated, these stresses canbe either detrimental or favourable with regard to the mechanical responseof the components.Process parameters such as layer thickness, orientation, raster angle,raster width, air gap, etc. affect the component’s final mechanical response(see Sood et. al. [16]). For instance, the anistropic nature of a part fab-ricated by FDM (see Song et. al. [17]) can be controlled by choosing themost suitable raster angle for each individual layer during the process, insuch a way that the component will survive under the required service con-ditions. This has been demonstrated by Casavola et. al. [18], wherein a5ectangular ABS specimen with a raster angle of ± ◦ produced the leastresidual stresses in the specimen. The optimal raster angle was determinedby repeated experimentation, with a different choice of the angle in each ex-periment. Developing a consistent theory for capturing the thermodynamicprocess and using the resulting constitutive relations in simulations can helpin developing standardized methods with optimized parameters for manufac-turing the component.Many of the theories currently available, concentrate on the analysis ofpolymers which are at a temperature below θ g , the reasoning being that, in-ternal reaction forces due to volume shrinkage are too low in the melt phaseto give rise to appreciable residual stresses (see Xinhua et. al. [19]; Wang et.al. [20]; Park et. al. [21], Macedo et. al. [22]). Although such an argumentseems reasonable, in a process like FDM, the mechanical and thermal histo-ries of the polymer melt have a huge impact on the mechanical response ofthe final component.Several computational frameworks have been developed, wherein somehave assumed incompressible viscous fluid models for the melt phase (seeXia and Lu et. al. [23]; Dabiri et. al. [24], Comminal et. al. [25]) andlinear elastic models (see Xinhua et. al. [19]; Wang et. al. [20]; Macedo et.al. [22]; Moumen et. al. [26]; Casavola et. al. [27]) or elasto-plastic models(see Armillotta et. al. [28]; Cattenone et. al. [29] ) for the solid phase.However, it is to be noted that polymer melts are generally modelled as vis-coelastic fluids (see Doi and Edwards [30]). This is supported by the factthat, residual stresses in injection molded components have been predictedmuch more accurately by viscoelastic models than by elastic or elasto-plastic6odels (see Kamal et. al. [31] and Zoetelief et. al. [32]). In this paper,the polymer melts are assumed to be amorphous, and we make an entropicassumption to simplify the analysis. Hence, the internal energy is solely afunction of temperature and therefore, at a constant temperature, the freeenergy accumulates only due to the reducing configurational entropy of themolecules, such as during stretching. When the load is removed, the meltagain assumes a high entropy state, which is also the stress free configurationof the melt.Once the amorphous melt has been laid on the substrate, it undergoesphase transformation due to cooling. Numerous “ad hoc” methods havebeen implemented by various authors to capture the phase change. One ofthe simplest methods used was to assume the melt to transform into a solidat a material point where the temperature reaches θ m (see Xia and Lu et.al. [33]) or θ g (see Macedo et. al. [22], see Armillotta et. al. [28]; Moumenet. al. [26] and Cattenone et. al. [29]). Generally, polymer melts undergophase transition through a range of temperature values. The heat transferbetween a polymer melt and the surroundings is usually proportional to thetemperature difference between them, and solidification is initiated at thematerial points where the germ nuclei are activated into a growth nuclei.An amorphous polymer melt normally transforms into an amorphous solid,unless the melt consists of specific germ nuclei which can cause crystalliza-tion. In the case of FDM, such germ nuclei can be induced in the melt whileextruding the melt out of the nozzle at high shear rates combined with a lowtemperature (see Northcutt et. al. [34]).In the transition regime, polymer molecules exist in both, the melt, and7he solid phases simultaneously, and usually the body is modelled as a con-strained mixture (see Rao and Rajagopal [35]; Rao and Rajagopal [36]; Kan-nan et. al. [37]; Kannan and Rajagopal [38]; and Kannan and Rajagopal[39]). While considering such a complex phase transition process, it is muchmore prudent to assume the primary thermodynamic functions to be depen-dent on the weight fraction of the phases. The weight fraction can be easilyrepresented as a linear function of temperature (see Liu et. al. [40]) or it canbe determined by statistical methods (see Avrami et. al. [41]).Considering this work to be a preliminary step in developing a consis-tent theory for FDM, we assume an amorphous polymer melt to transforminto a solid through a drastic increase in the viscosity at the glass transi-tion temperature ( θ g ). In such a phase change process, the melt transformscompletely into a solid when the temperature at all the material points fallbelow the glass transition value ( θ g ). This assumption limits the applicabil-ity of our theory to very particular processes having certain specific processparameters, such as:(1) The temperature of the nozzle is maintained above θ m throughout. (2) The intial temperature of the melt is also kept above θ m .These predefined conditions ensure that the melt does not undergo shearinduced crystallization. The proposed framework is based on the generaltheory for viscoelastic fluids developed by Rajagopal and Srinivasa [42] andon the theory of the crystallization of polymer melts developed by Rao andRajagopal [35] [43] [36]. The framework requires an appropriate assumptionof the Helmholtz potential and the rate of entropy production, and the evolu-tion of the “natural configuration” of the body is determined by maximizingthe rate of entropy production. 8 . Methodology Figure 1:
The various configurations of the melt.
Let us consider the melt to be subjected to thermal loading. Initially theamorphous polymer melt (atactic polymer) is at the reference configuration κ R (B) (refer to Fig.1). This configuration is stress free and the polymermelt is at the reference temperature θ o . As the melt is cooled from θ o , itundergoes deformation and occupies the current configuration κ t (B), at thecurrent temperature θ . The motion of the body from κ R (B) to κ t (B) is givenby x = χ κ R ( X , t ) . (1)The motion is assumed to be a diffeomorphism. The velocity and the defor-mation gradient are given by v = ∂χ κ R ∂t , (2)and F κ R = ∂χ κ R ∂ X . (3)9he deformation gradient mapping κ pm ( t ) (B) to κ t (B) is F κ pm ( t ) . Let α m ( θ ) I , where α m ( θ ) is a scalar valued funtion, represent pure thermal ex-pansion/contraction. When α m ( θ ) = 1, i.e., an isothermal process, Fig.1 willcoincide with the configurations employed by Rajagopal and Srinivasa [42].The isothermally unloaded configuration κ pm ( t ) (B) evolves due to thermal ex-pansion/contraction of the melt (see Rajagopal and Srinivasa [44]). The con-figuration κ ′ pm ( t ) is stress free, and since pure thermal expansion/contractiondoesn’t induce stresses in the material, the configuration κ pm ( t ) is also as-sumed to be stress free. The deformation gradient is multiplicatively decom-posed into a dissipative, thermal and an elastic part as shown F κ R = α m ( θ ) F κ pm ( t ) G m . (4)In general, G m need not be the gradient of a mapping.The right Cauchy-Green elastic stretch tensor and the left Cauchy-Greenelastic stretch tensor are given by C κ R = F Tκ R F κ R , and C κ pm ( t ) = F Tκ pm ( t ) F κ pm ( t ) , (5) B κ R = F κ R F Tκ R , and B κ pm ( t ) = F κ pm ( t ) F Tκ pm ( t ) . (6)The ‘right configurational stretch tensor’ associated with the reference config-uration ( κ R ) and the instantaneous natural configuration ( κ ′ pm ( t ) ) is definedas C κ R → κ ′ pm ( t ) = G Tm G m . (7)Therefore, Eq.(6) can be represented as B κ pm ( t ) = 1 α m ( θ ) F κ R C − κ R → κ ′ pm ( t ) F Tκ R . (8)10he velocity gradients and the corresponding symmetric and skew-symmetricparts are defined as L = ˙ F κ R F − κ R , (9) D = L + L T , W = L − L T , (10) L κ ′ pm ( t ) = ˙ G m G m − , (11) D κ ′ pm ( t ) = L κ ′ pm ( t ) + L Tκ ′ pm ( t ) , W κ ′ pm ( t ) = L κ ′ pm ( t ) − L Tκ ′ pm ( t ) . (12)The material time derivative of Eq.(8) will lead to˙ B κ pm ( t ) = − α m ( θ ) α m ( θ ) B κ pm ( t ) + LB κ pm ( t ) + B κ pm ( t ) L T + 1 α m ( θ ) F κ R ˙ C − κ R → κ ′ pm ( t ) F Tκ R . (13)The upper convected derivative is used to ensure frame invariance of the timederivatives, which, for any general second order tensor A is defined as ▽ A = ˙ A − LA − AL T . (14)On comparing Eq.(13) and Eq.(14), we arrive at ▽ B κ pm ( t ) = − α m ( θ ) α m ( θ ) B κ pm ( t ) + 1 α m ( θ ) F κ R ˙ C − κ R → κ ′ pm ( t ) F Tκ R . (15)Using Eq.(4), Eq.(7) and Eq.(15), we obtain ▽ B κ pm ( t ) = − F κ pm ( t ) D κ ′ pm ( t ) F Tκ pm ( t ) − α m ( θ ) α m ( θ ) B κ pm ( t ) . (16)11ote that, the evolution of B κ pm ( t ) would depend on the structure of theconstitutive equation chosen for α m ( θ ) (or more generally, on the material).In the event of an isothermal process, the scalar valued function α m ( θ ), willbecome unity, and Eq.(16) will reduce to the form derived by Rajagopal andSrinivasa [42] as given below ▽ B κ pm ( t ) = − F κ pm ( t ) D κ ′ pm ( t ) F Tκ pm ( t ) . (17)To capture the effects of volume relaxation we will be using normalisedinvariants, where unimodular tensors are defined such that − C κ pm ( t ) = det( F κ pm ( t ) ) − C κ pm ( t ) , and − B κ pm ( t ) = det( F κ pm ( t ) ) − B κ pm ( t ) . (18)The invariants that we employ are defined through I κ pm ( t ) = tr( − B κ pm ( t ) ) , II κ pm ( t ) = tr( − B κ pm ( t ) ) , and III κ pm ( t ) = det( F κ pm ( t ) ) . (19)However, we plot Lode invariants (see Chen et. al. [45]) in section 3, tointerpret the results better. This set of invariants are defined asK = tr( A ) √ , K = q tr( A d ) , and K = 13 sin − (cid:16) √ A d )tr( A d ) (cid:17) , (20)where A is any tensor and A d is the deviatoric part of A . When K = 0,the tensor is purely composed of the deviatoric components. Similarly, whenK = 0, the tensor is volumetric and K represents the mode of the tensorwhich varies from − π to π . 12 .2. Modeling2.2.1. Modeling the melt Helmholtz free energy per unit mass is defined with respect to the config-uration κ pm ( t ) (B), and is assumed to be a function of θ and the deformationgradient F κ pm ( t ) . We represent it as a sum of the contribution from the ther-mal interactions and the mechanical working, where the contribution due tomechanical working is assumed to be of the Neo-Hookean form. Further, werequire it to be frame indifferent and isotropic with respect to the configura-tion κ pm ( t ) (B), to arrive atΨ m ( θ, I κ pm ( t ) , III κ pm ( t ) ) = Ψ thm ( θ ) + Ψ mechm (I κ pm ( t ) , III κ pm ( t ) ) , (21)whereΨ thm ( θ ) = A m + (B m + C m )( θ − θ o ) − C m ( θ − θ o ) − C m θ ln (cid:18) θθ o (cid:19) , (22)andΨ mechm ( I κ pm ( t ) , III κ pm ( t ) ) = µ m θ ρ κ pm ( t ) θ o h tr( − B κ pm ( t ) ) − i + k m θ ρ κ pm ( t ) θ o h det( F κ pm ( t ) ) − i , (23)where θ o , A m , B m , C m , C m , µ m and k m are the reference material tempera-ture and the material constants respectively, and µ m , k m ≥
0, and ρ κ pm ( t ) isthe density of the melt in the configuration κ pm ( t ) (B).The rate of entropy production per unit volume due to mechanical work-ing is defined as a function of θ and L κ ′ pm ( t ) , and requiring the function to beframe indifferent and isotropic, we assume ξ mechm ( θ, D κ ′ pm ( t ) ) = η m ( θ )dev( D κ ′ pm ( t ) ) . dev( D κ ′ pm ( t ) ) + η m ( θ )tr( D κ ′ pm ( t ) ) , (24)13here η m and η m are the shear and bulk modulus of viscosity and η m , η m ≥ ξ mechm ≥ , (25)and the form given by Eq.(24), ensures that Eq.(25) is always met.The balance of energy can be written as ρ κ t ˙ ǫ m = T.D − div( q ) + ρ κ t r, (26)where T is the Cauchy stress tensor, ˙ ǫ m is the rate of change of internal energyper unit mass of the melt with respect to time, q is the heat flux and r is theradiation of heat per unit mass of the melt. We express the second law ofthermodynamics as an equality by introducing a rate of entropy productionterm, and it takes the form (see Green and Naghdi [46]) T.D − ρ κ t ˙Ψ m − ρ κ t η m ˙ θ − θ q . grad( θ ) = ρ κ t θζ m = ξ totm , ξ totm ≥ , (27)where η m is the entropy per unit mass of the melt, θ is the absolute temper-ature, ζ m is the rate of entropy production per unit mass of the melt and ξ totm is the total rate of entropy production per unit volume. Substituting Eq.(21)into Eq.(27), we can re-write Eq.(27) as T.D − ρ κ t (cid:26) ∂ Ψ thm ∂θ + µ m ρ κ pm ( t ) θ o h tr( − B κ pm ( t ) ) − i + k m ρ κ pm ( t ) θ o h det( F κ pm ( t ) ) − i + η m (cid:27) ˙ θ − ρ κ t (cid:26) µ m θ ρ κ pm ( t ) θ o tr( ˙ − B κ pm ( t ) ) − µ m θ ˙ ρ κ pm ( t ) ρ κ pm ( t ) θ o h tr( − B κ pm ( t ) ) − i + k m θρ κ pm ( t ) θ o h det( F κ pm ( t ) ) − i ˙det( F κ pm ( t ) ) − k m θ ˙ ρ κ pm ( t ) ρ κ pm ( t ) θ o h det( F κ pm ( t ) ) − i (cid:27) − θ q . grad( θ ) = ξ totm , ξ totm ≥ . (28)14he balance of mass for the melt undergoing deformation between theconfigurations κ R (B) and κ pm ( t ) (B) is ρ κ R = ρ κ pm ( t ) det( G m )( α m ( θ )) , (29)where ρ κ R is the density at the reference configuration, κ R (B). The balanceof mass for the melt between the configurations κ R (B) and κ t (B) is given by ρ κ R = ρ κ t det( F κ R ) , (30)where ρ κ t is the density at the current configuration, κ t (B). Taking the ma-terial derivative of Eq.(29), we arrive at˙ ρ κ pm ( t ) ρ κ pm ( t ) = − tr( D κ ′ pm ( t ) ) − α m ( θ ) α m ( θ ) . (31)Next, from Eq.(6) and Eq.(16)˙det( F κ pm ( t ) ) = det( F κ pm ( t ) ) (cid:26) tr( D − D κ ′ pm ( t ) ) − α m α m (cid:27) , (32)and from Eq.(16), Eq.(18) and Eq.(32), we obtaintr( ˙ − B κ pm ( t ) ) = 2dev( − B κ pm ( t ) ) . ( D − D κ ′ pm ( t ) ) (33)It follows from Eq.(4) and Eq.(28)-(32), that (cid:26) T − µ m θ det( F κ pm ( t ) ) θ o dev( − B κ pm ( t ) ) − k m θθ o h det( F κ pm ( t ) ) − i I (cid:27) . D − ρ κ t (cid:26) ∂ Ψ thm ∂θ + µ m ρ κ pm ( t ) θ o h tr( − B κ pm ( t ) ) − i + k m ρ κ pm ( t ) θ o h det( F κ pm ( t ) ) − i + η m (cid:27) ˙ θ + µ m θ det( F κ pm ( t ) ) θ o dev( − B κ pm ( t ) ) . D κ ′ pm ( t ) + (cid:26) k m θθ o h det( F κ pm ( t ) ) − i − µ m θ F κ pm ( t ) ) θ o h tr( − B κ pm ( t ) ) − i − k m θ F κ pm ( t ) ) θ o h det( F κ pm ( t ) ) − i (cid:27) tr( D κ ′ pm ( t ) ) + 3 ˙ α m ( θ ) α m ( θ ) (cid:21) − θ q . grad( θ ) = ξ totm , ξ totm ≥ . (34)One of the ways to satisfy Eq.(34) is to assume that the three equations givenbelow hold η m = − ∂ Ψ thm ∂θ − µ m ρ κ pm ( t ) θ o h tr( − B κ pm ( t ) ) − i − k m ρ κ pm ( t ) θ o h det( F κ pm ( t ) ) − i , (35) T = µ m θ det( F κ pm ( t ) ) θ o dev( − B κ pm ( t ) ) + k m θθ o h det( F κ pm ( t ) ) − i I , (36)and ξ totm = T . D κ ′ pm ( t ) − Ψ MECHm I . D κ ′ pm ( t ) + 3 ˙ α m ( θ ) α m ( θ ) (cid:26) − µ m θ θ o det( F κ pm ( t ) ) h tr( − B κ pm ( t ) ) − i + k m θ θ o det( F κ pm ( t ) ) h(cid:0) det( F κ pm ( t ) ) (cid:1) − i(cid:27) − θ q . grad( θ ) , ξ totm ≥ , (37)where Ψ MECHm = ρ κ t Ψ mechm . (38)The Helmholtz free energy is given byΨ m = ǫ m − θη m . (39)From Eq.(21) and Eq.(35), the internal energy is defined as ǫ m = A m − B m θ o + 12 C m ( θ − θ o ) + C m ( θ − θ o ) , (40)and the specific heat is ∂ǫ m ∂θ = C m θ + C m . (41)16hen the process is homothermal, Eq.(37) will reduce to ξ totm = T . D κ ′ pm ( t ) − Ψ MECHm I . D κ ′ pm ( t ) , ξ totm ≥ . (42)Therefore, in the current analysis, it seems prudent to suppose the following ξ mechm = T . D κ ′ pm ( t ) − Ψ MECHm I . D κ ′ pm ( t ) . (43)Eq.(43) can be re-written as ξ mechm ( D κ ′ pm ( t ) , θ ) = ξ mechm − T . D κ ′ pm ( t ) + Ψ MECHm I . D κ ′ pm ( t ) = 0 . (44)Now, we assume the natural configuration (i.e. D κ ′ pm ( t ) ) to evolve in such away that, for a fixed B κ ′ pm ( t ) and θ at each time instant, the rate of entropyproduction due to mechanical working is maximum, subject to the constraintin Eq.(44) (see Rajagopal and Srinivasa [47]). Therefore, to extremize ξ mechm constrained by Eq.(44), we write the augmented form asΦ aug = ξ mechm + λ ( ξ mechm − T . D κ ′ pm ( t ) + Ψ MECHm I . D κ ′ pm ( t ) ) , (45)where λ is the Lagrange multiplier. On taking the derivative with respect to D κ ′ pm ( t ) , we obtain ∂ Φ aug ∂ D κ ′ pm ( t ) = (1 + λ ) λ ∂ξ mechm ∂ D κ ′ pm ( t ) + (Ψ MECHm I − T ) = 0 . (46)At any time instant, assuming B κ pm ( t ) and θ to be fixed, D κ ′ pm ( t ) will take onlythose values which will satisfy Eq.(46). Substituting Eq.(24) into Eq.(46) andtaking the scalar product with D κ ′ pm ( t ) Φ aug ∂ D κ ′ pm ( t ) . D κ ′ pm ( t ) = 2(1 + λ ) λ (cid:20) η m ( θ )dev( D κ ′ pm ( t ) ) . D κ ′ pm ( t ) + η m ( θ )[tr( D κ ′ pm ( t ) )] (cid:21) + h Ψ MECHm tr( D κ ′ pm ( t ) ) − T . D κ ′ pm ( t ) i = 0 . (47)Comparing Eq.(47) with Eq.(44), we get (1+ λ ) λ = . Since the melt is assumedto be isotropic, we can assume that F κ pm ( t ) = V κ pm ( t ) , (48)Also the eigenvectors of B κ pm ( t ) and D κ ′ pm ( t ) are the same and hence thetensors commute (see Rajagopal and Srinivasa [42]). Therefore, from Eq.(16),Eq.(46) and Eq.(48), we arrive at the form for the upper convected derivativeof B κ pm ( t ) , which will determine the evolution of natural configuration of themelt, as ▽ B κ pm ( t ) = 2 (cid:18) Ψ MECHm − tr( T )3 η m ( θ ) (cid:19) I − η m ( θ ) dev( T ) − ˙ α m ( θ ) α m ( θ ) I ! B κ pm ( t ) . (49)
3. Application of the model
A prototypical problem to test the efficacy of the theory that has beendeveloped is considered, i.e., we assume a particular geometry to be fashionedby FDM and use the constitutive equations to determine the residual stressesinduced in the material during fabrication and the consequent dimensionalinstability of the geometry.
In FDM, a polymer filament is heated above the melting temperature andextruded out through a nozzle attached to a robotic arm. A stereolithogra-phy (STL) file, containing information of the 3D geometry to be made, is18hen fed into a processing system that controls the robotic arm. The poly-mer melt is laid layer by layer through the nozzle with the help of the arm toget the desired 3D geometry, and the melt is allowed to cool. Polymers usedin this process are usually thermoplastics like ABS (Acrylonitrile ButadieneStyrene), PLA (Polylactic Acid), PS (Polystyrene), PC (Polycarbonates),PEEK (Polyetheretherketone), PMMA (Polymethylmethacrylate) and elas-tomers (see Gonz´alez-Henr´ıquez et. al. [48]). Thermoplastics are preferredover thermosets, because the former’s melt viscosity enables smooth extru-sion through the nozzle, and at the same time helps to retain the shape oncethe melt is laid. Complex part geometries can also be made with the help ofsupport materials.In the current analysis, we assume each of the layers to be very long andbe made of a single raster. Further, we assume a rectangular cross section forthe layers, which finally adds up to each of the layers being in the shape of avery long ribbon. To lay such a layer, we need a nozzle with an appropriategeometry, such as a nozzle with a slot (see Loffer et. al. [49]) or a flat headnozzle which is capable of evening out the melt into the desired shape (seeKim et. al. [50]). Four layers of the polymer melt are assumed to be laid oneover the other consecutively, in the direction of the length of the ribbon. Allfour layers are assumed to have a sufficiently long rectangular cross-sections,which enables plane strain conditions to be enforced (refer to Fig.2).
The constitutive form for the scalar valued function α m ( θ ), which rep-resents the isotropic volume expansion/contraction, can be deduced from19q.(4). When t >
0, taking the determinant of Eq.(4) α m ( θ ) = ρ κ R ρ κ t det( F κ pm ( t ) )det( G m ) ! . (50)Dilatometry experiments that are conducted to measure the specific volumechanges are homothermal processes, and hence, the deformation gradients, G m and F κ pm ( t ) , shown in Fig.1, can be assumed to be identity transforma-tions. Thus, density will become solely a function of temperature and theisotropic volume expansion/contraction can be represented as α m ( θ ) = ρ κ R ρ κ t ! . (51)The density can be expressed in terms of specific volume, and therefore, α m ( θ ) is obtained as α m ( θ ) = v ( θ ) v ( θ o ) ! = vv o ! , (52)where v is the specific volume at the current temperature θ and v o is thespecific volume at the reference temperature θ o which can be represented bythe Tait equation (see Haynes [51]). At time t = 0 second, i.e. at referencetemperature (temperature of the reference configuration) θ = θ o α m ( θ o ) = 1 . (53)The coefficient of volumetric thermal expansion/contraction, ˆ α m ( θ ), is de-fined as ˆ α m ( θ ) = 1 v (cid:18) ∂v∂θ (cid:19) = ρ κ t (cid:18) ∂ ( ρ κt ) ∂θ (cid:19) . (54)20he relationship between isotropic volume expansion/contraction α m ( θ ) andthe coefficient of volumetric thermal expansion/contraction ˆ α m ( θ ) is obtainedby taking the derivative of Eq.(52) with respect to the current temperature, θ , as ˆ α m ( θ ) = 3 α m ( θ ) ∂α m ( θ ) ∂θ . (55)As the layers cool, the viscosity of the polymer melt shoots up at glasstransition temperature ( θ g ). The variation of the bulk and shear viscositiesof the melt are defined as η mn ( θ ) = η sn (cid:18) − tanh (cid:0) a ( θ − θ g ) (cid:1)(cid:19) η ln (cid:18) (cid:0) a ( θ − θ g ) (cid:1)(cid:19) , n = 1 , , (56)where η m is the shear viscosity and η m is the bulk viscosity, η s and η s are theshear and the bulk viscosities below the glass transition temperature ( θ g ), η l and η l are the shear and the bulk viscosities above the glass transition tem-perature ( θ g ) and a is a constant. The shear viscosity and the bulk viscosityabove the glass transition, i.e., η l and η l , are defined as η ln ( θ ) = η on e (cid:18) C (cid:0) θ − θo (cid:1)(cid:19) , n = 1 , , (57)where η o is the shear viscosity at the melting temperature ( θ m ), η o is thebulk viscosity at the melting temperature ( θ m ) and C is a constant.The dependent variables are assumed to be a function of the current co-ordinates and current time. They are the displacement vector ∼ u ( x , t ), thecurrent temperature ∼ θ ( x , t ) and the left Cauchy-Green stretch tensor as de-fined in Eq.(6) , i.e., ∼ B κ pm ( t ) ( x , t ). The bases in the local and the global21oordinate system are aligned with each other. Therefore, defining the de-pendent variables with respect to the local coordinate system is as good asdefining it with respect to the global co-ordinate system. The position vectorof a particle in the reference configuration is defined as X = x ( t ) − ∼ u ( x , t ) . (58)The velocity vector in the current configuration is derived by taking the totaltime derivative of Eq.(58) ∼ v ( x , t ) = (cid:16) I − ∂ ( ∼ u ( x , t )) ∂ x (cid:17) − ∂ ∼ u ( x , t ) ∂t , (59)and the velocity gradient is defined as ∼ L ( x , t ) = ∂ ∼ v ( x , t ) ∂ x . (60)The components of the displacement vector are (cid:2) ∼ u (cid:3) = u x ( x, y, t ) u y ( x, y, t ) , (61)where (cid:2) ∼ u (cid:3) represents the component form of the displacement vector and, u x ( x, y, t ) and u y ( x, y, t ) are the components in the x and y directions re-spectively. We represent the components of the velocity vector as givenbelow (cid:2) ∼ v (cid:3) = v x ( x, y, t ) v y ( x, y, t ) , (62)where v x = ∂u x ∂t (cid:0) ∂u y ∂y − (cid:1) − ∂u x ∂y ∂u y ∂t∂u x ∂x + ∂u y ∂y − ∂u x ∂x ∂u y ∂y + ∂u x ∂y ∂u y ∂x − , (63)22nd v y = ∂u y ∂t (cid:0) ∂u x ∂x − (cid:1) − ∂u x ∂t ∂u y ∂x∂u x ∂x + ∂u y ∂y − ∂u x ∂x ∂u y ∂y + ∂u x ∂y ∂u y ∂x − . (64)Consequently, the component form of the velocity gradient follows fromEq.(60) as (cid:2) ∼ L (cid:3) = ∂v x ( x,y,t ) ∂x ∂v x ( x,y,t ) ∂y ∂v y ( x,y,t ) ∂x ∂v y ( x,y,t ) ∂y
00 0 0 . (65)We represent the components of ∼ B κ pm ( t ) ( x , t ) and the stress tensor T (cid:16) ∼ θ, ∼ B κ pm ( t ) (cid:17) as (cid:2) ∼ B κ pm ( t ) (cid:3) = B xx ( x, y, t ) B xy ( x, y, t ) 0B xy ( x, y, t ) B yy ( x, y, t ) 00 0 B zz ( x, y, t ) , (66)and (cid:2) T (cid:3) = T xx ( x, y, t ) T xy ( x, y, t ) 0T xy ( x, y, t ) T yy ( x, y, t ) 00 0 T zz ( x, y, t ) , (67)where T xx ( x, y, t ), T yy ( x, y, t ), T xy ( x, y, t ) and T zz ( x, y, t ) are derived by sub-stituting θ ( x, y, t ) and the components of ∼ B κ pm ( t ) ( x , t ) into Eq.(36).There are seven coupled field equations that have to be solved for, i.e., twocomponents of the linear momentum equation, the energy equation [Eq.(26)]and four components of the evolution equation [Eq.(49)] represented in theEulerian form. The linear momentum balance is given by ∂ T xx ∂x + ∂ T xy ∂y = 0 , (68)23nd ∂ T xy ∂x + ∂ T yy ∂y = 0 . (69)We assume radiation to be absent in the current analysis and hence theenergy equation is written asT xx ∂v x ∂x + T yy ∂v y ∂y + T xy (cid:18) ∂v x ∂y + ∂v y ∂x (cid:19) + k (cid:16) ∂ θ∂x + ∂ θ∂y (cid:17) = ρ κ t ∂ǫ ( θ ) ∂θ (cid:16) ∂θ∂t + ∂θ∂x v x + ∂θ∂y v y (cid:17) . (70)The evolution equations are ∂ B xx ∂t + ∂ B xx ∂x v x + ∂ B xx ∂y v y − (cid:18) ∂v x ∂x B xx + ∂v x ∂y B xy (cid:19) = 23 η m ( θ ) " Ψ MECHm − (T xx + T yy + T zz )3 B xx − η m ( θ ) "(cid:16) xx − T yy − T zz (cid:17) B xx + T xy B xy − α m ( θ ) ∂α m ( θ ) ∂θ (cid:16) ∂θ∂t + ∂θ∂x v x + ∂θ∂y v y (cid:17) B xx , (71) ∂ B yy ∂t + ∂ B yy ∂x v x + ∂ B yy ∂y v y − (cid:18) ∂v y ∂x B xy + ∂v y ∂y B yy (cid:19) = 23 η m ( θ ) " Ψ MECHm − (T xx + T yy + T zz )3 B yy − η m ( θ ) "(cid:16) − T xx yy − T zz (cid:17) B yy + T xy B xy − α m ( θ ) ∂α m ( θ ) ∂θ (cid:16) ∂θ∂t + ∂θ∂x v x + ∂θ∂y v y (cid:17) B yy , (72) ∂ B zz ∂t + ∂ B zz ∂x v x + ∂ B zz ∂y v y = 23 η m ( θ ) " Ψ MECHm − (T xx + T yy + T zz )3 B zz − η m ( θ ) "(cid:16) − T xx − T yy zz (cid:17) B zz − α m ( θ ) ∂α m ( θ ) ∂θ (cid:16) ∂θ∂t + ∂θ∂x v x + ∂θ∂y v y (cid:17) B zz , (73)24nd ∂ B xy ∂t + ∂ B xy ∂x v x + ∂ B xy ∂y v y − (cid:18) ∂v x ∂x B xy + ∂v x ∂y B yy (cid:19) − (cid:18) ∂v y ∂x B xx + ∂v y ∂y B xy (cid:19) =23 η m ( θ ) " Ψ MECHm − (T xx + T yy + T zz )3 B xy − η m ( θ ) "(cid:16) xx − T yy − T zz (cid:17) B xy + T xy B yy − α m ( θ ) ∂α m ( θ ) ∂θ (cid:16) ∂θ∂t + ∂θ∂x v x + ∂θ∂y v y (cid:17) B xy , (74)where Ψ MECHm in Eq.(71-74) is derived by substituting θ ( x, y, t ) and the com-ponents of ∼ B κ pm ( t ) ( x , t ) into Eq.(38). We solve the seven equations for theseven unknowns: two components of the displacement field ( ∼ u ( x , t )), thetemperature field ( ∼ θ ( x , t )) and four components of the left Cauchy-Greenstretch tensor ( ∼ B κ pm ( t ) ( x , t )). The polymer melt is assumed to be of Polystyrene (PS). High molecularweight atactic PS does not undergo crystallization (see Chai et. al. [52]),and hence is a good candidate for our analysis. The thermal and mechanicalproperties of PS are given in Table 1. The mechanical properties providedin Table 1 are associated with the κ pm ( t ) configuration. However, in theliterature, the properties are associated with the κ t configuration. Therefore,we have assumed the values to lie in a ballpark range of the actual values.Since the values of shear modulus and bulk modulus given in Table 1 aretowards the higher end at temperatures above θ g , it leads to a very smallrelaxation time ( η o µ m = 10 − sec) of the melt. However, this seems to be areasonable approximation. 25 able 1: Properties of PS Material Properties Value Unit Ref.Thermal properties
Melting temperature, θ m
513 K [53]Glass transition temperature, θ g
373 K [53]Coefficient of heatconduction, k (averaged value) 0 . WmK [54]Coefficient of heat transferfrom the material to the 90 Wm K [55]environment, h e Coefficient of heat transferfrom the material to the 100 Wm K [55]substrate, h s Specific heat ( C p ):(Table 1a and 1b of [56])C m JK kg [56]C m JKkg [56]Constants in the Tait equation:A cm g [51]A . − o C − [51] Mechanical properties
Shear modulus ( µ m ) 10 Pa [57]26ulk modulus ( k m ) 3x10 Pa NAShear viscosity below glass transitiontemperature, θ g ( η s ) 10 Pas NABulk viscosity below glass transitiontemperature, θ g ( η s ) 2x10 Pas NAShear viscosity at the referencetemperature, θ o ( η o ) 10 Pas [58]Bulk viscosity at the referencetemperature, θ o ( η o ) 2x10 Pas NAC 22873 K [58] a The first layer of the melt is laid instantaneously on the substrate (whichis assumed to be at the ambient temperature θ sub = 27 o C throughout the pro-cess) at time t = 0 second. The surface that is in contact with the substrateis fixed along with the substrate acting as a heat sink. Free convection andtraction free conditions are assumed on the outer surfaces (refer to Fig.2a).The layer is then allowed to cool for one second, at the end of which, thesecond layer is laid instantaneously on top of the first layer. Once the secondlayer is laid, continuity is established between the layers, thus both the layersact as a single contiguous body (refer to Fig.2b), i.e., we assume a perfectinterface between the layers. The third and fourth layers are laid similarlyone after the other with an interval of one second each (refer to Fig.2c and27 a) The first layer laid at time t = 0 second and cooled till time t = 1 second. (b)
The second layer laid on top of the first layer at time t = 1 second and both layers cooledtogether till time t = 2 second. c) The third layer laid on top of the second layer at time t = 2 second and all the three layerscooled together till time t = 3 second. (d)
The fourth layer laid on top of the third layer at time t = 3 second and all the four layerscooled together till time t = 60 second.
Figure 2: Representative images showing the boundary conditions on the layers atdifferent time steps of the simulation. The dashed (dark red) line provides a breakand is intended to denote that the layers are very long. o C). Symmetry conditions have been imposed to reduce the computationaltime (refer to Fig.2a, Fig.2b, Fig.2c and Fig.2d). The first, second, third andfourth layers have an aspect ratio of 80:1,76:1, 72:1 and 68:1 respectively.Before each layer is laid, the layer is at 241 o C, which is 1 o C above themelting point of Polystyrene (PS). At this temperature, the isotropic vol-ume contraction/expansion ( α m ( ∼ θ )) is unity along with ∼ B κ pm ( t ) ( x , t ) beingidentity, i.e., (cid:2) α m ( θ o ) I (cid:3) = (cid:2) ∼ B κ pm ( t ) (cid:3) = at t = 0 second (75) Representing non-linear equations in the Lagrangian description and inte-grating them in the material frame can cause excessive mesh distortions. Thisdiscrepancy of the method is because, each of the nodes of the mesh, followthe corresponding material points during the motion. The Eulerian descrip-tion in which the nodes remain fixed while the material distorts, overcomesthe above discrepancy. However, it compromises the numerical accuracy ofthe solution while computing with coarse meshes. It also fails to preciselyidentify the interface between domains, leading to poor numerical accuracyof the interfacial fluxes.The disadvantages of both the descriptions can be avoided by using theArbitrary Lagrangian-Eulerian (ALE) method. ALE combines the best fea-30ures of both Lagrangian and Eulerian descriptions. Each node of the ALEmesh is free to either remain fixed or rezone itself (see Donea et. al. [59]),thus avoiding excessive mesh distortions and at the same time predicting thedomain interfaces with better accuracy. It also provides comparatively betternumerical accuracy of the solutions than the Eulerian approach.The non-linear equations, Eq.(68)-Eq.(74), have been solved using “Coef-ficient Form PDE” module in COMSOL Multiphysics TM . The coefficients ofthe general PDE in the module, are populated by extracting the respectivecoefficients of the governing equations by using a MATLAB code. It is tobe noted that the equations are provided to the module with respect to thespatial frame and are integrated in time by using the BDF method. Theorder of the scheme varies from 1 to 5 with free time stepping. Temperature in the melt starts reducing as heat flows out through thefree surfaces into the surroundings and through the fixed surface into thesubstrate (note that the rate of heat flow is higher into the substrate) andconsequently α m ( θ ) begins to fall. Due to the rapid cooling process, thetemperature reaches the ambient value in about 60 seconds. The temperaturedistribution in the layers is given in Fig.3 and the value of temperature ( θ )is higher towards the core than towards the surface, as is to be expected.Volumetric shrinkage in the first three seconds is not too evident. Thiscan be attributed to the re-heating of the relatively cooler layer each timecontinuity is established between consecutive layers. Once all the layers arelaid, and the entire system starts cooling, the volumetric shrinkage becomes31 a) Distribution of temperature in o C at 1 second. The layer has cooled from the intial temperature of241 o C. (b) Distribution of temperature in o C at 2 second. An increase in the temperature levels of the first layeris observed, which is attributed to the re-heating of the layer by the comparatively hotter second layer. c) Distribution of temperature in o C at 3 second. The temperature levels in the second layer increasesdue to re-heating by the comparatively hotter third layer. (d)
Distribution of temperature in o C at 60 second. The inset shows the longest layer having a length of0.02m, and the red box highlights the area of interest. The solid black line superimposed on the deformedconfiguration represents the initial geometry of the layers before warping and volumetric shrinkage has setin.
Figure 3: Banded contour plots of temperature in o C, visualised on the deformedconfiguration. The white arrows depict the local displacement vector, and itslength represents the magnitude. The relative length of arrows are in logarithmicscale. a) Distribution of the mean normal stress (Pa), K √ , at 1 second. The stress values are low because thetemperature values are well above θ g . (b) Distribution of the mean normal stress (Pa), K √ , at 2 second. Re-heating of the first layer ensures adelayed accumulation of stress in the material. c) Distribution of the mean normal stress (Pa), K √ , at 3 second. Appreciable stress values start devel-oping in the region where the temperature values are close to θ g . Re-heating delays the development ofany significant stress at the core of the geometry. (d) Distribution of the mean normal stress (Pa), K √ , at 60 second. The stress values have shot upsignificantly, and it is higher towards the core, in comparison to the values at the regions close to thestepped end. Figure 4: Banded contour plots of the mean normal stress (Pa), K √ , visualised onthe deformed configuration. a) Distribution of the norm of deviatoric part of stress (Pa), K , at 1 second. The deviatoric stressvalues are also insignificant due to the temperature values being sufficiently above θ g . (b) Distribution of the norm of deviatoric part of stress (Pa), K , at 2 second. The value remainsinsignificant due to the re-heating of the material. c) Distribution of the norm of deviatoric part of stress (Pa), K , at 3 second. A significant amount ofstress accumulation is observed in the region where the temperature values are close to θ g . (d) Distribution of the norm of deviatoric part of stress (Pa), K , at 60 second. The inset shows a regionhaving stress concentrations in the geometry. The point D inside this region is in a stress state given by:T xx = 1.38E8, T yy = 1.64E7, T zz = 9.94E7, T xy = 2.65E7. During service, a ductile failure is mostlikely to happen in the red zone containing the point D. Figure 5: Banded contour plots of the norm of deviatoric part of stress (Pa), K ,visualised on the deformed configuration. a) Distribution of the modes of deformation (K ) at 1 second. (b) Distribution of the modes of deformation (K ) at 2 second. c) Distribution of the modes of deformation (K ) at 3 second. (d) Distribution of the modes of deformation (K ) at 60 second. Point A is close to uniaxial tension(K ≈ . ≈
0) and point C is close to equibiaxial tension (K ≈− . Figure 6: Banded contour plots of the modes of deformation, (K invariant of thestress tensor), visualised on the deformed configuration. The dark red zones atthe stepped end is close to uniaxial tension and dark blue zones at the core of thegeometry is close to equibiaxial tension throughout the process. a) Distribution of the total rate of entropy production ( ξ totm ) at 1 second. (b) Distribution of the total rate of entropy production ( ξ totm ) at 2 second. c) Distribution of the total rate of entropy production ( ξ totm ) at 3 second. (d) Distribution of the total rate of entropy production ( ξ totm ) at 60 second. Figure 7: Banded contour plots of the total rate of entropy production ( ξ totm )visualised on the deformed configuration. u x components of the displacement field are positive and u y compo-nents are negative (the orientation of the basis are given in Fig.2). Towardsthe core, the u x components become negligible and only the u y componentssurvive. However, when the existing layers are re-heated each time a newlayer is laid, the core of the geometry starts expanding due to the inflow ofheat, which leads to negative u x components and positive u y components atthe core, although, the net effect of u y is in the negative direction through-out the material. Thus, the resultant displacement vectors near the cornerswould be directed downwards at an acute angle (clockwise) with respectto the x -direction, whereas the resultant displacement vectors at the corewould be initially directed downwards (refer to Figs.3a-7a), and later due tore-heating, as time proceeds it would be directed downwards at an obtuseangle (clockwise) with respect to the x -direction (refer to Figs.3d-7d).The thermal and mechanical interactions lead to a complicated distribu-tion of the residual stresses (refer to Fig.4 and Fig. 5) and consequently themodes of deformation (refer to Fig.6). The complexity is profound in thevicinity of the stepped geometry, whereas towards the opposite end, the dis-tribution becomes one dimensional, i.e., in the y -direction. The magnitudeof the stresses are negligible during the initial time steps (refer to Figs.4a-4c42nd Figs.5a-5c) due to the relatively lower volume shrikange. Re-heating ofthe layers, each time continuity is established, further adds to the relaxationof stresses. Magnitude of the stresses shoot up drastically when the tem-perature falls below θ g , which is expected since the viscosity of the polymerincreases sharply at θ g , and finally, the stresses freeze when the temperatureis sufficiently below θ g (refer to Fig.3d, Fig.4d and Fig.5d). Another rea-son for a large magnitude of stresses being induced in the material (refer toFig.4d and Fig.5d), is due to the high rate of cooling (refer to Table1), whichdoes not give enough time for the material to undergo stress relaxation.To understand the distribution of the modes, we consider a shear super-posed unequi-triaxial stress state of the material points. We mark points A,B and C on three different zones corresponding to three different modes ofdeformation, as shown in the Figs.6d. Note that the plane strain assumptionconstrains stress in the z -direction, i.e. T zz , to remain positive through-out the process. The stress state at point “A” is given by: T xx =0.5MPa,T yy =0.73MPa, T zz =27.7MPa and T xy =-0.19MPa. The magnitude of theT zz component is two orders higher than T xx , T yy and T xy components.Therefore, we can conclude that the mode of deformation of the materialpoint “A” is close to uniaxial tension in the z-direction with K ≈ .
52. Sim-ilarly, the mode of deformation of other points lying in the same zone will beclose to uniaxial tension in the z -direction.Next, the stress state at point “B” is given by: T xx =16.6MPa, T yy =-0.186MPa, T zz =34.8MPa and T xy =5.13MPa. The T xy component is quitesignificant, and therefore, the zone in which point “B” lies is close to apure shear mode of deformation with K ≈
0. Finally, the stress state at43oint “C” is given by: T xx =45.2MPa, T yy =-1.49MPa, T zz =47.5MPa andT xy =9.56MPa. T xx and T zz components are approximately same and pos-itive. Also, they are one order higher than that of T yy (which is negative)and T xy . Therefore, the mode of deformation at point “C” is close to equib-iaxial tension, and hence, the zone in which point “C” lies has K ≈ − . K √ ) to the norm of the deviatoric stresses(K ) is much less than 1 at the top corners of the layers and close to 1 atthe core. Therefore, the corners are dominated by the distortional stresses,whereas, the core is dominated by the mean normal stress, implying that, di-mensional instability is more at the corners than at the core. This is furthersupported by the bulging and the warping observed in the neighbourhoodof the stepping at the final stages of the simulation (refer to Figs.3d-7d).Warping can also be understood based on the direction of the componentsof the displacement fields at the top corners, where, u x is positive and u y isnegative. Therefore, the material at the topmost corner will get pushed out,leading to warping.During the entire process, the total rate of entropy production has to bemaintained positive. This needs to be checked a posteriori at every time stepof the process. The plots of total rate of entropy production at some chosentime steps have been given in Figs.7a-7d, and the values remain positive asexpected. 44 .6.2. Stress concentrations and delamination Figure 8: The stress components (T xy and T yy ) that are responsible for delam-ination (combination of mode I and mode II) of the layers, visualised along thelength of all the three interfaces (semi-log plot), after the melt has solidified (t =60 second). In a practical setting, the interfaces between the layers have discontinu-ities, and therefore, the inter-layer strength is usually the weakest, whichoften causes delamination of the layers for a component fabricated by FDM.Even though, in the current analysis, the interfaces are perfect, we assumethat the layers have a “tendency” to delaminate. Fig.8 helps us visualisethe stress components responsible for delamination along the length of eachof the interfaces. Only T yy and T xy components of the stress tensor causedelamination. T yy component tends to pull the layers apart (mode I delam-ination), whereas, T xy component tends to slide the layers in the plane ofthe interface (mode II delamination). It is observed that, T yy component45ecomes compressive in the region that is sufficiently away from the steppedend. Therefore, in this region, delamination can occur only due to the T xy component (mode II delamination). The tendency to delaminate is the high-est for the interface between the first and the second layer.The effect of mode I and mode II delamination tendencies, and the trac-tion free condition on the free edges, lead to the formation of stress concen-trations in the neighbourhood of the stepping (refer to Fig.5d). If delami-nation occurs, and the interface breaks apart near the stepping, the stressconcentrations would be relieved. The stress components are the highest (or-der of seven and eight) in the region of the stress concentration of the firstlayer. Therefore, yielding is most likely to occur in this region and lead to aductile failure termed as “intra-layer fracture”. Usually, inter-layer fracturedue to delamination and other discontinuities is the most common type offailure, compared to the intra-layer fracture caused by stress concentrations(see Macedo et. al. [22]). Therefore, a component fabricated by FDM, whenin service, will most probably fail due to the delaminated interfaces termedas“inter-layer fracture”.
4. Conclusion
A thermodynamic framework was developed to determine the differentialvolumetric shrinkage caused by drastic temperature gradients, the ensuingresidual stresses and consequently the dimensional instability of the geom-etry, by building on an existing theory which required the definition of twoprimitives: Helmholtz free energy and the rate of entropy production. Thetheory is most suitable for processes that use amorphous polymers or poly-46ers which do not crystallize much. The efficacy of the constitutive relationwas verified by considering a prototypical FDM process in which four layersof polystyrene melt were laid, such that, plane strain conditions could beenforced.As soon as the layers were laid, it underwent a complex thermo-mechanicalprocess which included rapid cooling, reheating, solidification at the glasstransition temperature ( θ g ) and redistribution of the residual stresses. Al-though, a high value of the shear modulus caused stresses in the range of kiloPascals in the melt phase, it is quite insignificant as compared to the stressesin the solid phase. This conforms with the reasoning put forward by Xinhuaet. al. [19], Wang et. al. [20], Park et. al. [21] and Macedo et. al. [22],who have confined the analysis to the solidified polymer since the stressesin the melt phase can be assumed to be negligible. The mean normal stressdistribution at the final time step (refer to Fig.4d) seems to follow a similartrend, in comparison to the mean normal stress distribution reported by Xiaand Lu et. al. [33] (refer to Fig.8 in [33]).The plane strain approximation caused the outer corners of the layers tobe close to uniaxial tension and the core to be close to equibiaxial tension.During service, there is a high probability for a pre-existing crack at thecore to open up in the x - z plane due to the equibiaxial tension. Further,in a practical setting, analogous to the current analysis, if the interfaces arehighly porous, a large amount of inter-layer fracture (a combination of mode Iand mode II) is expected at the stepped ends as compared to intra-layer frac-tures caused by stress concentrations. Another major dimensional instabilitythat occurred was the warping of the geometry. Warping can be reduced by47aking the corners in the stepping smoother. More smooth the corners are,less the degree of warping. The possibility of delamination, warping andthe magnitudes of the stress concentrations, can be reduced by choosing theright combination of rate of cooling, substrate temperature, re-heating andthe geometry for the PS layers, which reduces the magnitude of the “freezedin” stresses (see Macedo [22]). 48 eferencesReferences [1] K. Chaoui, A. Chudnovsky, A. Moet, Effect of residual stress on crackpropagation in mdpe pipes, Journal of Materials Science 22 (1987) 3873–3879.[2] L. E. Hornberger, K. L. Devries, The effects of residual stress on the me-chanical properties of glassy polymers, Polymer Engineering and Science27 (19) (October 1987) 1473–1478.[3] P. So, L. J. Broutman, Residual stresses in polymers and their effect onmechanical behavior, Polymer Engineering and Science 16 (12) (Decem-ber 1976) 785–791.[4] F. P. Melchels, M. A. Domingos, T. J. Klein, J. Malda, P. J. Bartolo,D. W. Hutmacher, Additive manufacturing of tissues and organs,Progress in Polymer Science 37 (8) (2012) 1079– 1104, topical Issue on Biorelated polymers. doi:https://doi.org/10.1016/j.progpolymsci.2011.11.007 .URL [5] A. Melocchi, F. Parietti, A. Maroni, A. Foppoli, A. Gazzaniga, L. Zema,Hot-melt extruded filaments based on pharmaceutical grade polymers for 3d printing by fused deposition modeling,International Journal of Pharmaceutics 509 (1) (2016) 255 – 263. doi:https://doi.org/10.1016/j.ijpharm.2016.05.036 .URL doi:10.1108/13552541011034825 .[7] E. L. Nyberg, A. L. Farris, B. P. Hung, M. Dias, J. R. Garcia, A. H.Dorafshar, W. L. Grayson, 3d-printing technologies for craniofacial re-habilitation, reconstruction, and regeneration, Annals of biomedical en-gineering 45 (1) (2017) 45–57. doi:10.1007/s10439-016-1668-5 .[8] D. Muller, H. Chim, A. Bader, M. Whiteman, J.-T. Schantz, Vascularguidance:microstructural scaffold patterning for inductive neovascular-ization, Stem Cells International.[9] X. Han, D. Yang, C. Yang, S. Spintzyk, L. Schei-deler, P. Li, D. Li, J. Geis-Gerstorfer, F. Rupp,Carbon fiber reinforced peek composites based on 3d-printing technology for orthopedic and dental applications,Journal of Clinical Medicine 8 (2). doi:10.3390/jcm8020240 .URL [10] M. Vaezi, S. Yang, Extrusion-based additive manufacturing of peek for biomedical applications,Virtual and Physical Prototyping 10 (3) (2015) 123–135. arXiv:https://doi.org/10.1080/17452759.2015.1097053 , doi:10.1080/17452759.2015.1097053 .URL https://doi.org/10.1080/17452759.2015.1097053 [11] G. Syuhada, G. Ramahdita, A. J. Rahyussalim, Y. Whulanza,Multi-material poly(lactic acid) scaffold fabricated via fused deposition modeling and direct hydroxyapatite injection as spacers in laminoplasty,50IP Conference Proceedings 1933 (1) (2018) 020008. arXiv:https://aip.scitation.org/doi/pdf/10.1063/1.5023942 , doi:10.1063/1.5023942 .URL https://aip.scitation.org/doi/abs/10.1063/1.5023942 [12] S. Harish, V. R. Devadath, Additive manufacturing and analysis of tibialinsert in total knee replacement implant, International Research Journalof Engineering and Technology 02 (2015) 633–638.[13] R. A. Borges, D. Choudhury, M. Zou,3d printed pcu/uhmwpe polymeric blend for artificial knee meniscus,Tribology International 122 (2018) 1 – 7. doi:https://doi.org/10.1016/j.triboint.2018.01.065 .URL [14] J. ten Kate, G. Smit, P. Breedveld,3d-printed upper limb prostheses: a review, Disability and Rehabilita-tion: Assistive Technology 12 (3) (2017) 300–314, pMID: 28152642. arXiv:https://doi.org/10.1080/17483107.2016.1253117 , doi:10.1080/17483107.2016.1253117 .URL https://doi.org/10.1080/17483107.2016.1253117 [15] L. Freitag, M. G¨ordes, P. Zarogoulidis, K. Darwiche, D. Franzen,F. Funke, W. Hohenforst-Schmidt, H. Dutau, Towards individual-ized tracheobronchial stents: Technical, practical and legal consider-ations., Respiration: Interventional Pulmonology 94 (2017) 442–456. doi:10.1159/000479164 . 5116] A. K. Sood, R. K. Ohdar, S. S. Mahapatra, Parametric ap-praisal of mechanical property of fused deposition modellingprocessed parts, Materials and Design 31 (2010) 287–295. doi:10.1016/j.matdes.2009.06.016 .[17] Y. Song, Y. Li, W. Song, K. Yee, K.-Y. Lee, V. Tagarielli,Measurements of the mechanical response of unidirectional 3d-printed pla,Materials & Design 123 (2017) 154 – 164. doi:https://doi.org/10.1016/j.matdes.2017.03.051 .URL [18] C. Casavola, A. Cazzato, V. Moramarco, G. Pappalettera,Residual stress measurement in fused deposition modelling parts,Polymer Testing 58 (2017) 249 – 255. doi:https://doi.org/10.1016/j.polymertesting.2017.01.003 .URL [19] L. Xinhua, L. Shengpeng, L. Zhou, Z. Xianhua, C. Xiaohu, W. Zhongbin,An investigation on distortion of pla thin-plate part in the fdm process,The International Journal of Advanced Manufacturing Technology79 (5) (2015) 1117–1126. doi:10.1007/s00170-015-6893-9 .URL https://doi.org/10.1007/s00170-015-6893-9 [20] T.-M. Wang, J.-T. Xi, Y. Jin, A model research for prototype warp deformation in the fdm process,The International Journal of Advanced Manufacturing Technology33 (11) (2007) 1087–1096. doi:10.1007/s00170-006-0556-9 .URL https://doi.org/10.1007/s00170-006-0556-9 doi:doi://10.26678/ABCM.COBEM2017.COB17-0124 .[23] H. Xia, J. Lu, S. Dabiri, G. Tryggvason, Fully resolved numerical simula-tions of fused deposition modeling. part i: fluid flow, Rapid PrototypingJournal 24 (2) (2018) 463–476.[24] S. Dabiri, S. Schmid, G. Tryggvason, Fully resolved numerical simu-lations of fused deposition modeling, in: Proceedings of the ASME2014 International Manufacturing & Science and Engineering Confer-ence, MSEC 2014, Detroit, 2014.[25] R. Comminal, M. P. Serdeczny, D. B. Pedersen, J. Spangenberg,Numerical modeling of the strand deposition flow in extrusion-based additive manufacturing,Additive Manufacturing 20 (2018) 68 – 76. doi:https://doi.org/10.1016/j.addma.2017.12.013 .URL [26] A. El Moumen, M. Tarfaoui, K. Lafdi, Modelling of the tem-perature and residual stress fields during 3d printing of poly-mer composites, Int J Adv Manuf Technol 104 (2019) 1661–1676. doi:https://doi.org/10.1007/s00170-019-03965-y .5327] C. Casavola, A. Cazzato, V. Moramarco, C. Pappalettere,Orthotropic mechanical properties of fused deposition modelling parts described by classical laminate theory,Materials & Design 90 (2016) 453 – 458. doi:https://doi.org/10.1016/j.matdes.2015.11.009 .URL [28] A. Armillotta, M. Bellotti, M. Cavallaro,Warpage of fdm parts: Experimental tests and analytic model,Robotics and Computer-Integrated Manufacturing 50 (2018) 140– 152. doi:https://doi.org/10.1016/j.rcim.2017.09.007 .URL [29] A. Cattenone, S. Morganti, G. Alaimo, F. Auricchio,Finite Element Analysis of Additive Manufacturing Based on Fused Deposition Modeling: Distortions Prediction and Comparison With Experimental Data,Journal of Manufacturing Science and Engineering 141 (1), 011010. arXiv:https://asmedigitalcollection.asme.org/manufacturingscience/article-pdf/141/1/011010/6278955/manu\_141\_01\_011010.pdf , doi:10.1115/1.4041626 .URL https://doi.org/10.1115/1.4041626 [30] M. Doi, S. F. Edwards, The Theory of Polymer Dynamics, Oxford Uni-versity Press, Walton Street, Oxford, 1986.[31] M. R. Kamal, R. A. Rai-Fook, J. R. Hernandez-Aguilar, Residual ther-mal stresses in injection moldings of thermoplastics: A theoretical andexperimental study, Polymer-Plastics Technology and Engineering 42 (5)(May 2002) 1098–1114.[32] W. F. Zoetelief, L. F. A. Douven, A. J. I. Housz, Residual thermal54tresses in injection molded products, Polymer Engineering and Science36 (14) (July 1996) 1886–1896.[33] H. Xia, J. Lu, G. Tryggvason, Fully resolved numerical simulations offused deposition modeling. part ii – solidification, residual stresses andmodeling of the nozzle, Rapid Prototyping Journal 24 (6) (2018) 973–987.[34] L. A. Northcutt, S. V. Orski, K. B. Migler, A. P. Kotula,Effect of processing conditions on crystallization kinetics during materials extrusion additive manufacturing,Polymer 154 (2018) 182 – 187. doi:https://doi.org/10.1016/j.polymer.2018.09.018 .URL [35] I. J. Rao, K. R. Rajagopal, Phenomenological modelling of polymer crys-tallization using the notion of multiple natural configurations, Interfacesand Free Boundaries 2 (2000) 73–94.[36] I. J. Rao, K. R. Rajagopal, A thermodynamic framework for the studyof crystallization in polymers, Zeitschrift fur Angewandte Mathematikund Physik 53 (2002) 365–406. doi:10.1007/s00033-002-8161-8 .[37] K. Kannan, I. J. Rao, K. R. Rajagopal,A thermomechanical framework for the glass transition phenomenon in certain polymers and its application to fiber spinning,Journal of Rheology 46 (4) (2002) 977–999. arXiv:https://sor.scitation.org/doi/pdf/10.1122/1.1485281 , doi:10.1122/1.1485281 .URL https://sor.scitation.org/doi/abs/10.1122/1.1485281 arXiv:https://doi.org/10.1122/1.1879042 , doi:10.1122/1.1879042 .URL https://doi.org/10.1122/1.1879042 [39] K. Kannan, K. R. Rajagopal, A thermomechanical framework for the transition of a miscoelastic liquid to a viscoelastic solid,Mathematics and Mechanics of Solids 9 (1) (2004) 37–59. arXiv:https://doi.org/10.1177/1081286503035198 , doi:10.1177/1081286503035198 .URL https://doi.org/10.1177/1081286503035198 [40] J. Liu, K. L. Anderson, N. Sridhar,Direct simulation of polymer fused deposition modeling (fdm) — an implementation of the multi-phase viscoelastic solver in openfoam,International Journal of Computational Methods 17 (01) (2020)1844002. arXiv:https://doi.org/10.1142/S0219876218440024 , doi:10.1142/S0219876218440024 .URL https://doi.org/10.1142/S0219876218440024 [41] M. Avrami, Kinetics of phase change. i general theory,The Journal of Chemical Physics 7 (12) (1939) 1103–1112. arXiv:https://doi.org/10.1063/1.1750380 , doi:10.1063/1.1750380 .URL https://doi.org/10.1063/1.1750380 [42] K. Rajagopal, A. Srinivasa, A thermodynamic frame work for rate type fluid models,Journal of Non-Newtonian Fluid Mechanics 88 (3) (2000) 207 – 227.56 oi:https://doi.org/10.1016/S0377-0257(99)00023-3 .URL [43] I. Rao, K. Rajagopal, A study of strain-induced crystallization of polymers,International Journal of Solids and Structures 38 (6) (2001) 1149 –1167. doi:https://doi.org/10.1016/S0020-7683(00)00079-2 .URL [44] K. R. Rajagopal, A. R. Srinivasa,On the thermomechanics of materials that have multiple natural configurations part i: Viscoelasticity and classical plasticity,Zeitschrift f¨ur angewandte Mathematik und Physik ZAMP 55 (5) (2004)861–893. doi:10.1007/s00033-004-4019-6 .URL https://doi.org/10.1007/s00033-004-4019-6 [45] M. Chen, Y. Tan, B. Wang, General invariant representations of the constitutive equations for isotropic nonlinearly elastic materials,International Journal of Solids and Structures 49 (2) (2012) 318 – 327. doi:https://doi.org/10.1016/j.ijsolstr.2011.10.008 .URL [46] A. E. Green, P. M. Naghdi, On thermodynamics and the nature of the second law,Proceedings of the Royal Society of London. A. Mathe-matical and Physical Sciences 357 (1690) (1977) 253–270. arXiv:https://royalsocietypublishing.org/doi/pdf/10.1098/rspa.1977.0166 , doi:10.1098/rspa.1977.0166 .URL https://royalsocietypublishing.org/doi/abs/10.1098/rspa.1977.0166 [47] K. R. Rajagopal, A. R. Srinivasa,On thermomechanical restrictions of continua, Proceedings of57he Royal Society of London. Series A: Mathematical, Phys-ical and Engineering Sciences 460 (2042) (2004) 631–651. arXiv:https://royalsocietypublishing.org/doi/pdf/10.1098/rspa.2002.1111 , doi:10.1098/rspa.2002.1111 .URL https://royalsocietypublishing.org/doi/abs/10.1098/rspa.2002.1111 [48] C. M. Gonz´alez-Henr´ıquez, M. A. Sarabia-Vallejos, J. Rodriguez-Hernandez, Polymers for additive manufacturing and 4d-printing: Materials, methodologies, and biomedical applications,Progress in Polymer Science 94 (2019) 57 – 116. doi:https://doi.org/10.1016/j.progpolymsci.2019.03.001 .URL [49] R. L¨offler, M. Koch, Innovative extruder concept for fast and efficient additive manufacturing,IFAC-PapersOnLine 52 (10) (2019) 242 – 247, 13th IFACWorkshop on Intelligent Manufacturing Systems IMS 2019. doi:https://doi.org/10.1016/j.ifacol.2019.10.071 .URL [50] T. Kim, R. Trangkanukulkij, W. Kim,Nozzle shape guided filler orientation in 3d printed photo-curable nanocomposites,Sci Rep 8 (2018) 57 – 116. doi:https://doi.org/10.1038/s41598-018-22107-0 .URL [51] W. M. Haynes, Hand book of Chemistry and Physics, CRC Press, 6000Broken Sound Parkway NW, Suite 300, 2014.[52] Y. Chai, J. A. Forrest, Using atomic force microscopy to probe crystallization in atactic polystyrenes,Macromolecular Chemistry and Physics 219 (3) (2018) 1700466.58 rXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/macp.201700466 , doi:10.1002/macp.201700466 .URL https://onlinelibrary.wiley.com/doi/abs/10.1002/macp.201700466 [53] D. van Krevelen, K. te Nijenhuis, Properties of Polymers, Fourth Edi-tion: Their Correlation with Chemical Structure; their Numerical Es-timation and Prediction from Additive Group Contributions, ElsevierScience, Linacre House, Jordan Hill, Oxford OX2 8DP, UK, 2009.[54] N. Sombatsompop, A. Wood, Measurement of thermal conductivity of polymers using an improved lee’s disc apparatus,Polymer Testing 16 (3) (1997) 203 – 223. doi:https://doi.org/10.1016/S0142-9418(96)00043-8 .URL [55] S. Costa, F. Duarte, J. Covas, Thermal conditions affecting heat transfer in fdm/ffe: a contribution towards the numerical modelling of the process,Virtual and Physical Prototyping 10 (1) (2015) 35–46. arXiv:https://doi.org/10.1080/17452759.2014.984042 , doi:10.1080/17452759.2014.984042 .URL https://doi.org/10.1080/17452759.2014.984042 [56] E. Marti, E. Kaisersberger, E. Moukhina,Heat capacity functions of polystyrene in glassy and in liquid amorphous state and glass transition,J Therm Anal Calorim 85 (2006) 505–525. doi:https://doi.org/10.1007/s10973-006-7745-5 .URL https://link.springer.com/article/10.1007/s10973-006-7745-5 [57] J. E. Mark, Polymer Data Handbook, Oxford University Press, NewYork 10016, 1999. 5958] A. A. Miller, Analysis of the melt viscosity and glass transition of polystyrene,Journal of Polymer Science Part A-2: Poly-mer Physics 6 (6) (1968) 1161–1175. arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/pol.1968.160060609 , doi:10.1002/pol.1968.160060609 .URL https://onlinelibrary.wiley.com/doi/abs/10.1002/pol.1968.160060609 [59] J. Donea, A. Huerta, J.-P. Ponthot, A. Rodr´ıguez-Ferran, Arbitrary Lagrangian–Eulerian Methods,American Cancer Society, 2004, Ch. 14. arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/0470091355.ecm009 , doi:10.1002/0470091355.ecm009 .URL https://onlinelibrary.wiley.com/doi/abs/10.1002/0470091355.ecm009https://onlinelibrary.wiley.com/doi/abs/10.1002/0470091355.ecm009