This paper describes the actual stage of the research work dedicated to the investigation, development, implementation and validation of a new-advanced shape optimization approach, particularly tailored for two-dimensional electric structures with complex geometry. The new proposed numerical approach is based on the efficiency of the Extended Finite Element Method (XFEM) and the flexibility of the Level Set Method, to handle moving interfaces without remeshing at each optimization step. The result is a powerful and extremely robust numerical shape optimization algorithm with a superior degree of automation in comparison with other methods reported in the literature. The new methodology is able to handle effortlessly the topological changes as well as to provide a high fidelity representation of the deformable interfaces. A test problem consisting in the minimization of the electric resistance of a corner type resistor is given as an example. The obtained results suggest that combining XFEM with the level set based boundary representation is a very promising approach for shape optimization of complex electric structures.
This paper describes the actual stage of the research work dedicated to the investigation, development, implementation and validation of a new shape optimization approach, particularly tailored for two-dimensional electric structures.
Design/methodology/approach
The proposed numerical approach is based on the efficiency of the Extended Finite Element Method and the flexibility of the Level Set Method, to handle moving interfaces without remeshing at each optimization step.
Findings
This new approach eliminates the conventional use of discrete elements and provides efficient, stable, accurate and faster computation schemes in comparison with other methods.
Research limitations/implications
This research is limited to shape optimization of two-dimensional electric structures, but however the work can be extended to 3D ones too.
Practical implications
The implementation of the proposed numerical approach for the shape optimization of a planar resistor is described.
Originality/value
The main value of the proposed approach is a powerful and robust numerical shape optimization algorithm that demonstrates outstanding suppleness of handling topological changes, fidelity of boundary representation and a high degree of automation.
Keywords: Shape Optimization, Level Set Method, Extended Finite Element Method, Genetic Algorithm, Electric Structures.
Category of the paper: Research paper
Authors information:
Vasile Topa1, Marius Purcar1, Calin Munteanu1, Laura Grindei1, Claudia Păcurar1,
Ovidiu Garvasuc1, Vasilis Chatziathanasiou2
1 Technical University of Cluj Napoca, Department of Electrical Engineering,
G. Baritiu street 26-28, 400020, Cluj Napoca, Romania
E-mail: [email protected]
2 Aristotle University of Thessaloniki, School of Engineering,
Department of Electrical and Computer Engineering, 54124, Thessaloniki, Greece
1. Introduction
The use of optimization tools in the early stage of the development process offers new potential in the process chain. The development process becomes faster and more efficient by using shape and topology optimization. This results in parts and devices that are lighter, better performing and more durable. Although significant progress has been made in the last decade in developing production quality CAD tools there is still a high need for more effective and fully automated optimization tools that are CAD integrated. This is driven by the ever increase in model complexity and the need to address and incorporate design and manufacturing requirements at the very early stages (Bendsøe and Sigmund, 2003; Haslinger and Makinen, 2003, etc.).
When shape or topology optimization using Finite Element Methods (FEM) is involved, a large number of mesh adaptations or even re-meshes are required during the optimization process. After each optimization step, a new geometry is obtained, and the finite elements (FE) mesh is possible to be distorted. When, during the optimization process excessive mesh distortion occurs, the accuracy of the numerical solution is essentially affected. As a consequence, a totally new mesh must be constructed. Therefore, developing robust algorithms for shape and topological optimization problems based on this approach still represent a real problem (Neittaanmaki et al., 1995; Haslinger and Neittaanmaki, 1997; Brandtstatter et al., 1998; Allaire, 2002, etc.).
Latest developments in mechanical engineering on modeling moving discontinuities without re-meshing provide challenging opportunities. This is due to the eXtended Finite Element Method (XFEM) that enables a local enrichment of the approximation space. From its conception (Belytschko and Black, 1999; Moes et al., 1999) the XFEM has been used to model a large number of applied mechanic problems such as growing cracks (Sukumar et al., 2000), solidification (Chessa et al., 2002), phase changes (Dolbow and Merle, 2002), material interfaces (Belytschko et al., 2003), two phase fluids (Chessa and Belytschko, 2003), dynamic crack propagation (Belytschko et al., 2003), dislocations and interfaces (Belytschko and Gracie, 2007).
On the other hand very powerful mathematical techniques and algorithms are available now in order to deal with moving boundary modeling problems, as for instance the Level Set Method (LSM), introduced by Osher and Sethian, 1988. Appling the LSM in place of performing geometrical operations, a convection equation is solved. In this way the new geometry, including the shape and topology changes are provided. As a consequence, the LSM becomes very popular in many scientific areas, such as computer graphics, image processing, computational fluid dynamics and structure optimization (Sussman et al., 1994; Sethian, 1999; Wang et al., 2003, etc.).
The combination of XFEM and LSM offers a very elegant and efficient way for solving shape and topology optimization problems, due to the fact that XFEM overcome the deficiency of using the classical FEM, as analysis tool. Literature with XFEM in level set framework, for shape and topological optimization problems is rather limited to structural optimization problems, such as the minimum mean compliance problem, the maximum natural frecquency problem and the minimum stress problem, but received more and more interest in the last years (Duysinx et al., 2006; Miegroet and Duysinx, 2007; Miegroet and Duysinx, 2009. etc.).
If the XFEM is recognized today in the scientific world as an important computational tool for complex applied mechanical problems, not the same can be said about the use of it for the study of the applied electromagnetic problems. So far, only in a few papers, XFEM was applied for solving electrical problems, as like e.g. electrostatic forces computation acting on micro-electro mechanical systems (Rochus et al., 2008; Rochus et al., 2009) and not at all for the optimal design of electromagnetic devices. As a consequence, we propose in this paper a shape optimization approach of complex electric resistors, based on XFEM and LSM.
The optimal design of the electrical resistance has a large impact on the performance of the electric circuits. Often, the designed pattern contains corners, which play an important role in the final value of the electric resistance and increase the power consumption and heat dissipation (Munteanu et al., 2001; Purcar et al., 2004). Such a case is an "L" resistor embedded in an insulator photoresistive substrate like, e.g. in Purcar et al., 2004. The shape optimization problem is solved by Purcar et al., using a numerical approach based on the combination of FEM with LSM. The FEM is used at each step of the optimization process in order to compute the current density distribution and the dissipated power on the given electric resistor. Then the LSM, represents the design structure through an embedded implicit function in order to "perturb the shape" and the finite element (FE) mesh is regenerated at each step of the optimization process.
In the present paper a similar resistor as in Purcar et al. 2004 is used for testing the proposed algorithm, based on XFEM and LSM. In the next section of the paper the XFEM is described in detail, being explained the corresponding terms due to the "enrichment", which it means that the displacement approximation is enriched (incorporated) by additional problem-specific functions. In section 3, the basis of the LSM, choice of the level function and computation of the velocity are described. In sections 4 and 5, details about the algorithm, a study case consisting in a two dimensional stationary electric field and the obtained results are presented. Conclusions and perspectives of the research are given in the last section.
The extended finite element method
Discontinuities, as for instance the moving material interfaces, play an important role in many types of optimisation problems. Scientific world has given more and more attention in the last years to the problems where it is necessary to model the motion of these kinds of discontinuities. Due to the fact that standard FEM are based on piecewise differentiable polynomial approximations, they are not well suited to problems where the solutions contain discontinuities, discontinuities in the gradient, singularities or boundary layers. Typically, FEM requires significant mesh refinement or meshes that conform to these features, to yield acceptable results.
In response to this deficiency of standard FEM, the extended finite element method is a novel approach tailored to simulate problems involving moving discontinuities. Initially the methodology was introduced for the analysis of crack propagation by the work of Belytschko and Black, 1999. The XFEM is a local partition of unity enriched finite element method (Melenk and Babuska, 1996), where by local understanding that only a region near the discontinuties such as cracks, holes, material interfaces are enriched. However this method has been quickly applied to several other types of problems, especially in the mechanical field.
In contrast with the FE meshes, where the mesh conforms to the interface, the XFEM uses a fixed mesh which does not need to conform to the interface. This is done by extending the standard FE approximation with extra basis functions that capture the behaviour of the solution near the interface. This is particularly useful for problems involving moving interfaces where the mesh would otherwise require regeneration at every time step, i.e. shape optimization problems. Consider the elliptic equation:
, (1)
in a domain Ω in two dimensions.
Embedded within Ω, there is an interface as in Figure 1. The coefficients α, κ and f may be discontinuous across and jump conditions are given on the interface. Inside the studied domain Ω, there are two sub-domains with different material properties. This type of problem arises in a broad spectrum of mathematical models and hence, a wide range of numerical methods have been devised to solve it. Often, the location of varies in time.
Figure 1
As a result, methods which are easily adapted to an arbitrary are important. In order to reduce the continuum problem described by equation (1) to a discrete system, the XFEM approach is used. The domain Ω could be meshed by an arbitrary FE mesh, but in this paper it is meshed with regular triangular elements, independent of the moving interfaces The XFEM approximation is:
, (2)
where:
Ni (x) are the classical shape functions associated to degree of freedom ui ;
Nj (x)· Ф(x, t) are the discontinuous shape functions constructed by multiplying a classical Nj (x) shape function with the enrichment function (presenting a switch value where the discontinuity lies);
ui and aj are the unenriched and the enriched degrees of freedom, respectively;
S is the set of nodes whose support is crossed by the moving interface.
Figure 2 shows how the elements placed in the vicinity of the mobile interface are enriched for a typical 2D discontinuity.
Figure 2
To include the interface's effect, enrichment functions are added to the standard finite element approximation for each element cut by the interface.
Three types of elements are defined:
fully enriched elements, all the nodes of the triangle elements are enriched with the enrichment function;
partially enriched elements, some nodes which belongs to the triangle elements are enriched with the enrichment function;
unenriched elements, no nodes are enriched. These are the classical FEM elements.
In other words, only the elements near the material discontinuity support extended shape functions, whereas the other elements remain unchanged.
With this assumption the gradient discontinuity at the interface Γ (t) is computed as:
(3)
The choice of enrichment function is based on the behaviour of the solution near the interface. These functions are chosen a priori by the knowledge of the physical problems at hand. In Table 1 are given the most typical enrichment functions.
Table 1
Type of problem
Displacement
Strain
Enrichment
Bi-material
Continuous
Discontinuous
Ramp
Crack
Discontinuous
-
Heaviside
Crack tip
Discontinuous
order
High gradient
order
In our applications the Ramp function is used:
(4)
This enrichment function yields only continuous solutions. The advantage is that it automatically satisfies the continuity condition [u] = 0 and does not require the use of Lagrange multipliers (as for example the Heaviside function).
The graphics of the discontinuous shape functions obtained by multiplying a classical Nj (x) shape function with an enrichment function are given in Figure 3. Case (a) match to a standard FE shape functions, while case (b) to a Ramp enrichment function. The discontinuity is visibly observed in Figure 3(b).
Figure 3
The change of the classical FE field approximation does not introduce a new form of the discretized finite element matrix equation, but leads to an enlarged problem to be solved:
. (5)
Kuu is the standard FEM element matrix, given by:
, (6)
with: B the standard discretized gradient operator and C the material law. The Kua , Kau and Kaa matrices are the new terms that arise from the addition of the enriched degrees of freedom and are computed in accordance with:
, (7)
where: is the enriched discretized gradient operator.
Below there are given the computed matrix coefficients if a Ramp function is used as an enrichment function and the material matrix C is supposed to be constant:
(a) The matrix coefficients (λ, β = i, j, k)
(8)
(b) The matrix coefficients (λ, β = i, j, k)
(9)
Evaluating the domain integral terms requires a numerical quadrature method. Elements away from the interface are evaluated using standard Gaussian quadrature.
Elements that are cut by the moving interface must be treated differently due to discontinuities in the coefficients and the enrichment function. The interface is first interpolated as a line segment (ab) and the element is divided into triangles that conform to the interface (Figure 4). In this case the integration term becomes:
(10)
The sub-elements are for integration only and do not introduce any extra degrees of freedom.
Figure 4
3. The level set method
The explicit representation of the interface that is used in the classical FEM forbids deep boundary or topological changes, such as creation of holes. This limitation is the main reason of the low performance generally associated to the shape optimization problem. In opposite, the LSM developed by Osher and Sethian, 1988 which consist of representing the interface with an implicit function, overcomes this kind of deep changes.
The LSM is a numerical technique first developed for tracking moving interfaces. It is based on the idea of representing implicitly the interfaces as a level set curve of a higher dimension function.
The moving interface , considered between the two sub-domains Ω1 and Ω2 with different conductivity (see Figure 1), is conventionally represented by a zero Level Function, and the sign convention, as follows [9]:
(11)
This function is approximated on a fixed mesh by a discrete function which is usually the signed distance function:
. (12)
The sign is positive (negative) if x is inside (outside) the interface as in Figure 5. The evolution of the interface is then embedded in the evolution equation for, which is given by [12]:
, (13)
Where: is the speed function defined on the interface in the outward normal direction to the interface.
Figure 5
Appling the XFEM framework, the Level Set is defined on the structured mesh and at each finite element node is associated a geometrical degree of freedom representing its Level Set function value. The Level Function is discretized on the whole design domain Ω, with the standard FE and in all cases the same mesh and shape functions are used as for the dependent variable [10]:
. (14)
The system of equations derived from equation (13) can be written as [12]:
, (15)
where: is the time variation vector of Φ in every node, {Φ} the vector of the unknown Φ at each time step, M the global mass matrix and C the global convection matrix.
Solving equation (13) for Φ = 0 on the boundary and in the vicinity of the boundary in Ω1 and Ω2 respectively, the new shape of the interface is found through interpolation of the LSM solution, on the FEM mesh (Equation 14).
4. Results
The test problem consists of a rather academic case, an "L" resistor embedded in insulator photoresist substrate, like in Purcar et al. (2004). All the dimensions of the electric resistor are given in meters as in Figure 6.
Figure 6
Inside the studied domain Ω, there are two sub-domains with different conductivity. The substrate conductivity (sub-domain Ω1) is = 1.0e-6 and the resistor conductivity (sub-domain Ω2) is = 5.72e+6 . The resistor terminals are marked with thick lines and are considered having a rigid position. The potentials of these terminals are u1 = 10 (V) and u2 = 0 (V), respectively. The interfaces between the two domains of different conductivities are set as moving boundaries and they will account for the resistor's shape optimization. Denote the union of the interfaces between Ω1 and Ω2 by. The optimization objective is to minimize the value of the resistor, equivalent to find the shape of the interface which assures this value.
The resistance is computed using a stationary electric field model, without charge distributions inside Ω2 and governed by a Laplace's equation:
, (16)
where u represents the electric potential distribution, and σ the electric conductivity.
The applied methodology should be synthesised as below and consists in the following steps:
The computational domain of the resistor is defined using a cross-section through the model in figure 6 right, by a plane perpendicular to the middle of the electrodes (1, 2). In order to accurately compute the boundary shape change rate equation (14) and (15) and the shape of the resistor at every optimization step, a uniformly spaced grid allover the computational domain Ω and the photoresit domains Ω2 and Ω3 is generated. The mesh of the domains Ω1,2,3 is merged in a big one of c.a. 2800 points and 5500 elements by keeping track of the boundary conditions and material properties of each of the underlying domains and used unchanged for all of the optimization steps.
Using the finite element method (FEM), the governing Laplace equation is solved in order to compute the current density distribution over the resistor surface and the boundary shape change rate (14) and (15). As it has been already mentioned, in most of the LSM applications has a meaning only on the boundary where is calculated and therefore the velocity extension on the computational domain is compulsory. Equations (14) and (15) give velocity at the nodes of the triangular elements with conductivity σ1 while at the other nodes is 0.
The velocity extension is not performed on the whole computational domain but only on a narrow band of triangles around the moving front [6] (see figure 9 left). After the initialization (and reinitialisation) of the level set function Φ with the signed distance function, figure 9 right, level set equation () is iterated over the narrow band (see figure 9 left). The usual number of triangles within the narrow band is c.a. 500 and considerably reduces the computational effort of level set equation ().
The micro scale time step ∆tLSM is accordingly computed such that it stays within the narrow band. Level Set equation is integrated implicitly within roughly 5 to 10 ≈ ΔtF/ tLSM micro scale time steps. The new zero level set solution gives the new position and shape of the resistor boundary for the next macro scale time step. As vmax is subject to changes at the next macro scale time step, the micro scale time step ∆tLSM and the narrow band are accordingly updated.
Figure 9: Velocity extension on the narrow band (left); sign distance to the shape to be optimized (right).
The maximum dissipated power for the initial geometry where the corners of the resistor play a negative role in the power dissipation, can be seen in figure 10 left. By imposing as dissipated power pimp= 2.77e+8 (W/m3) in equations (14) and (15), it can be clearly seen in figure 10 right the power dissipation tends to become uniform overall surface of the resistor while the optimization proceeds with a perturbation of the boundary under the speed law (14). As already mentioned this speed law described the balance between the maximal initial dissipated power below which the material can be eliminated and above which material should be added
(a)
(b)
(c)
Figure 10 Boundary evolution profiles during the optimization process. pmean -pimp on the initial geometry (left); after 10 optimization steps (right).
Beside the objective function (16) there are also two constraints imposed. Since the resistor's topology and the computation domain are reconstructed during the optimization process, the electrodes (1 and 2) are constrained to not move or change in size.
The optimisation is carried out over 20 optimisation steps. The maximum and minimum dissipated power variation over the optimisation steps is presented in figure 11.
Figure 11: Variation of pmax and pmin with the optimization step
At the end of the optimisation process, the geometry is changed in such way (see figure 12) that the dissipated power is almost constant in the entire resistor if refers the initial geometry.
Figure 12, indicates how the evolution of the boundary profile attenuates in the neighborhood of the optimal profile and level set march less than one layer of triangles. In order to stop the algorithm a condition as objective function is far more too stringent and therefore equation (16) is used.
The used structured finite element mesh in the XFEM approach, non-conform to the interface is given in the figure below, as well as the potential distribution for the initial structure.
Conclusions and perspectives
In this paper, a new method based on XFEM and LSM methods for shape optimization of 2D electric structures is presented. No significant problems have been encountered in the optimization procedure and quite good results were obtained in comparison with other results found in the literature, for the similar study case [11], [12]. The XFEM method has proven to be very useful: no remeshing process is needed in this application during the optimization process.
The ongoing work will extend the results already in hand: first, by computing the sensitivity analysis in order to use deterministic optimization methods, more faster as the stochastic ones; the next step will certainly focus on the development of a "real" Level Set framework in order to allow much more performance in the sense of geometrical modification.
Obviously, the coupling of the XFEM with the Level Set has proven to be a really promising method for the shape as well as for the topology optimization.
This new approach removes the conventional use of discrete elements and provides efficient, stable, more accurate and much faster computation schemes, eliminating the remeshing at each optimization step.
The XFEM and LSM exhibits really promising characteristics, as allows deep topological changes and a very flexible modeling of the geometry. However numerical applications pointed out some specific difficulties and problems to be handled in the implementation.
