In this paper we introduce a computational framework for optimizing 3D prestressed concrete structures. The procedure combines shape optimization of post-tensioning cables with topology optimization of concrete, enabling novel design configurations of 3D beams and plates. Prestressing cables are modeled as B-spline curves that are embedded within a concrete continuum and whose control points serve as design variables. Concurrently, the distribution of concrete is determined by density-based topology optimization. A special filter enforces concrete cover around the cables, creating a coupling between the two sets of design variables. The optimization problem formulation mimics the design intent of prestressed concrete structures, separating the load-balancing part with respect to permanent loads from stiffness maximization with respect to live loads. Results show that material savings exceeding 30% can be achieved compared to solidprestressedmembers,strengtheningtheargumentinfavorofusingprestressedconcreteforreducingembodiedCO 2. The complete implementation in MATLAB is provided as supplementary material. Keywords: topology optimization, shape optimization, prestressed concrete, embodied CO 2 1. Introduction Prestressed concrete is used extensively in new infras- tructure,e.g.roadandrailwaybridges,aswellasinbuild- ings with relatively large spans. The use of prestressing steel increases the durability of concrete structures due to thereduction,orevenelimination,oftensilestresses. Pre- stressing also reduces the sizes of concrete members that areutilizedmoreefficientlythantheirreinforcedconcrete counterparts,providinganeconomicadvantageaswellas the capability to design for larger spans. Optimization of concrete structures, and prestressed onesinparticular,isgaininginterestfortwomainreasons. Oneisthesignificantcontributionofcementproductionto global CO 2emissions (Andrew, 2019). Prestressed con- crete structures have an advantage w.r.t. embodied CO 2 (Milleretal.,2015;BroylesandHopper,2023)andformal optimizationcanfurtherreducethequantitiesofrawmate- rials and hence also the embodied CO 2. Another motiva- tionforoptimizingprestressedconcreteisthegrowingin- terestin3Dconcreteprinting,withprestressingbeingone of the methods that is being explored for joining printed partsandresistingbending(Vantyghemetal.,2020;Ooms et al., 2022; Li et al., 2024). These motivations call for formulating and developing shape and topology optimization procedures that are tai- lored specifically for prestressed concrete. The challenge isthattheshapeoftheprestressingcablesandthetopology of the concrete structure need to be optimized simultane- ously. Thecablestransferequivalentforcestotheconcrete, that depend on the shape and should counteract external loads–indicatingacouplingbetweenshapeandtopology. Moreover, in most practical scenarios, the cables are cov- ered by concrete – indicating another coupling between shape and topology .Simultaneous shape and topology optimization of pre- stressed structures is a relatively new research topic. It was approached by Amir and Shakour (2018a) where 2D beamswereoptimizedforstiffness,withasinglecablethat was parameterized by a set of linear segments. Later, the prestressing force was added as a design parameter to be minimizedandB-splineswereemployedformodelingthe cable,allowingtoimposecurvatureconstraints(Amirand Shakour, 2018b). More recently, Zhang et al. (2021) re- formulatedthedesignparameterizationusingsplinebases for the cable and the continuum concrete, and introduced stress constraints which are necessary in realistic applica- tions. Anotherrelatedworkpresentedtopologyoptimiza- tion of concrete girders with prestressed shape memory alloys,thatholdtheadvantageofapplyinglocalprestress- ing effects (Sung and Andrawes, 2023). Finally, Jha and Pathak(2023)presentedconcurrentshapeoptimizationof aconcretebeamanditspost-tensioningcable. Theauthors considered also frictional losses which have not been ac- counted for in previous studies. In the current work, we present a framework for op- timizing 3D concrete structures with 3D post-tensioning cables. The purpose is to allow complete design freedom in 3D that can generate innovative design concepts. Our formulation enables direct design of 3D beams, eliminat- ingmanualinterpretationof2Doptimization(Vantyghem et al., 2020; Ooms et al., 2022). Another exciting appli- cation is the optimization of post-tensioned plates, where research has shown that optimal cable layouts are curved in 3D (e.g., Sarkisian et al., 2018; Zelickman and Amir, 2021) and the interaction with topology optimization of concrete has yet to be explored. The remainder of the article is organized as follows. Modeling and design parametrization of both concrete and prestressing steel are described in Section 2. The Preprint submitted to Elsevier May 27, 2024 optimization problem and subsequent sensitivity analysis arepresentedinSections3and4. Severaldesignexamples are presented and discussed in Section 5 and concluding remarks are provided in Section 6. 2. Modeling and design parameterization In this section, we present the methodology adopted inthecurrentpaperformodelingprestressedconcreteele- mentsandprovidedetailsontheirdesignparameterization. The prestressed elements are modeled as a 3D continuum discretized using tri-linear cube finite elements, with em- bedded cables that are parameterized by B-splines. The optimizationisformulatedtosimultaneouslyoptimizethe topology of the concrete, that is controlled by pseudo- density design variables (Bendsøe, 1989; Bendsøe and Sigmund,2003),andtheshapeofthecables,thatisdeter- mined by the splines’ control points. Specifically, we use quadratic B-spline curves C(𝜉)=𝑁𝑐𝑝∑︁ 𝑖=1P𝑖𝑁𝑖,2(𝜉) where C(𝜉)is the quadratic three-dimensional B-spline curve;𝜉is the parameter along the curve; P𝑖and𝑁𝑖,2(𝜉) represent the control point 𝑖of the curve and its corre- sponding quadratic B-spline basis function, respectively; and𝑁𝑐𝑝isthenumberofcontrolpoints. AlltheB-splines curves in this paper are defined over open uniform knot vectors. The spatial coordinates of the curve’s control points serve as shape design variables in the optimization problem. An illustration of a continuum domain with a single embedded cable is depicted in Figure 1. 2.1. Pseudo-density representation Thedensitydesignvariablesinthispaperaresubjected to three consecutive filters to prevent numerical instabili- ties (Sigmund and Petersson, 1998) and meet specific ge- ometricalrequirements. Thefirstisadensityfilter(Bruns and Tortorelli, 2001; Bourdin, 2001) which is formulated asaPDE(LazarovandSigmund,2011). Initsdiscretized form, we solve  𝑙2 𝑜K𝑓+Mvol 𝑓+𝑙𝑠Msurf 𝑓 ˜𝝆=T𝝆 (1) where 𝝆is the vector of density design variables; ˜𝝆is the vector of their filtered values; 𝑙𝑜is the bulk length scale parameterthatisrelatedtothetraditionalfilterradius 𝑟min; and𝑙𝑠isthesurfacelengthscaleparameterthatguarantees consistentfilteringontheboundaries(Wallinetal.,2020). Further details on the implementation of Eq. (1) and the finite element assembly of K𝑓,Mvol 𝑓,Msurf 𝑓andTcan be foundinLazarovandSigmund(2011);Wallinetal.(2020). In the current context, we impose 𝑙𝑠=𝑙𝑜on faces whose materialdistributionisfree,and 𝑙𝑠=0onfacesthateither contain prescribed solids or belong to symmetry planes. As commonly done in topology optimization procedures, Eq. (1) is formulated based on the same finite element discretization that is used for solving the displacements. The second filter is applied to ensure an adequate con- crete cover for the prestressing cables. In the majority ofapplications involving prestressed elements, it is essential fortheprestressingcablestobefullyembeddedwithinthe concrete. In the numerical model, this is materialized by ensuring that the B-spline curve representing the cable is covered by finite elements with solid material. For this purpose, we adopt the filter suggested in Amir and Shak- our (2018a,b) and extend it to accommodate 3D designs. In the remainder of this paper, we refer to this filter as the cable-to-concretefilter. Byapplyingthefilter,concreteel- ements located within a specified radius around the cable are pushed to the material phase using a super-Gaussian function. The influence region of this filter is illustrated in Figure 2. Thecable-to-concretefilterisappliedtothefilteredden- sities ˜𝝆, while each ˜𝜌𝑖is updated as follows, ˆ𝜌𝑖=˜𝜌𝑖+(1−˜𝜌𝑖)𝑒−1 2𝑑𝑖 𝑟𝑐𝑓𝜇 (2) where ˆ𝜌𝑖represents the modified density of concrete el- ement𝑖;𝑑𝑖is the minimal distance between the concrete element𝑖and a cable segment; 𝑟𝑐𝑓is the radius of the cable-to-concrete filter; and 𝜇represents the sharpness of thesuper-Gaussianfunction. Noticethatthisfiltercouples thedensityandshapedesignvariablesthroughthedistance 𝑑𝑖. While the control points representing the cables are adjusted, the distance changes, affecting the resulting fil- tereddensities. Aplotofthesuper-Gaussianfunctionwith respecttothedistance 𝑑𝑖,for ˜𝜌=0.5,𝑟𝑐𝑓=4,andvarious values of𝜇is displayed in Figure 3a. As the sharpness parameter𝜇increases, the region surrounding the cable becomes more crisply defined. Notice that when the dis- tanceisbelow 𝑟𝑐𝑓,thedensityispushedto1,whilewhen thedistanceisgreaterthan 𝑟𝑐𝑓,thedensitytendstoremain the same. For further demonstration of the effect of the filter,aplotfordifferentvaluesof ˜𝜌ispresentedinFigure 3bandanillustrationofa2Dcross-sectionofacablewith its surrounding elements is presented in Figure 4. Finally, to ensure a clear separation between the void and material regions, the densities ˆ𝝆are subjected to a smooth Heaviside projection function (Guest et al., 2004; Xuetal.,2010). Specifically,thisprojectionisformulated following the ‘robust’ topology optimization formulation (Wang et al., 2011; Lazarov et al., 2016) that allows to consider worst-case scenarios. Herein, we employ the erodedlayoutforstiffnessevaluationandthedilatedlayout for volume evaluation, even though our stiffness measure is not strictly a compliance functional hence the eroded layout is not strictly the worst case. The element-level eroded and dilated densities are defined as (Wang et al., 2011) 𝜌ero𝑖=tanh(𝛽𝐻𝑆𝜂ero)+tanh(𝛽𝐻𝑆(ˆ𝜌𝑖−𝜂ero)) tanh(𝛽𝐻𝑆𝜂ero)+tanh(𝛽𝐻𝑆(1−𝜂ero))(3) 𝜌dil 𝑖=tanh(𝛽𝐻𝑆𝜂dil)+tanh(𝛽𝐻𝑆(ˆ𝜌𝑖−𝜂dil)) tanh(𝛽𝐻𝑆𝜂dil)+tanh(𝛽𝐻𝑆(1−𝜂dil))(4) where𝛽𝐻𝑆governs the sharpness of the smooth projec- tion;𝜂erois the projection threshold for the eroded layout (e.g.,𝜂ero=0.6);and𝜂dilistheprojectionthresholdforthe dilated layout (e.g., 𝜂dil=0.4). The intermediate layout, with𝜌intwhichisconsideredastheactualdesignintended 2 Figure 1: An illustrative model of a a continuum domain with a single embedded cable. The cable shape is outlined in red, its four control points are marked by blue dots and the control polygon is outlined by a black polyline. Figure 2: The influence region of the cable-to-concrete filter. 0 2 4 6 8 1 0 d i0 . 50 . 60 . 70 . 80 . 91^ ; iS u p e r - G a u s s i a n r c f = 4 ~ ; i = 0 : 57 = 2 7 = 4 7 = 6 7 = 8 7 = 1 (a) 0 2 4 6 8 1 0 d i00 . 20 . 40 . 60 . 81^ ; iS u p e r - G a u s s i a n ~ ; i = 0 : 2~ ; i = 0 : 5~ ; i = 0 : 8 (b) Figure3: Thesuper-Gaussianfunctionusedforthecable-to-concretefilter. Left: ˜𝜌=0.5,𝑟𝑐𝑓=4anddifferentvaluesof 𝜇;Right:𝑟𝑐𝑓=4,𝜇=8and different values of ˜𝜌. for manufacturing, is evaluated similarly using 𝜂int=0.5. Wenotethattheresponseaccordingtothislayoutdoesnot need to be evaluated during the optimization. However, thevolumeofthislayoutisevaluatedtoadjustthevolume constraint, as will be discussed in Section 5. After implementing all filters, the structural response is evaluated over the eroded layout with 𝝆ero, where the Young’s modulus of each concrete element is determined by a Modified SIMP interpolation scheme (Sigmund and Torquato, 1997), 𝐸𝑖(𝜌ero 𝑖)=𝐸𝑚𝑖𝑛+(𝐸𝑚𝑎𝑥−𝐸𝑚𝑖𝑛)𝜌ero 𝑖𝑝𝐸(5) where𝐸𝑚𝑖𝑛is a relatively small positive number in orderto avoid singularity of the stiffness matrix; 𝐸𝑚𝑎𝑥is the actual value of Young’s modulus for concrete; and 𝑝𝐸is a penalization factor that drives the design towards a 0-1 distribution. 2.2. Prestressing forces The prestressed elements in this paper are subjected to four different forces: prestressing forces, self-weight, external dead loads and external live loads. The nodal force vectors associated with these forces are denoted as f𝑝,f𝑠,f𝑑, and f𝑙, respectively. Dead and live loads are applied to the prestressed structural elements in specified regions based on the problem settings. The self-weight is 3 d ir c fFigure 4: An illustration of a 2D cross-section of a cable with its sur- rounding elements. Concrete elements that should be pushed to the materialphasearemarkedwithbluecrosses,whileelementsthatshould not be affected are marked with blue dots. evaluatedaccordingtotheerodedlayoutwhiletheweight of each finite element is applied to its associated degrees offreedom. Foraspecificfiniteelement 𝑖,thelocalnodal force vector due to self-weight is evaluated as follows, f𝑠,𝑖=𝛾𝑐𝑣𝑖𝜌ero 𝑖N𝑖, where𝛾𝑐isthespecificweightofconcrete; 𝑣𝑖isthevolume ofthefiniteelement;and N𝑖representstheshapefunctions of the finite element. Theprestressingforcesaremodeledasdistributedforces acting along the cable in its normal direction due to the cable’scurvature,aswellaspointloadsattheanchorages. These forces are denoted as FCurv 𝑝andFAnc 𝑝, respectively. A demonstration of is presented in Figure 5. To evaluate theforcesandsubsequentlycomputetheprestressingnodal force vector f𝑝, the following steps are performed: 1. Theintersectionsofthecablewiththefiniteelement grid are found. Thanks to the second-order B-spline representationofthecable,theintersectionpointscan bedeterminedanalyticallywith 𝜉beingtheunknown variable. Specifically, we solve the equations 𝑥(𝜉)= 𝑥𝐺𝑟𝑖𝑑,𝑦(𝜉)=𝑦𝐺𝑟𝑖𝑑, and𝑧(𝜉)=𝑧𝐺𝑟𝑖𝑑to find the intersection points with the 𝑦−𝑧,𝑥−𝑧, and𝑥−𝑦 planes, respectively. 2. For each cable segment between two intersection points, FCurv 𝑝is defined: FCurv 𝑝(𝜉)=𝑇𝑝˜𝜅ˆn, where ˜𝜅is the average curvature between the two intersection points, ˜𝜅=(𝜅(𝜉𝑠)+𝜅(𝜉𝑒))/2;𝜉𝑠and𝜉𝑒 arethevaluesof 𝜉attheintersectionpoints; 𝑇𝑝isthe prestressing force in the cable; and ˆnis the normal vector of the cable. The curvature and the normal vector as a function of 𝜉are evaluated as follows, 𝜅(𝜉)=|C′×C′′| |C′|3and ˆn(𝜉)=C′×(C′′×C′) |C′||C′′×C′|. A demonstration of FCurv 𝑝within a single finite ele- ment is presented in Figure 6. 3. Within each finite element, the force, FCurv 𝑝, is pro- jected to the finite element’s nodes using the ele- ment’s shape functions N𝑖, fCurv 𝑝,𝑖=∫𝜉𝑒 𝜉𝑠FCurv 𝑝N𝑖𝐽𝑑𝜉 with𝐽(𝜉)representing the Jacobian of the curve which is calculated by 𝐽(𝜉)=√︃ 𝑥(𝜉)′2+𝑦(𝜉)′2+𝑧(𝜉)′2 andNwith respect to 𝜉is defined as, N𝑖(𝜉)= N𝑖(C)=N𝑖(𝑥(𝜉),𝑦(𝜉),𝑧(𝜉)). 4. The loads at the anchors are applied as point loads at the cable’s ends. These loads are equal to the prestressing force and act in the direction tangential to the cable, FAnc 𝑝(𝜉ends)=𝑇𝑝t(𝜉ends), where 𝒕is the tangential direction of the cable and canbeevaluatedby C′/|C′|. Subsequently,thenodal prestressing forces at hosting elements of the cable’s endpoints are evaluated by fAnc 𝑝,𝑖=FAnc 𝑝N𝑖. 5. Finally, the prestressing force vector f𝑝is assembled as follows, f𝑝=𝑁𝐸∑︁ 𝑖=1 fAnc 𝑝,𝑖+fCurv 𝑝,𝑖 where𝑁𝐸is the number of finite elements. To facilitate the subsequent presentation and discus- sions regarding displacements, we denote the displace- ments caused by f𝑝,f𝑠,f𝑑, and f𝑙asu𝑝,u𝑠,u𝑑, and u𝑙, respectively. Additionally, the combined displacements resultingfrommultipleforcesarerepresentedbycombin- ing their subscripts. For instance, the displacements due tof𝑠+f𝑑aredenotedas u𝑠𝑑andthedisplacementsarising from f𝑠+f𝑑+f𝑝are denoted as u𝑠𝑑𝑝. We note that the stiffness of the cable and force losses along the cable are neglected. Nevertheless, the formu- lation is capable of accommodating immediate losses re- sulting from friction, anchorage slip and elastic shorten- ing, given their explicit dependence on the prestressing force,thecurvatureofthecableandconcretedeformation. Furthermore, long-term losses due to creep, shrinkage, andrelaxationcanbeevaluatedbasedoncableforcesand concrete deformation, and can be incorporated into the computation of prestressing forces. 4 Figure 5: A demonstration of the prestressing forces acting on a three-dimensional beam model. FCurv𝑝represents the prestressing forces due to the curvature of the cable while 𝐹Anc𝑝represents the prestressing forces at the anchors. Figure 6: A demonstration of 𝑭Curv𝑝within a single finite element. The blue dots indicate intersection points between the cable and the finite element. 3. Problem formulation The optimization problem is formulated to simultane- ously optimize the topology of concrete and the shape of the embedded cables. The topology is optimized follow- ing the density-based approach (Bendsøe, 1989; Bendsøe and Sigmund, 2003). The shape optimization is simul- taneously conducted by adjusting the positions of the B- splines’ control points. The objective of the optimization involvesminimizingcomplianceanddisplacement-related terms, while limiting the amount of concrete in the opti- mized design. The objective function in this paper has been uniquely developed to incorporate a commonly adopted concept in the design of prestressed elements. In prestressed con- cretedesign,theprestressingforcesbalancethestructural response due to self-weight and permanent dead loads, while live loads are subsequently applied to the already- balanceddesign. Toalignwiththisconcept,weformulate theobjectivefunctiontominimizetwoseparatemeasures. The first part of the objective is defined as 𝜙1=(l◦u𝑠𝑑𝑝)𝑇(l◦u𝑠𝑑𝑝), (6)and aims at minimizing the absolute displacements due to self-weight, permanent dead loads and prestressing – henceachievingabalancedstructuralresponse. InEq.(6), the displacements to be balanced are designated by the degreesoffreedomindicatedbythevector landthesymbol ◦indicates component-wise multiplication. Typically, l will refer to the primary loading surface of the structure, where also the live loads are applied. Without loss of generality,hereinthevector lcorrespondstoverticalDOF in the top surface of the prestressed structure. The second part of the objective is defined as 𝜙2=f𝑇 𝑙u𝑙 and aims at minimizing the compliance due to live loads. The key difference between the two parts of the objective isthat𝜙2drivestheconcrete’stopologytowardsthestiffest design under live loads, while 𝜙1is not directly related to stiffness. Byminimizing 𝜙1,thedesignshouldaccommo- dateseveralsetsofloadsandreachaneutraldisplacement field. Previousworkhasshownthatwhentheprestressing force is too large, the topology adapts itself by soften- ing parts of the structure to reach neutral displacements (Amir and Shakour, 2018a). One way to avoid this phe- 5 nomenonistominimizealsotheprestressingforce(Amir andShakour,2018b). Herein,wechoosetominimizetwo separate measures, a choice which is aligned with prac- tical preliminary design of prestressed structures. At the same time, a reasonable value for the prestressing force is determined prior to the optimization. We note that in principle, one can choose to apply only a certain portion ofthedeadloadsin 𝜙1andallowupwardsdeflections,and subsequently add the remainder of the dead loads to 𝜙2. This depends on design p