Topology Optimization with Matlab: Geometrically Non-Linear Optimum Solid Structures at Random Force Strengths

: This paper aims to investigate multiple large-strain topology-optimized structures, by interpreting their overlay as a probability density function. Such a strategy is suited to ﬁnding an optimum design of silicon electrodes subject to a random contact. Using this method, and prescribing a zero net-force constraint on the global system, the optimum structure is identiﬁed with a Schwarz P minimum-surface structure. Then, the optimum structure is subject to chemo-mechanically coupled cycling, in terms of an irreversible thermodynamic process, which shows the interplay between the mechanical and chemical ﬁelds. The Matlab-based optimization code is attached


Introduction
More than 60 years ago, Lucien Schmit [1] studied how to combine optimization techniques with the finite element method (FEM), and later, Bensøe and Kikuchi reported their seminal paper [2]. Since then, topology optimization techniques have been typically used in lightweight design [3][4][5][6][7], where only a percentage of a given volume is occupied by material, while meeting similar overall system stiffness requirements as the entire body. Technological advances in additive manufacturing in particular, make it possible to create structures that exceed conventional design [8]. However, optimum designs result predominantly from prescribed static boundary condition problems [9,10], even if the superordinate system is based on dynamic motion, e.g., a bicycle chain [11]. In such a system, a tensile force acts on the individual link of the chain. The optimum shape of a plate shows the positions for the bolted connections, see Figure 1. The resulting structure might be subject to a low volume constraint, so recent optimization algorithms also account for, e.g., buckling [12], fracture [13], or multi-material demands [14]. In addition, the optimum structure might be part of a periodic structure (the chain) [15], and other manufacturing constraints could be implemented by a projection-based approach [16]. However, while it is similar, it is also different from the silicon (Si) electrode particles in batteries, which are exposed to large deformations. According to [17], up to now, only a few structural designs, such as porous, core-shell, yolk-shell, and solid nanostructures, have been proposed for silicon anodes, some of them inspired by nature. However, recently, the topology optimization of electrodes has gained more and more attention [18,19]. For the anode design and the improvement of its electrochemical performance, it is crucial to identify the repeated volume expansions (up to about 300% [20]) and contractions during lithiation and de-lithiation, as the key failure mechanism of Si. Particle swelling leads to structural pulverization, electrical disconnection between the active materials and the current collector, continuous exfoliation, and consumption of the active Si and electrolyte of active Si particles [17,21,22]. All these failures can be assigned to the topology of the anode: (i) the increased consumption of the electrolyte (side reactions) could be related to a high surface area, (ii) a materials-dependent fracture occurrence above a specific particle size, and (iii) electrical conductivity issues due to point-to-point contact. As a result of the swelling, interactions with adjacent particles occur randomly, as shown in Figure 2. Most recently, contact has been added to the optimization, e.g., by [23][24][25][26]. Moreover, it was found that the contact stress is in the range of diffusion-induced stress under the free-expansion state [27]. A topology optimization problem is generally solved when finding the optimum material distribution subject to specific demands such as maximum stiffness, maximum conductivity, or others. The common algorithms rely on the finite element method, filtering methods, and the formulation of an optimality criterion [28]. Although numerous techniques exist in the scientific literature, e.g., evolutionary topology optimization (ETO [29]), smooth-edged material distribution for optimizing topology (SEMDOT [30]), floating projection topology optimization (FPTO [31]), etc., there are two popular methods: the Solid Isotropic Material with Penalization (SIMP) method and the Bidirectional Evolutionary Structural Optimization (BESO) method. The SIMP method is the most widely utilized approach, although the BESO method shows a fast convergence rate, with fewer iterations than the SIMP method [32]. In the literature, mainly the optimization problems of a linearelastic material model are solved [32,33] and some studies address non-linear problems, e.g., [34][35][36]. Recently, a two-dimensional BESO algorithm for non-linear problems has been reported [32], which claims one can hardly find convenient and compact published code.
In this paper, the aim is to elaborate on the BESO approach, in terms of (i) extending the code to 3D, (ii) making it comparable with the Matlab code reported in the SIMP approach [37], (iii) validating the comparability between both approaches for small deformations, (iv) using the BESO approach for random large force interactions as typical in silicon, and (v) defining an optimum structure at large strains by overlapping the single results. Since we are interested in large deformations, we need to recap in Section 2.1, briefly, some continuum mechanics-based essential equations of non-linear elasticity, cf. [38,39]. Then we proceed in Section 2.2 by solving the non-linear equations using FEM. Densitybased methods such as BESO and SIMP, and filtering techniques are briefly sketched in Section 2.3. In Section 2.4, we state the optimization problem at large strains. How the resulting structure must be modified to be used in a chemo-mechanically coupled inves-tigation is sketched in Section 2.5. Finally, we show and discuss the optimum results in Section 3. The 3D non-linear optimization code is provided in Appendix A.

Materials and Methods
In the following, we elaborate on the governing equations and the method for solving a non-linear optimization problem. In general, the structure of a software presupposes a certain convention and symbolic expressions, which are declared in the following.

Continuum Mechanics
Any problem in elasticity is usually based on a strain-displacement relation, a constitutive equation (stress-strain relation), a traction-stress relation, and the formulation of linear and angular balance laws. Since all equations rely on the formulation of the displacement of a material point, we embed it in the environment of continuous media first.

Displacement of Continuous Media
Let X refer to a point of a continuum in its reference ("material") configuration B 0 ⊂ R 3 , and a motion moves the point to its actual ("spatial") position x(X, t) at time t, in the current configuration B t ⊂ R 3 . Then, the corresponding displacement vector u, has the shortest distance between x and X, and is defined by their difference as From Equation (1), there follows an expression for the differential line element directly, where 1 ∈ R 3×3 is the identity matrix and F the deformation gradient. As usual, variables are denoted with capital letters when referring to the reference configuration. Moreover, the displacement has, in both configurations, the same values but different arguments, U(X, t) = u(x, t).

Green-Lagrangian Strain Tensor
Both of the differential line elements in Equation (2), dx and dX, are associated with a differential arc length, and their squares are coordinate-independent quantities. Thus, the difference between the squares of the differential line elements is considered as a measure of deformation, and defined as the Green-Lagrangian strain tensor. Together with Equation (2), there follows an expression for the Green-Lagrangian strain tensor which is based on the gradient of the displacement. As long as ∇u 1 is satisfied, we are in the regime of infinitesimal strain theory, and then the strain tensor composes additively, E ≈ ε = 1 2 (F T + F) − 1. Otherwise, the strain tensor composes multiplicatively, E = 1 2 (F T F − 1), and measures how much the right Cauchy-Green tensor C = F T F differs from 1.

Hooke's Law
In a small strain regime, Hooke's law for continuous media relates the Cauchy ("true") stresses σ ∈ R 3×3 and the strains ε ∈ R 3×3 by using a fourth-order elasticity tensor C ∈ R 3×3×3×3 , σ = Cε. Due to the inherent and material symmetries, and the assumption of isotropic media, the elasticity tensor can be reduced to only two independent Lamé parameters (λ, µ), and then the Cauchy stresses reads Alternatively, the Lamé parameters can be expressed by using the Young's modulus, E, and Poisson number, ν.

Voigt Notation
Due to the symmetry of the Cauchy stress tensor, σ = σ T , it can be condensed to an array with six entries. Then, using a superscript V, to indicate Voigt notation, Hooke's law is rewritten as Alternative arrangements of the off-diagonal stress-tensor entries are possible, e.g., as used by the authors in [37], or given by the Nye notation which is often adapted due to convention in other software coding. However, the elasticity tensor in Voigt notation is dependent on the elasticity constants (E, ν) and partitioned into four 3 × 3 blocks as where 0 = 0 1 is a zero matrix.

Cauchy Traction Vector
Now, the stresses (6) are related to infinitesimal acting forces d f per infinitesimal area of a region da, with outward normal n to the surface, both in the spatial configuration, and defined as the Cauchy traction vector t.

Nanson's Formula
The force-to-area ratio in Equation (8) refers to quantities in the current configuration. Nanson's formula links the area in the current configuration da, with its outward normal n, to the same area in the reference configuration dA, with outward normal N, by The determinant of the deformation gradient J, is the ratio of the volume elements dv and dV, from the actual and reference configurations.

Static Equilibrium
With Equation (9) at hand, the transformations of the area and volume elements from the actual to the reference configuration are given. Let ∂B t denote a surface within or on the boundary of B t . Then, the total force acting on this surface is given by f (8) = ∂B t t da (8) = ∂B t σn da (9) = ∂B 0 where the Cauchy traction vector (8), the Nanson formula (9), the divergence theorem (D), and the definition for the 1st Piola-Kirchhoff stress tensor P = JσF −T (11) have been used. P relates the forces f , in the spatial configuration, with areas in the material configuration. The forces f , can be assigned with a spatial body force density b, With Equations (10)- (12), the static equilibrium equation reads in differential form ∇ · P −B = 0(quasi-static balance of linear momentum) (13) with boundary conditions PN =T on ∂B T (Neumann boundary condition) (14) u =ū on ∂B u (Dirichlet boundary condition) (15) where ∂B 0 = ∂B u ∪ ∂B T with ∂B u ∩ ∂B T = ∅, and the quantities with a bar are prescribed. Moreover, the 2nd Piola-Kirchhoff stress tensor relates forces in the reference configuration to areas in the reference configuration.

Hyperelastic Material Model
For a hyperelastic material, the stress tensor (5) is derived from a strain energy density function Ψ, However, it is known that σ is variant to pure rotations, while ε is not. Therefore, it is not suited as a constitutive model for large deformations. In the case of finite deformations, instead of σ, often the 2nd Piola-Kirchhoff stress tensor S is used, In Equation (18), the most straightforward hyperelastic strain energy density is formulated, i.e., the Saint Venant-Kirchhoff model, which is an extension of Equation (17). This gives S = 2µE + λtr(E)1, which is similar to Equation (5). The limitations of the Saint Venant-Kirchhoff material model in large strain regimes are discussed in [40]. Alternatively, many other strain energy density formulations, such as the Ogden model [41,42], neo-Hookean model, Money-Rivlin model, Simo-Pister model [43], etc., could be adapted to match with a specific material behavior of interest.

Non-Linear Hookes Law
In the case of non-linear strain, the 2nd Piola-Kirchhoff stress tensor in Voigt notation has a similar form to Equation (6), and for infinitesimal deformations, the difference between the Cauchy and Piola-Kirchhoff stress tensors is marginal. The strain in Equation (19) has the components with indices as they are arranged for ε V in Equation (6).

Finite Element Analysis
Having formulated S in (18), the transformation of S to P in (16), and the equilibrium Equation (13) on the reference domain B 0 , it is now discretized into n el elements within the finite element analysis. In the following, the weak form of the problem (13) is formulated, ansatz functions are declared, and the Newton-Raphson method is used to solve the non-linear problem.

Weak Form of Problem
As usual, the multiplication of Equation (13) with a smooth test function vector λ, and integration over the domain B 0 , gives where (PI) denotes partial integration, the Frobenius inner product reads A : B = tr(A T B), the trace of a product tr(A T B) = tr(BA T ), and the test function vector λ is chosen so that it vanishes at the boundary ∂B u . The integral expression in Equation (21) is called the weak form of the problem.

Element ID, Connectivity, and Edof
As in the simplest case, the reference geometry is assumed to be a cuboid to solve Equation (21). This reference cuboid is subdivided into (n 1 ,n 2 ,n 3 ) equidistant elements, such that the total number of elements is n el = n 1 n 2 n 3 . The elements are numbered successively along the X 1 , then the X 2 , and lastly the X 3 direction, see Figure 3a. Each node, which is a corner stone of an element, has three degrees of freedom (dof). Therefore, one element has in total 8 · 3 = 24 dofs and the total number of dofs is n dof = 3(n 1 + 1)(n 2 + 1)(n 3 + 1). According to [37], the global node identification is ordered column-wise, with the rule: up-to-bottom, left-to-right, and back-to-front, see Figure 3b. The local node identification of the eight nodes within a single element does not necessarily follow the same convention, and is ordered in counter-clockwise direction instead, see Figure 3c. The global connectivity of each single local element, i.e., the mapping of the blue circled numbers to the red circled numbers, is stored for each element (index "el") within an array edof el = [node 1x node 1y node 1z node 2x . . . node 8z ] T ∈ N 24 (22) can be extracted for each element from the global displacement vectorū = [ū 1ū2 . . .ū n dof ]. The underbar in the expression•, indicates that a quantity • refers to the ordered nodal values.

Parametric Domain
Now, the solution field inbetween the nodal values must be interpolated. Let a local coordinate vector ξ = (ξ 1 , ξ 2 , ξ 3 ) T and the parametric domain B ξ = {ξ ∈ R 3 |ξ 1 , ξ 2 , ξ 3 ∈ [−1, 1]} be given. As we already performed the mapping from the current to the reference domain, B t → B 0 in Equation (10), by using the integration by substitution We must map each element to the reference domain B ξ , to perform Gaussian quadrature. In the following case, it is an affine transformation between the reference and the parametric domain, where h is the equidistant nodal spacing in the physical domain (see Figure 3a) and 2 is the length of the parametric domain in each direction. Without loss of generality, we set h = 2 in the following.

Ansatz Functions
The ansatz functions for both the displacement and the test function on each element are written in terms of the local coordinate vector ξ, as where it is necessary to note that the entire dependency of u el and λ el on ξ, is carried by the interpolation functions N i only. For computational reasons, often low polynomial degree functions are used, e.g., linear Lagrangian interpolation functions, which read and i ∈ {1, . . . , 8}. In accordance to the node identification of Section2.2.2, the sequence of s i In the case of derivatives, one should note that

Strain-Displacement Matrix
With the ansatz functions at hand, and by using Equation (26), the weak problem (21) 4 reads on each element A quantity called the "non-linear nodal strain-displacement matrix" is introduced, where F ij are the entries of the deformation gradient F el , and ∂ j ≡ ∂ ξ j is the derivative in the jth direction. The elemental strain-displacement matrix is composed by the nodal matrices as Alternatively, one can show that the elemental variation in the strain δE el = B el δu el relates the non-linear B el -matrix to the variation in the elemental displacements.

Residual of Equilibrium Equation
Now, we can elaborate on the weak form (21) on element basis, where the term in parentheses is identified as elemental residuum R el , which in turn is given by the difference between internal and external forces on that element, and both, the internal and external elemental forces read Since the nodal values λ i , of the test function are arbitrary, we can find a residual expression Practically, the assembly of the global residual vector R, is given by adding the local values on the corresponding global position as

Tangential Stiffness Matrix
The tangential stiffness matrix follows from a variation in the residuum, where G el stores the linear derivatives of the ansatz functions, and the components of the 1st Piola-Kirchhoff stress tensor S ij are stored in Again, similar to Equation (33), the assembly of the global tangential stiffness matrix reads K(edof el , edof el ) = K(edof el , edof el ) + K el ∀ el ∈ {1, . . . , n el } 2.2.8. Gaussian Quadrature Equation (38) has polynomial degree d = 3, so exact integration of some function f , dependent on a polynomial expression of ξ, by summation, is given if the number of Gauss points is n gp = d+1 2 = 2, at positions ξ gp = ± 1 √ 3 and weights w gp = 1.

Newton-Raphson Method
The nodal displacement is found iteratively by using the Newton-Raphson method, where it is the iteration number. To calculate the displacement update ∆ū it , we use a direct solver until ∆ū it ū it < 0.001 as the termination condition.

Compliance
The compliance of the system is defined in accordance with the strain energy as with tangential stiffness K, from Equation (41). Sometimes compliance is defined without the factor 1 2 [37]. In terms of an optimization problem, it relates a maximum stiffness demand by means of a minimum compliance.

Density Based Methods and Filtering Techniques
Although the term density, suggests that the method relies on a quantity extended in space, the density variable x i ∀i ∈ {1, . . . , n el }, is considered only at the element's center of mass. It is conceived as an effective quantity during the FEM calculation. The SIMP and BESO methods are based on the formulation that the density interpolates Young's modulus, see Equation (7). We first review the core idea of the SIMP method, then explain the filtering methods used in the optimization procedure, and show parallels with the BESO method.

SIMP Method
There are many different approaches to keep the binary system, consisting of material and void, from becoming singular in the stiffness tensor, by ascribing either a minimum density x min , or a minimum Young's modulus E min > 0, to the voids. Within the (modified) SIMP approach the Young's modulus is given by a convex combination between the solid material (E 0 ) and the void material (E min 0) which is a function of a (normalized) densityx, with penalization power p > 1. Here, it is important to distinguish between the (relative) element density x and the filtered (physical) densityx, as described in the following.

Density Filter
According to [44], a basic linear filter density function is defined as where the weights should fulfill the partition of unity, ∑ j w ij = 1 ∀i, and r min is an (user) input radius defining the neighborhood N i = (j : dist(i, j) ≤ r min ). In general, the weights influence the resulting structure, and since there is no overall optimum result, they are somehow arbitrary [45]. Other filter types could be defined by Gaussian weights [44], Heaviside filter, dilate and erode filter [46], or filter based on the geometric mean [45]. Nevertheless, what is known for sure, without filtering, the optimum result tends to become a scattered structure, often referred to as a checkerboard, with thin structural parts and/or many tiny holes [45].

Sensitivity Filter
A similar procedure to Equation (47) is applied to the sensitivities, which are the derivatives of the objective function, and reads which is strictly speaking not a density filter because no filtering is conducted regarding the density itself. The sensitivities are used together with the design domain,ṽ = [v 1 , . . . , v n ] T and a prescribed volume limitv, to solve the Karush-Kuhn-Tucker (KKT) condition, which leads to an optimality condition B i = 1, where should be fulfilled. The value of the Lagrange multiplier l i , satisfying the condition v(x(x new (l i ))), is the only unknown, and is found by the bisection method. With this, the update reads where m is a positive move-limit, and η is a numerical damping coefficient.

BESO Method
The BESO method relies also on the penalty method of the Young's modulus, but here, the sensitivity of the objective function is assigned with the total strain energy, α i , of the removed element, Then, the element sensitivity is given by the linear density filter (47), and reads To improve the convergence, [32], the sensitivity at the current iteration step (it) is averaged with that from the previous step, Introducing the two threshold constants, α th del for material removal and α th add for material addition, the update reads The threshold constants are found by the bisection method. Moreover, the target volume where c er is the evolutionary ratio determining the percentage of material to be removed. The evolutionary ratio is a number initialized in the previous iteration, and is zero when the target volume is reached. Then, only the topology alters.

Problem Formulation
Finally, the optimization problem to be solved reads where in the non-linear case, the displacementū is found iteratively, by using the Newton-Raphson method. The tilde in the expression•, indicates that a quantity • refers to the ordered elemental values. The optimum structure is found if a termination criterion is fulfilled, e.g., if a certain number of iterations is reached (failure) or the density/compliance difference is marginal between two iteration steps (success), as explained in the following.

Termination Criterion: Density Change
The first termination criterion is density-based, and states that an optimum structure is found if the change between the densities according to the L ∞ -norm is sufficiently small, ≈ 0.001. With this at hand, the updated structure would not change in the next iteration step regarding the previous one.

Termination Criterion: Compliance Change
Another termination criterion is based on the change in compliance, defined as a moving average | ∑ Here, the size of the moving window is somehow arbitrary, and found empirically.

Chemo-Mechanically Coupled Model
If the optimum structure is found, we need a chemo-mechanically coupled material model to test the particles' performance. In [47], we formulated a multi-field model for charging and discharging of lithium-ion battery electrodes. Briefly sketched, the model is based on a multiplicative split of the deformation gradient into an elastic (index e) and purely volumetric inelastic (index i) deformation. The inelastic part should depend on the lithium (Li) concentration c Li , and reads whereΩ is a normalized partial molar volume and c 0 a reference concentration. With this at hand, the total free energy density of Equation (18) is now, besides the right Cauchy-Green tensor C, dependent on another field variable c Li , where Ψ ela (C) is the elastic strain energy, e.g., Ψ ela (C) = Ψ E from (18) but now with concentration dependent elastic moduli, Ψ con (c Li ) is a double-well potential and Ψ int (∇c Li ) accounts for interfacial reactions. Moreover, besides the quasi-static balance of linear momentum, the Cahn-Hilliard equatioṅ is an extension of Fickean diffusion, i.e., accounts for the irreversible temporal evolution of the additional concentration field and follows from mass conservation; M, is called a mobility tensor. The Cahn-Hilliard equation is a fourth-order partial differential equation and needs, in the sense of finite element analysis at least, continuous differentiable ansatz functions. Therefore, the obtained optimum structure should be re-meshed by using quadratic B-Splines, and we follow the technique as proposed in [48]. Re-meshing is needed to evaluate the optimum structure in a post-process only, and thus does not affect the mesh during optimization.

Results
Here, we first compare the results obtained from linear and non-linear optimization. Then, we elaborate on the force limit, where the results are comparable (within a 3% deviation). With this, we show the simulation results obtained from a single force load. In the next step, we generate structures that result from random force strengths, subject to different constraints on the net force. Finally, all structures' (normalized) overlays give the most probable shape of a resistant structure to a random contact. This structure is subject to a periodic mechanical motion, and the induced chemo-mechanically coupled change in lithium concentration is studied.

Reference Parameter Definition
The optimization code listed in Appendix A, relies, according to [32], on the user input (n x , n y , n z ,v, c er , r min ) and some constants anchored inside the program (see Table 1). If there is no user input, the parameters from the calibration calculation from Section 3.2, are passed automatically. Moreover, the common linear density filter (47) is used in all numerical examples, but can be exchanged by the user if necessary.

Calibration
For the calibration of the optimization algorithm, we first check that the results obtained from the small strain BESO method are comparable with the small strain SIMP method, in terms of a binarized density at threshold x i = 0.5 ∀i ∈ {1, . . . , n el }. To this end, we choose the well-known cantilever beam, and set the number of elements to (50,16,4), see Figure 4a. Next, we declare a desired volume of 50% and the normalized force strength to be 10 −9 . The convergence plot is shown in Figure 4b and is very similar to the 2D case reported in [32]. The resulting structure from the BESO method is used as input structure for the SIMP optimization calculation, with similar parameters, i.e., we adapted E 0 = 1, p = 3, r min = 1.5,v = 0.5, and adjusted E min = 3.1 × 10 −2 (see Equations (45) and (46)) and tolx = 0.015. Since there is almost no difference (two elements change) between the SIMP and the BESO solution at small force strengths, we declare this structure as a reference solution. Then, we increase the force strength and define a deviation D, by element-wise density comparison. Finally, we count only those elements (dark cubes) that do not match the reference solution (see Figure 4c), i.e., the higher D, the more the structure deviates from the reference solution. As expected, the higher the force strength, the higher the deviation.

Central Single Force
Now, the reference geometry is given by a (16,16,16) cube, which is subject to a fixture by Dirichlet boundary conditions on the bottom (see blue dots in Figure 5a), a single force acting on top of the surface, and a prescribed volume, of 30% and 50%, respectively. The optimization converges very fast (see Figure 5c) and leads to a conical structure (see Figure 5b), which is almost independent of the force strength applied. Conical structures have already been investigated experimentally in the context of Li-ion battery anodes, by [49], and numerically by [50]. Nevertheless, this kind of structure is resistant to one-sided loadings only, and therefore, we proceed with all-round random contact in the following.

Forces with Random Strength and Zero Net Strength
Now, we consider again the (16,16,16) cube, which is fixed at the center of mass (see blue dots in Figure 6a). We apply forces with random strength, but overall, they are constrained to zero net force (see Figure 6b). Doing so makes the resulting structures look arbitrary (see Figure 6c). The observed island phenomenon is likely due to the applied boundary conditions in the center of the geometry. Due to the randomness of the force strength, the optimum structures are not comparable if repeating the simulation. Especially in Figure 6d, we have omitted the structure's walls, which were mostly filled with material. Therefore, we must consider a statistical interpretation of the internal structure, as explained in the following subsection.

Probability Density
Since the underlying domain is always kept the same, we can define a probability for each element, which indicates how often the optimized result in the specific domain is provided with material. The cube's core has the highest probability that it must be filled with solid material, which decreases towards the cube's corners. Instead of the prescribed volumev, we can now consult the probability quantile as the superior optimization constraint, which we had to declare for the post-processing. Regarding the color distribution in Figure 7, the most likely structure is a ball (yellow region), very unlikely a fully filled cube (dark blue), and a Schwarz P minimum surface delimits the transition area to an intermediate structure. The Schwarz P surface itself can be described analytically by the function cos(X 1 ) + cos(X 2 ) + cos(X 3 ) = 0 (63) In the following, we will mesh it as a solid and expose it to chemo-mechanical cycling.

Periodic Loading of Schwarz P Structure
The topology-optimized structure is now subjected to sinusoidal cyclic loading, and an initial value problem is solved using the chemo-mechanically coupled model. In Figure 8, the concentration distribution is shown, where the concentration is pushed outward in a compressive motion and vice versa. Since the constitutive law follows the thermodynamics of irreversible processes, it can be seen that after the first mechanical cycle, the concentration distribution does not return to the uniformly distributed original state, but an induced pattern persists.

Forces with Random Strength and Non-Zero Net Strength
Previously we assumed the constraint of a vanishing total force. Now, a net force on the object is prescribed, and then, the optimum topology changes to an eighth of a bee-shaped structure, with spikes in the diagonal corners, see Figure 9.

Discussion
The core question of this paper, was related to the optimal geometry of a silicon anode exposed to a random contact. It is known that silicon anodes swell to multiples of their original volume. Thus, depending on the charging behavior, the interaction between the particles is somewhat random, regarding the force strength and the local position of the occurrence. We have developed and provided a compact Matlab code to address the question of an optimal geometry that is as resistant as possible to this type of force action. The code is based on the two papers by [32,37]. The published code of the first reference is known for its particular simplicity of implementation of boundary conditions, and the second reference for its fast convergence rate, whereas the latter code had been published only in 2D. In addition, all studies we have performed using the new code, provide symmetric structures under symmetric loading, which is infrequent in the literature. Then, we subjected a cube to random loading but obtained structures challenging to interpret. However, we considered the normalized overlap as the probability density when overlapping the structures resulting from many simulations. As a result, we obtained the symmetric Schwarz P structure at a specific threshold, which has a minimum surface property. Such results sound promising in the context of electrodes, where structures with minimum surface area are desired. Furthermore, this method also yields what one would intuitively expect, e.g., a structure concentrated along the diagonal in the case of diagonal loading. Finally, the material model is chosen, so that it can be extended in the future by a crack density field, for a crack propagation study.

Conclusions
Here, we elaborated on large, geometrically non-linear deformations of battery anodes, intending to find an optimum (maximum stiff) electrode design when subject to random contact, in terms of random forces. The optimization technique relies on the BESO method, where the solution procedure includes non-linear finite element analysis, specific sensitivity filtering, and an improved weighting of the historical information, to ensure fast convergence. The existing 2D Matlab code reported in [32] was extended here to 3D, attached in the appendix, and applied to several design problems. Even at moderate volume fractions, the resulting topologies are hard to compare with each other. Thus, in such cases, the optimum structures are found by performing multiple simulations and interpreting their (normalized) overlap as a probability density.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: