We model two common fused lament fabrication mesostructures, square and hexagonal, using an orthotropic constitutive model and derive closed-form expressions for all nine e ective elastic constants. The periodic void shapes are modeled using three and four point hypotrochoid curves with a single shape parameter that controls the sharpness of the points. Using the complex variable method of elasticity, we derive the in-plane elastic constants ( Exx,Eyy,xy,Gxy) as well as out-of-plane antiplane shear constants (GzxandGzy). The remaining out-of-plane elastic constants ( Ezz,zx,zy) are derived by directly solving the linearelasticity equations. We compare our results by conducting unit cell simulations on both mesostructures and at various porosity values. The simulations match the closed-form expressions exactly forEzz,zx, andzy. For the remaining elastic constants, the simulation results match the closed-form expressions better for the square mesostructure than the hexagonal mesostructure. Di erences between simulation and closed-form expressions are less than 10% for porosity values less than 6% (hexagonal mesostructure) and 10% (square mesostructure) for any of the nine elastic constants. Keywords: fused lament fabrication, constitutive model, orthotropic, porosity, voids, homogenization 1 Introduction Fused lament fabrication is an extrusion-based additive manufacturing technique that uses a heated nozzle and motor to melt thermoplastic lament feedstock and extrude the molten lament in a pre-de ned set of print paths [1]. The technique is widely used for rapid prototyping [2]{[4], tooling [5], [6], production parts [7], and teaching [8]{[10]. However, widespread adoption of FFF for commercial production is still limited. Several limitations are theorized to be reasons driving down commercial demand: long print times [11], coarse surface nish [12], coarse geometric resolution [13], low interlayer bonding strength [14]{[16], and highly anisotropic behavior [17]{[21]. Nonetheless, there are certain applications where FFF is a better t; for example, the International Space Station operates a FFF machine because it can work in microgravity, the feedstock is relatively nonhazardous, and the printed parts are recyclable [22], [23]. In order to design parts speci cally for FFF, a detailed understanding of process parameters and how they a ect part behavior is necessary. Prior work has shown that a variety of parameters, such as layer height [24], air gap [12], raster angle and pattern [19], in ll pattern and density [25], [26], processing temperature [27]{[30], and even lament color [31] can all a ect the overall mechanics of FFF parts. In this work, we isolate and look at how the void shape and the overall porosity of two common FFF mesostructures a ect the overall e ective compliance of the printed material. Although the void shape and porosity are a result of many of the aforementioned process parameters, we use the void shape and porosity as starting points, abstracting away any upstream causes. However, we will brie y talk about process parameters that may a ect the compliance of a single extruded road compared to the lament feedstock. Cooresponding author: rchensix@alumni.stanford.edu ydsenesky@stanford.edu 1 1.1 Outline In Section 2, we will review prior literature on constitutive models for FFF materials as well as some relevant process parameters. In Section 3, we will de ne the mesostructures and material coordinate system that all subsequent analysis will use. In Sections 5.4, 5, and 6, we will derive all nine e ective elastic constants with the aid of the complex variable method of elasticity. Finally, in Section 7, we compare our closed-form solutions for the nine elastic constants with nite element analysis on unit cells and talk about ranges of porosity values where our model matches well with simulation. 2 Literature review 2.1 Existing constitutive models Methods for characterizing the apparent or e ective sti ness of FFF materials are generally borrowed from techniques used in the composites industry. The rule of mixtures, rst proposed by Voigt [32], is applicable for a variety of engineering material systems, including ber-reinforced composites [33], honeycombs [34], [35], and FFF materials [36], [37]. In fact, the rule of mixtures yields an exact solution for the e ective Young's modulus of multiphase unidirectional material systems (like the ones previously mentioned) along the longitudinal material axis, as long as the Poisson ratio of all material phases are equal. However, the rule of mixtures makes poor predictions for transverse loading. Reuss [38] proposed an inverse rule of mixtures for transverse loading, but this yields a trivial result of zero sti ness in the case of materials with voids. Modi cations to the rule of mixtures, such as adding empirical t factors, measuring the e ective area in di erent cross sectional planes, and re-formulating it as a power law have been used by several authors for FFF materials [35], [37], [39]. A comprehensive survey of empirical and semi-empirical models can be found in Choren's work [40]. Due to a lack of closed-form models speci cally for FFF materials, there has been signi cant interest in simulating the e ective elastic constants using homogenization theory. Rodriguez simulated the as-printed void by converting cross section photographs into nite element meshes [37], resulting in a monoclinic constitutive model (due to asymmetries in the as-printed void). By assuming an ideal void shape, simpler orthotropic models have been simulated [41], [42]. These models are expensive to compute as they require modeling the void geometry and solving a nite element problem for every void geometry and mesostructure variation. Nonetheless, they can be highly accurate if the cross section geometry is known. While present-day computers make the computational expense less of a concern, there is still a need for closed-form solutions, especially for performance-critical applications. For example, topology optimization that takes mesostructure geometry into consideration can become prohibitively expensive without a closed-form e ective sti ness model. Closed-form solutions for non-ellipsoidal voids would be impossible if not for the complex variable method of elasticity developed by Muskelishvili [43], which enabled mapping complex void shapes into the complex unit circle. This technique was used by several authors to derive e ective 2D elastic constants for elliptical [44], [45], polygonal [44], [46], [47], and arbitrary void shapes [48], [49]. In this work, we build on this line of work to develop an e ective orthotropic constitutive model for FFF materials. 2.2 Process parameters that a ect feedstock compliance Before developing our model, we mention some process parameters that may a ect the compliance of a single printed road, compared to the lament feedstock. Since all of the aforementioned methods require elastic constants of the bulk material (i.e. no voids) as an input, any process parameters that may change the bulk material properties are relevant. Rodriguez [27] and Bellini [50] performed tensile tests on both feedstock lament as well as single extruded roads. Bellini's results show almost no change in Young's modulus, while Rodriguez's results show a 4% decrease in Young's modulus. The slight decrease noted by Rodriguez is theorized to be caused by thermally- induced randomization of the polymer chains in the extruded road. Apart from these two studies, there is scant literature on experiments comparing the feedstock and a single extruded road; most studies compare the feedstock directly with printed tensile specimens, which already incorporate voids and a particular 2 (a) (b)zxyFigure 1: Unit cells of (a) hexagonal and (b) square mesostructures. mesostructure. In this study, we simply assume isotropic material properties Eandare provided as input parameters. While we believe using values derived from experiments on a single extruded road may yield marginally better results, the di erence seems to be small and will probably be dwarfed by more in uential factors, such as void shape and porosity. 3 Material frame We de ne a material aligned coordinate system, as shown in Figure 1. The z-axis is aligned with the direction of the extruded road. We will refer to properties in the xyplane as in-plane properties, and anything else as out-of-plane properties. We consider two mesostructure systems: the hexagonal mesostructure and the square mesostructure, which are named for the periodic arrangement of extruded roads (in works by Rodriguez, they are referred to as the skewed and aligned mesostructures [27], [36], [37]). When the extruded roads have perfectly circular cross sections, the voids formed by the gaps between adjacent roads form curved triangles and diamonds. In this article, we will interchangeably associate the hexagonal mesostructure with curved triangle voids and the square mesostructure with curved diamond voids. 4 Axial loading in z-direction For this section only, we assume that the voids can take any shape, as long as the overall structure is prismatic (i.e. the cross section does not change along the z-direction and the front/back faces are at). To nd the e ective elastic constants Ezz,zx, andzy, we apply a known displacement to the front/back faces while keeping the side face traction free. We also assume no body forces. The boundary value problem at 3 S/g3675S/g3676 S/g3677 zFigure 2: A general prismatic solid containing voids. For deriving elastic constants using axial loading in thez-direction, the voids can take any shape. hand is given by Equation 1. r= 0 in uz= 0 on S1 uz=u0 onS2 zx=zy= 0 on S1 zx=zy= 0 on S2 n= 0 on S3(1) where is the material domain and surfaces S1,S2, andS3are de ned in Figure 2. The solution to this problem is given by a constant stress eld zz=Eu0 L(2) whereEis the Young's modulus of the lament material and Lis the distance from S1toS2. All other stress components are zero. We de ne the e ective Young's modulus Exxas the constant of proportionality relating the average stress to the average strain (averaged over the volume). We also de ne a porosity ratio fas the ratio of the void volume to the total volume. The e ective Young's modulus can then be expressed as Ezz=avg zz avg zz =(1f)Eu0 L u0 L = (1f)E(3) The result is the venerable rule of mixtures for a two phase material system, where one phase is the void itself. Similarly, we can de ne the e ective major Poisson ratios zxandzyas zx=avg xx avg zz zy=avg yy avg zz(4) 4 Forzx, we get zx=zz E u0 L =(5) whereis the Poisson ratio of the lament material. Following the same approach for zyyieldszy=. These results are well known in linearelasticity theory and have been used to derive e ective elastic constants for other prismatic porous structures such as honeycombs [35]. Finally, we mention that the results of this section are exact for anyporosity ratio f2[0;1]. This will not be the case in subsequent sections. 5 In-plane loading in x- and y-directions To derive closed form expressions for the e ective elastic constants in the x-y plane, we break the problem into two parts. The rst part involves solving for the stress and displacement elds around a single void of a known geometry embedded in an in nite media. This part can be solved exactly using the complex variable method of elasticity [43], [51], [52]. Once the stress and displacement elds are known, we nd the e ective elastic constants by equating strain energy in the void-containing material to a homogeneous, void-free material [44], [46], [47]. For low porosities, we can assume that the voids are non-interacting. This assumption allows us to calculate the strain energy of a material containing multiple voids by summing up the strain energy contribution of each individual void. 5.1 Mapping function For fused lament fabrication, the hexagonal and square mesostructures produce periodic voids in the shape of curved triangles and diamonds, respectively. The shape of the voids can be approximated by hypotrochoid voids [53], [54]. Speci cally, the curved triangle is approximated by the conformal map mt=+ct 2(6) and the curved diamond is approximated by md=+cd 3(7) where=ei represents a coordinate in the complex -plane. These mapping functions transform the complex unit circle (given by !ei for 2[0;2]) into a hypotrochoid void. The shape parameters ct2[1 2;1 2] andcd2[1 3;1 3] control the sharpness of the corners of the voids. When ctandcdare zero, the voids are circular; when they take their maximum allowable values, the voids contain cusps at the corners; when they take their minimum allowable values, the voids are the same shape as the maximum case, but rotated by 60and 45, respectively. We can also write the void shapes in parametric form by taking the real and imaginary parts of Equations 6 and 7 to get xt= cos +ctcos 2 yt= sin ctsin 2 (8) for the curved triangle and xd= cos +cdcos 3 yd= sin cdsin 3 (9) for the curved diamond. A variety of voids are shown in Figure 3. 5 -1/2 -1/4 0 1/4 1/2 -1/3 -1/5 0 1/5 1/3ct cdFigure 3: Curved triangle (top row) and curved diamond voids (bottom row) constructed from various shape parameter values. 5.2 Stress and displacement eld around a single void We consider a single void embedded in an in nite media, subject to an applied uniaxial stress far away from the void. The direction of loading will determine which elastic constant we solve for, but for now we will formulate it for the general in-plane loading case. The 2D boundary value problem at hand becomes r= 0 in xx=0 xx forjrj!1 yy=0 yy forjrj!1 xy=0 xy forjrj!1 n= 0 on S(10) wherer2=x2+y2andSis the surface of the void. To solve this boundary value problem, we use the complex variable method of elasticity. Here, we give a brief description of the method. More details of the procedure can be found in textbooks by Muskhelishvili, England, and Barber [43], [51], [52]. The complex variable method of elasticity involves nding two holomorphic functions () and(), also known as the Kolosov-Muskhelishvili potentials. The potentials are related to stress and displacement elds by the following:  = 2 d =d dm=d+d =d dm=d!  =2 dm=d m d2 =d2 dm=dd dd2m=d2 (dm=d )2! +d d! u=1 2G  md =d dm=d!(11) where xx+yy yyxx+ 2ixy u=ux+iuy =3 1 +(12) 6 an overbar refers to the complex conjugate, and mis the mapping function. We note that xx,yy,xy, and  are all real-valued functions, but expressed in terms of the complex variable . We can rewrite Equation 12 as xx=<f g 2 yy= yy xy==fg 2(13) where<f:gand=f:grefer to taking the real and imaginary components, respectively. To nd and, we can break up the problem into two pieces: nding the uniform stress solution with no void, then nding the perturbation solution that makes the void traction free. Thus, = 0+ 1 =0+1(14) where 0-subscripts refer to the uniform solution and 1-subscripts refer to the perturbation solution. The uniform solution can be found by inspection by setting xx=0 xx,yy=0 yy, andxy=0 xyin Equation 12. Noting that the expressions for  and  in Equation 11 can be re-written as  = 2d dm+d dm  = 2 md2 dm2+d dm (15) we nd that 0=0 xx+0 yy 4m 0=0 yy0 xx+ 2i0 xy 2m(16) To nd 1and1, we rst de ne a traction function on the void surface s=ei 2S: (s) = (s) +m(s)d =d dm=d =1 s+(s) (17) This function can also be split up as = 0+ 1where the subscripts refer to the uniform solution and perturbation solution as before. Substituting in Equation 16 into 17 yields 0(s) =0 xx+0 yy 2m(s) +0 yy0 xx2i0 xy 2m(s) (18) However, the void must be traction free (i.e. = 0, so therefore, 1= 0=0 xx+0 yy 2m(s)0 yy0 xx2i0 xy 2m(s) (19) Given 1, we can nd 1and1by solving the following integro-di erential equations [52]: 1=1 2iI S 1(s) sds+1 2iI Sm(s)d 1(s)=ds dm(s)=ds(s)ds 1=1 2iI S 1(s) sds+1 2iI Sm(s)d 1=ds dm(s)=ds(s)ds(20) 7 To nd 1, we express 1and its derivative using a Laurent series with only the negative powers 1() =1X n=1an n d 1 d=1X n=1nan n+1 d 1 d=1X n=1nan n+1(21) Only the negative powers are included because we want the perturbation solution to disappear as jj!1 . In addition, even though Equation 21, in general, has an in nite number of terms, there will only be a nite number of terms in the nal solution due to simple form of the mapping functions given in Equations 6 and 7. After substituting Equation 21 along with a mapping function into the rst equation in Equation 20, we can nd the coecients anby comparing terms on both sides of the equation. Finally, after nding 1, we can nd1by substituting our result for 1into the second equation in Equation 20. The contour integrals in Equation 20 can be solved by using Cauchy's integral formula or the residue theorem. Once andare known, the stress, strain, and displacement elds can be calculated using Equation 11. 5.3 Strain energy equivalence Once the stress, strain, and displacement elds around a single void are known, we can use strain energy equivalence to derive e ective elastic constants. The strain energy of a single void under uniform stress far away from the void can be expressed using Eshelby's formula [55]: W=W0+W1 W0=1 2ZZZ V0 ij0 ijdV W1=1 2ZZ S0 ijnjuidS(22) whereWis strain energy, W0is the strain energy contribution due to the uniform applied stress, W1is the strain energy contribution due to the presense of a single void, Vis the volume of the in nite media, Sis the void surface, 0 ijand0 ijare the uniform stress and strain far away from the void, uiis the displacement, andnjis the unit normal vector along the surface of the void. We have to be very careful here, as ui() andnj() represent quantities in the x-y plane, but expressed as functions of . Since we're evaluating the line integral in Equation 22 along the void, uiandnjwithin the integrand are really functions of s=ei (or equivocally, just ). We can de ne the void surface in the x-y plane parametrically as a function of : x( ) =<fm()j!ei g y( ) ==fm()j!ei g(23) This allows us to express the components of the unit normal as nx=dy d r dx d 2 + dy d 2 ny=dx d r dx d 2 + dy d 2(24) In addition, the quantity dSbecomes dS=sdx d 2 +dy d 2 d dz (25) 8 Using these relations, the term W1in Equation 22 becomes W1=1 2Zt 0Z2 0 0 xxuxdy d 0 yyuydx d +0 xy uydy d uxdx d  d dz (26) For periodic mesostructures in FFF materials that contain multiple voids, we assume that the voids are small and do not interact with each other, such that we can express the strain energy contribution of all the voids as a sum of the individual voids. Under these assumptions, we can set Vto be the volume of a representative volume element (RVE) containing a single void. We de ne the porosity fas f=Vvoid VRVE=Avoid ARVE(27) whereArefers to area and Vrefers to volume. To nd the e ective elastic constants, we equate the strain energy in the RVE with a single void to the strain energy of an e ective homogeneous RVE (containing no void) subject to the same uniform loading condition to get W=W0+W1=0 ijSeff ijkl0 kl (28) whereSeff ijklis the e ective compliance matrix of the homogenized RVE. By applying stresses in particular directions, we can isolate and solve for the components Exx,Eyy,xy, andGxyembedded in the e ective compliance matrix. 5.4 Axial loading in x-direction The uniform loading far away from the void for this case is given by 0 xx=0,0 yy=0 xy= 0. Evaluating Equation 20 for the case of a curved triangle void, we nd the Kolosov-Muskhelishvili potentials to be =0 4+0 2ct0 42 =0 4+2+ 2 21 c2 t2ct1 4ct23(29) for the curved triangular void and =0 41 2(cd1)cd 43 =0 123c2 d2+ 3c3 d24+cd 3 +2+4 2 (3c2 d+4cd(3 +4))(30) for the curved diamond void. For the uniaxial loading case in the x-direction, Equation 22 simpli es to W=1 22 0 EARVEt+1 20Zt 0Z2 0uxdy d dzd =1 22 0 EAvoid ft+1 20tZ2 0uxdy d d (31) 1 22 0 EAvoid ft+1 20tZ2 0uxdy d d =1 22 0 EyyAvoid ft (32) The void areas for the curved triangular void and curved diamond void are given by At void= 12c2 t Ad void= 13c2 d (33) 9 Substituting uxfor the curved triangular void into Equation 32, we get Et xx E=12c2 t 2fc2 t+ 3f2c2 t+ 1(34) Doing the same for the curved diamond void, we get Ed xx E=(cd1) 3c2 d1 (cd1) (3c2 d1)f(3c3 d3c2 d+cd3)(35) By argument of symmetry, for the curved diamond, Exx=Eyy, so there is no need to repeat the exercise for loading in the y-direction. We will show in Section 5.7 that, for the curved triangle, Exx=Eyyas well. This is due to the fact that a regular triangular void (and all regular polygonal voids except for the square) is an isotropic shape [44]. Due to this, the \orthotropic" constitutive model we are deriving for the hexagonal mesostructure is in fact a transversely isotropic constitutive model. 5.5 Shear loading in y-z plane We repeat the process with 0 xy=0,0 xx=0 yy= 0. The Kolosov-Muskhelishvili potentials for the curved triangle void are =i0  =i0 4+ 2c2 t2ct+ 1 2ct3(36) and for the curved diamond void are =i0 (cd+ 1) =i0 4+ 43 cd+ 1 43c2 d+ (43)cd(37) Equation 28 simpli es to 1 22 0 GAvoid ft+1 20tZ2 0 uydy d uxdx d  d =1 22 0 GxyAvoid ft (38) Substituting in ux,uyfor the curved triangle void and rearranging yields Gt xy G=12c2 t f(+ 1)2c2 t+ 1(39) Doing the same for the curved diamond void yields Gd xy G=(cd+ 1) 13c2 d (cd+ 1) (13c2 d) +f(+ 1)(40) 5.6 Equibiaxial loading We repeat the process with 0 xx=0 yy=0,0 xy= 0. The Kolosov-Muskhelishvili potentials for the curved triangle void are =0 2ct0 22 =2 2c2 t+ 1 0 2ct3(41) 10 and for the curved diamond void are =0 2cd0 23 =3 3c2 d+ 1 0 3cd4(42) Equation 28 simpli es to 1 22 0(1) EAvoid ft+1 20tZ2 0 uxdy d uydx d  d =1 22 0(1xy) ExxAvoid ft (43) Solving for xyfor the curved triangle void yields t xy= 12c2 t (f+) f(2c2 t+ 3)2c2 t+ 1(44) and for the curved diamond void yields d xy= 3c3 d+ 3c2 d+cd1 f 3c3 d3c2 d+cd+ 1 f(3c3 d3c2 d+cd3)3c3 d+ 3c2 d+cd1(45) 5.7 Axial loading at an oblique angle At this point, we have derived all of the four in-plane elastic constants. However, it is worth further investigating the e ective Young's modulus Ewhen the mesostructures are uniaxially loaded at some arbitrary oblique angle (as measured positive counterclockwise from the positive x-axis). We rst transform our oblique angle uniaxial stress 0 =0into axis-aligned components: 0 xx=0cos2 0 yy=0sin2 0 xy=0cossin(46) The Kolosov-Muskhelishvili potentials for the curved triangle void are =0(cos(2) +isin(2)) 2+0 4ct0 42 =0 cos(2) 42c2 t2ct1 i sin(2) 4+ 2c2 t2ct+ 1 +i2 2c2 t+ 1 4ct23(47) and for the curved diamond void, they are =0 40((cd+ 1) cos(2)i(cd1) sin(2)) 2(c2 d1)cd0 43 =0 (cd+ 1) cos(2) 4+ 4+ 3 cd+ 1 2 (c2 d1) (3cd4) +0 (cd1) 2 3c3 d+ 3c2 d+cd+ 1 isin(2) 4+ 43 cd+ 1 2 (c2 d1) (3cd4)(48) The Young's modulus ratio for the curved triangle void is Et  E=12c2 t 2fc2 t+ 3f2c2 t+ 1(49) aod for the curved diamond void is Ed  E=3c4 d+ 4c2 d1 3c4 d+ 4c2 d1 +f(3c4 d2c2 d3)2fcdcos(4)(50) 11 −1.00 −0.75 −0.50 −0.25 0.00 0.25 0.50 0.75 1.00−1.00−0.75−0.50−0.250.000.250.500.751.00Eη E circle cusped triangle cusped diamondFigure 4: Polar plot of the Young's modulus ratio under oblique loading for a circular, cusped triangle (ct= 0:5), and cusped diamond ( cd= 0:33) voids at a porosity value of f= 0:1. Both the circle and the curved triangle void shapes are isotropic. We note that the results match those presented by Kachanov [44], and by setting = 0, we recover the results presented in Section 5.4. We also note that, as brie y mentioned in Section 5.4, the curved triangle void is an isotropic shape, and this is evident by the fact that the result in Equation 49 is the same as Equation 34. For the curved diamond void, the modulus ratio changes as a function of . As shown in Figure 4, the ratio is largest (for any given porosity value f) for the 45loading case and lowest for the 0/90cases. This implies that if the FFF mesostructure needs to be loaded in the x-y plane, the 45direction o ers the highest sti ness. 6 Shear loading in z-x and z-y planes In Section 5, we used the complex variable method of elasticity to solve the 2D elasticity problems around a single void. In this section, we will use the same method to solve the antiplane shear problem [52]. The antiplane shear problem assumes that all displacements in the x-y plane are zero, and the displacement in the z-direction is a function of x and y only, i.e. ux=uy= 0 uz=uz(x;y)(51) The stresses are given by xx=yy=zz=xy= 0 zx=G@uz @x zy=G@uz @y(52) 12 In the absence of body forces, the boundary value problem at hand becomes r2uz= 0 in zx=0 zx forjrj!1 zy=0 zy forjrj!1 n= 0 on S(53) The boundary condition on S, expanded out, is zxnx+zyny= 0 (54) Since both the real and imaginary parts of a holomorphic function independently satisfy Laplace's equation, and sinceuzis a real function, we can solve the boundary value problem by nding a holomorphic function h(m()) =hR+ihIde ned such that uz=hR (55) where=ei andm() =x+iyis the mapping function as de ned previously. Taking the derivative of h with respect to mgives us dh dm=@h @x@x @m+@h @y@y @m =@uz @xi@uz @y =zx Gizy G(56) which implies that zx=G<dh dm =Gdh dm+dh dm 2 zy=G=dh dm =iGdh dmdh dm 2(57) The normal components can be expressed as nx=dy ds =i 2dm dsdm ds ny=dx ds =1 2dm ds+dm ds(58) wheres=ei is an arc length variable along the void surface. Substituting Equations 57 and 58 into Equation 54 yields 0 =i 2Gdm dsdm dsdh dm+dh dm 2i 2Gdm ds+dm dsdh dmdh dm 2 =iG 2dh dmdm dsdh dmdm ds =iG 2d(hh) ds(59) 13 This means that, to satisfy the traction boundary condition on the void, we need hh= 2ihI= 0 (i.e. the imaginary part of hmust be zero). As before, we split the problem into nding a uniform solution and a perturbation solution due to the presence of a single void to get h=h0+h1 (60) where superscript (0) represents the uniform solution and (1) represents the perturbation. The uniform solution is given by h0() =0 zx Gm() +i0 zy Gm() (61) For the void to be traction free, =fh0(s) +h1(s)g= 0 (62) As before, we will represent h1using a Laurent series with negative power terms h1() =1X n=1an n(63) Substituting Equations 61 and 57 into Equation 62 yields 0 zx G=fm(s)g+0 zy G<fm(s)g+1X n=1=fanei g= 0 (64) Given a mapping function, we can solve for the unknown constants anby matching coecients in Equation 64. As before, only a nite number of coecients will be nonzero due to the simple form of the mapping function given by Equation 6 and 7. We start with nding Gzxby setting0 zx=0,0 zy= 0. For both the curved triangle and diamond void, his given by h=0 G +1  (65) The e ective modulus can be found by equating strain energies using Equation 22 1 22 0 GAvoid ft+1 20tZ2 0 uzdy d  =1 22 0 GzxAvoid ft (66) The e ective modulus ratios for the curved triangle and curved diamond are given by Gt zx G=12c2 t 2f2c2 t+ 1(67) Gd zx G=13c2 d 2f3c2 d+ 1(68) Due to symmetry in the curved diamond, Gd zx=Gd zy. We also observe that Gt zx=Gt zydue to thexyplane being isotropic for the curved triangle. 7 Unit cell simulations We simulate the e ective elastic constants by conducting volume average homogenization on unit cells using the procedure described by Sun [56]. Circular, cusped triangle ( ct= 0:5), and cusped diamond ( cd= 0:33) voids are used. For computational eciency, 2D simulations under plane stress conditions are used to nd Exx,Eyy,xy, andGxy. The remaining elastic constants are found using 3D simulations. Appropriate periodic boundary conditions are used in all cases. For material properties, E= 3500 MPa and = 0.33 are used, which represent the linearelastic material properties of polylactic acid (PLA), a commonly used FFF material. In the following plots, closed-form solutions are denoted by cfand simulation results are denoted bysimwhen appropriate. 14 0.0 0.2 0.4 0.6 0.8 1.0 Porosity f0.00.20.40.60.81.0Young's modulus ratio Ezz Ecircle (sim) triangle (sim) diamond (sim) cfFigure 5: Simulation results for Ezzmatch the closed-form solution for all porosity values and void shapes. 0.0 0.2 0.4 0.6 0.8 1.0 Porosity f0.00.10.20.30.40.5Poisson ratio νzx,νzycircle (sim) triangle (sim) diamond (sim) cf Figure 6: Simulation results for zxandzymatch the closed-form solution for all porosity values and void shapes. 15 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 Porosity f0.00.20.40.60.81.0Young's modulus ratio Exx E, Eyy Ecircle (sim) triangle (sim) diamond (sim) circle (cf) triangle (cf) diamond (cf)Figure 7: Young's modulus ExxandEyyas a function of porosity fand void shape. Results for Ezz,zx, andzyare shown in Figures 5 and 6. As expected, the numerical results match the closed-form solution exactly (apart from numerical error) for the entire range of porosity values. Results for ExxandEyyare shown in Figure 7, and the relative di erence between closed-form predictions and simulation results are shown in Figure 8. For the circular void, the closed-form solution is very accurate, di ering by 10% for porosity value around 35%. For the triangle and diamond void, 6% and 17% porosity lead to 10% di erence, respectively. This is a trend shown in all of the remaining elastic constants as well: the closed-form solution deviates from simulation faster than the diamond, which deviates faster than the circle. In addition, the sti ness ratio at a given porosity ratio is lowest for the triangle, followed by the diamond, with the circle having the highest sti ness ratio. Similar results have been observed by Jasiuk [47], who noted that the sti ness ratio is lower for polygonal holes with fewer sides. Results for xyare shown in Figure 9, and the relative di erence is shown in Figure 10. The results for the diamond void are interesting in that the closed-form solution predicts xyincreasing slightly from an initial value of 0 :33 and leveling o . This initial increase is captured by the simulation (see inset of Figure 9), but the simulations show that xydecreases towards 0 as porosity is increased for all void shapes. Intuitively, this makes sense; at high porosity values, there is a thin narrow band within the unit cell along the direction of loading that is under stress, while most of the the unit cell is relatively stress free, causing very little displacement transverse to the loading direction. Under 10% of di erence can be achieved by using porosity values less than 6%, 7%, and 15% for the triangle, diamond, and circular voids, respectively. However, we would like to note that in certain situations, such as uniaxial loading and bending of slender beams), the Poisson ratio plays a minimal role in the overall stress, strain, and displacement elds. In the aforementioned situations, having a less accurate Poisson ratio value may not matter. The results for Gxy,Gzx, andGzyare given in Figures 11 and 13, with relative di erences shown in Figures 12 and 14. The shear modulus results are similar to the Young's modulus result: the triangle void still has the lowest ratio, and the closed-form solution diverges from the simulation fastest for the triangle void. To achieve under 10% di erence for Gxy, porosity values less than 7%, 10%, and 11% should be used for the triangle, diamond, and circular voids, respectively. To achieve under 10% di erence for GzxandGzy, porosity values less than 8%, 14%, and 24% should be used for the triangle, diamond, and circular voids, respectively. 16 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 Porosity f01020304050Exx E, Eyy E relative difference (%)circle triangle diamondFigure 8: Relative di erence between closed-form solution and simulation for ExxandEyy. 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 Porosity f0.00.10.20.30.40.5Poisson ratio νxy circle (sim) triangle (sim) diamond (sim) circle (cf) triangle (cf) diamond (cf)0.00 0.05 0.10 f0.300.320.34νxy Figure 9: Poisson ratio xyas a function of porosity fand void shape. 17 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 Porosity f01020304050νxy relative difference (%) circle triangle diamondFigure 10: Relative di erence between closed-form solution and simulation for xy. 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 Porosity f0.00.20.40.60.81.0Shear modulus ratio Gxy Gcircle (sim) triangle (sim) diamond (sim) circle (cf) triangle (cf) diamond (cf) Figure 11: Shear modulus Gxyas a function of porosity fand void shape. 18 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 Porosity f01020304050Gxy G relative difference (%)circle triangle diamondFigure 12: Relative di erence between closed-form solution and simulation for Gxy. 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 Porosity f0.00.20.40.60.81.0Shear modulus ratio Gzx G, Gzy Gcircle (sim) triangle (sim) diamond (sim) circle (cf) triangle (cf) diamond (cf) Figure 13: Shear modulus GzxandGzyas a function of porosity fand void shape. 19 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 Porosity f01020304050Gzx G, Gzy G relative difference (%)circle triangle diamondFigure 14: Relative di erence between closed-form solution and simulation for GzxandGzy. One reason the hexagonal mesostructure performs worse than the square mesostructure is because the extruded roads in the hexagonal mesostructure become fully disconnected from neighbors at a lower porosity. We can conduct a geometric analysis to determine at what porosity the extruded roads become disconnected. For the cusped triangle void, the distance from the center to the void to one of the tips is 1 :5. The voids are arranged in a hexagonal array, with tips touching each other. The dimensions of the rectangular unit cell is 3p 39. Since there are the equivalent of four triangle voids per unit cell, the total void area in a unit cell (found using Equation 33) is 2 . This means that the maximum porosity for the hexagonal mesostructure is fmax t=2 3p 39 = 0:134(69) Similar analysis for the square mesostructure leads to fmax d=2=3 9 = 0:295(70) The close-pack nature of the hexagonal mesostructure essentially leads to worse sti ness-per-weight ratio than the square mesostructure. The only exception is loading in the x-direction, in which case they both perform the same for any porosity value. While these upper bounds on porosity assume ideal mesostructures and void shapes, real world parts printed with porosity values near (or exceeding) the upper bounds are expected to have many intersecting voids and bond defects{all of which lead to signi cantly degraded mechanical performance. The sharpness of the voids also has a signi cant e ect on the e ective sti ness, with rounded corners (lower absolute values of ctandcd) leading to higher modulus values. One method to round out the void shapes is to increase the negative gap between extruded roads, causing the material to over ow and ll in the voids. There is a limiting point: at very high negative gaps, the geometric distorsion can be extreme and the as-printed part geometry may di er signi cantly than the original input geometry [36], [57], [58]. 20 Table 1: Summary of e ective elastic constants. Elastic constant Curved triangle Curved diamond Ezz E1f 1f zx   zy   Exx E12c2 t 2fc2 t+3f2c2 t+1(cd1)(3c2 d1) (cd1)(3c2 d1)f(3c3 d3c2 d+cd3) Eyy E12c2 t 2fc2 t+3f2c2 t+1(cd1)(3c2 d1) (cd1)(3c2 d1)f(3c3 d3c2 d+cd3) xy(12c2 t)(f+) f(2c2 t+3)2c2 t+1(3c3 d+3c2 d+cd1)f(3c3 d3c2 d+cd+1) f(3c3 d3c2 d+cd3)3c3 d+3c2 d+cd1 Gxy G12c2 t f(+1)2c2 t+1(cd+1)(13c2 d) (cd+1)(13c2 d)+f(+1) Gzx G12c2 t 2f2c2 t+113c2 d 2f3c2 d+1 Gzy G12c2 t 2f2c2 t+113c2 d 2f3c2 d+1 8 Conclusion We have derived closed-form expressions for all nine e ective elastic constants for two mesostructures: hexag- onal and square. The e ective elastic constants depend on the void shape, porosity, and material properties of the bulk material. Simulation results for several void shapes match the closed-form expressions very well at low porosity values, with the largest error coming from the xyprediction. A summary of e ective moduli ratios as a function of porosity, shape parameters, and material Poisson ratio is given in Table 8. Funding Ruiqi Chen is supported by the National Defense Science and Engineering Graduate (NDSEG) as well as the Stanford Graduate Fellowship (SGF) programs. Con icts of Interest The authors declare no con ict of interest.