Combined Level-Set-XFEM-Density Topology Optimization of 4D Printed Structures undergoing Large Deformation

Advancement of additive manufacturing is driving a need for design tools that exploit the increasing fabrication freedom. Multi-material, 3D printing allows for the fabrication of components from multiple materials with different thermal, mechanical and “active” behavior that can be spatially arranged in 3D with a resolution on the order of tens of microns. This can be exploited to incorporate shape changing features into additively manufactured structures. 3D printing with a downstream shape change in response to an external stimulus such as temperature, humidity or light is referred to as 4D printing. In this paper, a design methodology to determine the material layout of 4D printed materials with internal, programmable strains is introduced to create active structures that undergo large deformation and assume a desired target displacement upon heat activation. A level set approach together with the extended ﬁnite element method (XFEM) is combined with density-based topology optimization to describe the evolving multi-material design problem in the optimization process. A ﬁnite


Introduction
Advanced additive manufacturing (AM) technologies like three-dimensional (3D) printing have improved vastly in recent years in terms of accuracy, material variety and reliability.Recently, the concept of 4D printing was introduced by [1], where external stimuli are used to trigger shape-changes after a structure is 3D printed.The shape-change over time is therefore seen as the fourth dimension.In [2] this concept was used to design transforming shape memory polymer structures utilizing the thermomechanical response of a glassy polymer within an elastomer matrix.A similar material response was studied in [3] to create active origami structures.Exploiting the shape memory behavior of a 3D printed glassy polymer, [4] used level set (LS) topology optimization to create printed active composites which assume a target shape after a thermomechanical training and activation cycle.As stated by [5], the drawback of these traditional approaches of 4D printing is the fact that a complex thermomechanical training and activation cycle is required through which a shape-change can be triggered.To alleviate this issue, [5] introduced a novel approach called "direct 4D printing" where no external training cycle is required before the shape change occurs upon printing.The so called "printing strain" used in this method for creating a shape change, is an inelastic eigenstrain that is programmed into the 3D printed structure during the printing process.This strain is compressive, and its magnitude can be designed by controlling the time and intensity of the UV-curing of photo polymers in the polyjet process.In the current work, two polymers with different magnitudes of eigenstrain and different glass transition temperatures are used to achieve a desired shape change upon release of the built-in eigenstrain.A component fabricated in this manner therefore consists of multiple well-bonded polymers with a high-fidelity geometry, and with each polymer containing a spatially-variable eigenstrain.Upon printing, the component is heated to release the built-in eigenstrain of the rubbery polymer (active material) by lowering the stiffness of the glassy polymer (passive material).This is achieved through heating of the printed structure beyond the glass transition temperature of the passive material causing a permanent shape change.Alternatively, a frontal polymerization process has been proposed by [6] and [7].This process uses the change in intensity of photo polymerization to create spatially varying material parameters leading to 4D printed self-folding origami structures.Shape changes are triggered by a solvent entering a loose polymer matrix which then causes isotropic swelling.In this work, the direct 4D printing method is employed as a convenient way to incorporate spatially varying eigenstrains during fabrication of self-actuating plate-like structures.For a more detailed study of the process of direct 4D printing, the interested reader is referred to [5].
In order to exploit spatially varying inelastic strains to obtain complex target displacements upon activation, an inverse design problem is formulated and solved.Given a desired target deformation, the optimal geometry and material layout of active 4D printed components are determined in a systematic manner.In previous works [8] and [9], this has only been done for 1D rod structures or for simple target shapes [5] where it was possible to pre-determine the optimal material layout intuitively.A similar approach based on classical origami designs is demonstrated in [10] for self-folding composite structures.In this paper, a new multi-material topology optimization formulation is used to determine the material arrangement of shape-changing 4D printed structures undergoing large deformation.The proposed computational design framework accurately captures the physical response of the structure by employing a hyperelastic thermomechanical model combined with a higher-order extended finite element method (XFEM) formulation.A hybrid approach combining LS-XFEM and densitybased topology optimization is introduced to solve the multi-material design problem.These topology optimization schemes have each demonstrated individually a wide range of applicability.An overview of LS and density topology optimization methods is given in [11], [12], and [13] respectively.Some applications of topology optimization can also be found in the field of design optimization of self-folding structures, for example: [14] used a density-based method to optimize the layout of monolithic liquid crystal elastomers in order to create folding liquid crystal elastomer actuators.However, a simplified linear elastic model operating in the small strain regime greatly limits the accuracy of this approach.A different simplified approach was taken by [15] where shape optimization is employed to determine the optimal layout of cuts to design active 3D origami and kirigami structures using 4D printing.And just recently, [16] applied the concept of moving morphable components (MMC) [17] in combination with a genetic algorithm (GA) for topology optimization of post-buckled 3D kirigami structures.Even though each of these approaches greatly reduces the number of design variables and hence the complexity of the optimization problem, only a sub-space of all possible designs is explored.To take advantage of the entire design space while considering a fully nonlinear structural response, a combined LS-XFEM and density topology optimization framework for designing 4D printed active structures is proposed in this paper.
The design, fabrication and activation process of 4D printed structures as proposed in this work is conceptually shown in Fig. 1, where target displacement matching is the objective.This can be achieved by either specifying the desired target shape, i.e. position of points of the structure in the deformed (activated) configuration, or by specifying the required displacement of certain points of the structure to achieve the desired deformation.In this work, the target displacement of target points on the domain boundary is specified to achieve target deformations, see Fig. 1 (a).Starting from an initially flat plate (b), a level set function (LSF) is used to define an initial arrangement of solid and void domains.Within the solid domain, a fictitious density field is used to interpolate the different properties of an active and a passive material.Both the LSF and the density distributions are discretized on the XFEM background mesh.The resulting parameter optimization problem is then solved by a nonlinear programming method.After a final design (i.e.geometry and material distribution) has been found (c), sub-domains of active and passive material are extracted in a post-processing step (d).Once the flat structure has been printed (e), the shape-change is activated through relaxation of the compressive printing strain in the active material in a 65.0 • C hot water bath.This is achieved by heating up the structure above the glass transition temperature of the passive, stiff material at which its stiffness is significantly reduced.At this state, the eigenstrains of the active material can relax leading to the desired change in shape.After a subsequent cool-down, the passive material stiffens and a permanently deformed structure is obtained (f).
The theoretical background of the physical model, the spatial discretization of state and optimization variables, and the optimization problem along with examples are presented in the remainder of this paper, which is organized as follows: Sec. 2 introduces the continuum model of the nonlinear thermomechanical phenomena.Sec. 3 provides details on the optimization approach and the XFEM model.Sec. 4 describes the validation of the thermomechanical model, and design optimization examples are presented in Sec. 5 along with experimental results.A summary of the work is presented in Sec. 6.

Hyperelastic Thermomechanical Model
The total design domain Ω 0 D is comprised of solid and void sub-domains (i.e.Ω 0 S and ), balance of linear momentum in the undeformed configuration is described by the static equilibrium.The weak form of the governing equation using a total Lagrangian formulation is stated as: where F = ∂x/∂X is the total deformation gradient tensor with x = u + X defining the relation between spatial coordinates in the undeformed (X) and deformed (x) configurations Ω 0 S and Ω S respectively.The displacements are denoted by u and δu are the admissible test functions.The first Piola-Kirchhoff stress tensor is denoted by P and B is the prescribed body force vector in the undeformed configuration.The density is denoted by ρ 0 and ū is the prescribed displacement field on Γ 0 ū, see Fig. 2 (a).The prescribed traction vector on Γ 0 T is T and the total design domain boundary Γ 0 is comprised of A hyperelastic Saint Venant-Kirchhoff constitutive model for isotropic compressible solids is used in the current work.The hyperelastic Saint Venant-Kirchhoff material model is suitable for large deformations at small strains, which is the case for the slender structures targeted in this work.The proposed design optimization framework can easily be extended to operate on nonlinear constitutive models for both compressible and incompressible materials.This is especially relevant for simulating rubbery, nearly-incompressible polymers and therefore should be treated in future work.The Saint Venant-Kirchhoff material model is defined as: where S is the second Piola-Kirchhoff stress tensor in the intermediate configuration (see Fig. 2 (b)) and E M is the Green-Lagrange strain tensor defined as the deformation of the intermediate configuration to the thermomechanical deformed configuration.The second order identity tensor is denoted by I.The Lamé constants are denoted by µ and λ which are related to the Young's modulus E and Poisson's ratio ν as follows: In the context of topology optimization using the density method, the Young's modulus is manipulated as a function of the design variables whereas the Poisson's ratio is kept constant.A more detailed discussion of the optimization problem is presented in Sec.3.1.
To model the inelastic printing strain built into the elastomer, a residual thermal strain model is used as suggested by [5].In a finite deformation case, special consideration with respect to the decomposition of the total deformation gradient is required.As indicated in Fig. 2 (b), a multiplicative decomposition of the total deformation gradient F into an inelastic thermal deformation gradient F T and an elastic mechanical deformation gradient F M is used.The multiplicative decomposition of the total deformation gradient can be stated as: For the thermal deformation gradient tensor, a linear thermal expansion model is assumed: where α k is the linear coefficient of thermal expansion (CTE) of either solid material (k = active or passive) representing the amount of inelastic expansion.The externally applied temperature is T which is a linear function of pseudo time of the deformation process and T 0 is the reference temperature.Using the previously defined relationships, the Green-Lagrange strain tensor is defined as: where the right Cauchy-Green deformation tensor is For integration of the governing equation stated in Eqn.(1) in the undeformed configuration, the second Piola-Kirchhoff stress tensor S is pulled-back from the intermediate configuration ΩS to the undeformed configuration Ω 0 S using the Piolatransform.This transformation is defined as: where J T is the Jacobian of the thermal deformation gradient tensor J T = det(F T ).The first Piola-Kirchhoff stress tensor in the undeformed configuration is finally computed as:

Follower Pressure Load
In order to enforce an end-stiffness constraint on the final, deformed structure, a follower pressure load is introduced and added to the weak form of the governing equation (1).To this end, a deformation dependent (i.e.follower) Neumann boundary condition is employed, which is formulated in the deformed configuration and mapped back to the undeformed configuration.Using Nanson's formula [18], the relationship between surface areas in the undeformed and deformed configuration is defined.The non-conservative follower pressure can simply be evaluated in the undeformed configuration as follows: where dA is an infinitesimal surface area in the undeformed configuration, J = det(F) is the determinant of the total deformation gradient and N is the surface normal in the undeformed configuration.The surface pressure scalar reformulated in the undeformed configuration is denoted by P. For simplicity it is assumed that P is not state dependent.A more detailed discussion of the end-stiffness constraint is presented in Sec.3.4.

Adaptive Load-Stepping Approach
To efficiently predict the nonlinear deformation of 4D printed self-deforming structures, an adaptive load-stepping scheme is used for the XFEM analysis.Schematically, the evolution of the temperature load T (t), applied follower load P(t) and displacement response u(t) as a function of pseudo time t is depicted in Fig. 3.An adaptive time-stepping scheme is used such that the initially set time step size ∆t 0 is reduced by a factor f in case a Newton-Raphson solve does not converge to an equilibrium solution for a set maximum number of nonlinear iterations.This is, for example, encountered when domains of intermediate (weak) material dominate the structural response during the topology optimization process.In such an event, the time step is adaptively reduced until a specified number of (converged) nonlinear solutions are obtained at the reduced time step.After that, the original time step size is gradually restored.The total number of time steps is adjusted accordingly to reach the final loading time t Load .The temperature load T is increased linearly between t 0 and t Load to achieve a maximum value T while the follower surface pressure load P is zero during this time.Conceptually, this is shown in Fig. 3 (a) and (b), respectively.The maximum value of T corresponds to the inelastic printing strain.Between t Load and t Pert the external temperature load is kept constant at T while a constant non-zero follower perturbation load P is applied.The magnitude of the non-zero follower perturbation load is chosen sufficiently small in order to achieve fast convergence of the nonlinear solver, ideally in a single iteration.This is valid for enforcing an end-stiffness constraint measuring the change in total strain energy based on a linear concept.MD-18-1343, Geiss, 5

Multi-Material Topology Optimization
A multi-material topology optimization approach is adopted to determine the geometry and the spatial material arrangement of 4D printed active structures.This approach builds on a LS-XFEM optimization framework previously used to study problems in structural mechanics [4,19,20] and fluid mechanics [21,22].This section provides an overview of the multi-material topology optimization approach, while more details regarding immersed boundary techniques used for design optimization are provided in [23].

Combined LS-Density Geometry and Material Description
The three-material problem depicted in Fig. 2 (a) is described by a combined LS-density approach.A nodally discretized LSF is used to distinguish between a solid and a void domain where negative LS values (φ i < 0) represent the solid domain Ω 0 S and positive LS values (φ i > 0) represent the void domain Ω 0 V .The zero LS iso-contour (φ i = 0) represents the phase boundary Γ 0 S,V between the solid and the void domain.Nodal LS values φ i are explicitly defined in terms of nodal design variables s φ j using a linear filtering scheme as proposed by [24].This linear filtering technique enhances convergence of the optimization problem by increasing the zone of influence of each design variable and is formulated as: where N n is the number of finite element (FE) nodes within the smoothing radius, r s , and |X i − X j | is the Euclidean distance between node i and j measured in the undeformed configuration (indicated by a superscript 0 ).Index i denotes the current node for which the LS value is computed and index j denotes each node within the smoothing radius contributing to the LS value φ i .
Combining the level set method (LSM) with the XFEM, yields a crisp solid-void interface which is naturally suited for AM reducing the need for post-processing.However, as mentioned in [24] and [12], a strong influence of the final solution on the initial guess is observed in LS-XFEM topology optimization.Moreover, since LS-XFEM topology optimization is solely driven by localized sensitivities along the interface, the appearance of new holes within the solid domain is not possible [4].One way to mitigate this issue and the dependency of the final design on the initial guess is to use topological derivatives, as suggested by [25].In the current work, a sufficiently large number of initial void inclusions is used to mitigate this dependency.
In addition, a density-based topology optimization approach is used within the solid domain Ω 0 S in order to distinguish between an active material (Ω 0 S A ) and a passive material (Ω 0 S P ).This combined approach allows for the description of the multi-material topology optimization problem at hand.Nodally discretized, fictitious density design variables 0 ≤ s ρ i ≤ 1 are used to track the material distribution within the solid phase.A standard solid isotropic material with penalization (SIMP) approach [26] is adopted.Considering an element-wise constant density interpolation, the material property within an element is defined as: where p e represents any elemental material property like the Young's modulus, density or CTE.Physical properties corresponding to the active material and passive material are denoted by p e min and p e max , respectively and β is the so called SIMP exponent used to achieve different interpolation behavior.Different SIMP exponents are used for interpolation of different material properties.Elementally averaged fictitious density values ρe j are obtained from nodal fictitious density design variables s ρ i as: where N e n is the number of nodes per element j.A linear elemental filter similar to the one stated in Eqn. ( 10) is then used to compute a smoothed fictitious elemental density ρe i from the elementally averaged fictitious densities ρe j .In addition, a smoothed Heaviside projection scheme proposed by [27] is used subsequently to mitigate the blurriness introduced by the linear filtering scheme.The projection function used to compute the fictitious projected elemental density ρe i is stated as: MD-18-1343, Geiss, 6 where γ is the projection sharpness parameter and η is the projection threshold parameter.The projection ( 13) is applied through a continuation approach where the projection threshold parameter is gradually increased in the optimization process.More details regarding this approach are found in the discussion of the results in Sec. 5.

XFEM Model
The XFEM is used in this work to discretize the weak form of the governing equation ( 1) on non-conforming subdomains.This immersed boundary method is utilized to distinguish between the solid (Ω 0 S ) and the void (Ω 0 V ) sub-domains along the interface (Γ 0 S,V ) defined by the zero LS iso-contour as depicted in Fig. 2 (a).The advantage of this discretization method is that a high spatial resolution of the interface described by the LSM is retained throughout the optimization process while operating on a fixed background mesh, thus avoiding the need for re-meshing.The XFEM approach, however, requires enriching the classical FE approximation spaces with additional shape functions [28] in order to avoid spurious coupling and load transfer between physically disconnected domains.Note that, in general, enrichment functions are also needed to capture weak and strong discontinuities across material phases.However, this requirement does not apply here, as the XFEM is applied to distinguish only between the solid (Ω 0 S ) and the void (Ω 0 V ) phase.Thus, the displacement field only exists within the solid sub-domain and any element within the void phase can be omitted.A generalized Heaviside enrichment approach [29] is used where the displacements in the solid phase are approximated by standard FE shape functions.The nodal displacements u i (X) of node i within the solid domain Ω 0 S are approximated as [22]: where H is the Heaviside function as a function of the LS value φ(X) defined as: The maximum number of enrichment levels is denoted by M, N k (X) is the elemental shape function and δ k mq is the Kronecker delta which selects the active enrichment level q for node k. δ k mq ensures that displacements at node k are only interpolated by a single set of degrees of freedom (DOFs) defined at node position X such that the partition of unity principle is satisfied [20].The Heaviside function H is used to only activate shape functions within the solid phase as displacement solutions in the void domain would be physically meaningless.For more details about the generalized Heaviside enrichment strategy employed in this work, the interested reader is referred to [23], [30], and [31].
For stabilization of the XFEM discretization, a geometric preconditioning scheme as proposed by [32] is used.This geometric preconditioner ensures that DOFs with vanishing zones of influence, which arise when the LS intersection is too close to a FE node, are re-scaled or eliminated in order to provide numerical stability.
Dirichlet boundary conditions applied to the solid sub-domain Ω 0 S are enforced weakly using Nitsche's method [33].The weak form of the governing equation ( 1) is augmented with the weak boundary condition residual contribution R D Γ stated as: where the jump operator [[•]] is defined as: The penalty factor to enforce the prescribed Dirichlet boundary condition ū is denoted by γ D .The first term in Eqn. ( 16) corresponds to the standard consistency term, the second term to the adjoint consistency, and the last term is the Nitsche penalty term which explicitly controls the accuracy at which the Dirichlet boundary condition is enforced.

Selective Structural Springs
To avoid ill-conditioning of the linear system due to rigid body modes (RBM) of disconnected solid sub-domains surrounded by the void phase (see Fig. 2 (a)), selective structural springs are introduced.This approach was initially proposed by [21] for fluid problems and is extended to nonlinear structural mechanics in this work.To identify disconnected solid domains "floating" within the void domain, an auxiliary indicator field modeled as a linear diffusion problem is introduced within the solid domain.Alternatively, globally applied structural springs could be used, as done by [19].The weak form of the governing equation for the auxiliary scalar indicator field θ is formulated as: where δ θ are the admissible test functions and κ denotes the thermal conductivity.The bulk heat transfer coefficient is denoted by h and θre f is the reference indicator value.The diffusion problem (18) along with the appropriate boundary conditions will lead to indicator values of close to 0.0 in domains that are connected to where the Dirichlet boundary condition of θ = 0.0 is applied and indicator values close to 1.0 in disconnected sub-domains.It should be noted that a non-physical "bulk" convection term in Eqn.(18) prevents ill-conditioning of the linear diffusion system even if disconnected sub-domains exist.After the auxiliary indicator field solution has been obtained, a smoothed Heaviside projection function is used to enforce a 0-1 indicator field.The weak form of the governing equation of the nonlinear structural problem ( 1) is augmented by the following selective spring stiffness term: where r is the relative spring stiffness ratio, and E k is the Young's modulus of either solid material (k = active or passive).
For more details regarding selective structural springs, the reader is referred to [34].In summary, the governing equation in weak form is composed of the following terms which have been introduced above, see ( 1), ( 9), ( 16), (18), and ( 19):

Formulation of Optimization Problem
The combined LS-XFEM and density approach introduced in Sec.3.1 describes the geometry and the material distribution in parameterized form.These parameters define the optimization variables s.The optimization problem considered here can be written as follows: where the nodal displacements u satisfy the discretized governing equation R = 0 and implicitly depend on the design variables s.The number of design variables is denoted by N s and the number of state variables is N u .The design variables include both the nodal LS values s φ i and the fictitious nodal SIMP densities s ρ i , such that The lower and upper bounds of the optimization variables are denoted by s L and s U , respectively.The number of inequality constraints g j is denoted by N g .In this work, a nested analysis and design approach (NAND) [35] is used where the displacements u are considered dependent variables of s and satisfy the governing equations for a given design.The advantage of this approach is that different solution algorithms can be utilized for solving the "forward" analysis problem and the optimization problem.
Displacement matching of active structures is the main objective of this work.The first part of the objective function is therefore formulated as a minimization of the squared difference between the nodal displacements u and the target displacements u tar at a specified target set Γ 0 tar ⊂ Γ 0 : MD-18-1343, Geiss, 8 It should be noted that even though a target displacement is specified a priori, no constraints with respect to the required geometry and material layout are imposed.The resulting, in general, non-intuitive geometry and material layout required to achieve a deformation that best matches the target displacements u tar is solely a result of solving the optimization problem of Eqn.(21).Selecting the target set Γ 0 tar , which in this work is defined as a subset of the design domain boundary Γ 0 , is in general non-trivial and highly dependent on the desired target deformation.Mechanical constraints like self-penetration or non-uniqueness of geometric mappings of planar structures onto the desired target deformation need to be considered in order to define a well-posed optimization problem.More details on the target displacements and the corresponding target sets considered in this work are provided in Sec. 5.
In addition to target displacement matching, a regularization term is added to the objective function in order to avoid the emergence of irregular geometric artifacts [11].Here, regularization is introduced through a perimeter penalty that is formulated as: where γ per is the perimeter penalty factor chosen such that smoothing of the interface geometry is obtained while not allowing Eqn.(23) to dominate the overall objective.
Besides the objective contributions defined in Eqn.(22) and Eqn. ( 23) two inequality constraints are imposed.The first one is a volume constraint bounding the maximum amount of solid phase allowed within the entire design domain: where γ v controls the maximum allowed solid volume Ω 0 S relative to the entire design domain volume To control the stiffness of the structure in the activated state, an end-stiffness constraint is enforced.Following the work initially proposed by [36], it is formulated as: where S is the strain energy of the system after activation (with no external load applied) and S is the strain energy after applying an additional external perturbation load.The limit in relative amount of change in strain energy with and without the final perturbation load is denoted by γ s .The end-stiffness constraint therefore requires a certain stiffness of the structure in the activated configuration in order to resist the perturbation pressure applied in the opposite direction of the desired deformation.

Design Sensitivity Analysis using the Adjoint Method
Due to a large number of design variables, the design sensitivities of objective and constraints are computed using the adjoint method.Since the mechanical model is static and conservative, the adjoint problem only needs to be solved at the end of the loading process [37].When the end-stiffness constraint is enforced, the design sensitivities for the displacement matching objective need to be evaluated at the load increment t Load while the sensitivities for the end-stiffness constraint are evaluated at t Pert .
Following the work of [24,38], the derivative of the objective z with respect to the vector of design variables s for a quasi-static case is: where the first term represents explicit dependencies while the second term represents the implicit sensitivities.Considering the two sets of nodal design variables, s φ and s ρ, and the filtering and projection relationships defined in Sec.3.1, the total MD-18-1343, Geiss, 9 derivative of ( 26) can be further expanded into: Satisfying the governing equation (20) for every design, i.e.R = 0, the derivative du/ds can be computed from: Solving Eqn.(28) for du/ds and combining it with Eqn. ( 26) yields: where the following adjoint problem can be identified: The adjoint solution is denoted by λ, which is used to finally compute the expression for the design sensitivities as: It should be noted that in the current work, the explicit contribution ∂z/∂s and the post-multiplication term ∂ R/∂s are obtained via finite differences on an elemental level.In a similar fashion as discussed above, the design sensitivities of the constraints with respect to the design variables can be obtained.For more details regarding the computation of design sensitivities with the XFEM, the interested reader is referred to [39].

Validation of the Hyperelastic Thermomechanical Model
The finite deformation thermomechanical model introduced in Sec. 2 is validated against experimental results and an analytical beam model to establish confidence in the XFEM model for design optimization.The test specimen for model     The model parameters and the material parameters of the XFEM model are given in Tables 1 and 2 respectively.Only one quarter of the domain is simulated and mechanical symmetry boundary conditions are applied along the X 1 and X 2 plane using the weak boundary condition formulation stated in Eqn.(16).The XFEM is used to describe the solid-void boundary of the bi-layer strips along the X 1 axis as highlighted in the insert in Fig. 4 (a).Note that the void domain is not shown in Fig. 4 (a) for clarity.This approach replicates a non-conforming mesh for analyzing the solid domain just as it is present during the subsequent design optimization process.The varying volume fraction of active material (Tango+ represented in blue in the top layer) and passive material (Vero represented in red in the bottom) is modeled by uniformly changing the material ratios respectively.The computational mesh consists of 1630 HEX20 XFEM elements of which 400 are intersected resulting in a total of 26,227 DOFs.HEX20 elements are 20 node hexahedral serendipity elements with a total of 60 DOFs.
The bending behavior of the bi-layer strips due to different thermal expansion of the layers is studied for five distinct volume fractions of active material.The curvature is measured at the mid-plane of the strips in the X 2 symmetry plane.A    [5] and physical experiments for different Tango+ volume fractions is presented in Fig. 5.It should be noted that due to a small sample size of 8 samples per volume fraction, error bars are not meaningful and therefore omitted when plotting the mean curvature in Fig. 5. Good agreement is achieved between the 3D XFEM model and the physical experiments, whereas the 1D model tends to overestimate the curvature, especially for large volume fractions of Tango+.This is due to assumptions made by the Timoshenko beam model which do not account for Poisson effect as well as transverse thermal expansion.Both of those phenomena are, however, present in the physical experiments where a double curvature (i.e.cylindrical bending) is observed for Tango+ volume fractions greater than 50%.These phenomena are accurately captured in the 3D XFEM model where the curvature along the X 1 direction is in fact reduced due to an increasing effect of curvature along the X 2 direction.

Design Optimization Examples
The proposed design methodology for finding the optimal design of 4D printed active structures is applied to four design problems.The LSF and the nodal densities for all design examples are discretized using trilinear shape functions.The material properties are approximated element-wise constant and are obtained using an elemental average of all nodal fictitious density values (12) together with the SIMP power law (11).In this work, β = 3.0 is used for interpolating the Young's modulus and CTE while β = 1.0 is used for interpolating the physical density.The values of the SIMP exponent β have been chosen such that a good convergence towards either active or passive material is achieved.
The projection ( 13) is applied through a continuation approach which represents a trade-off between convergence behavior and computational efficiency of the optimization process.Initially, the projection threshold parameter is set to η = 0.5 and the sharpness parameter is γ = 0.01.After a converged initial design is obtained, (e.g. after 100 design iterations), the sharpness parameter is increased to γ = 3.0 and the end-stiffness constraint is enforced.The projection sharpness parameter is then doubled every 100 design iterations.This procedure is repeated four times until γ = 48.0 and a sufficient approximation to a bi-material design within the solid domain is obtained.After the overall geometry has converged to an optimum (e.g. after 100 design iterations), the LS design variables stay unchanged while the material distribution within the solid domain is further optimized, to yield a bi-material design with a certain end-stiffness.For both sets of design variables (nodal LS values s φ i and nodal fictitious densities s ρ i ), a smoothing radius of r s = 4.0 mm is used.The optimization problem ( 21) is solved using the gradient-based nonlinear programming scheme Globally Convergent  For all examples, the initial design is a flat square plate composed of uniform, intermediate material of density s ρ = 0.5 in Ω 0 S with an initial seeding of square holes representing Ω 0 V , see Fig. 6 (b).The size of the cuboid inclusions is r i = [14.5,14.5, 14.5] mm for all examples.Different target displacements are prescribed to certain subsets of the domain boundary Γ 0 in order to achieve target deformations of varying complexities.A design with the desired mechanical response is found through optimizing the solid-void geometry as well as the active-passive material distribution within the solid domain by solving the optimization problem outlined in Sec.3.4.The material properties used for all design examples are listed in Tab. 2 and parameters specific for each design example are listed in Tab. 5.
For all subsequent design problems, quarter symmetry of the mechanical problem about the X 1 and the X 2 plane is exploited by analyzing only one quarter of the design domain and enforcing appropriate mechanical boundary conditions weakly via Nitsche's method (16).In addition, for the design problems of Sec.5.2, Sec.5.3, and Sec.5.4, design symmetry about the X 1 − X 2 diagonal is introduced by assigning a set of independent nodal LS values and nodal fictitious density values to only one eighth of the total design domain.The initial design used for all design optimization examples along with the symmetry boundary conditions is depicted in Fig. 6 (b).
The LSF of the initial design is shown in Fig. 6 (a) and computed as a signed-distance function from an array of cuboids defined as: where i is the number of individual cuboids, r i is the radius of the i-th cuboid, and Xi is the position of the center of the i-th cuboid in X 1 , X 2 and X 3 direction respectively.The roundness parameter is set to n = 100.0.The parameters of this arbitrary initial design have been determined through numerical studies to minimize the dependence of the final design on the initial geometry.
The plane creating the zero LS iso-contour is shown in gray (Fig. 6 (a)) and the corresponding solid-void material layout is depicted in Fig. 6 (b).It should be noted that LS design variables are only defined on the top surface of the initially flat plate, such that φ(X 1 , X 2 , X 3 ) = φ(X 1 , X 2 ).The LS values of all nodes which are not on the upper surface are dependent on the corresponding design variable hosted by the respective node on the top surface.This guarantees to always have a vertical cut in X 3 direction represented by the XFEM (see insert in Fig. 6 (b)) and it prevents the optimizer from locally reducing the plate thickness during the optimization process.
The target displacements are monitored at target sets Γ 0 tar defined at the tips of the structure (along the symmetry boundaries).To guarantee that these tips are mechanically connected to the base of the structure, the LS field along the symmetry boundaries is fixed to be negative (see Fig. 6 (a)).This means, strips along the symmetry boundaries are excluded from the LS-XFEM design domain and remain solid throughout the optimization process.However, regarding the density optimization, this domain is still considered design domain in which the material distribution can be altered.
The nonlinear thermomechanical model is discretized by quadratic HEX20 XFEM elements.While the in-plane discretization varies for each example, 10 HEX20 brick elements are used for discretization in thickness direction in order to accurately capture the bending behavior of the slender structures and to avoid shear locking exhibited by lower-order brick elements.An iterative Newton-Raphson scheme is used to solve the nonlinear problem that is considered as converged when a relative nonlinear residual norm drop greater 1.0 • 10 6 is achieved.Convergence is facilitated by the adaptive load-stepping MD-18-1343, Geiss, 13  [42,43].
To improve the computational efficiency of the optimization approach, the converged nonlinear solution of the state variables of the previous design is used as an initial guess for analyzing the current design.This approach is well suited for problems with static conservative mechanical models and incremental design changes where the displacement solution of the current design only differs slightly from the displacement behavior of the previous design [44].If the design change and the resulting change in the displacement response is too large, i.e. no equilibrium configuration is obtained using the converged solution of the previous design as an initial guess, the design is analyzed by simulating the entire load path.
In order to stabilize disconnected solid sub-domains throughout the optimization process, selective structural springs as introduced in Sec.3.3 are used.The thermal conductivity in this auxiliary diffusion problem is set to κ = 10.0, the bulk heat transfer coefficient is h = 0.01 and the reference indicator value is θre f = 1.0.A uniform initial indicator value of θ0 = 1.0 is applied to the entire domain and a relative spring stiffness ratio of r = 1.0 • 10 −6 is used.Adiabatic boundary conditions on the auxiliary indicator field are assumed on boundaries where no displacements are prescribed.To further improve the numerical stability of the proposed approach, a staggered solution algorithm is employed.First, the linear diffusion problem of the auxiliary indicator field ( 18) is solved and subsequently the nonlinear thermomechanical problem (1) is solved in a one-way coupled manner.
The final design for each example problem is fabricated and activated using the direct 4D printing method introduced by [5].As described in Sec. 4, for fabrication of the optimized shape-changing structures, a Stratasys Keshet J750 multimaterial 3D printer is used to deposit the Tango+ and Vero material accurately onto a build tray.The zero LS iso-contour is used to extract the solid-void boundary of the final design, while an iso-volume created along a fictitious density threshold of ρe = 0.5 is extracted to create distinct active and passive material domains.These post-processing steps are performed in ParaView [45], as it provides a convenient interface between the ExodusII mesh format and many computer-aided design (CAD) file types.In the current work, a stereo-lithography (STL) mesh file is extracted from the converged optimization result (provided in ExodusII format by the employed LS-XFEM optimization framework) for each material sub-domain.The STL file format is commonly supported by 3D printing software and therefore used to import the CAD data into the Stratasys printing software.In the Stratasys pre-processing environment, Tango+ is assigned to the active material domain represented by 0 ≤ ρe < 0.5, while Vero is assigned to the passive material domain corresponding to 0.5 ≤ ρe ≤ 1.The active structures are printed in a flat configuration and only deform upon release of the inelastic printing strain in the Tango+ material.This is triggered by submerging the initially flat, bi-material structures into a water bath of 65.0 • C which causes the stiffer Vero material to soften due to its glass transition temperature at around 53.0 • C.This allows the built-in compressive printing strain in Tango+ to be released resulting in a shape-change.The deformation is made permanent after the structure has cooled off to room temperature.A more detailed discussion of the direct 4D printing process as well as the activation steps can be found in [5].
It should be noted that all experimental results are shown for qualitative comparison only as a quantitative comparison is beyond the scope of this paper.The experimental verifications show the feasibility of the numerically determined designs and demonstrate the applicability of the proposed design framework to solving real-world design problems in a systematic manner.

Twisted Figure-Eight Example
The first example demonstrates the proposed computational design framework with the design of an initially flat plate that deforms into a twisted figure-eight like structure.The target displacement is tracked at the following target set Γ 0 tar = [X 1 , X 2 ]; at X 1 = [80.0,0.0, 0.0] mm the target displacement is set to u tar = [−80.0,0.0, +50.0] mm and at X 2 = [0.0,80.0, 0.0] mm the target displacement is set to u tar = [0.0,−80.0, −50.0] mm (see Fig. 6 (b)).This target displacement corresponds to a structure where two of the target points meet at the top center, and the other two target points meet at the bottom in the center of the full domain.The dimensions of the (quarter) design space and the mesh size used for this example are listed in Tab. 3. Representative for the first three design optimization problems, a total number of 47,200 design variables result from the problem setup shown in Fig. 6.The number of DOFs of the XFEM model reduces from initially about 126,000 to 30,500 at the lowest due to an expanding void phase which is excluded from the XFEM analysis.
The final design, obtained after about 624 design iterations is shown in Fig. 8 (a) in the activated, deformed stage (Tango+ is shown in blue, Vero is shown in red, void phase is not shown).A clear solid-void interface defined by the zero iso-contour of the LS field can be seen along with an almost binary material layout within the solid phase.A small fraction of elements with intermediate fictitious densities remains which is attributed to the finite transition regions between active and passive material.Using a yet higher projection, smaller filter radius or finer background mesh would help to further reduce those intermediate density domains, however, their effect on the performance of the final design is negligible.Fig. 8 (b) shows the corresponding printed and activated structure with good qualitative agreement compared to the XFEM prediction.Both the numerical XFEM model as well as the physical sample assume the twisted figure-eight target deformation well.In the printed sample, Tango+ is transparent and Vero is white.
MD-18-1343, Geiss, 15 A slight overlap of the tips of the four legs is seen in the numerical prediction.This is due to the lack of a self-contact formulation which is beyond the scope of this work.Intuitively, a perfectly symmetric design is expected in order to meet the twisted figure-eight target deformation.A slight asymmetry can, however, be observed in Fig. 8, for both the XFEM simulation and the physical experiment.This is attributed to a local minimum identified by the GCMMA during the design optimization process.The increasing projection parameter then locks the design into this local minimum as the design space becomes increasingly more non-convex with a higher projection parameter γ.
The evolution of the normalized objective z/z 0 , where z 0 is the objective value of the initial design, the volume constraint ( 24) and the end-stiffness constraint (25) are shown in Fig. 7 for this first design optimization example.A similar behavior is also observed for all other numerical examples.After initial oscillations in the objective, which are due to the violation of the volume constraint, a smooth and converging behavior is seen.The slight increase of the objective every 100 iterations is due to an increasing projection parameter γ which makes the design space increasingly more nonlinear.Therefore, the optimization algorithm needs a few design iterations to minimize the objective again until convergence is obtained.It is observed that even though initially the volume constraint is active, it becomes inactive within the first 20 design iterations.The motivation for the solid strips to become thinner is to reduce the effect of double curvature due to the isotropic inelastic printing strain which impedes on the primary bending behavior.Therefore, the volume constraint stays inactive for the remainder of the optimization process.As mentioned before, the end-stiffness constraint is not enabled during the first 100 design iterations.Once enabled, it is initially violated but quickly becomes inactive after about 50 design iterations.Using the smoothed Heaviside projection in combination with the end-stiffness constraint greatly benefits the final stiffness of the structure as intermediate material is mitigated with an increasing projection parameter γ.
Fig. 9 (a) shows the final design in the undeformed configuration after two distinct material phases have been extracted in a post-processing step.In the current and in future examples, Tango+ is printed transparent while Vero is printed in either white or magenta, see Fig. 9 (b).Overall, this example shows a first application of the proposed method and it demonstrates that desired target deformations can be achieved upon activation of an initially flat, plate-like structure.

Cylinder Gripper Example
A similar design study as performed in Sec.5.1 is repeated for finding the optimal design of an initially flat structure which deforms into a gripper enclosing a cylinder as its target.The target displacement of the gripper is monitored at the target set Γ 0 tar which is comprised of A smooth evolution of objective and constraints is observed as discussed before.The convergence history plots are omitted here for brevity.Fig. 10 shows the final design (a) in the deformed configuration next to (b) the experimental result for the cylinder gripper.Qualitatively, the structure fabricated and activated through direct 4D printing assumes the anticipated target deformation well.From the predicted XFEM simulation it can clearly be seen that the optimizer took advantage of the possibility to use double curvature in order to achieve the cylindrical target deformation.Moreover, a nonuniform curvature in the outward pointing legs is used to best match the target deformation at the tips of the gripper.The optimal curvature is controlled by the amount of stiff, passive Vero material along the top of the gripper (shown in red).Pure Vero is also placed as fillets at the tips of the gripper to control the bending behavior and to also contribute towards a higher stiffness of the overall structure.
In Fig. 11 the corresponding initial, flat configuration of the optimal design is shown.As before, distinct material phases were extracted from the final fictitious density field in a post-processing step; see Fig. 11 (a).Fig. 11 (b) shows the printed physical specimen for the self-deforming cylinder gripper.This example successfully demonstrates that the proposed design method is not limited to finding the optimal geometry of simple target deformations, but also handles more sophisticated ones.Non-uniform double curvature is used in the optimized gripper design to obtain the desired deformation upon activation.

Four-Legged Gripper Example
The design of an active four-legged gripper is the aim of the third design optimization example.The target displacement is formulated such that gripping of an object is simulated.In order to achieve this, a target displacement of u tar = [−80.0,0.0, −50.0] mm is defined at X 1 = [0.0,80.0, 0.0] mm and a similar target displacement of u tar = [0.0,−80.0, −50.0] mm is prescribed at X 2 = [80.0,0.0, 0.0] mm.This target displacement describes a structure where all four target points coincide in the center of the design domain below its initial, flat configuration.As in previous examples, the XFEM model parameters listed in Tab. 3 are used and the parameters specific to this four-legged gripper example are found in Tab. 5.
The optimal design layout in the deformed configuration is depicted in Fig. 12 showing both (a) the XFEM prediction and (b) the physical sample of the four-legged gripper.Due to the increased complexity of the target deformation, this design problem exhibits design phenomena previously not observed in Sec.5.1 and Sec.5.2.The final design of the gripper consists of four legs which deform to match the target deformation in a non-intuitive manner.Since an end-stiffness constraint is enforced, convex-concave curvature sections are created to increase the stiffness of the four legs.These interesting features demonstrate how geometry features triggered by a mechanical stiffness constraint yield a meaningful, complex active gripper structure.Due to the finite width of the legs, self-penetration of all four gripper legs is observed at the bottom center of the structure.This is caused by inability of the XFEM model to account for contact.In the physical experiment, self-penetration is avoided by overlapping of the four gripper legs at different X 3 depths.However, in order to more accurately account for a behavior like this, contact should be included in future work.Overall, good qualitative agreement between the XFEM prediction and the direct 4D printed four-legged gripper is observed.
Fig. 13 shows the corresponding initial, undeformed design for the four-legged gripper.The distinct domains of convexconcave curvature can be identified in both the XFEM model (Fig. 13 (a)) as well as the printed specimen (Fig. 13 (b)) where the bi-layer beam composition is locally inverted at about the half-way point of each gripper leg.

Elevated Plane Example
The last example used to demonstrate the capabilities of the proposed multi-material topology optimization framework is the design of a self-elevating plane.For this design example, the domain size is increased to 120.0 × 120.0 × 1.0 mm for the quarter domain.This increase in in-plane domain size improves the ability to capture double curvature with the given eigenstrain provided by the active material.A mesh size of 48 × 48 × 10 elements is used to discretize one quarter of the design domain, leading to initially 276,000 free DOFs.To achieve self-elevation of the domain in the center of the design space by 40.0 mm in the X 3 direction, a target displacement of u tar = [0.0,0.0, +40.0] mm is defined at the center of the domain spanned by X 1 = [0.0≤ X 1 ≤ 15.0, 0.0 ≤ X 2 ≤ 15.0, 0.0] mm.In addition, the out-of-plane displacement is fixed at  The final design of the self-elevating plane in the deformed configuration is shown in Fig. 14.Both, (a) the XFEM prediction and (b) the direct 4D printed physical specimen are depicted.The target deformation is achieved well in both simulation and physical experiment.A cross-shaped structure attached to the platform at the center is created.Each member of the cross exhibits a non-uniform active-passive material layout to yield distinct domains of concave-convex bending upon activation.This change in curvature is necessary to elevate the center of the structure in X 3 direction while introducing no out-of-plane displacement and rotation at the tips.
Fig. 15 shows the undeformed structures of the elevated plane corresponding to (a) the XFEM simulation and (b) the printed sample before activation.A thresholding scheme is employed to eliminate any remaining intermediate densities and to clearly identify active and passive material domains required for direct 4D printing.

Conclusions
A topology optimization approach for designing direct 4D printed, shape-changing structures undergoing large deformations was proposed.A combined LS-XFEM and density-based topology optimization approach was introduced in order to describe the multi-material optimization problem.The LS-XFEM was employed to describe the solid-void domains in a crisp manner, while the SIMP method was used within the solid domain to distinguish between active and passive material sub-domains.The hyperelastic thermomechanical model was discretized by quadratic displacement elements using an XFEM formulation.The inelastic printing strain introduced through the direct 4D printing process was modeled as a residual isotropic eigenstrain.Accurate prediction of the behavior of self-deforming structures by the proposed large strain thermomechanical XFEM model was verified by comparing numerical results against physical experiments and an analytical beam model.
The capabilities of the proposed design optimization methodology were demonstrated through four example problems.The objective was to match a given target displacement subject to a volume constraint and an end-stiffness constraint in order to assure structural integrity of the final design.Optimal designs matching target displacements of a twisted figure-eight design, a cylinder gripper, a four-legged gripper, and a self-elevating plane were successfully obtained using the proposed framework.Geometries and material arrangements of increasing complexity were created that take advantage of mechanical phenomena such as double curvature, locally concave-convex curvatures, and domains of uniform passive material in order to meet the end-stiffness constraint.
After having demonstrated the design capabilities of the proposed framework, there are a few shortcomings that need to be addressed in future work.These include the simplifications made with regards to structural instabilities, such as buckling and snap-through, which are a common phenomenon in slender structures.These structural instabilities have been avoided by using sufficiently thick structures in the current work, but need to be considered in general, when large compressive stresses are present in the structure initiated by the active material.Previous works in the field of modeling of slender structures incorporating eigenstrains include [44] and [46] where especially [46] and [47] consider fully nonlinear behavior including instabilities.Future work should study the influence of structural instabilities on the optimal design of 4D printed structures and develop a systematic optimization approach to either leverage them as a desired design feature or to avoid them in order to increase design robustness.Moreover, the third example experienced a significant amount of self-penetration in the activation stage due to the lack of a self-contact formulation within the employed XFEM framework.Incorporating an XFEM contact formulation with multi-material topology optimization as studied by [20] into the proposed design framework is another topic which needs to MD-18-1343, Geiss, 19 be addressed in future work.
This paper has presented an initial demonstration of a systematic design approach for 4D printed structures undergoing large deformations.In addition, design problems with more complex target displacements should be studied in future work.
Physical specimen for all design examples were fabricated using direct 4D printing.Upon activation, qualitatively good agreement between the XFEM prediction and the physical response was seen.Discrepancies between the numerical results and the experiments are mainly attributed to minor anisotropy of the material properties with respect to the print direction, ambient effects like gravity, and viscosity of the water bath on the soft structures during the activation process.None of those physical effects have been accounted for in this work.However, in order to achieve better agreement between XFEM simulations and 4D printed specimen, all of those effects should be accounted for in future studies.paper are those of the authors and do not necessarily reflect the views of the sponsoring organization.

Fig. 1 .
Fig. 1.Conceptual steps of the design process from (a) the definition of a target displacement to (f) an activated structure using direct 4D printing.

Fig. 2 .
Fig. 2. (a) Design domain decomposition into solid (Ω 0 S ) and void (Ω 0 V ) sub-domains using LS-XFEM along the interface Γ 0 S,V .Within the solid domain, a further distinction is made between an active material (Ω 0 S A ) and a passive material (Ω 0 S P ).(b) Multiplicative decomposition of the total deformation gradient into thermal and mechanical deformation gradient.

Fig. 3 .
Fig. 3. Adaptive load-stepping technique used to solve the nonlinear problem.(a) Temperature load profile, (b) follower load profile and (c) displacement response over pseudo time.

Fig. 4 .
Fig. 4. Bi-layer strips with different active material ratios (in %) used for XFEM model validation.The solid-void interface of the strips along the X 1 axis is cut by the XFEM and tetrahedralized for volume integration.(a) Simulation results and (b) experimental results.

K
Poisson Ratio of Tango+, Vero ν T = ν V = 0.4 validation are rectangular, bi-layer strips with different volume ratios of active versus passive material, leading to different curvature values.The bi-layer composite is made out of Tango+ as the active elastomer in the top layer (Ω 0 S A ) and Vero as the passive glassy polymer in the bottom layer (Ω 0 S P ) printed on a Stratasys Keshet J750 multi-material 3D printer, with a printing layer thickness of 27 µm.The XFEM predictions of the bi-layer strips are shown in Fig. 4 (a) while the printed and activated specimen are shown in Fig. 4 (b).

Fig. 5 .
Fig. 5. Comparison of mean curvatures obtained by simulations and experiments for different Tango+ volume ratios.
curvature values obtained by the 3D XFEM model, a 1D Timoshenko beam model as outlined in

Fig. 6 .
Fig. 6.Initial design of a quarter of the self-deforming structure with (a) the initial LSF and (b) the mechanical boundary conditions.

Table 4 . 2 Lower
Optimization problem parameters.Perimeter Penaltyγ per = 5.0 • 10 −Asymptotes (GCMMA)[40] without inner iterations.The respective parameters of the optimization problem are listed in Tab. 4. The optimization problem is considered converged once a relative residual norm drop of the Karush -Kuhn -Tucker (KKT) conditions[41] greater 1.0 • 10 10 is achieved and all constraints are satisfied.

Fig. 7 .
Fig. 7. Evolution of (a) normalized objective, (b) volume constraint and (c) end-stiffness constraint for the twisted figure-eight example.

Fig. 8 .
Fig. 8. Final design of the twisted figure-eight example in deformed configuration.(a) XFEM prediction and (b) activated 4D printed sample.

Fig. 10 .
Fig. 10.Final design of the cylinder gripper in deformed configuration.(a) XFEM prediction and (b) activated experiment where Tango+ is transparent and Vero is white.
30.0, 75.0 ≤ X 2 ≤ 80.0, −0.5] mm and X 2 = [75.0≤ X 1 ≤ 80.0, 0 ≤ X 2 ≤ 30.0, −0.5] mm spanning the tips of the gripper.The desired deformation of the tips of the structure is described by a surface of a cylinder aligned with the X 3 axis, with radius R tar = 50.0mm and a depth of X 3 tar = −45.0mm.The XFEM model parameters listed in Tab. 3 are used.Parameters specific to this example are found in Tab. 5.

Fig. 12 .
Fig. 12. Final design of the four-legged gripper example in deformed configuration.(a) XFEM prediction, (b) activated experiment.In the printed specimen, Tango+ is transparent and Vero is magenta.

Fig. 14 .
Fig. 14.Final design of the self-elevating plane in the deformed configuration.(a) XFEM prediction and (b) direct 4D printed specimen.In the printed sample, Tango+ is transparent and Vero is magenta.

Table 2 .
Material parameters of Tango+ and Vero.

Table 3 .
XFEM model parameters for design optimization problems.

Table 5 .
Parameters used for each design example.Sec.2.2 where the time step reduction factor is set to f = 0.25 and the maximum number of nonlinear iterations is set to 40.Again, the importance of the adaptive load-stepping scheme with respect to the convergence of the nonlinear problem during the optimization process should be emphasized.Especially for designs with large amounts of weak material, i.e. material with intermediate densities, adaptively reducing the load step is crucial in order to facilitate convergence of the Newton-Raphson scheme.The linearized sub-systems are solved using the Multifrontal Massively Parallel Solver (MUMPS)