Topology Optimization and 3D printing of Large Deformation Compliant Mechanisms for Straining Biological Tissues
TTopology Optimization and 3D-printing of Large DeformationCompliant Mechanisms for Straining Biological Tissues
P. Kumar ∗ , , C. Schmidleithner † , N. B. Larsen † , and O. Sigmund ∗ ∗ Department of Mechanical Engineering, Solid Mechanics, Technical University of Denmark, 2800 Kgs.Lyngby, Denmark † Department of Health Technology, Technical University of Denmark, 2800 Kgs. Lyngby, Denmark
Abstract:
This paper presents a synthesis approach in a density-based topology optimizationsetting to design large deformation compliant mechanisms for inducing desired strains in bi-ological tissues. The modelling is based on geometrical nonlinearity together with a suitablychosen hypereleastic material model, wherein the mechanical equilibrium equations are solvedusing the total Lagrangian finite element formulation. An objective based on least-square errorwith respect to target strains is formulated and minimized with the given set of constraints andthe appropriate surroundings of the tissues. To circumvent numerical instabilities arising due tolarge deformation in low stiffness design regions during topology optimization, a strain-energybased interpolation scheme is employed. The approach uses an extended robust formulation i.e.the eroded, intermediate and dilated projections for the design description as well as variation intissue stiffness. Efficacy of the synthesis approach is demonstrated by designing various compli-ant mechanisms for providing different target strains in biological tissue constructs. Optimizedcompliant mechanisms are 3D-printed and their performances are recorded in a simplified ex-periment and compared with simulation results obtained by a commercial software.
Keywords:
Topology Optimization; Biological Tissue; Compliant Mechanisms; 3D printing;Stereolithography; Flexible Poles Method
Development of new drugs is challenged by the limited predictive accuracy of current simplecell models on safety and efficacy in the human body (Mordwinkin et al., 2013). Functionalmini-organ models with higher predictive value are increasingly used in the pharmaceuticalindustry to meet this challenge (Ikeda et al., 2017). These mini-organ models can be derivedfrom a healthy/diseased person’s tissues, adult stem cells (which can be differentiated into theparticular type of tissues in vitro relatively faster and in few steps), human embryonic stemcells (hESCs) or by inducing pluripotency in human adult cells (hiPSCs) (Duelen and Sam-paolesi, 2017). To facilitate maturation of differentiated tissue cells, local static- and dynamic-mechanical forces are essential to induce the required strains (Vining and Mooney, 2017). Ingeneral, uni-axial stretching up to 15 −
20% in mini-organs (e.g., skeletal- and cardio-myocytes)is needed (Vandenburgh et al., 1995) for achieving alignment and proper contractile behavior.However, current available culture systems are limited in providing the needed strains (Cooket al., 2016; Riehl et al., 2012) and, in general, are designed by trial and error approaches. Corresponding author: [email protected], [email protected] Ethical issues a r X i v : . [ c s . C E ] A ug he motive herein is to provide a systematic approach, e.g., topology optimization, while con-sidering geometric nonlinearity to design largely deformable compliant mechanisms which canfurnish the required strains in the biological tissues in response to external stimuli. The non-linearity stems from the large desired strain the biological tissues.A compliant mechanism (CM) performs its tasks using motion obtained from the elastic de-formations of its members. These mechanisms, in general, have monolithic designs and theiramount of elastic deformations (small/large) depend upon the applications they are designedfor. By virtue of their geometrical features, CMs offer numerous advantages over their classicalrigid body counterparts, such as low manufacturing and assembly cost, less frictional lossesdue to absence of joints, low wear and tear, high precision and repeatability, to name a few(Frecker et al., 1997; Sigmund, 1997). Thus, the use of such mechanisms designed by topologyoptimization is continuously increasing in a wide variety of applications e.g. path generation(Kumar et al., 2019b; Pedersen et al., 2001; Saxena and Ananthasuresh, 2001), displacementdeliminator (Saxena, 2013), MEMS (Ananthasuresh et al., 1994; Jonsmann et al., 1999), inbiomedical/ biomechanics/drug-discovery (Frecker et al., 2005; Kollimada et al., 2017; Kumaret al., 2019a). To add to the list of their ever expanding applications, herein, we propose anapproach using topology optimization in the nonlinear continuum regime to design compliant(micro-)mechanisms which can induce programmable strain up to 20% in biological tissues.Topology optimization (TO) relocates material in an optimum fashion within a prescribed designdomain by extremizing desired objective(s) with a given set of constraints (Sigmund and Maute,2013). In a general structural TO framework, finite elements (FEs) are employed to describethe design domain, and each FE is assigned a design (density) variable ρ . ρ i = 1 implies, i th FE is in solid phase, whereas ρ i = 0 represents its void state. In a gradient-based TO,FEs with 0 < ρ < < ρ <
1) and impose length scale in the final solutions, penalization and a robust formulationwith Heaviside projection filter (Wang et al., 2011) is adopted in this paper.In large deformation TO as considered here, FEs with low stiffness are prone to undergo ex-cessive distortion/deformation and thus, cause numerical instabilities which may jeopardize theprogress of optimization. To circumvent these numerical instabilities, Buhl et al. (2000) modifiedthe Newton-Raphson convergence criterion in their approach by excluding the internal nodalforces originating from low stiffness FEs. Bruns and Tortorelli (2001) treated such instabilitiesvia removing and reintroducing low stiffness FEs and others suggested only to remove elementsin low density regions (Cho and Kwak, 2006). Yoon and Kim (2005) proposed a connectiv-ity parameterization approach by employing fictitious springs to connect FEs. Klarbring andStr¨omberg (2013); Lahuerta et al. (2013) and Luo et al. (2015) employed special hyperelasticmaterial laws for low stiffness FEs in their approaches. An approach based on scaling of theelement deformation was proposed by van Dijk et al. (2014). Wang et al. (2014) proposed amethod based on an energy interpolation scheme that models low density elements as linear.Based on our extensive experience from other applications, we find that treating low densityFEs as linear (Wang et al., 2014), provides stable and reliable convergence and is hence usedhere.In general, a bioreactor for muscle tissue maturation is expected to induce uni-axial cyclic strainup to 20% for cellular alignment and proper contractile behavior in addition to supporting therequired auxotonic resistance (Vandenburgh et al., 2008). Numerous bioreactors for static uni-axial straining via auxotonic resistance have been reported. In those, a bioreactor with flexiblepoles approach (Vandenburgh et al., 2008) requires minimum assembly, provides a method forhigh throughput solutions supporting the formation of 3D microsized tissue constructs, all ina low cost format. The mechanical challenge of stably attaching slippery microtissues to the2exible poles is solved by making the tissue form by self-aggregation of muscle cells seeded in acompliant gel matrix between and around the poles (Hansen et al., 2010; Vandenburgh et al.,2008). This is a well-established procedure that was also recently demonstrated in fully 3Dprinted devices (Christensen et al., 2020). However, there is no medium for inducing cyclicstrains mechanically in the tissue constructs. Indeed, this shortcoming of the approach moti-vated us to formulate the presented design problem. Compliant (micro-)mechanisms are en-visioned to be combined in a design with flexible poles environment for inducing the requiredstrains in the constructs mechanically in response to external mechanical stimuli applied on abase platform (see Sec. 4.2). The approach is conceptualized to support many pairs of flexiblepoles requiring different straining in their respective biological tissues.Development (amount and location) of stress/strain in a structure significantly depends uponits geometrical features and the loading conditions. In TO, one can find approaches that imposestress constraints (Duysinx and Bendsøe, 1998; Luo et al., 2013; da Silva et al., 2019), and also,on stress isolation at predefined regions within the design domain (Li and Wang, 2014; Luoet al., 2017; Picelli et al., 2018). Our aim here is not to isolate strain (Picelli et al., 2018) inthe design domain but rather to achieve a target strain level in the biological tissue substratesusing optimized compliant mechanisms. A least square objective (see Sec. 3.1) based on targetstrains is minimized with a given set of resource constraints. Geometric nonlinearity is consid-ered wherein an energy interpolation scheme (Wang et al., 2014) is exploited to handle numericalissues arising due to large deformations in the low stiffness regions during the topology opti-mization. Here, the robust formulation (Wang et al., 2011) is employed in a large deformationdensity-based TO setting, wherein the maximum value of the three objectives, evaluated forthe dilated, intermediate and eroded designs, is minimized. Further, the robust formulationis suitably modified to also accommodate varying biological tissue geometries with the flexiblepoles surrounding (Vandenburgh et al., 2008). Prototypes of the designs are printed and theirperformances are compared with corresponding ABAQUS analyses.In summary, the contributions of this paper are: • Formulation of the topology optimization problem for designing compliant mechanismswhich can induce desired strains in biological tissues in response to external loading, • Conceptualization of an objective based on least-square error which is minimized to achievethe optimum material layout for the compliant actuators, • Illustration of the robust approach in a large deformation topology optimization settingfor designing CMs to strain tissue constructs, which is further modified to cater for tissueconstruct variations, • Demonstration of the approach by synthesizing various compliant mechanisms for strain-ing biological tissue with and without flexible poles environment, • Comparison of performances of the optimized mechanism designs with a commercial soft-ware, ABAQUS, • Realization of the optimized mechanisms by 3D-printing and performing a simplified ex-periment, and also, comparing the experimental result with corresponding ABAQUS anal-ysis result.The layout of this paper is arranged as follows. Section 2 describes the problem definitionwith optimization formulation. For the sake of completeness, we present a brief description ofnonlinear finite element analysis and the energy interpolation formulation (Wang et al., 2014).3he objective based on the target strains and its sensitivity analysis with energy interpolationmethod are presented in Section 3. Section 4 presents two set of numerical examples withdiscussions. Optimized mechanisms with the flexible poles setting are fabricated and an exper-iment is performed, and their performances are compared with respective ABAQUS analyses.Lastly, in Section 5, conclusions are drawn.
A schematic diagram of the design problem is illustrated in Fig. 1. The aim is to obtain anoptimized mechanism using TO in Ω which can induce the target strains (cid:15) ∗ (up to 10-20%) inthe biological tissue Ω BT in response to external actuation F in (Fig. 1). A density-based topology Ω BT Bio-tissue F in Ω Fixed k s ǫ ∗ Target strain
Figure 1: A schematic representation of the micro-actuator design problem. Ω , Ω BT indicate the design domainand region for a biological tissue, respectively. F in denotes the actuating force and k s is the input spring stiffness. optimization approach (Bendsøe and Sigmund, 2003; Sigmund, 2001) is adopted wherein eachFE is assigned a design variable ρ e ∈ [0 , ρ e of an FE is defined using the smooth Heaviside projectionfilter (Wang et al., 2011) as¯ ρ e ( ˜ ρ e ( ρ e )) = tanh( βη ) + tanh( β ( ˜ ρ e − η ))tanh( βη ) + tanh( β (1 − ˜ ρ e )) , (1)where β ∈ [1 , ∞ ) defines the steepness of the projection filter and η ∈ [0 ,
1] is a thresholdparameter for ¯ ρ e , and ˜ ρ e is the filtered form of ρ e . Ideally, β → ∞ , for a discrete (0 − β is increased from an initialvalue β int = 1 to a maximum value β max using a continuation strategy (see Sec. 4 for specificvalues).The filtered variable ˜ ρ e is evaluated as˜ ρ e = (cid:88) i ∈ n e w ( x ) v i ρ i (cid:88) i ∈ n e w ( x ) v i , (2)where n e = { i, || x ci − x ce || ≤ r min } with x ci and x ce as center coordinates of the i th and e th elementsrespectively, r min is the filter radius and || . || denotes distance in the Euclidean space. v i denotesthe volume of i th element, and w ( x ), a linearly decaying weighting function, is defined as w ( x ) = r min − || x ci − x ce || . (3)4e use the modified SIMP (Solid Isotropic Material with Penalization) interpolation scheme torelate the physical density ¯ ρ e with the Young’s modulus of the given material as E e ( ¯ ρ e ) = E v + ( E s − E v )( ¯ ρ e ) p , ¯ ρ e ∈ [0 , , (4)where E s and E v are Young’s moduli of the actual and void material, respectively. E v is setto E s × − and the penalty parameter p = 3 is chosen, which guides topology optimizationresults towards close to “0-1” solutions. To avoid checkerboards, one-node connected hinges, mesh-dependencies, and other artifacts,Sigmund (2009) formulated the design problem in a robust way wherein he employed two pro-jection filters. Wang et al. (2011) modified the formulation using the smooth Heaviside pro-jection filters, which is considered herein with suitable modifications for the large deformationcontinuum setting. The formulation in (Wang et al., 2011) considered a set of three designsi.e. dilated, intermediate and eroded continua for a problem and minimized the worst objectiveobtained out of three designs. Dilate and erode are morphological-based image operators. Asper (Sigmund, 2007), they can be used in a TO setting for controlling the feature sizes andensure robustness. Designs obtained using these operation in association with suitable filteringare called dilated and eroded structures. The dilate operation corresponds to under-etching ofthe fabricated designs, whereas erode operation corresponds to over-etching. The intermediatedesigns indicate the desired or correctly etched structures, i.e., the blueprint structures.The dilated ¯ ρ d , intermediate ¯ ρ i and eroded ¯ ρ e design vectors are obtained via Eq. 1 using thethreshold 0 . − ∆ η, . . η , respectively (see Sec. 4 for specific ∆ η ). The optimizationproblem is formulated in the nonlinear continuum setting as a min/max problem (Wang et al.,2011) to also accommodate different geometries of the tissue construct, which can be writtenas: min ρ : max k : (cid:0) f k ( u k , ¯ ρ d ( ρ )) , f k ( u k , ¯ ρ i ( ρ )) , f k ( u k , ¯ ρ e ( ρ )) (cid:1) s.t. : R k ( u k , ¯ ρ d ( ρ )) = , k = 1 , , · · · , N BT : R k ( u k , ¯ ρ i ( ρ )) = , k = 1 , , · · · , N BT : R k ( u k , ¯ ρ e ( ρ )) = , k = 1 , , · · · , N BT : V f ( ¯ ρ d ( ρ )) − V ∗ d ≤ ≤ ρ ≤ , (5) where R k ( u k , ¯ ρ d ( ρ )), R k ( u k , ¯ ρ i ( ρ )) and R k ( u k , ¯ ρ e ( ρ )) are the residual terms (Eq. 8) for the di-lated, intermediate and eroded designs, respectively, f k is the formulated objective (see Sec. 3.1)and N BT is the number of variations of the tissue construct. V f ( ¯ ρ d ( ρ )) and V ∗ d are the volumefraction and its upper limit for the dilated design, respectively. Volume of the dilated design isupdated after every specific number of optimization iteration so that the volume fraction of theintermediate design becomes equal to the prescribed volume V ∗ i at the end of the optimizationprocess when the volume constraint is active (Wang et al., 2011). The formulation furnishesthree material distributions and 3 × N BT load cases but only one design variable field for adesign, one can select as per the manufacturing and material limits. However, the intermediatedesign is for the intended blueprint realization. For the sake of completeness, we present the used nonlinear finite element approach (Bathe,2006; Wriggers, 2008; Zienkiewicz and Taylor, 2005) in brief here. The total Lagrangian finite5lement formulation is considered in this approach, and all desired strains refer to the originalcoordinate system. The deformation gradient F is defined as F = I + ∇ u , (6)where I is the identity tensor and ∇ u indicates gradient of the displacement field with respectto reference configuration X ∈ Ω . In terms of F , the right Cauchy-Green deformation tensor C equals to F (cid:62) F , which is used further to find the Green-Lagrangian strain tensor E as E = 12 ( C − I ) . (7)In view of the standard FE method, the weak form of a mechanical equilibrium equation provides(Zienkiewicz and Taylor, 2005) R ( u , ρ ) = F int ( u , ρ ) − F ext = , (8)where R ( u , ρ ) is a residual term. The internal force at element level F int e is evaluated as F int e = (cid:90) Ω e B (cid:62) TL ( u ) S e ( u , ρ ) d Ω e , (9)where B TL ( u ) and S e are the total Lagrangian strain-displacement matrix (Bathe, 2006) andthe second Piola-Kirchhoff stress tensor of an FE Ω e , respectively. F ext is assumed to be aconstant force here. Eq. 8 is solved herein using the Newton-Raphson (N-R) iterative solver.The second Piola-Kirchhoff stress tensor is evaluated as S = 2 ∂W∂ C using strain-energy functionW and tensor C . W is taken here as (Zienkiewicz and Taylor, 2005) W = G J C − − ln J ) + κ J − , (10)where G = E ν ) is the shear modulus, however one evaluates κ = E − ν ) and κ = E ν )(1 − ν ) for 2D plane stress and 2D plane strain, respectively. J C is the first principal invariant ofthe right Cauchy-Green tensor C and J = det F . Further, E denotes Young’s modulus and ν indicates the Poisson’s ratio. Using the fundamentals of nonlinear continuum mechanics(Holzapfel, 2001), one finds the second Piola Kirchhoff stress S = G ( I − C − ) + κ ( J − J C − for the given strain-energy function (Eq. 10). To handle numerical instabilities in the TO due to large deformation, the energy interpolationscheme (Wang et al., 2014) which ensures smooth deformation, is employed herein. Mathemat-ically, the interpolated strain-energy ¯ W e for an FE is written as¯ W e ( u e ) = (cid:2) W e ( γ e u e ) − W Le ( γ e u e ) + W Le ( u e ) (cid:3) E e ( ρ e ) , (11)where W e ( . ) indicates the strain-energy density function of the actual material at unit Young’smodulus, W Le ( . ) is the strain-energy for small deformation at unit Young’s modulus, and E e is the Young’s modulus of an element (Eq. 4). Further, γ e , defined using a smooth Heavisideprojection filter, is used to determine behavior of a FE as γ e = tanh( β η ) + tanh( β ( ¯ ρ p e − η ))tanh( β η ) + tanh( β (1 − η )) , (12)where η is the threshold value. We use β = 500 and η = 0 .
01 for all examples as suggestedby (Wang et al., 2014). Italic font is used to write the field quantities, whereas the discrete quantities are written using normal font. Objective Formulation and Sensitivity Analysis
This section describes formulation of the objective and its sensitivity analysis.
Let (cid:15) ∗ xx , (cid:15) ∗ yy and (cid:15) ∗ xy indicate the target strains in x − , y − and xy (shear) − directions respectively.Then, the error objective f k can be formulated in a continuum setting as f k = 1 A (cid:90) Ω BT (cid:32) w ( k (cid:15) e xx − (cid:15) ∗ xx ) + w ( k (cid:15) e yy − (cid:15) ∗ yy ) + w ( k (cid:15) e xy − (cid:15) ∗ xy ) w ( (cid:15) ∗ xx ) + w ( (cid:15) ∗ yy ) + w ( (cid:15) ∗ xy ) (cid:33) dA, (13)and in its corresponding FE setting evaluated in element centroids f k = 1 N be N be (cid:88) i =1 (cid:32) w ( k (cid:15) i xx − (cid:15) ∗ xx ) + w ( k (cid:15) i yy − (cid:15) ∗ yy ) + w ( k (cid:15) i xy − (cid:15) ∗ xy ) w ( (cid:15) ∗ xx ) + w ( (cid:15) ∗ yy ) + w ( (cid:15) ∗ xy ) (cid:33) , (14)where w , w , and w are user defined weighing factors depending upon the desired axis ofstraining. Ω BT is the design region for the biological tissue and A is the associated area.Further, N be is the total number of FEs used to represent Ω BT . In the Voigt notation for a 2Dcase, k (cid:15) i xx , k (cid:15) i yy and k (cid:15) i xy are the first, second and third entries of the Green-Lagrange strain k E e which is evaluated at the center, of the i th FE associated to the k th biological tissue. Note thatstrains are not very accurately modelled at the center points when using the standard finiteelement approaches. However, we are minimizing a function which is an integral over a largearea and hence, errors are expected to be very small. A gradient-based approach is employed to solve the optimization problem (Eq. 5). Sensitivitiesof the objective and the constraints with respect to the design vector ρ are evaluated using theadjoint-variable method. The augmented performance function L , defined using the objectiveand the equilibrium equation (Eq. 8), can be written as L = f k ( u k , ¯ ρ ) + λ (cid:62) R k ( u k , ¯ ρ ) , (15)where λ is the Lagrange multiplier vector. Differentiating Eq. 15 with respect to ¯ ρ yields d L d ¯ ρ = ∂f k ∂ ¯ ρ + (cid:20) ∂f k ∂ u k + λ (cid:62) ∂ R k ∂ u k (cid:21)(cid:124) (cid:123)(cid:122) (cid:125) Term 1 ∂ u k ∂ ¯ ρ + λ (cid:62) ∂ R k ∂ ¯ ρ . (16)One chooses λ such that Term 1 vanishes (Bendsøe and Sigmund, 2003), i.e., λ (cid:62) K T = − ∂f k ∂ u k , (17)where K T = ∂ R k ∂ u k (Zienkiewicz and Taylor, 2005) is the tangent stiffness matrix at the equilib-rium state (Buhl et al., 2000). A procedure to evaluate ∂f k ∂ u k is mentioned in Appendix A. ∂ R k ∂ ¯ ρ can be evaluated in view of Eq. 8 as ∂ R k ∂ ¯ ρ = (cid:18) ∂ F int ∂ ¯ ρ + ∂ F int ∂ γ ∂ γ ∂ ¯ ρ (cid:19) . (18) Employed here to represent the stress, strain and material tangent tensors for the FE analysis. ∂ L ∂ ρ as (Wang et al., 2014) ∂ L ∂ ρ = ∂ L ∂ ¯ ρ ∂ ¯ ρ ∂ ˜ ρ ∂ ˜ ρ ∂ ρ . (19)One finds ∂ ¯ ρ ∂ ˜ ρ and ∂ ˜ ρ ∂ ρ using Eq. 1 and Eq. 2, respectively. The method of moving asymptotes(MMA) is used (Svanberg, 1987) to update the design vector. This section presents two sets of numerical examples to demonstrate efficacy and robustness ofthe compliant micro-actuator design optimization approach. In the first set (SET-1), concep-tualized academic examples with N BT = 1 (Eq. 5) are solved to show controllability of strainsin different directions, whereas the second set (SET-2) with N BT = 2 pertains to a practicalapplication, i.e., in conjunction with the flexible poles approach (Vandenburgh et al., 2008). Inaddition, the numerical results are appraised with discussions, 3D-printed prototypes, an ex-perimental setup and results, as well as ABAQUS analyses results. For the examples presented FixedA BCD L x L y a b F in k s dx dyFixedFixed Figure 2: Design domain with a central biological tissue (blue in color). L x × L y = 10 mm ×
10 mm and a = 2 mmand b = 2 mm. k s is the stiffness of the input spring and F in is the actuating force. The bottom, right and topsides of the domain are fixed. in both sets, some parameters are set common as follows. The dilated, intermediate and erodeddesigns are evaluated using ∆ η = 0 .
05. External move limit, i.e., change in design variablesper MMA iteration, is set to 0.1. The material definition given in Eq. 10 is used with Poisson’sratio ν = 0 .
45 and plane strain conditions. The volume fraction for the intermediate design isset to V ∗ i = 0 .
25 and the volume of the dilated design is updated every 25 th MMA iteration.
In this section, the mechanisms are designed in a general setting wherein a biological tissue(blue in color) is assumed to be placed in the middle of the design domain (Fig. 2).8igure 2 indicates the design domain specifications. Length and width of the design domainABCD are L x = 10 mm and L y = 10 mm, respectively. The central blue domain, length a = 2 mm and width b = 2 mm, denotes a biological tissue. The symmetric half of the designdomain is parameterized using N ex × N ey = 200 ×
100 quad-FEs, where N ex and N ey denotenumber of the FEs in x − and y − directions respectively. The bottom, right and top sides of thedomain are fixed, whereas the center of the left side of the domain is used to apply a strain-based actuator controlled by input spring with stiffness k s = 10 000 N / m and blocking force F in , as depicted in Fig. 2. The out-of-plane thickness for the mechanism domain and biologicaltissue is set to 2 mm. Young’s moduli for the biological tissue and remaining region (mechanismdomain) are taken as 0 . β is changed from 1 to β max = 32 using a continuation scheme whereinit is doubled each 60 th MMA iteration and once it reaches maximum value β max , it remainsso for the remaining optimization iterations. The filter radius is set to 5 . × max( L x N ex , L y N ey ).Maximum number of optimization iteration is set to 400. The objective is evaluated within theblack edged rectangle (Fig. 2), where dx= 3 L x N ex and dy = 3 L y N ey are taken. This is done toavoid inclusion of high localized strains, appearing on the corners of the biological tissue, in theobjectives evaluation.The binary nature of the optimized mechanisms is measured by a gray scale indicator M nd which is defined as (Sigmund, 2007) M nd = (cid:80) N e i =1 ρ i (1 − ρ i ) N e × , (20)where N e is the total number of FEs used to describe the design domain. The root mean square(RMS) errors in actual strain can be evaluated asErr x = (cid:113) f k | w =1 , w = w =0 × , Err y = (cid:113) f k | w =1 , w = w =0 × , Err xy = (cid:113) f k | w =1 , w = w =0 × , (21)where Err x , Err y and Err xy are the RMS errors in the x − , y − and shear directions, respectively. In this example, we seek a mechanism which can induce 20% axial (in the x − direction) strainin the biological tissue (Fig. 2) when it is actuated by a force F in = − . x − direction . (cid:15) ∗ xx = 0 . , w = 1 and w = w = 0 are used to evaluate the objectives.Figure 3a indicates the full optimized mechanism with the central biological tissue for theintermediate design. The gray scale indicator M nd for the optimized dilated, intermediate anderoded mechanisms are evaluated to be 0.56%, 0.47% and 0.46%, respectively.The strain distribution for x − direction is depicted in Fig. 3b. One notices, the strain distributionin the tissue (within the black edged rectangle) is close to uniform with value approximatelyequal to 0 . x error in the strain is calculated using Symmetric about a horizontal line Corresponding to an unloaded actuator with displacement 0 . Suitably transferred from the symmetric half results (a) Intermediate design (b) Actual strain distribution(c) Objectives convergence plot (d) Volume fraction convergence plotFigure 3: Solutions for inducing 20% strain in x − direction in the central tissue. (a) Optimized intermediatedesign, M nd = 0 . In this example, a mechanism which can provide a desired bi-axial straining in the biologicaltissue, is designed. The design domain and optimization specifications for this example aresame as those used for
Example 1 . The target strains (cid:15) ∗ xx = 0 .
125 and (cid:15) ∗ yy = 0 .
075 are set. Theactuating force for this case is F in = − . (cid:15) ∗ xx = 0 . (cid:15) ∗ yy = 0 . (cid:15) ∗ xy = 0 . w = w = w = 1. Figure 4 indicates the optimized intermediatedesign of the mechanism which can provide bi-axial strains in the central tissue. The actualstrain distribution for the x − and y − axes are reported in Fig. 4b and Fig. 4c respectively. Therecorded errors in strain distribution for x − , y − and xy − (shear) directions are Err x = 11 . y = 11 .
18% and Err xy = 10 . p − norm, if these errorsare critical. One notices that the strains at the edges of the biological tissue are higher thanthe desired ones (Fig. 4b and Fig. 4c). However, within the black-edged rectangle, the strain-distributions (Fig. 4b, Fig. 4c and Fig. 4d), by and large, are very close to their respective targetstrains. M nd = 2 .
91% is obtained for the intermediate design (Fig. 4a).10 a) Intermediate design (b) Strain in x − direction (c) Strain in y − direction (d) Shear strainFigure 4: Solutions for inducing 12.5%, 7.5% and 0.0% strains in x − , y − and shear directions in the biologicaltissue. (a) Optimized intermediate design, M nd = 2 .
91% (b) Strain distribution in x − direction, (c) Straindistribution in y − direction. The maximum and minimum observed strains in y − direction are 0.11 and 0.035,respectively. (d) Shear strain distribution. This example considers inducing strains in axial, transverse and shear directions in the centralbiological tissue (Fig. 2). The design domain specifications and optimization parameters are thesame as those used for
Example 1 . The full design domain is considered for the optimizationinstead of the symmetric half in order to accommodate the prescribed non-zero shear strain. Thetarget strains are set as (cid:15) ∗ xx = 0 . (cid:15) ∗ yy = − .
15, and (cid:15) ∗ xy = 0 .
20 and the objective is evaluatedwith w = w = w = 1 . The actuating force in the x − direction is set to F in = − . ×
200 bi-linear FEs to describe the design domain. (a) Intermediate design (b) Strain in x − direction (c) Strain in y − direction (d) Shear strainFigure 5: Solutions for inducing 30.0%, -15.0% and 20.0% strains in x − , y − and shear directions in the biologicaltissue. (a) Optimized intermediate design, M nd = 1 .
13% (b) Strain distribution in x − direction, (c) Straindistribution in y − direction. (d) Shear strain distribution. The optimized intermediate design of the mechanism is displayed in Fig. 5a, which inducesthe desired strains in axial (Fig. 5b), transverse (Fig. 5c) and shear (Fig. 5d) directions in thetissue. The gray scale indicator M nd is 1 . x = 19 . y = 13 .
30% and Err xy = 17 . The robustness of the synthesis approach using the numerical examples is illustrated in Sec. 4.1.This section presents a set of compliant actuators designed in accordance with the flexible poles method (Vandenburgh et al., 2008). Herein, it is envisioned that the presented approach shallnot only support one pair of flexible poles (one bioreactor) but can also provide support tomany other similar bioreactors requiring different induced strains in biological tissues having11 ǫ ∗ ǫ ∗ ǫ ∗ ǫ ∗ Base plateMechanism domainsFlexible polesBiological tissues n ǫ ∗ : Target strain in the n th tissue ∆ : Amount of actuationFixed part of the base plate Figure 6: A schematic diagram for the design problem in a region
A BD CF EM M M M P P P P P P P P B B ∆ Symmetry lineS S L x L y Figure 7: Design domain for designing a compliant mechanisms in conjunction with the flexible poles approach(Vandenburgh et al., 2008). The left and bottom sides of the base plate are with roller supports. The right sideof the base plate is used for applying actuation of amount ∆ = 0 . L x . Figure has same color scheme as Fig. 6. either same or different geometries (see Fig. 6). In Fig. 6, different regions are indicated. Foreach pair of flexible poles, we seek the optimized compliant mechanisms which can provide aspecific desired strain n (cid:15) ∗ xx in the n th biological tissue when the base plate is strained uniformlyby ∆.Figure 7 indicates the design domain specifications for designing compliant micro-actuator mech-anisms with two flexible poles, a sample of biological tissue construct and a base plate. Table 1depicts the dimensions, material parameters (Young’s moduli) and thicknesses for these domains(Pless, 2019). The mechanisms are designed for achieving 5%, 10%, 15% and 20% straining intheir respective biological tissue and named herein as compliant bio-mechanism CBM I, CBMII, CBM III and CBM IV, respectively. The color scheme of Fig. 7 is used further to show theresults wherein black color is used for the optimized mechanisms.Designing these robust compliant mechanisms poses unique challenges because: (i) the biological12 egion Name Young’s Modulus Length Height Thickness
ABCD Base Plate 18 .
68 MPa 3 mm 1 mm 2 mmM M M M Mechanism Domain 18 .
68 MPa 2 mm 1 mm 2 mmP P P P Flexible Pillar 18 .
68 MPa 0 . . . P P P Flexible Pillar 18 .
68 MPa 0 . . . B P P Biological Tissue Construct 0 .
20 MPa 1 mm 0 . . Table 1: Nomenclature, Material and dimensional specifications for the flexible poles design domain (Fig. 7)(a) Intermediate design (b) Strain distribution (c) ABAQUS strain distribution(d) Objectives convergence plot (e) Volume constraints convergence plotFigure 8: Full solutions to CBM I (5% desired straining). (a) Optimized intermediate design, M nd = 1 . tissue is comparatively very soft, (ii) the tissue construct is connected to the micro-mechanismvia flexible poles and (iii) it is essential to consider geometric nonlinearity as already mentionedbefore. To account for the challenges and also to permit different tissue construct geometries,we use the extended robust formulation with N BT = 2 (see Sec. 2.1). The first tissue constructcorresponds to the actual (minimum) given size, whereas the width of the second tissue constructBT2 is assumed to be two times that of the first one ensuring tissues with different (higher)stiffness.The symmetric half of the design domain ABEF (Fig. 7) is discretized using N ex × N ey = 200 ×
280 quad-FEs, and thereafter, FEs representing different regions are detected. FEs describingthe base plate, poles and biological tissue are attributed by ρ = 1 throughout the optimizationprocess. Likewise, ρ = 0 is designated to all FEs associated to the void regions. Roller boundaryconditions permitting movement of the bottom and left sides of the plate along x − and y − axes, Not directly actuated by the mechanism a) Dilated design (b) Intermediate design (c) Eroded design(d) Strain distribution (e) ABAQUS strain distributionFigure 9: Full solutions to CBM II (10% desired straining). (a) Optimized dilated design, M nd = 1 .
62% (b)Optimized intermediate design, M nd = 1 . M nd = 1 . respectively, are considered, whereas the right side of the plate is actuated by an amount∆ = 0 . L x , i.e., 10% of uniform straining in the base plate (Fig. 7). Filter radius is set to10 × max( L x N ex , LyN ey ). The Heaviside projection filter parameter β is altered from 1 to 128. Itis double at each 25 th MMA iteration till it reaches to 128 and thereafter, it remains 128 forthe remaining optimization iterations. The maximum number of optimization iterations is setto 250. The optimized results are shown with the actual tissue construct.Figure 8a, Fig. 9b, Fig. 10b, and Fig. 11a show the optimized intermediate designs for CBM I,CBM II, CBM III and CBM IV respectively. Figure 9a and Fig. 10a, and Fig. 9c and Fig. 10cdepict optimized dilated and eroded designs for the CBM II and CBM III. One notices theoptimized topologies CBM III and CBM IV are identical, however the latter one has compara-tively some thin slender sections facilitating more deformation and thus, help providing close to20% strain in the tissue construct (Fig. 11b and Fig. 11c). The objectives convergence plots forCBM I and CBM IV are indicated by Fig. 8d and Fig. 11d respectively and their correspondingvolume fraction convergence plots are depicted via Fig. 8e and Fig. 11e.As β changes, the approximating function in Eq. 1 alters and thus, as per β continuation, jumpsin the convergence curves of the objectives and volume constraints can be noticed. In addition,with updates in the volume of the dilated design, the convergence curves may also get altered.Convergence plots are smooth after 200 MMA optimization iterations and volume constraint issatisfied and remains active in the end of the optimization (Figs. 8d, 8e, 11d, and 11d).One can notice uniform distributions of actual strains, close to their desired ones, within therespective tissues of CBM I, CBM II, CBM III and CBM IV (Fig. 8b, Fig. 9d, Fig. 10d andFig. 11b). The actual strain distribution is demonstrated with respect to the optimized inter-14 a) Dilated design (b) Intermediate design (c) Eroded design(d) Strain distribution (e) ABAQUS strain distributionFigure 10: Full solutions to CBM III (15% desired straining). (a) Optimized dilated design, M nd = 1 . M nd = 1 . M nd = 1 . mediate designs. The strain errors Err x (Eq. 21) 12 . . .
4% and 14 .
8% are noticedfor CBM I, CBM II, CBM III and CBM IV, respectively.To demonstrate the accuracy of the obtained results, extracted and smoothed intermediatedesigns are also analyzed in ABAQUS with the same boundary conditions, actuating forces andmaterial properties as those used by the formulation. Figure 8b, Fig. 9d, Fig. 10d, and Fig. 11billustrate the actual strain distributions and their respective results obtained by ABAQUSanalyses are depicted by Fig. 8c, Fig. 9e, Fig. 10e, and Fig. 11c respectively. One can see thatthe strain distributions obtained by the presented approach and ABAQUS analyses are in closeagreement with each other.
A stereolithography-based 3D printing process (Zhang and Larsen, 2017) is employed to fabri-cate the optimized mechanisms CBM I, CBM II, CBM III and CBM IV with their flexible polesand a base plate wherein a hydrogel material based on poly (ethylene glycol) diacrylate is usedfor printing of the final prototypes. The adopted printing steps are as follows:1. The optimized mechanisms are converted into 3D CAD (Computer-aided design) mod-els with their flexible poles. These models are assembled on a base plate. A circularhole with diameter 500 µ m is extruded from the one end of the base-plate to facilitateactuation/stretching (Fig. 12a). 15 a) Intermediate design (b) Strain distribution (c) ABAQUS strain distribution(d) (e)Figure 11: Full solutions to CBM IV (20% desired straining). (a) Optimized intermediate design, M nd =1 . Base plateCBM I CBM II CBM III CBM IVSupporting structures A circular hole (a) CAD model
Image acquisitionActuation setupActuating direction (b) Experimental setup
Printed structureBeakerNeedleActuating direction (c) Sample setupFigure 12: (a) 3D CAD model, (b) Overview of the experimental setup and camera position and (c) Close-up ofsample fixation
2. Cylindrical supporting structures, radii 150 µ m and heights 750 µ m, are placed on a squaregrid of center-to-center spacing 500 µ m, i.e., with a 200 µ m gap to ease the printing process(Fig. 12a).3. The assembled 3D CAD design is sliced into 20 µ m layers using the open-source softwareSlic3r for printing process.4. A photo-curable resin consisting of 50 %v/v poly (ethylene glycol) diacrylate, M n ase plateActuating side Fixed sideCBM IV CBM III CBM II CBM I L Figure 13: Photo of 3D printed mechanisms, immersed in water, with a base plate of length L = 10 .
6. After printing, structures are washed in water and swollen to equilibrium. Support struc-tures are detached from the print with the help of a needle.In the experimental setup (Fig. 12), the cover-glass with attached printed structures is fixedon the bottom of a rectangular transparent beaker filled with water (Fig 12c). The prototypeis aligned by 3-axes stages (RB13M, Thorlabs, Inc.). A 500 µ m outer-diameter needle is fixedto the actuation stage and its tip is bent by 90 o and subsequently inserted into the hole onthe actuating side. Actuation is carried out by moving the needle in steps of 100 µ m in theactuation direction until failure (∆). An image for every position is acquired by a camera(UI-3880LE-C-HQ, IDS Imaging) placed orthogonally to the actuated platforms (Fig 12b).At the preliminary stage (presented herein), the experiment is performed without biologicaltissue constructs (Fig. 12). Figure 13 indicates prototypes of the mechanisms with their flexiblepoles and a base plate having length L . The base plate is stretched up to 8% of its lengthi.e. ∆ = 0 . L (Fig. 14a and Fig. 14b). The deformed profiles of the mechanisms are imagedand displayed in Fig. 14a. An animated image sequence of the gradual profile deformation asa function of stretching is included in the supplementary material (Movie 1). A correspondingABAQUS model is developed and analyzed. Figure 14b illustrates the deformed profiles of themechanisms with respective flexible poles obtained via the ABAQUS analysis. The printedmechanisms perform as they are expected, i.e, they could move apart their respective flexiblepoles and thus, can induce strains in the respective tissue constructs when base plate is actu-ated. The tiny dimensions and compliant materials used do not allow for functional modelingusing externally generated phantoms, such as suspended rubber bands. Functional testing withmicrotissues suspended between sets of poles is outside the scope of the current work. However,we have recently demonstrated such stable and reproducible muscle tissue generation in fully3D printed flexible pole devices of similar dimensions assisted by an integral microreservoir forcell seeding (Christensen et al., 2020), which we will combine with the presented TO approachin future work. The recorded axial stretching in the flexible poles pertaining to CBM I, CBM II,CBM III and CBM IV via the experiment are ∆ E = 0 . l , ∆ E = 0 . l , ∆ E = 0 . l and∆ E = 0 . l , respectively. Corresponding stretching obtained from the ABAQUS analysis are∆ A = 0 . l , ∆ A = 0 . l , ∆ A = 0 . l and ∆ A = 0 . l , where l = 1 mm. It can be notedthat experimental result show: (i) higher deformation for CBM III and CBM IV, (ii) closeagreement for CBM II and (iii) lower deformation for CBM I with respect to its ABAQUS anal-ysis. There could be many reasons for such discrepancies, e.g., geometry variations, descriptionof material properties including the Poisson’s ratio, boundary conditions, out-of-plane bendingof the 3D-printed flexible poles, and they need further and deeper investigations which are outof scope of the current manuscript and left for our future endeavors.17 E1 ∆ E2 ∆ E3 ∆ E4 L + ∆ (a) Deformed profiles: Experiment ∆ A1 ∆ A2 ∆ A3 ∆ A4 L + ∆ (b) Deformed profiles: ABAQUSFigure 14: Deformed profiles of the mechanisms. (a) Experimental deformed profiles with ∆ E = 0 . l , ∆ E =0 . l , ∆ E = 0 . l , and ∆ E = 0 . l b Deformed profiles obtained by ABAQUS analysis wherein ∆ A =0 . l , ∆ A = 0 . l , ∆ A = 0 . l , and ∆ A = 0 . l . Here, l = 1 mm i.e. the initial gap between each pair offlexible poles and ∆ = 0 . L . This paper presents a method using topology optimization to design large deformation compliantmechanisms for inducing the desired strains in biological tissues. An objective based on leastsquare error is formulated using the given target strains and minimized. To cater for largedeformation of the tissues, geometric nonlinearity with a suitable neo-Hookean material model isconsidered. The versatility of the presented approach is demonstrated by designing mechanismswhich can induce strains in the biological tissues in both axial and bi-axial directions.The mechanism design problem is conceptualized in the flexible poles environment, and variouscompliant mechanisms are successfully designed for inducing different strain levels in theirrespective biological tissues. The robust formulation is extended to accommodate differentgeometries of the tissue constructs. A base plate is used to actuate the mechanisms whichrender specific movements in their associated flexible poles and thus, help inducing the targetstrains in the tissues. Actual strain distributions in the tissues by the optimized mechanismsusing the approach closely resemble those determined using their respective ABAQUS analyses.The optimized mechanisms with their flexible poles and a base plate are 3D-printed usingpoly(ethylene glycol) diacrylate material and a simplified experiment is performed. With re-spect to its corresponding ABAQUS analysis, we observe good qualitative agreement with somediscrepancies in the stretches developed by their flexible poles. These discrepancies could haveresulted from geometry variations, human errors, boundary conditions, material properties, and,a subject for our near future study. In addition, extension to a 3D setting with flexible polesenvironment is one of the prime directions for future work.
Acknowledgment
All authors acknowledge support from Independent Research Fund Denmark, grant 7017-00366B. O. Sigmund acknowledges the support from the Villum Investigator project InnoTopprovided by the Villum Foundation. The authors also acknowledge Prof. Krister Svanberg forproviding MATLAB codes of the MMA optimizer.18 ppendix A Evaluating the derivative ∂f k ∂ u In view of Eq. 14, one finds the derivative ∂f k ∂ u e as ∂f∂ u e = 2 N be N be (cid:88) e =1 w ( (cid:15) e xx − (cid:15) ∗ xx ) ∂(cid:15) e xx ∂ u e + w ( (cid:15) e yy − (cid:15) ∗ yy ) ∂(cid:15) e yy ∂ u e w ( (cid:15) ∗ xx ) + w ( (cid:15) ∗ yy ) + w ( (cid:15) ∗ xy ) (A.1)Therefore, one needs ∂(cid:15) e xx ∂ u e and ∂(cid:15) e yy ∂ u e and they can be extracted from the derivative ∂ E e ∂ u e . Now,using Eq. 7 and Eq. 6, we have E = 12 (cid:16) ∇ u + ( ∇ u ) (cid:62) + ∇ u ( ∇ u ) (cid:62) (cid:17) , (A.2)In view of FE setting, the displacement vector u of an element in terms of its nodal displacements u AI and bi-linear shape functions N A can be written as : u = (cid:88) A N A ( ζ ) u AI = N A ( ζ ) u AI . (A.3)Now, Eq. (A.2) yields using Eq. (A.3) as E IJ = 12 (cid:18) ∂N A ∂X J u AI + ∂N A ∂X I u AJ + ∂N A ∂X I ∂N B ∂X J u AK u BK (cid:19) . (A.4)One finds derivative of E IJ with respect to u AI as ∂E IJ ∂u AK = 12 ∂N A ∂X J (cid:18) δ IK + 2 ∂N B ∂X I u BK (cid:19) , (A.5)and hence, ∂(cid:15) e xx ∂ u e and ∂(cid:15) e yy ∂ u e . References