SHAPE OPTIMIZATION USING FINITE EL EMENT ANALYSIS IN EDDY CURRENT TESTING AND ELECTRO - THERMAL COUPLED PROBLEM S By Victor Uthayakumar Karthik A DISSERTATION Submitted to Michigan State University i n partial fulfillment of the requirements for the degree of Electrical Engineering Doctor of Philosophy 2015 ABSTRACT SHAPE OPTIMIZATION USING FINITE ELEMENT ANALYSIS IN EDDY CURRENT TESTING AND ELECTRO - THERMAL COUPLED PROBLEM S By Victor Uthayakumar Karthik The inversion of electro - heat problems is important in areas such as electrical machine design, the metallurgical processes of mixing, and hyperthermia treatment in oncology. One of the important computations involves synthesizing the electromagnetic arrangement of coils to accomplish a desired heat distribution to achieve mixing and reduce machine heat , or to burn cancerous tissue. Two finite element problems need to be solved, first for the magnetic fields and the joule heat that the associated eddy currents generate , and then, based on these heat sources, the second finite element problem for heat distribution. This two part problem needs to be iterated on to obtain the desired thermal distribution by optimiz ing the shape of the current source . This represents a heavy computational l oad associated with long wait - times before results are ready. The graphics processing unit (GPU) has recently been demonstrated to enhance the efficiency of finite element field computations and to cut down solution times. In our work , given the heavy comp utational load from the two - part problem and the attendant optimization, we use the GPU to perform the electro - heat optimization by the genetic algorithm launched on several parallel threads to ach ieve computational efficiencies. To avoid the complexitie s of a two - physics problem, we first develop the shaping algorithms on a single physics problem for nondestructive evaluation. This shape optimizing concept is developed for defect characterization in a steel plate. When a steel plate in a ground vehicle i s found to be defective, it is usually taken out of service for repair without determining if the defect warrants withdrawal . W e investigate and estab lish a procedure s o that a decision to withdraw a vehicle may be justifiable and ensure s the safety of the vehicle and its passengers. To test o ur algorithms and elements and optimization with true device design. Flip - teaching was introduced to tackle the challenges of time. The traditional order of a) delivering theory b) programming ancillary tools (mesh generators, solvers) is flipped to do real design . After the algorithms we re understood from use in NDE and the f lip - t eaching experience , we successfully develop ed them for the two - physi cs system with reduced computational time with the speedup (CPU Time to GPU time ratio) of 28 and increased accuracy established through the problem of shaping a current carrying conductor so as to yield a constant temperature along a line . Finally we appl ied the electro - thermal software to hyp erthermia treatment planning by a numerical model of a human thigh with a tumor treated by current carrying conductors to be shaped to produce the desired temperature at the tumor. The bio - heat equation s under steady state condition are solved and the heat removal due to blood perfusion is also taken into account to determine the shape yielding the desired temperature profile . An efficient methodology for multi - physics systems has been developed with applications in fl ip - teaching, NDE and hyperthermia. iv T his dissertation is lovingly dedicated to my mother, Sakunthaladevi Uthayakumar and my sustained me throughout my life . v ACKNOWLEDGEMENT S It has been three and a half years since joining the PhD program at Michigan State University. As this journey comes to an end, I would like to express my gratitude to several people without whom successful completion of my degree woul d not have been possible. First and foremost, I thank my God for giving me the good health and guidance that were necessary to complete this program. I would like to express my sincere gratitude to my advisor Prof. S. Ratnajeevan. H. Hoole for the conti nuous intellectual and financial support of my PhD studies and related research, for his patience, motivation and immense knowledge. His guidance and prompt response to my questions helped me in all the time of the research and writing of this thesis. I co uld not have imagined having a better advisor and mentor for my PhD program. Besides my advisor, I would like to thank the rest of my PhD committee members: Prof. Lalita Udpa, Prof William. F. Punch, Prof. Prem J. Chahal and Dr. Paramasothy Jayakumar, no t only for their insightful comments and encouragement, but also for the hard questions which were the incentives for me to widen my research and look at it from various perspectives. Development Center (TARDEC) for funding our research under Contract N umber W911NF - 11 - D - 0001 and W56HZV - 07 - 2 - 001 . I am thankful to the MSU Graduate School and College of Engineering vi for awarding summer graduate excellence fellowship awards twice and offeri ng graduate teaching assistantship for 3 semesters. This support was instrumental in facilitating the completing my degree. I also would like to thank the members in my research group, S. Sivasuthan, T. Mathialakan and Mohammad R. Al Rawashdeh for sharin g their research knowledge and the friendships they provided me. Last but not least, I would like to thank my family my mother, my wife, my sister and my brothers and my friends for sustaining me emotionally and my Church for supporting me spiritual ly throughout writing this thesis and my life in general. vii TABLE OF CONTENTS LIST OF TABLES ................................ ................................ ................................ ................... ix LIST OF FIGURES ................................ ................................ ................................ .................. x 1 INTRODUCTION ................................ ................................ ................................ ................. 1 1.1 Motivation ................................ ................................ ................................ ........................ 1 1.1.1 The Electro - Thermal Coupled Problem ................................ ................................ .... 1 1.1.2 Reconstructing and Classifying Damage in a Steel Plate ................................ ......... 2 1.2 Inverse Problems for Design ................................ ................................ ............................ 3 1.3 Shape Optimization of Two - Physics Systems: Gradient and Zeroth Order Methods ...... 6 1.4 Method of Optimization The Genetic Algorithm ................................ ........................ 10 1.4.1 Decision to use the GA ................................ ................................ ........................... 10 1.4.2 The Genetic Algorithm ................................ ................................ ........................... 12 1.4.2.1 Introduction ................................ ................................ ................................ ..... 12 1.4.2.2 Basic Concepts Biological Background ................................ ....................... 12 1.4.2.3 Working Principle ................................ ................................ ............................ 13 1.4.2.4 Representing a Solution ................................ ................................ ................... 14 1.4.2.5 Assigning a fitness score ................................ ................................ ................. 15 1.4.2.6 Selection ................................ ................................ ................................ .......... 17 1.4.2.7 Crossover ................................ ................................ ................................ ......... 18 1.4.2.8 Mutation ................................ ................................ ................................ ........... 19 1.4.2.9 Real - coded GAs ................................ ................................ ............................... 20 1.5 Computer Aided Design ................................ ................................ ................................ . 22 1.5.1 In pursuit of Accuracy ................................ ................................ ............................ 22 1.5.2 Methods of Appr oximate Solution ................................ ................................ .......... 25 1.5.3 The Finite Element Method ................................ ................................ .................... 26 1.5.4 Two Dimensional, Linear, Triangular Finite Elements ................................ .......... 33 2 EDDY CURRENT BASED DEFECT DETECTION AND CHARACTERIZATION . 40 2.1 Single Physics Problem ................................ ................................ ................................ .. 40 2.2 Introduction to Eddy Current Testing ................................ ................................ ............ 41 2.2.1 Eddy Current ................................ ................................ ................................ ........... 41 2.2.2 Basis of eddy current testing ................................ ................................ ................... 41 2.3 Detection of Defect ................................ ................................ ................................ ........ 42 2.3.1 Field Changes ................................ ................................ ................................ .......... 42 2.3.2 Main Equations ................................ ................................ ................................ ....... 45 2.3.2.1 Maxwell Equations ................................ ................................ .......................... 45 2.3.2.2 Eddy Current Equation ................................ ................................ .................... 47 2.3.2.3 Finite Element Computation for Eddy Current Problem ................................ . 49 2.3.3 Problem Statement ................................ ................................ ................................ .. 52 2.3.4 Experimental Work in Defect Det ection ................................ ................................ . 56 2.4 Defect Characterization ................................ ................................ ................................ .. 58 viii 2.4.1 Design Parameters ................................ ................................ ................................ .. 58 2.4.2 Numerical Model ................................ ................................ ................................ .... 61 2.4.3 Constraints Handling ................................ ................................ .............................. 62 2.5 Results and Discussion ................................ ................................ ................................ ... 62 2.6 Experimental Work in Defect Characterization ................................ ............................. 69 2.6.1 Future Work ................................ ................................ ................................ ............ 74 2.7 Conclusion ................................ ................................ ................................ ...................... 75 3 FLIP - TEACHING TO TEACH FINITE ELEMENT OPTIMIZATION ..................... 76 3.1 Introduction - Flip Teaching ................................ ................................ .......................... 76 3.2 Teaching Engineering Optimization: Real Design ................................ ......................... 78 3.3 A Useful Tin Problem ................................ ................................ ................................ .... 81 3.4 Design of complex devices: The challenge ................................ ................................ .... 83 3.5 The Assignments ................................ ................................ ................................ ............ 85 3.6 Conclusions ................................ ................................ ................................ .................... 86 4 EL ECTRO - THERMAL PROBLEM ................................ ................................ ................ 87 4.1 Introduction ................................ ................................ ................................ .................... 87 4.2 Context and Novelty ................................ ................................ ................................ ....... 91 4.3 Finite Element Computation for the Electro - Thermal Problem ................................ ..... 92 4.4 GPU Computation for Genetic Algorithms for Electro - heat Problems ......................... 94 4.5 Test Problem: Shaping an Electro - heated Conductor ................................ .................... 97 4.6 Results and Discussion Thereof ................................ ................................ ................... 100 4.7 Conclusion s ................................ ................................ ................................ .................. 103 5 HYPERTHERMIA TREATMENT PLANNING ................................ ........................... 105 5.1 Introduction ................................ ................................ ................................ .................. 105 5.2 Problem Statement ................................ ................................ ................................ ....... 107 5.3 Main Equations ................................ ................................ ................................ ............ 109 5.3.1 The Finite Element Equation ................................ ................................ ................ 109 5.3.2 A note on electrical and thermal conductivity changes ................................ ........ 114 5.3.2.1 Electrical conductivity ................................ ................................ ................... 114 5.3. 2.2 Thermal conductivity ................................ ................................ ..................... 115 5.4 The Algorithm for the Inverse Method ................................ ................................ ........ 115 5.5 Simulation Results ................................ ................................ ................................ ........ 117 5.6 Conclusions ................................ ................................ ................................ .................. 121 6 SUMMARY AND CONCLUSIONS ................................ ................................ ................ 123 APPENDICES ................................ ................................ ................................ ....................... 126 Appendix A: First - Order Interpolations Using Triangular Coordinates ................................ . 127 Appendix B: Publications Raised from This Thesis ................................ ............................... 128 BIBLIOGRAPHY ................................ ................................ ................................ ................. 131 ix LIST OF TABLES Table 1 - 1 The Local Matrices of the Six Elements of Figure 1 - 16 ................................ ........... 36 Table 1 - 2 The Six Local Matrices Added into the 4x4 Matrix ................................ ................. 38 Table 2 - 1 Calculating the efficiency of reconstruction ................................ ............................. 65 Table 2 - 2 Real and binary solutions time need to compute ................................ ...................... 67 Table 2 - 3 Real and binary best fitness score achieved ................................ .............................. 67 Table 5 - 1 Geometrical dimensions of the model ................................ ................................ .... 109 Table 5 - 2 Physical properties of the numerical model ................................ ............................ 117 Tabl e 5 - 3 Physical parameters of blood ................................ ................................ .................. 11 8 Table 5 - 4 Reading at the measuring line ................................ ................................ ................. 121 x LIST OF FIGURES Figure 1 - 1 Typical Forward Problems: Analysis ................................ ................................ ........ 3 Figure 1 - 2 Pole - Faces to be shaped ................................ ................................ ............................. 4 Figure 1 - 3 Minimal Problem of Figure 1 - 2 u sing Symmetry ................................ ..................... 4 Figure 1 - 4 Jagged Pole Face of Right Half of Recording Head ................................ .................. 5 Figure 1 - 5 Genetic Algorithm Speed Compared with S imulated Annealing ............................ 11 Figure 1 - 6 A typical chromosomal representation ................................ ................................ .... 14 Figure 1 - 7 Flowchart of the working principle of a GA ................................ ........................... 16 Figure 1 - 8 A random population of six cans ................................ ................................ ............. 16 Figure 1 - 9 Tournaments played between six population members ................................ ........... 18 Figure 1 - 10 Changing design parameters in GA ................................ ................................ ....... 18 Figure 1 - 11 An illustration of the single - point crossover operator ................................ ........... 19 Figure 1 - 12 An illustration of the mutation operation ................................ .............................. 20 Figure 1 - 13 Probability distributions used in SBX crossover with different index n ............... 21 Figure 1 - 14 The Solution Process ................................ ................................ ............................. 23 Figure 1 - 15 A Triangle with Three Nodes ................................ ................................ ................ 28 Figure 1 - 16 Heat Flow Problem ................................ ................................ ................................ 33 Figure 1 - 17 Triangle with internal generic numbers ................................ ................................ . 37 Figure 2 - 1 Eddy current flow on a sur face with (a) and without (b) crack ............................... 41 Figure 2 - 2 Penetration of Eddy Currents in the Presence of Defects ................................ ........ 43 Figure 2 - 3 Defect M odel: The model has defects in different angle, depth and length ............ 44 Figure 2 - 4 Numerical Model ................................ ................................ ................................ ..... 53 xi Figure 2 - 5 Equipotential line s for the defects at various positions and orientations ............... 54 Figure 2 - 6 with defect location ................................ ................................ ............... 55 Figure 2 - 7 Rectangular defect on steel plate ................................ ................................ ............. 56 Figure 2 - 8 Lab setup ................................ ................................ ................................ ................. 57 Figure 2 - 9 Close - up - look of the defect ................................ ................................ ..................... 57 Figure 2 - 10 Voltage induced on pickup coil ................................ ................................ ............. 58 Figure 2 - 11 Defect Model ................................ ................................ ................................ ......... 59 Figure 2 - 12 Design Cycle fo r the computation process ................................ ............................ 60 Figure 2 - 13 Numerical Model for defect characterization ................................ ........................ 61 Figure 2 - 14 Optimum shape of the Reconst ructed Defect ................................ ........................ 63 Figure 2 - 15 Calculating efficiency of reconstruction ................................ ................................ 64 Figure 2 - 16 Solutions in Equipotential lines for the Nu merical Model ................................ .... 64 Figure 2 - 17 Time required for real and binary solutions ................................ .......................... 68 Figure 2 - 18 True and Reconstructed shape of defects ................................ .............................. 69 Figure 2 - 19 Planar coil and GMR sensor ................................ ................................ .................. 71 Figure 2 - 20 Induced eddy current in the plate ................................ ................................ .......... 71 Figure 2 - 21 Experimental set - up ................................ ................................ ............................... 72 Figure 2 - 22 Scanned Images for the aluminum plate with rivet ................................ ............... 72 Figure 2 - 23 Readings for the steel plate with defect ................................ ................................ . 73 Figure 2 - 24 2D model to match the experimental setup ................................ ........................... 74 Figure 3 - 1 Pole - faced to be shaped ................................ ................................ ........................... 78 Figure 3 - 2 Minimal problem using symmetry ................................ ................................ .......... 79 Figure 3 - 3 Fringing in initial geometry ................................ ................................ ..................... 79 xii Figure 3 - 4 The shaped pole face ................................ ................................ ............................... 80 Figure 3 - 5 Stretched Spring ................................ ................................ ................................ ...... 80 Figure 3 - 6 Obje ct function for spring (on the plane) ................................ ................. 81 Figure 3 - 7 Object function ................................ ................................ ................................ ........ 83 Figure 3 - 8 A simplified 5 - parameter mesh genera tor ................................ ............................... 84 Figure 3 - 9 The final optimum shape ................................ ................................ ......................... 85 Figure 4 - 1 Finite Element Analysis and Optimization of Coupled Magneto - Thermal Prob lems. ................................ ................................ ................................ ................................ ................... 88 Figure 4 - 2 The Design Cycle for the Coupled Magneto - Thermal Problem. ............................ 90 Figure 4 - 3 The Parallelized Process of the GPU ................................ ................................ ....... 96 Figure 4 - 4 Numerical Model for Coupled Electro - heat Problem ................................ .............. 97 Figure 4 - 5 The Parameterized Geometry ................................ ................................ .................. 99 Figure 4 - 6 Equi - Temperature Lines at different iterlations ................................ .................... 101 Figure 4 - 7 Optimal Shape by the Genetic Algorithm ................................ ............................. 102 Figure 4 - 8 Temperature Distribution: Desired, Initial and Optimized ................................ .... 102 Figure 4 - 9 Speedup: GA Optimization GPU Time: CPU Time ................................ .............. 103 Figure 5 - 1 Cross section of the human thigh with a tumor ................................ ..................... 107 Figure 5 - 2 Numerical model of the problem ................................ ................................ ........... 108 Figure 5 - 3 Design parameters ................................ ................................ ................................ . 116 Figure 5 - 4 Mesh for the numerical model ................................ ................................ ............... 118 Figure 5 - 5 Equi - Temperature lines for the initial shape ................................ ......................... 120 Figure 5 - 6 Optimum shape of the conductors to treat tumor ................................ .................. 120 Figure 5 - 7 Initial, optimum and desired tem perature ................................ .............................. 122 xiii Figure A - 1 Triangular Coordinates ................................ ................................ ......................... 127 1 1 IN T RODUCTION 1.1 Motivation 1.1.1 The Electro - Thermal Coupled Problem In electromagnetic field problems, optimizati on methods have usually been successfully developed and efficiently applied on a single processor . Most of these methods only deal with the single - physics problem. However, real problems are more complex and often coupled. Coupled field problems often are where one field system with the shapes to be optimized influences another field system in which the objects of design are defined. A good example is electro - heating [1] . The electrical system provid ing the joule heating has to be shaped so as to produce a particular heat distribution in the thermal system . An example application in i ndustry is metal forming , where molten metal is heated through heavy currents to produce the forces that make the molten metal subject to the extrusion or turning forces we want [2], [3] . A second example is hyperthermia treatment for oncology where exterior electrodes attempt to burn interior cancerous tissue [4] . The inversion of these two electro - heat problems involves synthesizing the electromagnetic arrangement of coils so as to accomplish a desired heat distribution. To this end two finite element problems need to be so lved, first for the magnetic fields and the joule heat that the associated eddy currents generate and then, based on these heat sources, the second problem for heat distribution. This two part problem needs to be iterated on to obtain the desired thermal d istribution by optimization. This problem represents a heavy computational load associated with long wait - times before results are ready. The graphics processing unit (GPU) has recently been demonstrated to enhance the efficiency of finite element field co mputations and cut down solution times [5] . In this work , given the 2 heavy computational load from the two - part problem and optimization, we use the GPU to perform the electro - heat optimization by the naturally parallel Genetic Algorithm (GA) [6] to achieve computational efficiencies far b etter than those reported for a finite element problem on a single processor . The methodology established has direct applications in electrical machine design, metallurgical mixing processes, and the treatment of cancer by hyperthermia. In these three practical applications, we need to compute the coil arrangement that would accomplish a desired heat distribution to achieve mixing, reduce machine heat or to burn cancerous tissue. Particularly the above - mentioned application will alleviate human sufferi ng through use in hyperthermia treatment planning . 1.1.2 Reconstruct ing and Classifying Damage in a Steel Plate Our main object is to optimize shape in two - physics problem s . To avoid the complexities and verify the design optimization algorithm , our algorithm i s experimented first on a single - physics problem to identify and characterize the internal defects in a steel plate using Non - Destructive Evaluation (NDE) Methods. When a steel plate in an army ground vehicle is found to be defective, it is usually taken o ut of service for repairing without determining if the defect warrants withdrawal or not, whereas the defect might be minor and withdrawal wasteful. On the other hand the defect might be invisible, yet warranting replacement. The methodology at present exa mines the response of the hull under test to an excitatory signal from an eddy current probe [7] . Kn owing the response when there is no defect, if the response is different because of the defect, the test object is presently flagged as defective and the vehicle is sent for repairs without assessing if the defect is serious enough for removal of the vehic le from service. In our work, w e characterize defects by parametrically describing them and then 3 optimizing the parameters to match computed fie lds to measurements. Therefore knowing the shape of the defect , any decision to withdraw may be tho ught - out, jus tifiable and ensuring of the safety of the vehicle and its passengers. We also backup our computational work with laboratory experiments. 1.2 Inverse Problems for Design The direct problem with which the finite element method started has a device governed by a particular differential equation, the Poisson equation in our case (1 - 1) as was solved by the finite element method by Zienkiewicz and C heung in their classic paper [8] . Once we have which may be electric potential, pressure , magnetic vector potential etc. depending on the system we may compute performanc e descriptions like inductance, force, etc. as shown in Figure 1 - 1 . That is, from the system description, we compute performance. This is analysis . Figure 1 - 1 Typica l Forward Problems : Analysis The inverse problem the more practically realistic problem, which is synthesis goes from the right hand side of Figure 1 - 1 to its left. That is, wanting a performance, computing t he system description from it. Thus the computational design assignment may be this: compute the size and other descriptions of a motor that can produce so much torque. In industry this 4 was done by the cycle of design - make - test - redesign. This required an e xpert to redesign and took long. In time by the 1970s instead of making and testing, it was analyzed, solving the direct problem by the finite element method. It was left to engineers dealing with stress analysis and fluids to couple optimization with the finite element method [9] [11] , and the second half of the 1970s and 1980s would be the time for true synthesis solving for geometric shape and material values from design criteria. The earliest persons to automate this cycle in magnetics were Marrocco and Pironneau in 1978 [10] . Figure 1 - 2 Pole - Faces to be shaped Figure 1 - 3 Minimal Problem of Figure 1 - 2 us ing Symmetry 5 They attempted to optimize the shape of the magnetic pole of a recording head so that the fringing effect at the edges of a pole could be countered so as to realize the object of constant flux density B in a recording head. Figure 1 - 2 gives a similar problem with a repeating alternating pole system from electric machinery where a constant flux density is required on top of each pole to facilitate alternating waveform generation. The minimal boundary va lue problem for analysis is shown in Figure 1 - 3 . Marrocco and Pironneau [10] located their work in the latter Université Pierre - et - Marie - Curie (also known as UPMC and Paris VI), optimizing structural and fluid systems. That is, their work may be seen as parallel to the 1976 work of Arora and Hang [11] who established finite element optimization in a journal. They approached this problem by defining an object function F consisting of the square of the difference be tween the computed and desired flux densities. Thus the problem is one of optimizing i.e., minimizing F which is a function of parameters defining the geometry and which are computed in device synthesis so as to minimize F: (1 - 2) Figure 1 - 4 Jagged Pole Face of Right Half of Recording Head 6 Their results are shown in Figure 1 - 4 . The nonsmooth jagged conto ur in Figure 1 - 4 - b that they realized is practically not a manufactural shape. This they addressed by manually smoothening the pole face as in Figure 1 - 4 - c . However, Marrocco and Pironneau, aware of the problem of jagged contours, have further addressed it to permit smoothening by allowing nodes to move only along prescribed paths. Although the comprehensive 198 4 book by Pironneau o n optimization [12] deals with the theory of imposing constraints and applies them to many systems, the results presented from magnetics are the same as presented by Marrocco and Pironneau much earlier without cons traints. Pironneau, a widely experienced pioneer scientist in optimization and finite element analysis, particularly in fluids , was not focused on magnetics. He with Marrocco solved this problem to demonstrate the broad applicability of their methods of fi nite element optimization and then moved on. 1.3 Shape Optimization of Two - Physics Systems: Gradient and Zeroth Order Methods In the optimization of a single physics problem as in magnetostatics [9], [12] [16] or eddy current analysis [17 ] where the design variables are defined in the vector {p} , we construct one finite element mesh, solve for the magnetic vector potential , and then change {p}. The method by which we change {p} depends on the method of optimization we employ [9] . In coupled problems like the electro - heat problem under discussion, two different meshes are often required [18] . For example, at a copper - air boundar y in magnetics, both regions are nonmagnetic and therefore have the same permeability, the permeability of free space . However, for the thermal problem, they need to be model ed as two different regions because air has little thermal conductivity whereas copper is highly conductive. 7 Moreover, the optimization process too imposes huge difficulties depending on the me thod employed. In the simpler zeroth order methods, only the value of , given {p}, needs to be computed. This takes simply a finite element solution for a mesh constructed for the present value of {p}. In the more powerful first order methods, however, d efining symbols , the gradient vecotor (1 - 3) needs to be computed [9] . This may be by finite differences that is, to get the derivative of F with respect to , we need to evaluate F( ), then adjust by a very small amount , redo the finite element solu tion (which means a new mesh have to be generated for the changed geometry), and then evaluate F, which would give us (1 - 4) thereby leading to the derivate by finite difference, . This is computationally expensive because of having to be done for all m parameters , which means the first finite element solution at {p} and evaluation of F({p}) has to be followed by m finite element solutions with only a particular adjusted by and then the evaluation of m values of . That is, m+1 finite element meshes a nd solutions are required at great cost. Furthermore, even after all that work, it is known that the derivative by this finite difference process has poor accuracy [19] . The less approximate way to obtain accurate derivatives is from the derivative information inherent to the finite element solution through the finite element trial function; that is, although t he solution for the magnetic vector potential is explicitly in terms of the values of at the finite element nodes, we are really solving for as given by the trial function 8 which is expressed through the nodal values of [20] . That is, in solving the finite element Dirichlet matrix equation (1 - 5) although we explicitly solve for the nodal values {A}, it is really for the trial function A (x,y) that we are solving. Both the Dirichlet matrix [P] and right hand side {Q}in (1 - 5) are expressed as known functions of the vertex coordinates of the finite elements of the mesh, permeability, and current density. After solving for {A} in (1 - 5) , differentiating the equation with respect to we obtain [12 14, 17, 19] (1 - 6) where {A} has already been solved for, and and are computable. There is further computational efficiency to be reaped [19] by using the Cholesky factorization method of splitting [P] into its lower triangular and upper triangular Cholesky factors [L] and [U] (1 - 7) which for symmetric [P] as in both magnetic and thermal field problems by finite elements , reduces to (1 - 8) so that solving (1 - 4), in its new form (1 - 9) reduces to first solving (1 - 10) for {z} where (1 - 11) 9 and then, having found {z}, solving (1 - 11) for {A}. The computational efficiency lies in the solution is in finding [L] by solving (1 - 8). Thereafter the forward elimination in solving (1 - 10) and back substitution in solving (1 - 11) are trivial [15] . That is, once [L] and [U] = are in hand, solving (1 - 6) with t he same coefficient matrix [P] as (1 - 5) for the m gradients is trivial since only forward elimination and back substitution are required. Be that as it may, while solving (1 - 6) is trivial, forming (1 - 6) is not because computing and to form (1 - 6) is, in terms of programming, an arduous task that is not easily amenable to building up as general purpose software. As a particular changes, some vertices of a few triangles will move [15], [20] and the analyst needs to keep track of whether one, two or all three of the vertices of a triangle move by . Very complex coding is required that is problem - specific rather than general purpose. In a coupled pro blem, computing the derivatives of the finite element equations for temperature with respect to parameters in the magnetic problem is prohibitively complex although there are ways around it [21] . Programming the problem specific computations of and is ill - advised because of the complexities. Therefore zeroth order optimization methods, which are more slowly convergent than gradient methods, are the best route to go for optimizing coupled electro - heat problem s. Simkin and Trowbridge [22] aver that simulated annealing and the evolution strategy (a variant of the genetic algorithm [23] ) take many more function evolutions. Although they are computationally intens ive , they are far easier to implement, especially as general purpose 10 software. Indeed commercial codes that need to be general purpose use zeroth order methods not necessarily because they are superior to gradient methods, but because once a finite element analysis program is developed by a company, giving it optimization capabilities only takes coupling it with an optimization package to which object function evaluations can be fed whereas feeding both the object function and its gradient ( as would be re quired when gradient methods are in use ) would take extensive code development. In the context of single physics problems, Haupt [ 24] advises that the genetic algorithm is best for many discrete parameters and the gradient methods for where there are but a few continuous parameters. We rationalize this position on the grounds that gradient computation though difficult is more manag eable when there are fewer parameters to optimize with respect to. Indeed, we have gone up to 30 continuous parameters using gradient methods without problems. But we are now dealing with multi - physics electro - heat problems to which these considerations b ased on single physics systems do not apply. 1.4 Method of Optimization The Genetic Algorithm 1.4.1 Decision to use the GA Having settled on a zeroth order method of optimization that does not need gradient information about the object function, we take cognizanc e that most zeroth order methods are statistical so that several - fold more object function computations need to be made. The Random Search method is too random in its searching and we need something more systematic in its searching of the domain space so a s to not lead to excessive computation times in relation to gradient methods. 11 The practicable alternatives are simulated annealing and the genetic algorithm, both good methods widely used in industry and a part of commercial software [25] . Going by the literature Preis, Magele and Biro [26] , staunch advocates of the zeroth order ev olution strategy, merely say it is competitive with its higher order deterministic counterparts (which we take to mean the same in time at best), but claim that we agree with because search methods will ne ver see mesh - induced artificial local minima as a problem [15] . The weight of evidence seemed to favor the genetic algorithm over simulated annealing but not firmly so. So we did a quick study, the results of which, shown in Figure 1 - 5 , support the genetic algorithm. This study was done on the coupled electro - heat problem described in greater detail in chapter 3 . We note that in results from other disciplines besides finite element optim ization, Manikas and Cain [27] search of rec ent literature also confirms this finding [28] . We therefore decided to work exclusively with the genetic algorithm. Figure 1 - 5 Genetic Algorithm Speed Compared with Simulated Annealing 12 1.4.2 The Genetic Algorithm 1.4.2.1 Introduction The Genetic Algorithm (GA ) is a zeroth order search optimization tool, which work s differently when compared to other search optimization methods. Because of its broad applicability, eas e of use , computational inexpensive ness and the global perspective it receive s [6] , GA ha s been increasingly applied to various search and optimization problems. The g enetic algorithm concept was first introduced by John Holland of University of Michigan, Ann Arbor. From that time he and his students contributed much to develop this field. Most of the initial works can be found in several conference proceedings and journals [29 - 30] . However t here are several text books that now exist for GA [28 34] as it matured . In this section, the basic GA working principle is described. 1.4.2.2 Basic Concepts Biological Background Man made many things gaining inspiration from nature such as the crane ( the machine) from the crane ( the bird), the submarine from fishes, the aircraft from bird s , the micro - processor from the brain processor, and artificial neural networks from biological neural networks. Likewise GA was developed by inspiration from the reproduction process seen in nature. The nam e, GA borrow s its working principle procedures are motivated by principle s of natural genetics and natural selection. It is interesting to see terminologies in GA which originated from a biological backgr ound. Every organism consists of cells. A cell consists of a set of chromosomes. Chromosomes consist of genes which are the block s for DNA. The complete set of chromosomes is called a genome. GAs borrow the same terminology too. Moreover when two organism s mate they share their genes. T he result ant offspring will have both parents genes. A chil d may be more father - like 13 or mother - like. In GA this process is called crossover. The evolutionary process can be expedited by improving the variety of the genes. T h is may change in bit error. Likewise in GA there migh t be bit error in the offspring. T his is called mutation. There is important inspiration from nature i n the survival of the fittest concept. Powerful, wealthy, good looking, educated people get more opportunities than the others to bear offspring . But this is random. Approximately 10% of couples do not have children biologically or they do not want to have them . In GA the chromosomes with better fitness score have mor e opportunity to create the next generation. 1.4.2.3 Working Principle GAs begin with a set of solutions (represented by chromosomes) called the population. Solutions from a population are taken and used to form a new population. This is motivated by the possibil ity that the new population will be better than the old. Solutions are selected according to their fitness to form new solutions ; the more suitable they are the more the chances they have to reproduce. ators called crossover and mutation. These operators are used to create a new population. Let us consider the can problem which is used by Deb [36] to explain the GA working principle . Consider a cylindrical can which has diameter and height . This can need s to have a volume o f at least 300 ml. The object of the design is to minimize the cost of the can material. For t his we first write the object function corresponding to the area of tin r equired. Minimize (1 - 12) Subject to (1 - 13) 14 Variable bounds are (1 - 14) (1 - 15) Here is the cost of the material per square cm, and the upper and lower boundaries of and are [ ], [ ]. 1.4.2.4 Representing a Solution Binary en coding is the most common way to represent infor mation contained in a chromosome . In GAs, it was first used because of its relative simplicity. Let us say that a five bit code is used to represent each design parameter. Then design parameter s and will be represented by a string of length 10. The fo llowing string represents a can of and ( Figure 1 - 6 [36] ). The cost of the can is marked as 23 currency units. Figure 1 - 6 A typical chromosomal representation The l ower and upper limits are considered to be zero and 31 respectively. Choosing the lower and upper limits like this allows GAs to consider only integer values in the range (0, 31) . Bu t, GAs are not restricted to using only integer values in the above range, in fact GAs can be assigned to use any other integer or non - integer by changing the string length and the upper and lower limits by the following mapping function : (1 - 16) 15 w here is the string length used to code the - th parameter and is the decoded value of the string . The a bov e mapping function in (1 - 1 6) allow s us to take any arbitrary level of precision to be achieved by using a long en ough string and also to take positive and negative values. 1.4.2.5 Assigning a fitness score Once a solution is created by genetic operators, it is necessary to evaluate it, pa rticularly the object function and constraint function. For example, the object function of the above can will be (1 - 17) a ssuming 4 unit. Since the object is minimizing , the solution with lower fitness score will be the better solution. T he fitness score is defined in terms of the o bject function . In most cases the fitness score is made equal to the object function. In some cases fitness score might be defined as (1 - 18) Here is to be minimized and the fitness score has to be maximized for the genetic algorithm . The w orking principle of a GA is shown by a flowchart in Figure 1 - 7 . In Figure 1 - 7 , first we randoml y generate hundreds of solutions ( and ) as vectors (each called a chromosome) and this set is termed the initial population. With parallels to evolution, a new generation is to be created based on the best of this population. Then the fitness sco re for each {p} is calculated and checked as to whether there is a score at 1 or close enough for our purposes. This computation involves computing the object function according to (1 - 1 2) and evaluating with the constraint function (1 - 13). If the termination criterion is not satisfied, the 16 population of solutions is changed by three operators called selection, crossover and mutation , and a new population is created. Figure 1 - 7 Flowchart of the work ing principle of a GA Figure 1 - 8 A random population of six can s Figure 1 - 8 [36] shows the random populat ion of six cans. The fitness score of each can is marked on the relevant can. T he t wo cans which have fitness score of 11 and 9, though they have a good fitness score, do not satisfy the constraint function (1 - 13). This means they do 17 not have 300 ml volume . Therefore these cans are penalized by ad ding the extra artificial cost , so as a result their fitness score will bad. 1.4.2.6 Selection Chromosomes are selected from the population to be parents to crossover. The problem is how to select these chromosomes . According to Darwin's evolution theory , the best ones should survive and create new offspring. There are many methods on how to select the best chromosomes; for example roulette wheel selection, Boltzman selection, tournament selection, rank selection, s teady state selection, elitism and some others [37] . In the following we illustrate the tournament selection and elitism selection. As the name suggests, tournaments are played between two solutions and the better solution is chosen and placed in a new popul ation set. If done systematically, each solution can be made to participate in exactly two tournaments. The best solution in a tournament will win both times; therefore there will be two copies of the best solution in the new population. In this way, the w orst solution will lose both times; therefore it will be eliminated from the new popula tion. It has been shown that tournament selection has better convergence and computational time than the other selection methods [37] . Figure 1 - 9 [36] shows the six different tournaments played between population members. When cans with a cost of 23 units and 30 units compete in a random tournament, 23 units is chosen. It is i nteresting to note that the better solutions have made themselves to have more than one cop y and the worst solutions have been eliminated from the population. 18 In elitism selection, the new generation is based on elitism which first copies the best half of the chromosomes to the new population without any changes ( Figure 1 - 10 ). Elitism can very rapidly increase the performance of GA because it prevents losing the best found solutions . The remaining half of the popul ation will be replaced by offspring of the best half after crossover, and mutation as shown in Figure 1 - 10 . Figure 1 - 9 Tournaments played between six population mem bers Figure 1 - 10 Changing design parameters in GA 1.4.2.7 Crossover Crossover is a genetic operator that combines two chromosome s to produce a new chromosome. The idea behind crossover is that the new chromosome may be better than the 19 parents if it takes the best characteristics from each of the parents. Like selection, there exist a number of crossover and mutation operators in the GA literature [38] . In s ingle point crossover, we randomly select one crossover point and then copy everything before this point from the first parent and then everything after the crossover point from the second parent. Let us illustrate the crossover operator by two solutions from the new population created by the selection operator. Figure 1 - 11 An illustration of the single - point crossover operator The cans and their string s are shown in Figure 1 - 11 . The c rossover point is taken at the third string at random and these contents after that point are swapped between the two strings. The created children solutions are likely to be good s trings. However, every crossover may not create better solutions. But, if bad solutions are created, they will be eliminated in the next iteration and will have a short life. 1.4.2.8 Mutation After a crossover is performed, mutation takes place. Mutation is an i mportant part of the genetic search: it helps to prevent the population from stagnating at local optima. The mutation operator simply inverts the value of a randomly chosen gene of a chromosome. Figure 1 - 12 shows h ow a string which has gone through selection and crossover operators has been mutated to another string. 4 th bit of t his string is mutated by a small mutation probability hoping that a slight change in the string would keep the diversity of the population. 20 Figure 1 - 12 An illustration of the mutation operation T here are indeed many variants in operators of the genetic algorithm [39] , but here we are interested in using any robust GA approach that works. 1.4.2.9 Real - coded GAs We have discussed GAs where a solution is represented by a binary string. Howev er, this is not the case always and variables taking r eal values can be used directly without conversion to binary form . In real - code d G A s , selection operators work the same as with binary - coded GA. But, using efficient crossover and mutation operator s is different [37, 38] . The real - coded GAs eliminate the difficulties of arbitrary precision in decision variables and the Hamming cliff problem [6] associated with binary string representat ion of real numbers. A drawback of encoding variables as binary strings is the presence of Hamming cliffs: large Hamming distances between the codes of adjacent integers. For instance, 01111 and 10000 are integer representations of 15 and 16, respectively , and have a Hamming distance of 5 because 5 bits are different . For the genetic algorithm to change the representation from 15 to 16, it must alter all bits simultaneously. Such Hamming cliffs present a problem for the algorithm, as both mutation and cros sover cannot overcome them easily. There are three properties of binary single - point crossover operator s which have been observed and the SBX operator for real - coded GAs has been developed [40] . These properties are i) both child solutions lie on either inside (contraction crossover) or outside (expanding 21 crossover) of the r egion bound by the parents solutions , ii) the distance from one child from a parent is exactly the same as the distance from the o ther child to the other parent , and iii) the overall probabilities of contracting and expanding crossovers are the same. In o rder to implement this crossover operator for any two parent solutions p 1 and p 2 , the spread factor has been defined as the ratio of the spread of created child solutions c 1 and c 2 (1 - 19) If , it will use the contracting crossover and if , it will use the expanding crossover. It is possible to calculate the probability of creating a pair of child solutions having a certain . That probability distribution has been approximated by a polynomial probability d istribution as follows: (1 - 20) Figure 1 - 13 Probability distributions u sed in SBX crossover with different index n In (1 - 20), n i s any non - negative integer. If n is large, there is a higher probability for creating solutions near the parents. If n is small, it allows us to create distan t solutions from the parents. 22 The f ollow ing procedure is used to create child solutions c 1 and c 2 from the parent solutions p 1 and p 2 using the probability distribution (1 - 20) [40] . Create a random number u between 0 and 1. Find a for which the cumulative probability (1 - 21) Knowing the value of , the child solutions are calcul ated as (1 - 22) This is how SBX crossover for the real - coded GA works. The main difference between binary - coded GA and real - coded GA is that , for real - coded GA, an explicit probability distribution is used. Kalyanmoy Deb has liberally made his well - commented generic code available on the web [42] which we adapted and applied for the first time to finite element coupled problem optimization [43] . 1.5 Computer Aided Design 1.5.1 In pursuit of Accuracy rical models constitute the core of computer - aided design . Classical methods have been referred to as exact methods whereas numerical methods are called approximate methods. Classical methods are limited by the mathematics by which the solutions are expres sed [44] . 23 Figure 1 - 14 The Solution Process As shown in Figure 1 - 14 [44 ] , for every problem, first a model is chosen and the problem is solved based on that model instead of the actual problem. When the solution is not reasonable, the model will be reexamined with a view to improving it. For example in the mechanics of proje ctiles , It work s most of the time. But when it comes to space travel, the solutions may not match measurements and will be modified to incorporate y. Likewise in electromagnetics, in classical electrodynamics modeling, say o f a transmission line, the assumption is made about symmetry in the direction of the line so that the problem is tw o dimensional. Thereafter it might be further assume d that the t ransmission line is circular. These assumptions permit the reduction of the problem to a simple one - dimensional di fferential equation of the form, (1 - 23) where is the current in the line. The solution of the model is as exact onl y as much as the model is exact; (1 - 24) w here is the reference point for . To get a more accurate solution, the model need s to be improved. For example the model of the earth above which this infinitely long line 24 runs might be modelled by considering the images of the line. This makes the problem two dimensional. With superpo sition, this too can be handled by our model and the corresponding solution method. Should this also prove inadequate for purposes of the required accuracy, we go further and allow for nearby transmission towers and their equipment like insulator strings, and the sag of the wire from tower to tower. At this point the model breaks down as the problem is no longer two dimensional and truly three dimensional. Although many three dimensional problems have classical solutions, they rely upon some form of neatnes s such as spherical shapes, homogen e ous systems and so on. This complex transmission line has no known classical solution. At this point a numerical solution is a must. Numerical solutions by definition are approximate insofar as derivatives and integrals assume linear variation over small discrete distances. An example of the approximation is for the derivative and integral [44] : (1 - 25) (1 - 26) Although a linear change of is assumed over the interval this so called approximation is increasingly more accurate for smaller such interv als. That is, while classical solutions are only as accurate as the models on which they are based, the numerical approximations for derivatives and integrals can be made as accurate as we want by our resorting to finer subdivisions of the relevant interva ls over the domain. This means the accuracy of n umerical methods depend s on the number of subdivisions. This involves data preparation (pre - processing), solving the problems, facing an excessive amount 25 of data corresponding to the numerous points of discr etization, make meaning of numbers and post processing. 1.5.2 Methods of Approximate Solution N umerical methods can be divided into two classes based on whether solving integral equations or differential equations [44] . Each method has advantages and disadvantages based on problem type. The B oundary Element Method (BEM) [45] is a popular integral method. For example in electrostatics, the Coulomb equation becomes the basis of the boundary element method. (1 - 27) Here (1 - 27) is discretized and the potential is expressed in terms of the charge distribution . But, on material boundaries there is a surface charge density according to the derivatives in the normal direction n at t he interface separating regions 1 and 2: (1 - 28) Equation (1 - 28) states that difference in normal flux density across the interface is the surface charge density. Therefore the integral Poisson equ ation relates at each point to the known charge distribution and distribution too may be unknown as the charges rearrange themselves because of Coulomb forces. Thus every discrete point has a n unknown related to the unknowns at all other discrete points. The resulting matrix equation is fully populated and not necessarily symmetric. Thus although, unlike in differential methods relating every point only to others in its neighborhood 26 thereby me asuring rates of change, every point in the domain has a relationship with all other points in the domain of solution when solving integral equations . However, here with integral methods, we have unknowns only where there are conductors and material interf aces that is we have fewer unknowns. Yet the method is not too popular especially because of the non - symmetric nature of the matrix associated with the solution equation, besides the fully populated nature of the matrix equation [44] . The best known differential methods are the finite difference and finite element methods. The finite difference method is best suited to homogeneous problems and therefore is very popular in solving the wave equation solving whi ch usually involves modeling in free homogeneous space. It is not so popular in the many problems where inhomogeneous materials are present. Here this thesis will focus on the finite element method because of its wi de applications and we use the finite ele ment method to solve the coupled problem. 1.5.3 The Finite Element Method The finite element method is one of the most powerful numerical methods available for solving partial differential equations which apply over complex shapes. Very often, in engineering sc ience, it is difficult to solve a partial differential equation, which applies over a complicated shape. Unlike integral methods and the finite difference method, it always yields symmetric sparse matrices which allow the use of the most powerful and elega nt matrix solvers available for use . It allows us to use a fine mesh in parts where fine detail has to be incorporated or high accuracy is required, while a sparse mesh may be used in other parts. 27 Inhomogeneit ies are not a problem and are handled by ensuri ng that an element is always ent irely in a homogeneous region. To describe the finite element method it is good to note that it has certain ing redients and to spell them out: Ingredient 1 A trial or test function on the domain of solution (the x - y plane in our 2 - dimensional demonstration) consisting of many interpolation functions with i = n. The test function would be a weighted sum of all the s and the s are functions that we prudently choose so as to best model the solution. This means having some knowledge of the solution. In this form method and was not quite successful since we do not usually know the form of the solution. So our postulates based on our choice of s is that the solution is of the form (1 - 29) Our task then would be to determine the numbers so as to yield the best fit of the actual solution that can be made by the trial functions chosen. Ingredient 2 A mesh o f finite elements of the domain involv e s a tessellation of the area of solution. Commonly the domain of solution is subdivided into triangles or quadrilaterals. In this lies the brilliance of the method. Once we divide the domain into small regions, in place of the complicated functions, our s, we may go for simple functions such as linear or quadratic polynomials. In simple terms, to approximate a squiggly function y(x) on a large domain, we 28 may need to add several s of different forms; but if we divide the domain into small bits called fin ite elements , a straight line on each bit or element would suffice. When the elements are smaller, our model will be more accurate . Thus on a triangular linear finite element, we have (1 - 30) where the functions 1, x and y take the place of th e three s we need for the small triangular domain. As we divide the domain into smaller and smaller triangles, the modeling of the actual shape of the unknown would become more and more accurate. This linear model is an approximation of the potential bec ause may not vary linearly and may have a complex relationship with x and y. However, x and y are linear first order functions. This is the key to determining a, b and c in a triangle. Figure 1 - 15 A Tr iangle with Three Nodes Usually from one triangle to its neighbor we would need continuity conditions so that the trial functions match on the edge common to two adjacent elements. For this a, b and c of one triangle need to have a relationship with the a, b and c of any adjacent triangle. This relationship being complicated, instead of a, b and c we deal with the three vertex potentials , i =1 , 2 & 3. So at each vertex of a triangle ( , wh ere values are known because we are solving the Poisson 29 equation in a specific domain divided into a specific mesh of known vertices connected into a triangular mesh, applying the l inear approximation ( Figure 1 - 15 ): Or (1 - 31) (1 - 32) = (1 - 33) = (1 - 34) where is twice the area of the triangle provide d the nodes 1, 2 and 3 go counter - clockwise ( and the negative of that if the nodes go clockwise); the numbers , and are the numbers 1, 2, and 3 in cyclic order (that is, the integers , and respectively are 1,2,3 if is 1, or 2, 30 3, 1 i f is 2 or 3, 1, 2 if is 3; or put another way mod 3 +1 and mod 3 +1 where the mod 3 or modulo 3 function gives the remainder when divided by 3 ) ; and (1 - 35) and (1 - 36) Ingredient 3 A functional (a function of a function, in this case of ) , minimizing which is the same as satisfying the differential equation [44] . There is an alternative approach called the Galerki n formulation relying on the theory of function spaces going into minimizing the residual when our trial function is put into the differential equation being solved but we will not get into that here. A limitation of the approach we describe, the variation al calculus approach, is that to take it we must have a functional corresponding to the solution of the differential equation wh ich is sought. Suffice it to say that the Poisson equation does have a functional corresponding to it: (1 - 37) We can justify this intuitively. Since the electric field strength and the flux density , we see that the first term is the stored energy density ½ times elemental volume dR. Since by definition the potential is the work done in bringing a unit charge to where it is, the second term is the work done in bringing the charge , the charge density times elemental volume, to where it is. The essence of the functional is that it says that the difference between the two co - energies ought to be a minimum when the Poisson equation is satisfied. But let us see this more precisely in mathematical terms. 31 Let the potential take a small excursion to , where is a constant an d is a small arbitrary function. The modified functional then is (1 - 38) Subtracting (1 - 37) from (1 - 38), the chang e is (1 - 3 9 ) where for small changes, the terms of order 2 (designated ) in may be dropped for being negligibly small but are retained to make a point we will encounter soon [44] . We note that is from which is strictly positive. To proceed we need to use a vector identity expressing the divergence of a vector scaled by a scalar s: (1 - 40) toge ther with the divergence theorem: (1 - 41) where the surface S bounds the domain R. Integrating the vector identity just encountered and applying the divergence theorem to it or (1 - 42) Now set s = and = (1 - 43) We note further that the elemental surface vector points in the normal direction so that where is a unit vector in the normal direction. Further since 32 and therefore , and since the direction x is arbitrary according to how we orient the axes, we may wr ite . So the above equation upon rearranging becomes (1 - 44) Using this in the expression for the change in the Lagrangian functional (1 - 45) Now the boundary S has either Dirichlet or Neumann conditions on each point of it. Where the Dirichlet condition applies since the potential is specified and cannot vary. On those parts of S where the Neumann condition is specified, Therefore the surface integral vanishes because the integrand is zero at every point on S, thereby leading to (1 - 46) Now we are in a position to say a) If the functional is at an extremum, whether a minimum or a maximum, the left hand side is zero since changes are flat at the point of extremum. Further the second order terms may be neglected so that (1 - 47) But since the change in potential is arbitrary, this is possible only if the other factor of the integrand or (1 - 48) meaning that the Poisson equation is satisfied at the extremum. 33 b) We need one more consideration to show that th e extremum is a minimum. At that extremum since the residual of the Poisson equation , we must have (1 - 49) which is a strictly positive quantity as we have noted so that any changes in about the extremum of will m ake rise , meaning that extremum has to be a minimum. 1.5.4 Two Dimensional , Linear , Triangular Finite Elements We used two dimensi onal triangular finite elements in this thesis ; therefore this thesis will presen t this model through an example. Figure 1 - 16 [44] shows a two dimensional system governed by the Poisson equation with zero Neumann conditions along the upper and lower boundaries and Dirichlet c onditions on the left and right boundaries , where is the thermal conductivity and is the heat distribution. Demonstrating through this example, th is heat flow problem ( Figure 1 - 16 ) has and in the shaded region and and elsewhere. Although this is a problem governing the flow of heat, the equations are similar to those developed for an electrostatic system with and . Let u s see how T will be solved everywhere u sing hand calculations. Figure 1 - 16 Heat Flow Problem 34 By linear approximations or trial functions, we have (1 - 50) The f irst order approximation on a triangle leads to (1 - 51) Now (1 - 52) The functional (1 - 53) w here is a constant that may be pull ed out of the integral sign so that the integral of dR is simply the triang le area . Likewise, by assuming that q is a constant at q 0 in a triangle, that too comes out of the integral sign so that the integral of the linearly varying T is the average of T at the three vertices times a measure o f the integration area which is ; this is a variant of - dimension for 2 - dimensions. Thus we have, substituting for b and c, while noting that (1 - 54) (1 - 55) With not ation for defined by correspondence for terms, and noting that a scalar b is its own transpose, we have (1 - 56) 35 where the brackets have been dispensed with in the absence of there being multiplicative compatibility. Similarly (1 - 5 7 ) Therefore (1 - 58) where, in corresponding notation, the local matrices (1 - 59) and (1 - 60 ) Now let us c onsider the first triangle of the mesh as shown in Figure 1 - 17 . Where we start numbering the nodes (whether 8 - 3 - 7 or 3 - 7 - 8 or 7 - 8 - 3) does not matter but that the order sho uld be counter - clockwise is important. 36 Table 1 - 1 The Local Matrices of the Six Elements of Figure 1 - 16 Ele. 1 PMatrix QMatrix Verts 8 3 7 8 4 - 2 - 2 1/6 3 - 2 2 0 1/6 7 - 2 0 2 1/6 Ele. 2 PMatrix QMatrix Verts 1 7 3 1 4 - 2 - 2 1/6 7 - 2 2 0 1/6 3 - 2 0 2 1/6 Ele. 3 PMatrix QMatrix Verts 3 4 1 3 1 - 0.5 - .5 0 4 - 0.5 0.5 0 0 1 - 0.5 0 0.5 0 Ele. 4 PMatrix QMatrix Verts 2 1 4 2 1 - 0.5 - .5 0 1 - 0.5 0.5 0 0 4 - 0.5 0 0.5 0 Ele. 5 PMatrix QMatrix Verts 4 6 2 4 1 - 0.5 - .5 0 6 - 0.5 0.5 0 0 2 - 0.5 0 0.5 0 Ele. 6 PMatrix QMatrix Verts 5 2 6 5 1 - 0.5 - .5 0 2 - 0.5 0.5 0 0 6 - 0.5 0 0.5 0 37 Figure 1 - 17 Triangle with internal generic numbers For this triangle (1 - 61) Therefore area A = ½. We could have said Area = ½ x base x height = ½ but this is not general enough an approach for triangles whose sides are not parallel to the Carte sian axes. We have (1 - 62) (1 - 63) Observe that only relative coordinates are important since we are dealing with differences in coordinate values which fact makes the actual values of the coordinates irrelevant. Now 38 = = 4x0.5 = (1 - 64) = (1 - 65) Table 1 - 2 The Six Local Matrices Added into the 4x4 Matrix The 4x4 Matrix [ ] of size 4x1 0 4 0.5 0.5 0 0 0 0 0 - 0.5 0 0 0 - 2 - 0. 5 0 0 0 0 0 0 0 0 0 200+1/6 0 0 0 0 0 0 0 0 - 0.5 0 0 0 0 0 1 0.5 0.5 0 0 0 0 0 0 0 0 0 - 0.5 - 0.5 0 0 0 0 0 0 0 0 - 2 - 0.5 0 0 0 0 0 0 0 0 0 2 2 1 0 0 0 0 0 - 0.5 0 - 0.5 0 1/6 200+1/6 0 0 0 0 0 0 0 0 0 0 0 0 0 - 0.5 - 0.5 0 0 0 - 0.5 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 39 The local matrices of the six elements of are shown in Figure 1 - 16 . The 4 rows and columns of the matrix [P] are filled up as shown in Table 1 - 2 . For clarity, each position of the global matrix [ ] is divided into 6, each sub - location corresponding to the contribution from one of the six elements. Similarly the column vector shows its contributions from th e 6 elements as it is built up. The numbers 200 in the right hand side column vector come from the two additional columns of [P] corresponding to the known T values at 100. The numbers zero when they need to be added are not refl ected. Thus we have from Table 1 - 2 , (1 - 66) Solving (1 - 67) So finally the solution gives the temperature at nodes 1, 2, 3 and 4. 40 2 EDDY CURRENT BASED DEFECT DETECTION AND CHARACTERIZATION 2.1 Single Physics Problem Our intention is to study and establish the most efficient, GPU based methodologies fo r shape optimization in two - physics system s . To avoid the complexities of a two - physics problem and learn the issues , we begin our study by examining shaping algorithms on a single physics problem for nondestructive evaluation with the view to extending th e principles to two - physics systems . Eddy C urrent T esting (ECT) is a widely accepted, cheap and portable method for the detection of cracks and other defects in conductive materials [46] . In NDE problems, the test measurements have to be matched with the computed results where as in inverse design problems , the desired performance results have to be matched with computed performance results ( I nverse design problems have been discussed in 1.2 ) . Therefore when an object function is defi ned for NDE problems, has to be replaced by in (1 - 2). Otherwise, NDE problem s and design problems are mathematically similar [7] . This then is the foundation for studying NDE shape optimization and then extending it to two - physics shape optimization. In particular we will study defect characterization where we will optimize the s hape of a postulated defect to make the measured and computed fields match. corrosion when a sample is exposed to a very corrosive environment for 3 months [47] while other researchers have established that eddy current testing can identify corrosion to depths down to 15 mm [48], [49] . For s ubmarines and aircraft , engineers employ eddy current probes for testing thin (up to 1 inch) conducting systems [50] , and ultrasound for thicker submarine 41 hulls [51] . For deeper penetration, pulsed eddy currents [52], [53] or remote field eddy currents are used [54] . ECT Technology is therefore perfect for the purpose at hand. 2.2 Introduction to Eddy Current Testing 2.2.1 Eddy Current Eddy currents are found in any conducting material which is subject to time varying magnetic fields. If there is a defect present inside the conducting material, the flow of eddy current will be interrupted. Due to this interruption and changes in magnetic field, one can get information about the defect inside the material. A thorough analysis about eddy current s and their behavior can be found in [55] . 2.2.2 Basis of eddy current testing When an AC current flows in a coil which is in close proximity to a conducting surface , the magnetic field from the coil will induce eddy currents in that surface ( Figure 2 - 1 .a [56] ). The eddy currents will affect the loading on the coil and thus its impedan ce. Figure 2 - 1 Eddy cu rrent flow on a surface with (a) and without (b ) crack 42 Assume there is a crack on the surface of the material ( Figure 2 - 1 .b ). This will interr upt or at least reduce the eddy current flow. Therefore it will decrease the loading of the coil and increasing the effective coil impedance. Further, the voltage induced in an exterior sensing coil will be changed. This is th e basis of eddy current testin g; by monitoring the voltage across the coil in such an arrangement we can detect the defects inside the material. There are some factors which will affect the eddy current response in a probe. The conductivity of material has a direct effect on the eddy current flow. If the conductivity is high, the flow of the eddy current also will be high. The material permeability ( ) also has a direct effect on eddy current flow. This may be described as the ease with which a material can be magnetized. For non - ferrous metals such as copper, brass, and aluminum , is the same as that of air. For ferrous metals the value of will be several hundred times that of air, and this has a very significant influence on eddy current flow. Eddy current flow is greatly affected by the frequency we cho o se for the exciting system . The effect of frequency on eddy current flow is discussed in se ction 2.3.4 . The proximity of the probe to the material also has to be considered. The closer a probe coil is to the surface , the greater will be the effect on that coil. Other than these factors, the geometry of the defect pla ys an important role in eddy current test ing . This will be discussed in the following section. 2.3 Detection of Defect 2.3.1 Field Changes Our goal in this chapter is cha racterizing the shape of the defect using eddy current . But, before we move into defect charact erization, a study on defect detection has to be done. 43 Defects can be in many different orientations. Some of these defect orientations make it hard for the defect to be detected, and therefore such defects will be impossible to characterize by our invers e algorithm. Therefore in this section we analyze the defect orientations and its relation to detectability. So this study will be useful for the defect characterization later in this chapter. Even though ECT is a good and effective method to detect a def ect, the reality is that, thinking qualitatively about the three situations in Figure 2 - 2 , the depth for detection depends on defect orientation. In the first situation , when a defect is parallel to the surface, th e defect will be detectable through the interruption of the eddies. In the second situation, it is perpendicular to the surface but not reaching the surface; we would expect the skin - effect to allow the currents to flow along the surface. So it will be har d to detect the defect in this second case. Whereas in the third situation if the defect is perpendicular to the surface without a gap for the eddy current to flow along the surface, the defect would interrupt the flow of eddy currents which would find it difficult to go deep because of the skin effect. Therefore this will be easily detectable. On the other hand if the defect is thin and long, a 3D model will show no current interruption. Figure 2 - 2 Penetra tion of Eddy Currents in the Presence of Defects 44 This interruption , if it occurs, would make the presence of defects easily detectable. Different types of defect models with the ECT coil are shown in Figure 2 - 3 [57] . Defects a, b, c, d and e in Figure 2 - 3 show the difference in angle, depth and length. The defects a , b, and c are at different angle s , whereas c and e are at different depth s . At the same time, though defects d and e are at same depth, they vary in length . Moreover, if the axis of the external exciting ECT coil is parallel to the surface, our considerat ions need to be entirely different. This study therefore through a series of finite element analyses of parameterized geometries of a single defect, seeks to identify the limits of eddy current testing in NDE. This study is confined to steel plates for a rmy ground vehicle armor [57] . A qualitative relationship between the depth of a detectable defect, the frequency of excitation and the geometric parameters is sought to help engineers ch oose whether or not to use ECT. Figure 2 - 3 Defect Model: The model has defects in different angle , depth and length 45 A defect in a material can be varied in shape and location. The defect detection method should concentrate on the characteristics of the defect. These effects are analyzed in this work for varying size, shape and location of a defect in materials. An iterative approach is presented that repeatedly emp loys the finite element technique for modeling the forward problem to calculate the effect caused by the defects in a steel plate. We calculate the magnetic flux density for the known defects from the forward problem using the finite element method [48] and examine the extent to which the exterior field is altered by the defect to see if the presence of the defect can be discerned through the changes in the measurable external field. 2.3.2 Main Eq uations 2.3.2.1 Maxwell Equations The electromagnetic field equations are the basic equations to deal with for design analysis purposes. Though there are many basic laws in electromagnetics such as by Coulomb, Lenz, and Ampere, the general Maxwell equations are m ore relevant and appropriate to discuss here as they are directly involve d with our work. The four Maxwell equations [58], [59] are; (2 - 1) (2 - 2 ) (2 - 3 ) (2 - 4) 46 where is electric f lux density , is the electric charge densi ty , is the magnetic flux density, is the magnetic field intensity, is the current density , and is the electric field intensity . Here is related to by the constitutive relationship (2 - 5 ) where is the material ele ctric conductivity. is related to by the constitutive relationship (2 - 6 ) where is the permeability. In the same way is related to by the constitutive relationship (2 - 7 ) where is the permittivity . In addition to the se equations, based on (2 - 2) the magnetic vector potential is defined by (2 - 8) Now for the vector to be unique and therefore determinable, i ts divergence has to be defined. This we do by the Coulomb gauge: (2 - 9) The subject o f computational electromagnetics may be divided into low and high frequency electromagnetics. It is a natural separation in view of the governing equations and the specialized nature of the division between low and high frequency electromagnetics. Typical low frequency devices are electrical machines, electronic devices, transmission lines, and magnetic recording heads whereas high frequency devices are waveguides, resonant cavities, and radiating devices such as antenna e . In terms of equations, the differe nce between these two 47 systems is that in low frequency devices, the displacement current is negligible while in high frequency devices it is not. Therefore , for low frequency devices, (2 - 3) will reduce to (2 - 10 ) In contrast, for high frequen cy devices , although is small , its rate of change is high, making of (2 - 3) directly applicable. However in many problems in lossless media, since the medi um has zero conductivity , no conduction current may exist in keeping with (2 - 5) . 2.3.2.2 Eddy Current Equation In magnetostatics, an electric field causes current and the current causes a magnetic field. If time variation is introduced, the resulting magnetic field does affect the initial electric field via (2 - 4). In the same way the electr ic field affect s the magnetic field with the extra displacement current term (2 - 3). Let us consider (2 - 4 ), the non - divergent may be modelled by a vector potential as in (2 - 9): (2 - 11) The two vectors and have the same curl, which is possible only if: (2 - 1 2 ) I t is realized that the term is the externally imposed electric field driving the current and is the induced electric field. Combining (2 - 5) and (2 - 10 ), and substituting (2 - 6 ), (2 - 8 ) and (2 - 11), 48 (2 - 13) s o that (2 - 14) a nd (2 - 1 5 ) Setting imposed current density , (2 - 1 6 ) and u sing the vector identity for the curl of the curl and (2 - 9), (2 - 1 7 ) we finally have , using (2 - 1 6 ), (2 - 1 7 ) and using phasor representation, where differentiation with respect to time is the equivalent of pre - multiplication by (2 - 18) This reduces to our final e quation (2 - 1 9 ) For eddy current problem, (2 - 1 9 ) is the key equation to be solved to find the magnetic vector potential for a given imposed current density . Once we find , the magnetic flux density c an be calculated using (2 - 8) and the electric field intensity can be calculated using (2 - 1 2 ) according to the requirement of the problem. Calculating from (2 - 1 9 ) using the finite element method is described in detail in the following section. 49 2.3.2.3 Finit e Element Computation for E ddy C urrent P roblem Magnetic fields in a ferromagnetic material can be generated by placing an AC (Alternative Current) coil on top of the material. For AC magnetization, the magnetic vector potential , and the exciting curren t density , are related by as described in (2 - 1 9 ) . The corresponding functional can be written , similar to the way we derived (1 - 53) ; (2 - 20 ) Since the first two part s of (2 - 16 ) are similar to the terms we have already derived in ( 1 - 58 ) , the last part of (2 - 16 ) which is is computed as follows; 50 (2 - 21 ) wh ere in a triangle element has be en written in triangular coordinates and , But , (2 - 22 ) By applying (2 - 22 ) i n (2 - 21 ), and rearranging (2 - 23 ) Now, substituting (2 - 23 ) i n (2 - 20 ), 51 (2 - 2 4 ) w here, in corresponding notation, the local matrices (2 - 2 5 ) and (2 - 2 6 ) Finite element analysis [20] provides the solution to ( 2 - 1 9 ) by applying certain boundary conditions. The local matrices of elements will be added to the corresponding positon of the global matrix to be solved for . This leads to the fi nite element matrix equation ( 2 - 2 7 ) Equation (2 - 2 7 ) is solved in the same way as explained in section 1.5.4 , calculating local matrices and adding them to the global matrix equation to be solved for . 52 2.3.3 Problem Statement The numerical example in Figure 2 - 4 [57] is used to study the detectability of d efects in different orientations . The coil ( with , current density , and ) excites the magnetic field in the steel plate (with and current density ). The conductor is surrounded b y air (with and current density ). Any defect in the material should affect the magnetic f lux density between the AC coil and the material. The magnetic flux density in the direction is calculated at 10 points as sh own in Figure 2 - 4 labeled as the measuring line. First, for the steel plate with no defect is calculated at the measuring line and named . Here is the number of a particular po int on the measuring line. Since our study is about whether the defect is detectable or not, we defined a single defect with parameterized geometries as shown in Figure 2 - 4 . By changing the par ameters of the defect, we change the position of the defect as shown in Figure 2 - 5 . and calculate at the measuring points, which we name . Then the effect on flux density as described by (2 - 28) (2 - 28) was calculated for each defect and tabulated. By considering the maximum of , we e stablish the defect detecting quantity which is measured. is the maximum value of flux density ratio between flux changes caused by the defect and flux without the defect, where a defect with length is at depth from the material surface and rotated by angle clock - wise e parallel to the steel surface ( Figure 2 - 4 ). 53 Figure 2 - 4 Numerical Model Figure 2 - 6 shows how the ratio varies with defect locations. When is high, there is a higher chance that the defect will be detected. Figure 2 - 6 .a sho ws the results for a horizontal defect when the depth d increases , decreases. So that means when the depth of the defect is higher, it would be harder to detect as to be expected. Figure 2 - 6 .b shows when the defect is vertical how varies with depth. Here also when d increases decreases. Here corresponds to the right - most part of Figure 2 - 2 . Anyway when we compare a horizontal defect and vertical defect, the horizontal defect has a higher chance to be detected, for the given value 54 This is to be expected because a horizontal defect will interrupt flux more since the flux flow is alon g the surface. Fig.6.c shows how varies with angle : goes to a maximum when . So we could say that when the angle of the defect is , there is a higher chance to be detected. a.No defect b. and c. and d. and Figure 2 - 5 Equipotential lines for the defects at various positions and orientations 55 a. varies with d for horizontal defect b. varies with d for vertical defect c. Figure 2 - 6 with defect location 56 2.3.4 Experimental Work in Defect Detection Figure 2 - 7 Rectangular defect on steel plate Experimental work was done to verify the computational work of defect detection. The coil was moved along the x axis of the steel plate ( Figure 2 - 9 ) and the voltage induced was measured as shown in Figure 2 - 8 . We created a rect angular defect on a steel plate as shown in Figure 2 - 7 . A close - up look of the defect is shown in Figure 2 - 9 . Alternati ng sinusoidal current was excited in the coil with 300 mV at 10 0 0 Hz using the waveform generator. A pick up coil was wound on top of the AC coil to measure the voltag e induced. Since the voltage induced is very small, we used the lock - in amplifier to measure voltage. We moved the coil along the surface of the steel and measured the induced voltage at the measuring points indicated i n Figure 2 - 7 is proportional to the change in flux ( 2 - 29 ). ( 2 - 2 9 ) The results were tabulated and plotted in Figure 2 - 10 . From the plot, when the coil is just above the defect, the voltage induced is at its maximum. The test was repeated with different input voltages and different frequencies and we obtained the same behavior as in Figure 2 - 10 . In this way we established that we could detect the defects and therefore moved to the defect 57 characterization part of this study. The frequency of excitation also plays an important role in detecting the defect. As in equation (2 - 25 ), when frequency goes high, will be high and can be measured easily. But the eddy current will flow close to the surface due to the skin effect. Therefore we cannot detect defects which are deeply embedded inside the plate. Figure 2 - 8 Lab setup Figure 2 - 9 Close - up - look of the defect 58 Figure 2 - 10 Voltage induced on pickup coil 2.4 Defect Characteri zation 2.4.1 Design Parameters After detecting the defect, we investigated more on defect characterization. It is important to know the size and character of the defect after establishing that there is a defect inside. So in our work we investigate and establish a procedure for defect characterization so that a decision to withdraw a defective part may be thought - out, and be justifiable. The methodology at present [7] examines the response of an army ground vehicle hull under test to an excitatory signal from an eddy current probe. By knowing the response when there is no defect, if the response is different b ecause of the defect, the test object is presently flagged as defective and the plate is sent for repairs without assessing if the defect is serious enough for removal from service. In our work, we extend that methodology to defect characterization. An ite rative approach is presented that repeatedly employs the finite element technique for modeling the forward problem to characterize the shape of defects in a steel plate. 59 Figure 2 - 11 Defect Model We can cal culate the magnetic flux density for known defects from the forward problem. But in our inverse problem we need to know the characteristics of the defect for that field configuration. In design optimization, the problem geometry is defined in terms of desi gn parameters contained in a vector ( Figure 2 - 11 ). An object function is defined as the sum of the squares of the difference between computed and measured (defect) perf ormance values at measurement points i, ( 2 - 30 ) is a function of defect shape. By minimizing the object function F with respective to the parameters by any of the op timization methods, the characteristics of the defect can be estimated. The computational process in the defect identification system is shown in Figure 2 - 12 . It needs to solve the design parameters . First, th e mesh needs to be generated for the given design parameter. Mesh generation is a very important part of finite element analysis based design optimization. We need to use a parameter based mesh generator, because in e very iteration , the 60 mesh has to be gene rated automatically when design parameters change. Therefore we use a special script - based parametric mesh generator [60] using as backend the single problem mesh generator Triangle [61] to get th e corresponding finite element solution for the magnetic vector potential . Triangle had to be modified for seamless, nonstop optimizations by a colleague in our research group [60], [62] . His work is used here. After computing by finite elements , the magnetic f lux density ( ) is computed and the object function is evaluated. When the object function is minimum , the parameters will be found. If is not minimum , the design parameters will be chang ed using the optimization method being employed . Figure 2 - 12 Design Cycle for the computation process For the known defect which is known as true profile, was calculated first. Then faking that we do not know the shape of the defect, we generated hundreds of different shapes of defects using our algorithm . was calculated each time. Comparing and the object function is evaluated each time and optimized . When is minimum, the reconstructed shape of the defect is generated. 61 2.4.2 Numerical Model The numerical model in Figure 2 - 13 [57] is used to validate the proposed algorithm. The coil (with , and current density ) excites the magnetic field in the steel plate (with and current density ). The conductor is surrounded by air (with and current density ). The magnetic f lux density in the direction B y using 10 points in the interval as shown in Figure 2 - 13 labeled as the measuring line. Figure 2 - 13 Numerical Model for defect characterization On each node on the defect, the vertical displacements are selected a s design parameters Figure 2 - 11 . In our numerical model we have 8 geometric parameters contained in the vector ( h 1 , h 2 , h 3 , h 4 , h 5 , h 6 , h 7 , h 8 ). The measuring line located at y = 4.5 cm, is sampled into 10 equa lly spaced points and tolerance boundaries on value are set to . 62 Each design variable is represented by 10 bits. For testing we took a particular defect as = {2.0. 1.2, 2.4, 2.2, 2.7, 3.0, 3.2, 3.5} cm and computed the field ( ). Now our algorithm has to reconstruct ). The d esign parameters are changed and the magnetic flux density is calculated along the measuring line. The o bject function is evaluated by comparing with . 2.4.3 Constraints Handling In inverse problem design optimization, getting a practically manufacturable shape is important. An erratic undulating shape with shar p edges arose when Pironneau optimized a pole face to achieve a constant magnetic flux density [12] . Since our design optimization is in defect characterization, there is no need to impose constraints to get a smoot h manufactural shape. But we have to get a single defect with a realistic shape. This problem was overcome by imposing constraints [63] . So as to maintain a realistic shape with a single defect , node 8 is on top of node 1, node 7 on top of node 2, node 6 on top of node 3 and node 5 on top of node 4 as in Figure 2 - 14 , we had imposed the constraints as h 8 >h 1 , h 7 >h 2 , h 6 >h 3 and h 5 >h 4 ; Therefo re we could get a single and rea li stic defect as shown in Figure 2 - 16 . 2.5 Results and Discussion Several simulations of the defect characterization problem were run and the results were tabulated. The best fitness score was achieved when the populat ion size is 200 and for 200 iterations. Figure 2 - 14 shows the optimum shape of the defect after 200 iterations for a 63 population size of 200. To measure the error in the reconstruction of the defect, using e , the le ngth from the centroid of the true profile to each point was calculated : (2 - 31 ) w here is length from the centroid to the true profile, and is the length from the centroid to the reconstructed profile as show n in Figure 2 - 15 and is the number of coordinates in the profile. The c alculation process for our numerical model is shown in Table 2 - 1 . This calculation gives error = 12% , which means an efficiency of 8 8 % in the average reconstruction of the defect. For the reconstructed profile, the f inite element solution of the magneti c vector potential is shown in Figure 2 - 16 . Figure 2 - 14 Optimum shape of the Reconstructed Defect 64 Figure 2 - 15 C alculating efficiency of reconstruction Figure 2 - 16 Solutions in Equipotential lines for the Numerical Model 65 Table 2 - 1 Cal c ulating the efficiency of reconstr uction true profile Reconstructed Profile No x y x y 1 8.00 2.00 1.589222 8.00 1.89 1.628872 0.0249 2 9.00 1.20 1.416201 9.00 1.3 5 1.276959 0.0983 3 10.00 2.40 0.515388 10.00 2.11 0.649788 0.2608 4 11.00 2.20 1.534805 11.00 2.38 1.506992 0.0181 5 11.00 2.70 1.510174 11.00 2.61 1.502406 0.0051 6 10.00 3.00 0.689656 10.00 2.65 0.515388 0.2527 7 9.00 3.20 0.840015 9.00 2.81 0.5755 22 0.3149 8 8.00 3.50 1.789029 8.00 3.57 1.828121 0.0219 It is important to know the efficiency of the algorithm s we develop and the efficiency of our coding before we move into more complex and time consuming problems like the coupled problem which is our ultimate goal . To study the efficiency of our inverse algorithm, several tests are required and also different numerical models have to be checked besides the numerical test model we have tested here . First we started our test with the optimization met hod we use, the genetic algorit h m. We have already discussed (in section 1.4.2 ) that in GA, solutions can be represented by a binary string or real numbers. There are advantages and disadvantages to whether we use binary or rea l numbers to represent the solution. Therefore we did the tests with both binary and real numbers for the same numerical model problem we discussed. Test s were carried out for the different population size s and different number of iterations and measured t he time taken to compute the problem and tabulated these in Table 66 2 - 2 . In the same way for the different population sizes and different number of iterations, the best fitness score is measured and tabulated in Table 2 - 3 . From Table 2 - 2 and Figure 2 - 17 , we can see that solutions th at are represented by real numbers give solutions faster than t he solutions th at are represented in binary numbers . Especially as the number of iterations and the number populations are increased the time g ap between real and binary solution times is increased. The reason would be that , in the binary based genetic alg orithm, every time the solution is evaluated, the encoding and decoding process es ha ve to be done; w here as in the real number based genetic algorithm implementation , this process in unnecessary as the solution is already a real number. When comparing the b est fitness score achieved by real and binary solution s in Table 2 - 3 , the binary solutions give slightly better solutions than the real number solutions . But if we consider 2 digit approximations, both real and bin ary solutions give almost the same fitness score. Therefore we conclude that the real number based genetic algorithm is more efficient than the binary implementation in our problem . So we decided to work with real number based genetic algorithm for all fu ture problems. After testing the genetic algorithm with real and binary solutions, we also tested with different shape s of the defect . For the known defect which we term the true profile, faking that we do not know the shape to measure the convergence , m any different shapes were generated using the optimization algorithm we use and matched using the objective function. When the objective function is minimum, the reconstructed profile is generated and matched with the true profile. We tested this process w ith two different defects , Defect - 1 and Defect - 2 to check the efficiency of the algorithm as shown in Figure 2 - 18 . We get 87.5 % of average reconstruction for D efect - 1 , and 92.7 % of average reconstruction for D efe ct - 2. 67 Table 2 - 2 Real and binary solutions time need to compute Table 2 - 3 Real and binary best fitness score achieved Population Size 30 iterations 40 iterations 50 iterations Time Taken(s) Time Taken(s) Time Taken(s) Real Binary Real Binary Real Binary 10 27 4.33 311.97 346.58 408.29 413.16 503.43 20 540.68 661.79 695.28 891.06 866.47 1105.42 30 861.79 1024.65 1095.82 1258.17 1588.32 1749.10 40 1319.36 1431.34 1666.79 1950.36 2024.03 2390.25 50 1583.73 1769.60 1994.77 2302.99 2323.93 2891.94 60 1757.78 21 31.94 2199.26 2534.36 2656.53 3303.85 Population Size 30 iterations 40 iterations 50 iterations Best Fitness Score Best Fitness Score Best Fitness Sco re Real Binary Real Binary Real Binary 10 0.0123 0.0019 0.0091 0.0020 0.0091 0.0015 20 0.0084 0.0012 0.0084 0.0008 0.0084 0.0008 30 0.0040 0.0009 0.0040 0.0009 0.0040 0.0009 40 0.0105 0.0017 0.0105 0.0017 0.0105 0.0017 50 0.0111 0.0017 0.0068 0.0017 0.0068 0.0017 60 0.0127 0.0017 0.0127 0.0016 0.0127 0.0000 68 a. 30 iterations b. 40 iterations c. 50 iterations Figure 2 - 17 Time required for real and binary solutions 69 a. Defect - 1 b. Close - up l ook of defect - 1 c. Defect - 2 d. Close - up look of defect - 2 Figure 2 - 18 True and Reconstructed shape of defects 2.6 Experimental Work in D efect Characterization The f requency of excitation is a very important fact or when it comes to experimental eddy current work which is where real testing occurs . When the frequency is high , like in the kHz range, the signal received is strong enough to measure by a pick - up coil. But, the depth of penetration was low for the high frequency signal. Therefore the defects which are not close to the surface of the steel would not be detected. On the other hand, low frequency signal s ha ve 70 good depth of penetration, but the signals are too weak to measure by a pick - up coil. Even when w e use a lock - in amplifier to amplify them, the signals were found to be weak of the order of 10 - 6 V. Therefore we used a n EC - GMR (Eddy Current Giant Magnetoresistive Sensor) to measure the signal. EC - GMR sensor s have been used to enhance the sensitivity of EC testing. This sensor measure s the magnitude of the magnetic flux directly [64], [65] . The EC - GMR sensor offers high sensitivity over a wide range of frequencies from DC to MHz [66] . An e xperimental study of automatic crack dete ction under a s teel fastener using a low frequency EC - GMR sensor was done by Yang et al. [67] . We use the same experimental set up to measure the magnitude of the magnetic flux , and then use that measurement for the inverse problem solution of defect characterization . The EC - GMR detection system contains a planar multi - line coil and set of GMR sensors m ounted on the line of symmetry as shown in Figure 2 - 19 [68] . The multi - line coil carrying current is employed as a sheet current source to induce the eddy current on the material we investigate. If there is a crack or defect in the material, the induced currents in uniform and linear patterns are distorted as shown in Figure 2 - 20 [68] . Figure 2 - 20 - a shows the induced eddy current patterns when there is no discontinuity. Figure 2 - 20 - b shows the patterns when there is a rivet and Figure 2 - 20 - c shows when there is a rivet with a defect. 71 Figure 2 - 19 Planar coil and GMR sensor a. No discontinuity b. A rivet c. A river with defect Figure 2 - 20 Induced eddy current in the plate We used this ex perimental set - up for our work as shown in Figure 2 - 21 . We started the experiment with an aluminum plat e with a rivet to study and investigate the selection of frequency. Figure 2 - 22 - a show s the aluminum plate with a rivet. This aluminum plate is covered by another 5 mm thick aluminum plate like a second layer , so it will look like a plate with an internal defect. Figure 2 - 22 - b and c show the scan images for the frequencies 1000 Hz and 100 Hz respectively. The testing material (aluminum plate) was scanned in the direction from - 30 mm to 30 mm distance considering the origin is approximately at the defect . There 72 are 32 GMR sensors fixed in the line of symmetry scanning the magnitude of the magnetic flux every mm distance. We found that at 1000 Hz, we get the sc anned image with the le a s t noise and magnitude of flux density with the range of 2 mT 12 mT. Figure 2 - 21 Experimental set - up a. Aluminum plate with rivet b. Scan ned image at 1000 Hz c. Scan ned image at 100 Hz Figure 2 - 22 Scan ned Images for the aluminum plate with rivet 73 Later we moved to a steel plate which has a n internal defect as shown in Figure 2 - 23 - a. In this exp eriment, we used only 3 sensors as our current computational work is limited to 2D. H ere the steel plate was scanned in the direction from - 1 0 mm to 3 0 mm distance considering that the origin is approximately at the defect. From the readings, the midd le GMR sensor readings , which was named S ensor - 2 readings , were filtered out and plotted in the direction as in Figure 2 - 23 - b . a. Steel plate with internal defect b. Reading of Sensor no.2 Figure 2 - 23 Readings for the steel plate with defect Then we worked on our computational model that wou ld match the experimental model. In the experimental work a multi - line coil is considered to be a current sheet. But in our 2D computation al mod el, the current sheet will be a single line and this is impossible to compute using a 2D finite element method. Therefore the current she et is assumed to have thickness and we designed our model as shown Figure 2 - 24 - a. The solution for the magnetic vector potential is shown in Figure 2 - 24 - b. 74 a. 2D Model b. Equi - potential lines for the model Figure 2 - 24 2D model to match the experimental setup The solution (magnitude of magnetic f lux density) of the experimental model and computational model did not match to progress further. There can be two reasons for this mismatching. One reason is that our app roximating the current sheet to be a thick infinite layer is not a good assumption. The second possible reason is that the GMR sensor is moving with the coil in the experimental model, but it is not in the computational model. In a 3D computational model , the experimental setup including the current sheet can be modeled to get better results. 2.6.1 Future Work Matching experimental and computational work in defect charac terization work is in progress as we are developing the 3D mesh generator [62] and validating the solutions. Once we match both in the objective function , will be ta ken from the experimental setup. T hen using the c omputational setup, the defect characterization will be effected . 75 2.7 Conclusion Our goal in this chapter was characterizing the shape of the defect. But, before move into defect characterization, a study on the defect detection has been done. Defects in ma ny different orientations and the detectability of these defects were tested. According to our investigation s when the depth of the defect increases, it is hard to detect. While that may be as to be expected, w e also found that , when we compare horizontal and vertical defects , a horizontal defect has a higher chance of being detected. Our results further showed that when a defect is at an angle of 45 0 , it has the highest chance of being detected. This chapter also presents a finite element technique for s olving inverse problems in magnetostatic NDE [69] . Defect shape reconstructing using the genetic algorithm optimization method is presented and validated using a numerical model. We also imposed constraints in the system to get a realistic single defect reconstruction that is smooth. The e fficiency of the genetic algorithm was checked with real and binary solutions. We co nclude that for our problem, real numbers based genetic algorithms take less computation time. Therefore we decided to work with only real numbers based genetic algorithms in all future problems in this thesis . We also did many tests for the different shap e of the defects and made sure of the accuracy of our inverse algorithm. Experimental work was done to measure the magnetic flux , so in our computational set up, experimental readings ( can be given as input to find the defect characterist ic. Shaping algorithms o n a single physics problem give confidence to work with the two - p hysics problem. So we moved to the electro - thermal coupled problem [43] which is discussed in the chapter 3 . 76 3 FLIP - TEACHING TO TEACH FINITE ELEMENT OPTIMIZATION 3.1 Introduction - Flip Teaching To test o ur algorithms and seek novel use for them , they were at Michigan State University that teaches the finite element method and optimization with true device design [70] . Flip - teaching was introduced to tackle the challenges of limited clas s time. The traditional order of a) delivering theory and next b) develop ing the ancillary programming tools (mesh generators, solvers and equipotential plotter ) is flipped to do real design . Giving students the ancillary tools created the efficiencies to cover finite elements and optimization, usually taking two courses, in the same single course [70] . To do real design, pre - constructed meshes are given to students described by design parameters and ancillary programs ( mesh generator, solvers and equipotential plotters) are also given to them to do the homework and final project. In the last three weeks of the semester, e ach student mastered one optimization method and presented his or her results and taught that speci fic method to others. The problem now is that to teach this process of design we need to teach two courses, one on fi nite element analysis and the other on optimization because in the world of design the one is not too useful without the other. Taking mu ltiple courses needs specialist commitment from the student such as at graduate level. Therefore we need a single course to address the curricular crunch. And within that course we need to teach finite element s and optimization and pack in programming exer cises. We do that by flip teaching. Flip teaching was resorted to as a solution. Flip teaching [71] , in alternative p hraseology is reverse - teaching [72] . By these terms what is broadly meant is the instructor recording lectures on video for the students to watch at home, thereby freeing valuable class time for face - to - face (F2F) discussions between 77 the instructor and class. It flips the order of homework and F2F teaching. While Bergman and Sams [73] and just 2 journal articles in the Web of Sc ience [71], [72] define flip teaching, the popular literature is cl earer in defining flip teaching [16], [53], [74] . words [74] : concepts and assign them as homework. That frees up p recious class time to work directly with students on projects, exercises or problem sets the stuff that students would traditionally do at home. Now instead, of struggling alone, students can do the most important work with a teacher or peers who can hel In the normal two - semester order of developing code, students would write their finite element code which includes the matrix solver in the first course. Thereafter in the optimization part they would first develop ancillary code like a mesh generator as a function call, plotting routines for object functions and then codes pertaining to optimization routines using the finite element code and the ancillary codes. The ancillary codes takes a lot of F2F time with the instructor for students to develop th emselves, and that time does not inculcate much by way of optimization theory. In this course the term flip teaching is extended to the instructor preparing such relevant ancillary computer codes for the students as do not directly apply to the task of opt imization. Students use and familiarize themselves with the ancillary code by themselves and then become free to switch away from the mesh generation concepts not directly relevant to the pedagogy of the course, to focus during F2F teaching on those concep ts directly relevant to the material of the course like the optimization methods and the synthesis process. 78 3.2 Teaching Engineering Optimization: Real Design T o illustrate using a bench mark design problem that is used to test commercial optimization progra ms [75] , Figure 3 - 1 , and Figure 3 - 2 describe the problem at the i nitial design. Figure 3 - 2 gives the minimal boundary value problem based on symmetry that w ould yield a unique solution [20] . The target is a uniform 1 T vertical flux density along the dotted line shown i n Figure 3 - 2 . The design question is this: what should be the shape of the pole face shown in Figure 3 - 2 , to overcome the effects of fringing and make the flux density cons tant along the measuring line ? The initial flux flow of Figure 3 - 3 for the starting configuration shows that fringing diminishes the flux at the corner from its value at the center - line; and to raise it , the reluc tance at the left edge needs to be diminished which is done by reducing the length of the reluctance of air from the corner to the steel above. Figure 3 - 1 Pole - faced to be shaped 79 Figure 3 - 2 Minimal problem using symmetry For optimization through automatic computation of the shape, we would designate measuring points (9 in this case, shown as dots in Figure 3 - 4 ) w here the vertical 1 T flux density required is computed [70] . Thus we define the object function (3 - 1) where the descriptions {p} would involve the heights of the pole - face a t various points going from left to right, the current in the coil and, as necessary, the permeability of the steel. Figure 3 - 3 Fringing in initial geometry 80 Figure 3 - 4 The shaped pole face The teaching of optimization is kept to a minimum genetic algorithm , simulated annealing, [9], [25] with only a week devoted to the lectures by the instructor but a total of three weeks thereafter for a synthesis project on enginee ring optimization making up a total of 4 weeks in the semester for optimization. Figure 3 - 5 Stretched Spring 81 The stretched spring problem of Figure 3 - 5 (with va lues shown) is an ideal problem whose energy - based object function is given by the work stored in the spring minus the work expended by the two forces, (3 - 2) This function is from Vanderplaats who has shown its properties to bring out the stren gths of the various methods [9] . The equal - F lines plotted by us are shown in Figure 3 - 6 . This object function allows us to see how the different methods work with different starting points. Venkataraman [25] gives several programs with his book so that with the MATLAB toolbox on optimization [76] , students have an invaluable aid to problem solving Figure 3 - 6 Object function for spring (on the plane) 3.3 A Useful Tin Problem The trivial tin probl em from a high school calculus book [77] was also used as a tool for students to play with. The exercise is simp ly to determine the radius r and height h of a 82 cylindrical tin of volume V such as to minimize the area A of the tin sheet to be used. Thus the object function becomes, the sum of the areas for the two lids and the curved surface (3 - 3) subject to the equality constraint that volume (3 - 4) This problem may be discussed as a three dimensional optimization problem subject to the equality constraint added to the object function with a Lagrange mult iplier [9, 25] (3 - 5) Alternatively it may be cast as a one - dimensional problem in by substituting for h in A using the constraint to obtain the one - dimensional object function (3 - 6 ) and sol ved by setting: (3 - 7 ) For volume V = 10 cm 3 , we obtain r = 1.1675 cm and h = 2.34 cm. Knowing the solution, the former three - dimensional approach helps students to deal with constraints through pe nalty functions wit h understanding. The theory is [9], [25] that as the Lagrange multiplier gets to infinity the total augmented object function settles down to the minimum of the unconstrained object function, but because of stability issues has to be increased progressively. As a very short assignment students may bypass the stability probl em, understand how works and verify the approach by random search. Using random numbers to generate r and h within a window with a large = 100,000, and searching for lower values of , we get convergence to the correct answers r = 1.1683 cm h = 2.3322 cm or something 83 similar fairly quickly. But with smaller or negative r and h converge to th e highest values of the window. An understanding also emerges from this tin problem of how the shape of the object function determines convergence. When plotted without the constraint, Figure 3 - 7 . a shows long vertical strips around the same h value wh ere solutions are possible with very slight variations in the object function. Figure 3 - 7 . b and Figure 3 - 7 . c show the effect of the constraint confining the solution. The large blank area can numerically accommodate many solutions where the object function varies but slightly with small gradien t, and emphasizes the numerical stability of methods which may be judged knowing the exact solution when students play with the various tools [25, 72] . In this case, the random method works very well. a. b. c. Constraint h 2 only Figure 3 - 7 Object function 3.4 Design of complex devices: The challenge The second difficult class assignment is automatic mesh generation. In finite element optimization, we start with a design {p} and iteratively change it towards the optimum. After 84 {p} changes, there corresponds a new geometry for which a new mesh has to be generated and the problem solved to evaluate the object function again. For practicable, seamless, no nstop optimization, given {p}, the mesh needs to be constructed without stopping the optimization iterations. While many commercial finite element mesh generators exist, they are incapable of being used as a function call that is, given the parameter set {p}, the function must return for the nex t solution the new mesh data [78] . After years of focused research on optimization for electromagnetic device synthesis as well as NDE, there is to date no general purpose mesh generator capable of this reported in the literature although there are commercial codes with this capability [60] . Figure 3 - 8 A simplified 5 - parameter mesh generator Thus a problem - specific (as opposed to general purpose) mesh generator was prepared for the design problem of Figure 3 - 4 . Figure 3 - 8 shows a 5 - parameter mesh generator where a 85 different but equally valid line of symmetry is used. The steel is shown in white and the coil and air in different shades of dark. The mesh is crude but suffices for the problem at hand. T he five parameters are the heights of the white pole - face at x = 0, 2.5, 5.0, 7.5, and 10 as measured up from y = 5.0 (these scales being relative). Mesh description {p} then consists of each height as a component of {p}. Each height is divided into 5 equal pieces. As a height is changed from one iteration to the next, the length above that pole - face going up to cm where the measuring points are) is appropriately adjusted and divided into three equal pieces as seen in Figure 3 - 8 . The mesh generator, handled as a MATLAB function, takes these 5 heights and the permeability of the steel and current in the coil as input and returns the mesh data. 3.5 The Assignments Figure 3 - 9 The final optimum shape 86 - hour course, an assignment was given with a mix of analysis and programming problems. It is really the latter that helped students understand the numerical method. Give n the difficulties and the immense workload, the programming assignment could be done by two students together. The eleventh week was devoted to lectures on optimization. During the last three weeks of class each student was asked to take over one particul ar method of optimization, become familiar with it, develop MATLAB code for it, do the the pole - face shaping with it and present it to the whole class through a 20 minute lecture. The experience helped even the instructor hear of methods (such as the taxi - and used. The optimum, yielding and in the two triangles where the design goal was measured is shown in Figure 3 - 9 obtained by simulated annealing [25] . Gradient methods, random algorithm were pres ented and thoroughly discussed. 3.6 Conclusions Finite elements and mathematical optimization have been taught together in a semester to make up a true course on engineering opt imization involving shape synthesis. Flip teaching was effectively used to give pre - prepared programs dealing with the ancillary, mundane tasks associated with finite element optimization. Students learned by themselves how to use these programs. As a resu lt the more detailed and complex programming assignments were taken up immediately for face - to - face teaching. The end of the course saw each student armed with a suite of programs for engineering optimization. Through this course our algorithms were tested and this give s confidence to work on coupled problem optimization. 87 4 ELECTRO - THERMAL PROBLEM 4.1 Introduction In electromagnetic field problems, optimization methods have been successfully developed and applied efficiently. But , most of these methods only de al with the direct problem. However , real problems are more complex and often coupled, with two or more physical systems interacting. Heavy currents always lead to heating through the joule effect. This heat is often undesirable as in electrical machinery like alternators where the heat not only diminishes the efficiency of the generator but also can damage the insulation [79] . In other cases this heat can be beneficial as in a) the metallurgical industry where the heat is used to melt the ore and mix it through electromechanical forces [3], [80], [81] or b) hyperthermia treatment in oncology where cancerous tissue is burnt off albeit with lower currents, achieving the heating by stronger eddy currents through a higher 1 kHz frequency [4], [82], [83] . Whatever the situation, it is often desirable to accomplish a particular thermal distribution whether to save an alternator from overheating or to accomplish the necessary melting of the ore or to burn cancerous tissue without hurting healthy ti ssue. As shown in Figure 4 - 1 , the design process involves setting the parameters {p} that describe the electro - heat system (consisting of geometric dimensions, currents in magnitude and phase, and material values), solving the eddy cur rent problem for the magnetic vector potential [84] : ( 4 - 1) 88 where is the magnetic permeability, is the electrical conductivity, the electrical field strength and is the externally imposed electric field driving the current [79] (see section 2.3.2 for the origin of this equation) . The frequency is relatively low (50 Hz to 1 kHz) so that the current density has only the conduction term and no displacement term After finding [18] - the finite element computation for finding for the eddy cu rrent problem is described in section 2.3.2.3 - we com pute the joule heating density from ( 4 - 2) and ( 4 - 3) Figure 4 - 1 Finite Element Analysis and Optimization of Coupled Magneto - Thermal Problems. Once we have the heat source distribution q, the second problem of finding the resulting temperature is addressed by solving [85] ( 4 - 4) 89 where is the thermal conductivity [18] . Here we assume that is isotropic ( meaning that hea t flow in all directions is at the equivalent rate) . We also assume is only from Joule heating. Therefore r adiation and convection effects a t the boundaries are neglected. W e avoid more exact details because our purpose here is to establish the feasibil ity of the numerical methods we use . For s olving ( 4 - 4), the electro - thermal finite element computational work is explained in section 4.3 . Since the problem began with defining the parameters of system description {p}, we note that T = T({p}) since the computed T will depend on the values of {p}. When a particular temperature distribution is desired, the problem is one of finding that {p} which will yield T({p}) = ( 4 - 5). This is recognized as inverting ( 4 - 5) to find {p} and therefore it is referred to as the inv erse problem which is now well understood in the literature, particularly when we are dealing with one branch of physics like electromagnetics [10], [11], [13] [16], [86] . In multi - branch , coupled physics problems like the electro - heat problem under discussion, {p} is often defined in the electromagnetic system and F in the t hermal system [18] . Further, when dealing with numerical methods such as the finite element method, T is given at the nodes (although technically T(x,y) may be derived from the finite element trial functions using the nodal values) and , rather than being a function of x and y, is more conveniently defined at measuring points i, numbering say m. The design desideratum then may be cast as an object function F to be minimized with respect to the parameters {p} ( 4 - 6) 90 where is as computed at the same m points i where is defined. The optimization process, by whatever method [9] , keeps adjusting {p} until F is minimized, making come as close to as it can as close as it can rather than exactly to because our design goal as expressed in ( 4 - 6) may not always be realistic and achievable. At that point {p} would represent our best design. The computational process in inverse problem solution as it relates to Figure 4 - 1 is shown in Figure 4 - 2 . It requires solving for the vector of design parameters{p}. We first, generate the mesh from the latest parameter set {p} and get the corresponding finite element solution for . Then from , we compute the rate of heating Q . Thereafter we find the finite element solution for the temperature T. From T we evaluate the object function F . The method of optimization used will dictate how the parameter - set of device description {p} is to be changed depending on the computed F . Figure 4 - 2 The Design Cycle for the Coupled Magneto - Thermal Problem. 91 4.2 Context and Novelty To set this work in context, we note that this two - part optimization problem has been worked on before and we need to state what is novel here. In our work the design vector {p} includes geometric dimensions . This is what is new [43] . That is, we do shape optimization by an accurate finite element eddy current problem followed by another accurate finite element temperature problem and this two part finite element problem is cyclically iterated on by the inherently parallel genetic algorithm. The inherent parallelism of the genetic algorithm was recognized as far back as 1994 by Henderson [87] and done on parallel computers. In 2009 Wong and Wong [88] and Robil liard, Marion - Poty, and Fonlupt [89] took that parallelization work on single field problems to the Graphics Processing Unit, working still with single field systems. We apply their work to the coupled field problem where solutions can take very long and reducing solution times is critical in creating practicably quick engineering design systems without sacrificing accuracy. The coupled electro - heat optimization problem has indeed been tackled before as we found in a Web of Science Core Collection literature search, but with key differences. Pham and Hoole [18] have used two finite element solutions and done shape optimization. But as pointed out below, that process is difficult to build gener al purpose software with because of difficulties with the gradient optimization; that deficiency is being rectified here. Siauve et al . [90] solve eddy current and thermal problems sequentiall y b ut , what they optimize is not shape but the antenna currents to obtain a specific absorption rate (SAR) for hyperthermia treatment. It is not sh ape optimization that they do. Thus the same mesh suffices for every optimization iteration and they do not need to tackle the problem of changing meshes in every iteration and along 9 2 every line search within an iteration as we do [9], [25] . Likewise in the electro - heat optimization work of Nabaei et al. [91] too, the current distribution in furnace transformers is being optimized using an impedance model on MATLAB and there is no shape optimization. Similarly Battistetti et al. [92] avoid a two - part finite element solution by using an analytical solution for the el ectromagnetic part and a finite difference solution for the thermal part. Di Barba et al . [93] similarly use circuit models where they find the axial position of coils, avoiding accurate field comp utation. What we present therefore is very complex, new material. 4.3 Finite Element Computation for the Electro - Thermal Problem Finite element compu tation for the electro - thermal problem is a little more complex than the example thermal problem derived in se ction 1.5.4 . The difference is in section 1.5.4 , where is given. But here we have to find in terms of . This thesis demonstrate s this procedure by the hand calculation below. By combining = ( 4 - 1) and ( 4 - 3), ( 4 - 7) For the thermal problem governed by , ( 4 - 4), the functional will be ( 4 - 8 ) Let u s consider , substituting ( 4 - 7) using the first order interpolation for 93 By applying , ( 4 - 9 ) Now substi tuting ( 4 - 9) i n ( 4 - 8), ( 4 - 10 ) where, in corresponding notation, the local matrices ( 4 - 11 ) a nd ( 4 - 1 2 ) 94 Finite element analysis [20] provides the solution to ( 4 - 4 ) by applying certain boundary conditions. The local matrices of elements will be added to the corresponding positon of the global matrix to be solved for . This leads to the finite element matrix equation ( 4 - 13 ) Equation ( 4 - 13) is solved in the same way as explained in section 1.5.4 , calculating local matrices and adding them to the global matrix t o be solved for . 4.4 GPU Computation for Genetic Algorithms for Electro - heat Problems Naturally the computational work in two - physics electro - heat problems in finite element GA optimization is far beyond that for a single finite element solution. We hav e a 2 - part coupled problem, and we have to solve the magnetic field for A and then the thermal problem where we solve for the temperature. For realistic problems this has to be done several times indeed tens of thousands of times in searching the solut ion space for the minimum object function. Wait times can be excessive, making optimization practicably infeasible. To cut down solution time, parallel processing needs to be resorted to [94] [97] . From the 1990s multiprocessor computers have been tried out. Typically with n proc essors (or computing elements), solution time could be cut down by almost a factor of (n - 1) that is (n - 1) rather than n because one processor is reserved for controlling the other (n - 1), and almost a factor of (n - 1) rather than exactly (n - 1) because of t he additional operations of waiting while one processor accesses the data being changed by another [94] [98] . Although much of this 95 parallelization work was moved to cheaper engineering workstations by the late 1990s [98] , the restrictions on processors remained. Recently the graphics processing unit, endowed with much computing prowess to handle graphics operatio ns, has been exploited to launch a computational kernel as several parallel threads [99] . This is ideally suited for object function evaluation as the kernel so that multiple threads can perform the finite element analyses and evaluation of for each {p} in parallel. (CUDA) [100] are today available on practically every PC as a standard and are increasingly exploited with more and more applications being ported thereto. Significantly the number of parallel threads is not limited as on a shared memory supercomputer. Wong and Wong [88] , Robilliard, Marion - Poty and Fonlupt [89] and Fukuda and Nakatani [101] have shown that the genetic algorithm with its inherently parallel structure may be efficiently implemented on the GPU to optimize magnetic systems. We extend that here to coupled problems. Cecka et al . [102] have also created and analyzed multiple approaches in assembling and solving sparse linear systems on unstructured meshes. The GPU coprocessor using single - precision arithmetic achieves speedups of 30 or so in comparison with a well optimized double - precision single core implementation [15] . We see that this is far better than the factor of just below 7 possible on a very expensive 8 processor supercomputer. So this is the way we will go, using the GPU to process the GA algorithm in parallel. 96 W e fork the fitness value computations as in Figure 4 - 3 . This fitne ss value calculation is the time consuming part as it involves both mesh generation, and finite element calculation. This forms a kernel that will be launched in parallel threads . Therefore we divide the GPU threads and blocks of the same number as the pop ulation size and compute fitness values simultaneously ( Figure 4 - 3 ) . Since we have 65,536 blocks ( and 512 threads in a general GPU, we can go up to a population size of using the one GPU card on a PC. Since we do not need such a large population size for effective optimization, this is not restrictive. All the finite element calculation parts were programmed on the GPU in the CUDA C language. So when, given each {p}, we launch the fitn ess computation kernel as several threads (one for each {p}). The fitness score for all chromosomes is thereupon calculated at the same time in parallel. But for each chromosome, the finite element calculation will be done sequentially ( Figure 4 - 3 ) and not parallelized along the lines of Cecka, Lew, and Darve [102] because that would be attempting to fork within a fork (for the ability to fork an alr e ady forked thread is discussed in a companion paper [103] ) . Figure 4 - 3 The Parallelized Process of the GPU 97 4.5 Test Problem: Shaping an Electro - heated Conductor The test prob lem chosen [43] is a simple one on which the method can be demonstrated and its feasibil ity established. Shown in Figure 4 - 4 a ) is a rectangular conductor which is heated by a current through it. The equi - temperature profiles would be circle - like around the conductor. But we want a constant temperature along two lines paralle l to the pre - shaping rectangular which are to be shaped ( Figure 4 - 4 ). The question is this: how should that edge be reshaped to accomplish a constant temperature along the lines on ei ther side of the conductor? This is the same problem that has been solved by the gradient method which, as noted [18] , needs an alternative computational process because of the difficulties in constructing general purpose software yielding gradient information for the coup led problem and the mesh induced fictitious minima which cause problems [21] and need special mesh generators to address [104] . a) Electrically Heated Conductor: The Actual Geometry Figure 4 - 4 Numerical Model for Coupled Electro - heat Problem 98 Figure 4 - 4 b) Symmetric Quarter: Boundary Value Problem (nominal values shown for magnetic permeability electrica l conductivity and thermal conductivity k) Figure 4 - 4 b ) presents the associated boundary value problem formed from a quarter of the minimal system for analysis consisting of a square conductor (with S/m, W/m/ . A current density has a relatively low frequency of which has been kept deliberately low to avoid a very fine mesh , our purpose here being to investigate and establis h methodology rather than to solve large problems in their full complexity. The top edge of the conductor has to be shaped to get a constant temperature profile of 60 o C, at at 10 equally spaced measuring points in the interval as shown in Figure 4 - 4 along the measuring line . W e define the object 99 function from ( 4 - 6), and for testing we set the desired temperature as 60 o C. Therefore and the will be ( 4 - 14 ) Figure 4 - 5 The Parameterized Geometry An erratic undulating shape with sharp edges arose when Pironneau optimized a pole face to achieve a constant magnetic flux density [12] and this was overcome by the others through constraints [63] . Haslinger and Neittaanmaki [105] suggest Bezier curves to keep the shapes smooth with just a few variables to be optimized, while Preis, Magele and Biro [26] have suggested fourth order polynomials which when we tried gave us smooth but undulating shapes. As such we follow a radiant of the method by Subramaniam et al . [63] and extend their principle, so as to maintain a non - undulating shape by imposing the constraints 100 h 1 >h 2 > h 3 >h 4 >h 5 >h 6 >h 7 ( 4 - 15 ) to ensure a smooth shape. P enalties were imposed by adding a penalty term to the object function F in ( 4 - 14 ) whenever it fails to satisfy the conditions for constraints [9 ] . Tolerance boundaries of each were set to ( 4 - 16 ) The parameterized problem - specific mesh is shown in Figure 4 - 5 where the device descriptive parameter set {p} consists of the 7 heights . The numerical model was uniformly meshed with 234 nodes and 408 elements. This was deliberately kept crude to control debugging; after succeeding with the method and establishing that it works as a method and as programmed on the GPU. In the process of optimization , as these heights change , the mesh connections remain the same but the element sizes and shapes change. For the specific example shown in Figure 4 - 5 , the heights are divid ed into six pieces so the locations of the seven equally spaced seven nodes along the height would change smoothly [104] . Accordingly the length from above the edge of the conductor being shaped to the vertical boundary will be adjus ted and divided into 11 equal lengths as shown in Figure 4 - 5 . 4.6 Results and Discussion Thereof Figure 4 - 6 shows temperature profile at different iteration s . By the images we can see that when the iteration number increases, we get the desired constant temperature along the measuring line and move towards the optimum shape. Figure 4 - 7 shows the optimum shape of 101 the conductor and temper ature profile after 40 iterations for a population size of 512. The corner of the conductor rising toward the line where the constant 60 temperature is desired is as to be expected. For as seen in Figure 4 - 8 [43] (which shows the design goals being accomplished), the constant 60 temperature is perfectly matched. The lower graph giving the initial temperature shows that the temperature drops above the corner of the conductor. Therefore to addr ess this, the corner has to rise close to the line of measurement to heat the line above the corner and Figure 4 - 7 shows that this is what the optimiz ation process has accomplished. Iteration - 1 Iteration - 2 Iteration - 1 0 Iteration - 2 1 Figure 4 - 6 Equi - Temperature Lines at different iterlations 102 Figure 4 - 7 Optimal Shape by the Genetic Algorithm Figure 4 - 8 Temperature Distribution: Desired, Initial and Optimized Significant speedup was accomplished with GPU computation as seen in Figure 4 - 9 . No gains in speedup beyond a factor of 28 were seen after a population size of 500. The meandering 103 nature of the gain after that may be attributed to the happenstance inherent to a statistical method like the genetic algorithm. Figure 4 - 9 Speedup: GA Optimization GPU Time: CPU Time 4.7 Conclusions Shape optimization for the electro - heat problem using GA has been presented and validated using a simple geometry and neglecting radiation and convection. Real - coded object functions give f aster more accurate solutions. This is the first two - stage finite element solution of a magnetic field problem and then a thermal problem repeated in optimization iterations for finding both shape and currents and is amenable to implementation as a general purpose software tool. This avoids the need to use circuit models or analytical solutions to obviate the difficulties of optimization in a two - stage finite element solution process. The procedure provides for shape optimization whereas the extant limited literature on two - stage electro - heat problem optimization shows only the current magnitude and phase being optimized without change of shape. 104 This problem was computed using parallel GPU computing techniques whereby speedups of 28 were demonstrated. This is comparable to the speedup of 30 recently demonstrated in the literature for a single finite element solution. Yet we have demonstrated a speedup of 148 for a single finite element matrix equation solution [62] . A companion paper shows how such an immense speedup was achieved [99] . 105 5 HYPERTHERMIA TREATMENT PLAN N ING 5.1 Introduction In chapter 4 , the electro - thermal problem was explained including how it would be solved and the shape of the current source optimized to get the desired temperature at the location s of interest . In this chapter , this electro - thermal problem is extended on its application side , p articularly to the treatment of tumors using hyperthe rmia . Historically, the treatment of cancer ous tumor s with hyperthermia can be tr aced back to 3000 B. C. when soldering sticks were inserted in tumors. [106] was introduced in the 19 th century. This produced whole body hyperthermia which resulted in tumor regression. In recent years, using hyperthermia or related forms of therapy has increased tremendously. Currently hyperthermia is an experimental treatment and usually applied to late stage patients. There are many heating methods such as whole body heating using wax, hot air, hot water suits, infrared or partial body heating using radiofrequency (RF), microwave, ultrasound and hot blood or fluid perfusion. Clinical and experimental results show a promising future for hyperthermia. However, the main problem is the gene ration and control of heat in tumors. An extremely important issue is to control the temperature distribution in the treated area to avoid excessive temperature s in the normal tissues surrounding the tumor [107] . There are many studies on the treatment of cancer using hyperthermia which demonstrate that this aspect is still important and more research is needed in this matter [108] [111] . Successful applications of electro - thermal stimulation with the aid of a low frequency field to treat the tumors have been reported [107] . There is evidence that hyperthermia can 106 significantly reduce/treat the tumor, but it is not clear how it works and how it should be applied [112] . Some clinical studies have demonstrated the efficiency of thermal therapy in suspending tumor growth [113] . Therefore numerical modeling to distribute the optimal temperature in hyperthermia can be helpful in identifying the better treatments. Hyper thermia treatment involves heating the tumor tissue to a temperature greater than 42 o C without exceeding the normal physiol ogic temperature which is lower than 44 o C 45 o C. Therefore the working temperature margin is very small. If the temperature at the tumor is lower than 42 o C, there is no therapeutic effect. On the other hand, if the temperature is greater than 44 o C 45 o C, both healthy and tumor ous cells are damaged [114] . The blood vessels in tumor cells have a greater diameter than in healthy cells, and therefore occupy a greater volume. The temperature at a tumor is greater than at the surrounding tissue during hyperthermia treatment. This is caused by the fact that the healthy cells ha ve usually greater conductivity than the cancerous cells. In other w ords, we may say tha t tumor tissues are more sensitive to heat. The temperature rise in tumors and tissues is determined by the energy deposited and the physiological responses of the patient. When electromagnetic methods are used, the energy deposition is a complex function of frequency, intensity, the polarization of the applied fields, geometry and size of the app licator, and the geometry and size of the tumor. The final temperature elevations are not only dependent on the energy deposition but also on blood perfusion which carries away heat and affects thermal conditions in tissues. Generally, it is not easy to obtain an accurate temperature distribution over the entire treatment region during clinical hyperthermia treatment . Furthermore, to ensure that the temperature is w ithin the desired range , the clinician usually monitors the temperature every few seconds. Thus it is 107 desirable to develop a mathematical model that can determine design parameters to optimize the temperature distribution in the target region before treatm ent. In this way the treatment efficiency can be assessed more preci sely [115] . 5.2 Problem Statement Figure 5 - 1 Cross section of the human thigh with a tumor Let us consider a c ross section of the human thigh, shown in Figure 5 - 1 [109] . It is assumed that the human thigh and bone inside have an ellipsoidal shape and that the tumor has a circular form. Kurgan and Gas [109] solve d this problem with a numerical model. In their model, the tumor inside the human thigh was heated by external RF hyperthermia with the frequency of 100 MHz . A s imilar model was used by Tsuda and Kuroda [115] . H ere they 108 optimize the electrode conf iguration, particularly driving voltages for radio frequency of 10 MHz using the gradient method. In our work, this model is extended to a low frequency field such as at 1 kHz to 5 kHz . Here we optimize the shape of the current carrying conductor to get the desire d temperature distribution at the tumor. The n umerical model of the problem to be solved is shown in Figure 5 - 2 . [116] . The g eometrical dimensions of the numerical model are described in Table 5 - 1 . The dimensions are taken from Gabriel et al . [117] . Figure 5 - 2 Numerical model of the problem 109 Table 5 - 1 Geometrical dimensions of the model Human bo dy - length of semi axes Tumor radius Skin thickness Bone - length of semi axes Current carry ing conductor Square shaped Near the human thigh a coil with the exciting current is placed. The exciting current in the coil generates a sinusoidal electromagnetic field which induces eddy currents in the human body. These currents act as sources of Joule heat and after some transient time a temperature distribution in the bod y is established. This final, steady state temperature distribution is what we shall seek to solve for. 5.3 Main Equations 5.3.1 The Finite Element Equation As discussed in chapter 2, for a given current , magnetic vector potential will be found by solving the eddy current problem ( 5 - 1) The frequency is relatively low ( 1 k Hz to 5 kHz) so that the current density has only a conduction term and no displacement term . Solving for using the finite element 110 method is explained in section 2.3.2.2 . After finding we compute the joule heating density from ( 5 - 2) and ( 5 - 3) Once we have the heat source distribution q, the second problem is finding the resulting temperature . In chapter 3 , for the electro - thermal problem, we solved the heat equation [85] . ( 5 - 4) But when it comes to electric currents in a human body, no clear consensus exists for an appropriate mathematical model for the evaluation of the temperature field dist ribution in biological tissues. An extremely important work is that by Pennes [118] in modeling the heat transfer in biological tissues. The equation he deri ved is named the bio - heat equation and this can be derived from the c lassical Fourier law of heat conduction. Pennes model is based on the assumption of the energy exchange between the blood vessels and the surrounding tumor tissues. It may provide better information about the temperature distributions in the whole body. The Pennes model states that the total heat exchange between the tissue surrounding a vessel and the blood flowing in it is proportional to the volumetric heat flow and the temperature difference between the blood and the tissue. The expression of the bio - heat equation in a body with uniform material properties in transient analysis is given by [119], [120] . ( 5 - 5) 111 where is the body temperature, is the thermal conductivity, is the density, is the specific heat, is the blood vessel temperature, is the blood density, is the blood specific heat, is the blood perfusion rate, is the heat generation by the external heat source which is responsible for the changing of temperature inside the body according to the e quations ( 5 - 3), ( 5 - 2) a nd ( 5 - 1), and is the metabolic heat generation rate. Metabolic heat production is caused by the chemical factors in the body. The s pecific dynamic action of food is often mentioned , especially protein, that results in a rise of metabolism; and a high environmental temperature that, by raising temperatures of the tissues, increases the velocity of reactions and thus increases heat production . In our problem, we only consider steady state analysis. Therefore the transient part of the bio - heat e quation will be neglected. For the rest of the bio - heat equation, the functional will be ( 5 - 6) In ( 5 - 6 ), the first two parts we have already calculated in the finite element computation for the electro - thermal problem in section 4.3 . Therefore in this section, heat generation due to the metabolic effects and hea t removal because of the blood circulation will be solved. Let us consider heat generation due to the metabolic effects, ( 5 - 7) 112 Let us now consider heat removal because of blood circulation, By applying , ( 5 - 8) By substituting ( 5 - 7) and ( 5 - 8) i n ( 5 - 6), The f ull e quation using (3 - 10) then is, 113 By arranging the same coefficients at one place , ( 5 - 9 ) where, in corresponding notation, the local mat rices ( 5 - 10 ) 114 a nd ( 5 - 1 1 ) Finite element analysis [20] provides the solution to ( 5 - 6 ) by applying certain boundary conditions such as the Dirichlet and Neumann conditions [20] . The local matrices of elements will be added to the corresponding positon of the global matrix to be solved for . This leads to the finite element matrix equation ( 5 - 12 ) Equation ( 5 - 12) is solved in the s ame way as explained in section 1.5.4 , calculating local matrices and adding them to the global matrix equation to be solved for . 5.3.2 A note on electrical and thermal conductivity changes 5.3.2.1 Electrical conductivity In general whe n temperature increases, electrical conductivity reduces. The e lectrical conductivity of most materials change s with temperature . If the temperature does not vary too much, we typically use the approximation [121] . ( 5 - 13) where is the temperature coefficient of resistivity, is the fixed reference temperature (20 o C), and is the electrical conductivity at . Isaac Chang [122] has done a study with temperature - dependent electrical conductivity and constan t electrical conductivity 115 for the electro - thermal coupled problem using fin ite element analysis. His results show that in temperatures below 45 o C, the change in electrical conductivity is less than 10% [122] . In hyperthermia treatment, the temperature does not exceed 43 - 45 o C [123] . Therefore i n our computation, we approximate electrical conductivity as a constant. 5.3.2.2 Thermal conductivity The effect of temperature on thermal conductivity is different for metals and non - metals. In metals, the conductivity is primarily due to free electrons. The t hermal conductivity in metal s is proportional to a multipl e of the absolute temperature and electrical conductivity [124] . ( 5 - 1 4 ) In metals decreases when the temperature increa ses. Thus the product of ( 5 - 1 4 ) , sta ys approximately constant. in nonmetals is mainly due to lattice vibrations ( phonons ). Except for high quality crystals at low temperatures, the phonon mean free path is not reduced s ignificantly at higher temperatures. Thus, the thermal conductivity of nonmetals is approximately constant at low temperatures. At low temperatures well below the Debye temperature , the thermal conductivity decreases, as does the heat capacity. Therefore the thermal conductivity in both metal s and non - metal s is approximate ly constant at low temperatures. S o in our hyperthermia problem too we consider therm al conductivity as a constant. 5.4 The Algorithm for the Inverse Method In this problem we define the shape of the conductor as a design parameter {p} . In our 2D model, we have two conductors, both carrying current in opposite directions to each other. For 116 practical implication purposes, we have the same design parameter s for both conductor s as shown in Figure 5 - 3 . The design desideratum then may be cast as an object function F to be minimized with respect to the par ameters {p} ( 5 - 1 5 ) The o ptimization process, using the genetic algorithm , keeps adjusting {p} until F is minimized, making come as close to as it can as close as it can rather than exactly to because our design goal as expressed in ( 5 - 1 5 ) may not always be realistic and achievable. At that point {p} would represent our best design. In out model, since the tumor need s to be treated, the desired temperature is set to be = 42.5 o C. Figure 5 - 3 Design parameters 117 The computational process of this hyperthermia treatment problem is similar to the problem of electro - thermal pro blem which has been described in Figure 4 - 4 . But, instead of the heat equation, we have to solve the bio - heat equation. 5.5 Simulation Results The m ain purpose of this work is to shape the heating coil, so as to get the treat ment temperatures at tumor s without damaging the healthy cells. In analyz ing with our model, the human body, tumor and conductor are consid ered as homogeneous media with these typical material properties. For simplicity we assumed a constant value for the blood perfusion rate in various biological tissues. Electrical and thermal properties are taken from [117] and presented in Figure 5 - 2 . The p hysical parameters of the blood are given in Table 5 - 3 Table 5 - 2 Physical properties of the numerical model Material 2 Air 1 0.0257 0 Tumor 1 0.07 0.56 0 Muscle 1 0.15 0.22 0 skin 1 0.1 0.22 0 Bone 1 0.7 0.013 0 Coil - 1 1 401 30 Coil - 2 1 401 - 30 118 Table 5 - 3 Physical parameters of blood Tissue Muscle/skin 300 1020 3640 310.15 0.0004 Tumor 480 1020 3640 310.15 0.005 bone 120 1020 3640 310.15 0.00001 Figure 5 - 4 Mesh for the numerical mo del The exciting current density in one conductor which is carrying current i s set to 30 and the other conductor is set to - 30 . The exciting frequency is kept relatively low in the 60 Hz 5 kHz range to avoid a very fine mesh. The param eterized mesh for the initial shape of the numerical model is shown in Figure 5 - 4 . The numerical model is meshed with 1625 nodes 119 and 3195 elements. The r adiation and convection effects at the boundaries are neglect ed, taking the boundary temperature to be the room temperature at 20 and the magnetic vector potential at boundary. The simulation result for the numerical model is shown in Figure 5 - 5 . This shows the equi - temperature lines for the human thigh with tumor. The t emperature is measured on muscle and inside the tumor along the axis where the measuring line is marked ( in green ) . The measuring line is divided into 10 points where readings were taken an d tabulated as in Table 5 - 4 . For hyperthermia treatment, the temperature at a tumor should be greater than 42 o C and should not be greater than 44 o C. Therefore the desired temperature is set to be 42.5 o C and to achi eve that the edges of the conductors are shaped as shown in Figure 5 - 5 . In the optimization process using the genetic algorithm, multiple solutions were created and checked as to whether it is minimizing the obje ctive function ( 5 - 1 5 ). The tolerance boundaries for each were set to for the conductor at the top and for the conductor at the bottom. As such we follow , like follow, Subramaniam et al . [63] and extend their principle, so as to maintain a non - undulating shape by imposing the constraints ( 5 - 1 6 ) to ensure a smooth shape. The pe nalties were imposed by adding a penalty term to the object function in ( 5 - 1 5 ) whenever it fails to satisfy the conditions for constraints. The optimization process using GA was experimented with different population size s and number s of iterations. Figure 5 - 6 shows the optimum shape of the conductors for the population size of 160 and 40 120 iterations which minimize the object function to a value of 0.97. When the conductor is at its optimum shape, the temperat ure at the measuring line is measured and tabulated as the last column of Table 5 - 4 Figure 5 - 5 Equi - Temperature lines for the initial shape Figure 5 - 6 Optimum shape of the conductors to treat tumor 121 Table 5 - 4 Reading at the measuring line Measuring point x(cm) y (cm) Initial Shape Optimum Shape 1 - 11 0 53.92 43.08 2 - 10.75 0 53.37 42.99 3 - 10.5 0 52.70 42.91 4 - 10.25 0 51.92 42.82 5 - 10 0 51.13 42.74 6 - 9.75 0 50.91 42.65 7 - 9.5 0 50.69 42.57 8 - 9.25 0 50.48 42.49 9 - 9 0 50.26 42.41 10 - 8.75 0 49.92 42.33 5.6 Conclusions The t emperatu re inside the tum or for the initial and optimum shapes of the numerical model is plotted with the measuring line points shown in Figure 5 - 7 . For comparison, the desired temperature is also plotted in the figure. The temperature for the optimum shape is in a very good range to treat the tumor as it is not beyond 44 o C and not below 42 o C. It is impossible to keep the temperature everywhere in a tumor the same. Our results show that we can keep it in the temperature range which would tr eat the tumor eff ectively . Theoretical studies of temperature distributions with magnetic induction methods of achieving hyperthermia have been presented. The most appropriate application for the electro - thermal coupled problem is to be in hyperthermia treatment planning. Our numerical model studies show that the proposed inverse algorithm can optimize the heating conditions in hyperthermia 122 treatment. In treatment planning, once the tumor is located, this system synthesizes the coils to give 42 o C - 44 o C at the tumor. It is very difficult to measure the temperature inside the tumor from the medical point of view. Therefore this system is useful to burn the tumor for late stage patients who have no other choice of treatment. Figure 5 - 7 Initial, optimum and desired temperature 123 6 SUMMARY AND CONCLUSIONS The predominant contribution of this thesis is in exploiting the capabilities of the graphics processing unit (GPU) to analyze and synthesize two - physics s ystems, where given the desired design criterion, we determine the geometric shape of the system. We show that GPU parallelization gives us a speedup of 28, an important feature in solving two physics inverse problems in reasonable time. The software syste m and its shape synthesis capabilities have been successfully applied a) to hyperthermia treatment planning; b) to defect characterization in the non - destructive testing of army ground vehicles; and c) to realize curricular efficiencies using flip teaching for a course on finite elements and optimization. Expanding on the above paragraph, the goal of this this thesis has been to solve the electro - thermal coupled problem to g et a desired solution by shape optimi zation with reduced computational time and wi th more accuracy. That has been accomplished. One of the important applications of the electro - thermal problem is in hyperthermia treatment planning; shaping the co nductor to accomplish a desired heat distribution to burn cancerous tissue. Two finite eleme nt problems need to be solved, first for the magnetic fields and the joule heat that the associated eddy currents generate , and then, based on these heat sources, the second finite element problem for heat distribution. This two part problem needs to be it erated on to obtain the desired thermal distribution by optimizing the shape of the current source. From the studies and experiment, we foun d that the genetic algorithm would be a suitable optimization method for our problem. Particularly the solutions rep resented by real numbers give better and faster results. So we decided to use real number based genetic algorithm for optimizing the problem. 124 To avoid the complexities of a two - physics problem, we first develop the shaping algorithms on a single physics problem for nondestructive evaluation. This shape optimizing concept is developed for defect characterization in a steel plate. But, before we move d into defect characterization, a study on de fect detection ha d to be done. Defects can be in many different orientations. Some of these defect orientations make it hard for the defect to be detected, and therefore such defects will be impossible to charac terize by our inverse algorithm . According to our investigation s when the depth of the defect increases, it i s hard to detect. While that may be as to be expected, w e also found that , when we compare horizontal and vertical defects , a horizontal defect has a higher chance of being detected. Our results further showed that when a defect is at an angle of 45 0 to th e surface, it has the highest chance of being detected. After studying defect detectability, we move d into defect characterization. D efect shape reconstructing using the genetic algorithm optimization method has been presented and validated using a numeric al model. We also imposed constraints in the system to get a realistic single defect reconstruction that is smooth. We also did many tests for the different shape of the defects and got reconstruction that is matching more than 85 % of true profile. To test o ur algorithms and elements and optimization with true device design. Flip - teaching was introduced to tackle the challenges of time. The traditional order of a) delivering theory b) p rogramming ancillary tools (mesh generators, solvers) was flipped to do real design . Flip teaching was effectively used to give pre - prepared programs dealing with the ancillary, mundane tasks associated with finite element optimization. Students learned by themselves how to use these programs ahead if class . 125 After the algorithms were understood from use in NDE and the f lip - teaching experience, we successfully developed them for the two - physics system with reduced computational time with the speedup (CPU Time to GPU time ratio) of 28 and increased accuracy established through the problem of shaping a current carrying conductor so as to yield a desired temperature profile along a line . This is the first two - stage finite element solution of a magnetic field problem and then a thermal problem repeated in optimization iterations for finding both shape and currents and is amenable to implementation as a general purpose software tool. Finally we applied the electro - thermal software to hyperthermia treatment plan ning by a numerical model of a human thigh with a tumor treated by current carrying conductors to be shaped to produce the desired temperature at the tumor. Theoretical studies of temperature distributions with magnetic induction methods of achieving h yper thermia have been presented . Numerical model studies show that the proposed inverse algorithm can optimize the heating conditions in hyperthermia treatment. In treatment planning, once the tumor is located, this system synthesizes the coils to give 42 o C - 44 o C at the tumor. An efficient methodology for multi - physics systems has been developed with applications in flip - teaching, NDE and hyperthermia. 126 APPENDI CES 127 Appendix A: First - Order Interpolations Using Triangular Coordinates For two - dimensional analysis we have decided to use the triangle as our basic element shape. Considering the triangle of Figure A - 1 [20] , the symmetry of the triangle results upon our changing to triangular coordinates defined by: (A - 1) (A - 2) - order functions of distance since they vary with the first power of altitude in (A - 1). Figure A - 1 Triangular Coordinates The b est and natural first - order interpolation to use within the triangle is (A - 2) 128 Appendix B: Publications Raised from This T hesis Journals 1. Victor U. Karthik, Sivamayam Sivasuthan, Arunasalam Rahunanthan , Ravi S. more accurate, parallelized inversion for shape optimization in electroheat problems on a graphics processing unit (GPU) with the real - COMPEL - Int. J. Comput. Math. Electr. Electron. Eng. , vol. 34, no. 1, pp. 344 356, Jan. 2015. 2. Depth Limits of Eddy Current Testing for Defects: A Computational Investi gation and Smooth - Manf., vol. 8, no. 2, pp. 450 457, Apr. 2015. 3. - teaching engineering optimization, electr omagnetic product design, and nondestructive evaluation Comput Appl Eng Educ , vol. 23, no. 3, pp. 374 382, May 2015. 4. S. R. H. Hoole, S. Sivasuthan, V. U. Karthik, A. Rahunanthan, R. S. Thyagarajan, and P. ic device optimization: The forking of already parallelized Applied Computational Electromagnetics Society Journal , vol. 29, no. 9, pp. 677 684, 2014 5. S. Ratnajeevan H. Hoole , Victor U. Karthik, S. Sivasuthan , A. Rah unanthan, R. - destructive Evaluation: A Review in Magnetics, and Future Directions in GPU - based, Element - by - Elemen International Journal of Applied Ele ctromagnetics and Mechanics , vol. 47(3), pp. 607 - 627, 2015 6. S. Sivasuthan, V. U. Karthik, A. Rahunanthan, P. Jayakumar, R. S. Thyagarajan, L. - Based, Parameterized Finite Element Mesh for IETE Te chnical Review , vol. 32, no. 2, pp. 94 103, Mar. 2015. 7. S. Sivasuthan, V. U. Karthik, A. Rahunanthan, P. Jayakumar, R. S. Thyagarajan, L. Defect Characterization: Element - by - J Nondestruct Eval , vol. 34, no. 2, pp. 1 7, Mar. 2015. DOI: 10.1007/s10921 - 015 - 0282 - z 8. Sivamayam Sivasuthan, Victor U. Karthik, Arunasalam Rahunanthan, Ravi S. PU Computations Revue roumaine des sciences techniques Série Électrotechnique et Énergétique Vol. 60, No. 3, pp 241 - 251, 2015 . 129 Conferences 1. Karthik, V.U., T. Mathialakan, P. Jayakumar, R.S. T Electromagnetics (ACES), 2015 31st International Review of Progress in, 1 2, 2015. 2. Reconstructing and Classifying Dama ge in a 2D Steel Plate Using Non - Destructive Evaluation (NDE) Methods Rochester, MI, April 2014 . 3. S. Sivasuthan, V. U. Karthik, A. Rahunanthan, P. Jayakumar, Ravi Thyagarajan, Lalita putation: Why Element by Element Conjugate Gradients", 16th Biennial IEEE Conference on Electromagnetic Field Computation (CEFC), Annecy, France, May 25 - 28, 2014. 4. S. Sivasuthan, V. U. Karthik, A. Rahunanthan, P. Jayakumar, Ravi Thyagarajan, Lalita Udpa and - based, Parameterized Mesh Generator Library for Coupled Gradient Design and NDE", 16th Biennial IEEE Conference on Electromagnetic Field Computation (CEFC), Annecy, France, May 25 - 28, 2014. 5. Sivamayam Sivasuthan, Victor U. Karthik, Arunasalam Rahunanthan, Ravi S. Element Method in Optimization: The Forking of Already Parallelized Threads on Graphics Processing Units to Realize Speedup" The 30th Internati onal Review of Progress in Applied Computational Electromagnetics, Jacksonville, Florida, March 23th - 27th, 2014. 6. Victor U. Karthik Genetic Algorithm on NVIDIA GPU Architecture for Synthe 40th Annual Review of Progress in Quantitative Nondestructive Evaluation , edited by Dale E. Chimenti, Leonard J. Bond, and Donald O. Thompson, AIP Conference Proceedings 2014, vol. 1581,pp. 1991 - 1998 , American Institute of Physics, Melville, NY, 2014. 7. S. Sivasuthan, Victor U. Karthik 40th Annual Review of Progress in Quantitative Nondestructive Evaluation , edited by Dale E. Chimen ti, Leonard J. Bond, and Donald O. Thompson, AIP Conference Proceedings 2014, vol. 1581, pp. 1967 - 1974 , American Institute of Physics, Melville, NY, 2014. 8. S. Sivasuthan, V. U. Karthik electrical engineering optimization: Parallelization on graphical processing units realizing high speed - ICPAM - LAE Satellite Conference , Port Moresby, Papua New Guinea, 2013. 9. Victor U. Karthik, S. Sivasuthan, A. Rahunanthan, P. Jayakumar, R. Thyagarajan, S. 130 Ground Vehicle Systems Engineering and Technology Symposium, Troy, MI, August 2013. 10. Lansing, 2012 (Digest) 131 B IBLIOGR APHY 132 BIBLIOGRAPHY [1] IEEE Trans. Magn. , vol. 24, no. 6, pp. 2521 2523, Nov. 1988. [2] IEEE Trans. Ind. Appl. , vol. IA - 9, no. 4, pp. 395 399, Jul. 1973. [3] A. A. Tseng, J. J. Comput. Phys. , vol. 102, no. 1, pp. 1 17, Sep. 1992. [4] ances in external Cancer Treat. Res. , vol. 93, pp. 213 245, 1998. [5] ACM SIGGRAPH 2003 Papers , New York, NY, USA, 2003, pp. 917 924. [6] D. E. Goldberg, Genetic Algorithms in Search, Optimization, and Machine Learning, First edition. Reading, Mass: Addison - Wesley Professional, 1989. [7] tion of inverse IEEE Trans. Magn. , vol. 34, no. 5, pp. 2924 2927, Sep. 1998. [8] The Engineer , pp. 507 5 10, Sep. 1965. [9] G. N. Vanderplaats, Numerical Optimization Techniques for Engineering Design: With Applications , 1st edition. New York: Mcgraw - Hill College, 1984. [10] Comput. Methods Appl. Mech. Eng. , vol. 15, no. 3, pp. 277 308, Sep. 1978. [11] Int. J. Numer. Methods Eng. , vol. 10, no . 4, pp. 747 766, 1976. [12] O. Pironneau, Optimal Shape Design for Elliptic Systems. , New York: Springer - Verlag, 1984. 133 [13] local jacobian derivative for direct finit e element computation of electromagnetic force, J. Magn. Magn. Mater. , vol. 26, no. 1 3, pp. 337 339, Mar. 1982. [14] S. Gitosusastro, J. - L. Coulomb, and J. - calculations and optimization proce IEEE Trans. Magn. , vol. 25, no. 4, pp. 2834 2839, Jul. 1989. [15] S. R. H. Hoole, S. Subramaniam, R. Saldanha, J. - L. Coulomb, and J. - C. Sabonnadiere, materials, IEEE Trans. Magn. , vol. 27, no. 3, pp. 3433 3443, May 1991. [16] S. R. H. Hoole, Ed., Finite Elements, Electromagnetics and Design, First edition. Amsterdam, New York: Elsevier Science, 1995. [17] Jacek Starz COMPEL - Int. J. Comput. Math. Electr. Electron. Eng. , vol. 17, no. 4, pp. 448 459, Aug. 1998. [18] on of coupled magneto - IEEE Trans. Magn. , vol. 31, no. 3, pp. 1988 1991, May 1995. [19] - and three - J. Appl. Phys. , vol. 57 , no. 8, pp. 3850 3852, 1985. [20] S. R. H. Hoole, Computer - aided Analysis and Design of Electromagnetic Devices . New York: Elsevier Science Ltd, 1988. [21] techniques in IEEE Trans. Magn. , vol. 28, no. 4, pp. 1948 1960, Jul. 1992. [22] IEEE Trans. Magn. , vol. 27, no. 5, pp. 4016 4019, Sep. 1991. [ 23] R. L. Haupt and S. E. Haupt, Practical Genetic Algorithms. John Wiley & Sons, Hoboken, New Jersey, 2004. [24] - based optimization algorithms for IEEE Trans. Magn. , vo l. 31, no. 3, pp. 1932 1935, May 1995. [25] P. Venkataraman, Applied Optimization with MATLAB Programming , 2 edition. Hoboken, N.J: Wiley, 2009. [26] electromagnetic d IEEE Trans. Magn. , vol. 26, no. 5, pp. 2181 2183, Sep. 1990. 134 [27] \ iComput. Sci. Eng. Res., Vol 1, pp 97 - 101, May 1996. [28] Eur. J. Oper. Res. , vol. 114, no. 3, pp. 589 601, May 1999. [29] SIAM J. Comput. , vol. 2, no. 2, pp. 88 105, Jun. 1973. [30] J. H. Holland, Adaptation in natural and artificial systems: an introductory analysis with applications to biology, control, and artificial intelligence . Ann Arbor: University of Michigan Press, 1975. [31] Z. Michalewicz, Genetic Algorithms + Data Structures = Evolution Programs . Springer Science & Business Media, 1996. [32] M. Mitchell, An Introduction to Genetic Algorithms , Reprint edition. Cambridge, Mass.: A Bradford Book, 1998. [33] M. Gen and R. Cheng, Genetic Algorithms and Engineering Design, First edition. New York: Wiley - Interscience, 1997. [34] K. Deb, Multi - Objective Optimization Using Evolutionary Algorithms, First edition. Chiches ter; New York: Wiley, 2009. [35] K. Deb, Optimization for Engineering Design: Algorithms and Examples, 2nd ed. PHI, New Delhi, 2012. [36] Adhan A , vol. 24, no. 4, pp. 293 315. [37] D. E. Goldberg and K. Deb, Kaufmann Publishers, pp. 69 93. [38] W. Spears, An Analysis of Multi - Point Crossover . 1991. [39] K. Deb, A. Pratap, genetic algorithm: NSGA - IEEE Trans. Evol. Comput. , vol. 6, no. 2, pp. 182 197, Apr. 2002. [40] - coded Genetic Algorithms with Simulated Binary Crosso vol. 9, no. 6, pp. 431 454, 1995. [41] -- Coded Genetic Algorithms and Interval - Morgan Kaufmann, 1992. 135 [42] Available: http://www.egr.msu.edu/~kdeb/codes.shtml. [Accessed: 23 - Feb - 2015]. [43] Victor U. Karthik, Sivamayam Sivasuthan, Arunasalam Rahunanthan, Ravi S. Thyagarajan, Paramsothy Jayakumar, Lalita Udpa, and S more accurate, parallelized inversion for shape optimization in electroheat problems on a graphics processing unit (GPU) with the real - COMPEL - Int. J. Comput. Math. Electr. Electron. Eng. , vol. 34 , no. 1, pp. 344 356, Jan. 2015. [44] P. R. P. Hoole, K. Piraphaharan, and S. R. H. Hoole, Electromagnetics engineering handbook: analysis and design of electrical and electronic devices and systems . Southamption: WIT Press, 2013. [45] BOUNDARY ELEMENTS: T heory and Applications. By: John T. Katsikadelis, Elsevier B. V, 2002. . [46] J. García - Martín, J. Gómez - Gil, and E. Vázquez - - Destructive Techniques Sensors , vol. 11, no. 3, pp. 2525 2565, Feb. 2011. [47] M. Ala - Insight - Non - Destr. Test. Cond. Monit. , vol. 54, no. 2, pp. 72 75, Feb. 2012. [48] rence in the detection limits of flaws in the depths of multi - layered and continuous aluminum plates using low - NDT E Int. , vol. 41, no. 2, pp. 108 111, Mar. 2008. [49] ction of Deep Corrosion using \ iProc. Seventh Joint DoD/FAA/NASA Conference on Aging Aicraft, New Orleans, 2003. [50] flux leakage pipeline Int. J. Appl. Electromagn. Mech. , vol. 25, no. 1 4, pp. 357 362, 2007. [51] - Przeglad Elektrotechniczny , vol. 86, no. 1, pp. 249 254, Jan. 201 0. [52] Eddy - Current Based Giant Magnetoresistive System for the Inspection of Aircraft IEEE Trans. Magn. , vol. 46, no. 3, pp. 910 917, Mar. 2010. [53] Z. Liu, Meas. Sci. Technol. , vol. 19, no. 8, p. 085701, Aug. 2008. 136 [54] lement study of the remote field IEEE Trans. Magn. , vol. 24, no. 1, pp. 435 438, Jan. 1988. [55] R. L. Stoll, The analysis of eddy currents . Oxford: Clarendon Press, 1974. [56] Scribd . [Online]. Available: https://www.scribd.com/doc/48834808/An - Introduction - to - Eddy - Current - Theory - and - technology. [Accessed: 19 - Feb - 2015]. [57] Depth Limits of Eddy Current Testing for Defects: A Computational Investigation and Smooth - SAE Int. J. Mater. Manuf. , vol. 8, no. 2, pp. 450 457, Apr. 2015. [58] J. C. Maxwell and Physics, Treatise on Electricity and Magnetism, Vol. 1 , 3rd edition. New York: Dover Publications, 1954. [59] Philos. Trans. R. Soc. Lond. , vol. 155, pp. 459 512, Jan. 1865. [60] S. Sivasuthan, V. U. Karthik, A. Rahunanthan , P. Jayakumar, R. S. Thyagarajan, L. - Based, Parameterized Finite Element Mesh for Design IETE Tech. Rev. , vol. 0, no. 0, pp. 1 10, 2014. [61] enerator and Delaunay Towards Geometric Engineering, Springer Berlin Heidelberg, 1996, pp. 203 222. [62] PhD Thesis, Michigan State University, 2015. [63] IEEE Trans. Magn. , vol. 30, no. 5, pp. 345 5 3458, Sep. 1994. [64] - based eddy - IEEE Trans. Magn. , vol. 37, no. 5, pp. 3831 3838, Sep. 2001. [65] d Proc. SPIE - Int. Soc. Opt. Eng. , 2005. [66] 927 930. [67] w frequency EC - GMR detection of cracks at ferromagnetic fastener 487, Denver, Colorado, 2013. 137 [68] - GMR array with rotating current excitation for AIP Conference Proceedings , 2015, vol. 1650, pp. 368 376. [69] Depth Limits of Eddy Current Testing for Defects: A Comput ational Investigation and Smooth - Warrendale, PA, SAE Technical Paper 2015 - 01 - 0595, Apr. 2015. [70] - teachin g engineering optimization, electromagnetic product design, and nondestructive evaluation in a \ iComput. Appl. Eng. Educ., , Jun. 2014. [71] Int. J. Emerg. Technol. Learn. IJET , vol. 5, no. 2, pp. 45 48, Jun. 2010. [72] Arquitetura Rev , vol. 8, pp. 76 87., 2012. [73] Jonathan Bergmann and /resources/product?ID=2285. [Accessed: 24 - Jul - 2015]. [74] Opinionator . [Online]. Available: http://opinionator.blogs.nytimes.com/2013/ 10/23/in - flipped - classrooms - a - method - for - mastery/. [Accessed: 24 - Jul - 2015]. [75] J. Mater. Process. Technol. , vol. 181, no. 1 3, pp. 136 141, Jan. 2007. [76] T. A. Davis, MATLAB Primer, Eighth Edition, Taylor & Francis, 2010 . [77] e: http://www.slader.com/textbook/9780534393397 - stewart - calculus - 5th - edition/. [Accessed: 24 - Jul - 2015]. [78] IEEE Trans. Educ. , vol. E - 29, no. 1, pp. 21 26, Feb. 19 86. [79] M. V. . Chari and S. J. Salon, Numerical Methods in Electromagnetism . San Diego, CA: Academic Press, 2000. [80] IEEE Tran s. Magn. , vol. 27, no. 5, pp. 4004 4007, Sep. 1991. [81] IEEE Trans. Magn. , vol. 19, no. 6, pp. 2566 2572, Nov. 1983. 138 [82] Hartmut Brauer, Marek Ziolkowski, Uwe Tenner, Jens Haueisen, and Hannes Nowak, COMPEL - Int. J. Comput. Math. Electr. Electron. Eng. , vol. 20, no. 2, pp. 595 606, Jun. 2001. [83] on of hyperthermia - SAR distributions in 3 - IEEE Trans. Magn. , vol. 26, no. 2, pp. 1011 1014, Mar. 1990. [84] J. Starzynski and S. Wincenciak, Benchmark Problems for Optimal Shape Design for 2D Eddy Currents . 1998. [85] Iliana Marinova and Valentin Matee COMPEL - Int. J. Comput. Math. Electr. Electron. Eng. , vol. 31, no. 3, pp. 996 1006, May 2012. [86] COMPEL - Int. J. Comput. Math. Electr. Electron. Eng. , vol. 22, no. 3, pp. 535 548, Sep. 2003. [87] Comput. Syst. Eng. , vol. 5, no. 4 6, pp. 441 453, Aug. 1994. [88] M. L. Wong and Sarker, and B. - T. Zhang, Eds.Intelligent and Evolutionary Systems, Springer Berlin Heidelberg, 2009, pp. 197 216. [89] D. Robilliard, V. Marion - Genet. Program. Evolvable Mach. , vol. 10, no. 4, pp. 447 471, Dec. 2009. [90] N. Siauve, L. Nicolas, C. Vollaire, A. Nicolas, and J. A. Vasconc 3 - IEEE Trans. Magn. , vol. 40, no. 2, pp. 1264 1267, Mar. 2004. [91] Distribution in Parallel Windings of Furna IEEE Trans. Magn. , vol. 46, no. 2, pp. 626 629, Feb. 2010. [92] design of an inductor for transverse flux heating using a co mbined evolutionary simplex COMPEL - Int. J. Comput. Math. Electr. Electron. Eng. , vol. 20, no. 2, pp. 507 522, Jun. 2001. [93] systems for induction heating: m COMPEL - Int. J. Comput. Math. Electr. Electron. Eng. , vol. 22, no. 1, pp. 111 122, Mar. 2003. 139 [94] simulation and possibilities for extension Int. J. Appl. Electromagn. Mater. , vol. 2, no. 1, pp. 99 108, 1991. [95] IEEE Trans. Magn. , vol. 26, no. 2, pp. 837 840, Mar. 1990. [96] S. - IEEE Trans. Magn. , vol. 26, no. 4, pp. 1252 1255, Jul. 1990. [97] fi J. Appl. Phys. , vol. 67, no. 9, pp. 5818 5820, May 1990. [98] IEEE Trans. Magn. , vol. 33, no. 2 , pp. 1966 1969, Mar. 1997. [99] S. R. H. Hoole, V. U. Karthik, S. Sivasuthan, A. Rahunanthan, R. S. Thyagarajan, and P. in magnetics, and future directions in GPU - ba sed, element - by - element coupled optimization Int. J. Appl. Electromagn. Mech. [100] [Online]. Available: http://www.nvidia.com/object/cuda_home_new.html. [Accessed: 23 - Feb - 20 15]. [101] Co - Optimization Using Genetic Algorithm and GPU - IEEE Trans. Magn. , vol. 48, no. 11, pp. 3895 3898, Nov. 2012. [102] C. Cecka, A. J. Lew, and E. Dar Int. J. Numer. Methods Eng. , vol. 85, no. 5, pp. 640 669, 2011. [103] S. R. H. Hoole, S. Sivasuthan, V. U. Karthik, A. Rahunanthan, R. S. Thyagarajan, and P. Optimization: The Forking of Already Parallelized pp. 677 684, 2014 [104] Parametric Geo \ iJ. Mater. Process. Technol., vol. 161, pp. 368 373, 2005. [105] J. Haslinger and P. Neittaanmäki, Finite Element Approximation for Optimal Shape, Material and Topology Design , 2 edition. C 140 [106] - Aug - 2012. [Online]. Available: http://www.cancerresearchuk.org/about - cancer/cancers - in - general/cancer - questions/coleys - toxins - cancer - treatment. [Accessed: 27 - Jun - 2015]. [107] C. - IEEE Trans. Instrum. Meas. , vol. 37, no. 4, pp. 547 551, Dec. 1988. [108] fizyczne [109] Przeglad Elektrotechniczny , vol. 87, no. 12b, pp. 103 106, 2011. [110] ues in External RF Przeglad Elektrotechniczny , vol. 86, no. 01, pp. 100 102, 2010. [111] temperature distributions produced by interstitial arrays of sleeved - sl ot antennas for IEEE Trans. Microw. Theory Tech. , vol. 51, no. 12, pp. 2418 2426, Dec. 2003. [112] B. Stea, K. Rossman, J. Kittelson, A. Shetter, A. Hamilton, and J. R. Cassady, ermoradiotherapy for supratentorial malignant Int. J. Radiat. Oncol. Biol. Phys. , vol. 30, no. 3, pp. 591 600, Oct. 1994. [113] i Int. J. Hyperthermia , vol. 7, no. 2, pp. 317 341, Jan. 1991. [114] optimization for multiple source applicat Int. J. Hyperth. Off. J. Eur. Soc. Hyperthermic Oncol. North Am. Hyperth. Group , vol. 15, no. 4, pp. 291 308, Aug. 1999. [115] in RF - IEEE Trans. Biomed. Eng. , vol. 43, no. 10, pp. 1029 1037, Oct. 1996. [116] - dimensional computer simulation in deep regional hyperthermia using the finite - difference time - IEEE Trans. Microw. Theory Tech. , vol. 38, no. 2, pp. 2 04 211, Feb. 1990. [117] Phys. Med. Biol. , vol. 41, no. 11, p. 2251, Nov. 1996. [118] of Tissue and Arterial Blood Temperatures in the Resting J. Appl. Physiol. , vol. 1, no. 2, pp. 93 122, Aug. 1948. 141 [119] H. Kroeze, J. B. van de Kamer, A. a. C. de Leeuw, M. Kikuchi, and J. J. W. Lagendijk, Int. J. Hyperth. Off. J. Eur. Soc. Hyperthermic Oncol. North Am. Hyperth. Group , vol. 19, no. 1, pp. 58 73, Feb. 2003. [120] temperature regulatory system -- IEEE Trans. Biomed. Eng. , vol. 23, no. 6, pp. 434 444, Nov. 1976. [121] M. R. Ward, Electrical Engineering Science . McGraw - Hill, 1971. [122] obes using Temperature - Biomed. Eng. OnLine , vol. 2, p. 12, May 2003. [123] Int. J. Hyperth. Off. J. Eur. Soc. Hypert hermic Oncol. North Am. Hyperth. Group , vol. 4, no. 1, pp. 17 23, Feb. 1988. [124] W. Jones and N. H. March, Theoretical Solid State Physics . Courier Corporation, 1985.