Universiteit Gent Faculteit Ingenieurswetenschappen. Promotor: Prof. dr. ir. L. DUPRE Begeleider: ir. G. CREVECOEUR

Universiteit Gent Faculteit Ingenieurswetenschappen

Vakgroep Elektrische Energie, Systemen en Automatisering Voorzitter: Prof. dr. ir. J. MELKEBEEK

‘Space Mapping’-optimalisatie van het magnetisch circuit van elektrische machines met inbegrip van lokale materiaaldegradatie door Stijn UITTERHAEGEN

Promotor: Prof. dr. ir. L. DUPRE Begeleider: ir. G. CREVECOEUR

Scriptie ingediend tot het behalen van de academische graad van burgerlijk werktuigkundig-elektrotechnisch ingenieur

Academiejaar 2005–2006

Universiteit Gent Faculteit Ingenieurswetenschappen



Promotor: Prof. dr. ir. L. DUPRE Begeleider: ir. G. CREVECOEUR



Voorwoord Het is intussen vele maanden geleden dat ik aan het maken van deze masterproef begonnen ben. Het zijn zeer leerrijke maanden geworden, waarin fantastische successen en vreselijke tegenslagen elkaar afwisselden. De voorliggende tekst is het resultaat van deze onderneming. Deze tekst had er totaal anders uitgezien zonder de bijdrages van een aantal personen. Daarom wil ik hier mijn dankbaarheid uitspreken voor de steun die ik de voorbije maanden heb genoten. In de eerste plaats wil ik mijn dank betuigen aan de promotor van deze scriptie, Prof. Luc Dupré. De eerste voorwaarde om een project als deze masterproef tot een goed einde te brengen is een welomlijnde, relevante en uitdagende opdrachtsbeschrijving. Hieraan was vanaf de eerste dag voldaan. Op geen enkel moment was er twijfel over wat er van me verwacht werd. Steeds kon ik met problemen bij hem terecht en geen enkele vraag bleef onbeantwoord. Ook de voorzitter van de vakgroep Elektrische Energie, Systemen en Automatisering, Prof. Jan Melkebeek, wil ik danken. Zonder zijn ondersteuning had deze tekst niet tot stand kunnen komen. Mijn dank gaat in het bijzonder ook uit naar ir. Guillaume Crevecoeur, die vele uren heeft opgeofferd om me in te leiden in de complexe materie van de optimalisatie-algoritmes, die zijn expertise in space mapping heeft ingezet en deze geavanceerde methode tot het geven van oplossingen heeft weten te overtuigen, die me steeds met woord en daad heeft bijgestaan om alle praktische problemen te overwinnen en die me steeds weer met prikkelende vragen en diepgaande overwegingen tot nieuwe inzichten heeft gebracht. Verder wil ik ir. Lode Vandenbossche bedanken. Toen ik in oktober met goede moed de koe bij de horens vatte en aan het opmeten van materiaaleigenschappen begon, hebben zijn aanwijzingen me behoed voor alle vergissingen en problemen die gewoonlijk inherent zijn aan het begin van een masterproef. De interessante gesprekken waarin hij zijn uitgebreide kennis in verband met magnetische materialen met me deelde hebben me inzichten opgeleverd die gedurende het hele jaar hun nut hebben bewezen. Ook dr. ir. Peter Sergeant verdient een dankwoord voor zijn waardevolle hulp bij het uitvoeren van de eindige-elementenberekeningen. Ten slotte wil ik alle professoren en lesgevers bedanken die me de voorbije jaren een uitstekende opleiding hebben verschaft en het me zodoende mogelijk gemaakt hebben deze opdracht tot een goed einde te brengen.

De auteur geeft de toelating deze scriptie voor consultatie beschikbaar te stellen en delen van de scriptie te kopiëren voor persoonlijk gebruik. Elk ander gebruik valt onder de beperkingen van het auteursrecht, in het bijzonder met betrekking tot de verplichting de bron uitdrukkelijk te vermelden bij het aanhalen van resultaten uit deze scriptie. Stijn Uitterhaegen 30 mei 2006 i




Promotor: Prof. dr. ir. L. DUPRE Begeleider: ir. G. CREVECOEUR Faculteit Ingenieurswetenschappen Universiteit Gent


Samenvatting Productieprocessen zoals snijden, uitgevoerd op laminaties van elektrisch staal, hebben lokaal invloed op de magnetische eigenschappen van deze laminaties. Aangezien het ontwerp van elektrische machines dit effect niet nauwkeurig in rekening neemt, kan het resultaat, of de productie, suboptimaal zijn. Daarom is er nood aan een procedure die elektrische machines optimaliseert met een grote nauwkeurigheid en een aanvaardbare rekentijd, rekening houdend met de lokale materiaaldegradatie. Wij hebben een procedure ontworpen die aan deze vereisten beantwoordt en we hebben de invloed van de materiaaldegradatie op het ontwerp van een geschakelde reluctantiemotor gekwantificeerd. Deze tekst stelt het gebruik van de spacemappingtechniek voor om een elektrische machine te optimaliseren, evenals de constructie van een ruw, semi-analytisch en een fijn, numeriek model van de geschakelde reluctantiemotor en de optimalisatie van de motor, rekening houdend met de lokale materiaaldegradatie. We bespreken ook het gebruik van een experimentele opstelling om de invloed van de machineproductie op de magnetische eigenschappen van het elektrisch staal te bepalen. We stellen vast dat de space-mappingtechniek, gebruik makend van een hybride methode en trust-regionstrategieën, kan gebruikt worden om een snelle en nauwkeurige optimalisatie van complexe systemen uit te voeren. Uit de resultaten van de oplossing van het inverse probleem blijkt dat het mogelijk is de invloed van de lokale materiaaldegradatie op het ontwerp van de machine te kwantificeren. Vooral de optimale afmetingen van stator- en rotorjuk ondervinden een relatief grote invloed van deze materiaaldegradatie.

Trefwoorden: optimalisatie, space mapping, elektromagnetische berekeningen, reluctantiemotor, materiaaldegradatie

Space mapping optimization of the magnetic circuit of electrical machines including local material degradation Stijn Uitterhaegen Supervisor(s): Luc Dupré, Guillaume Crevecoeur Abstract— Production processes like cutting, performed on electrical steel laminations, influence their magnetic properties locally. Since the design of electrical machines does not take this effect into account accurately, the result or the production may be suboptimal. The need exists therefore to develop a numerical procedure which is capable of optimizing electrical devices, taking into account the local material degradation and featuring high accuracy and acceptable calculation time. We developed a procedure meeting these requirements and we quantified the influence of the material degradation on the design of a switched reluctance motor (SRM). This paper presents the use of the space mapping technique to optimize an electrical machine, the construction of a coarse, semi-analytical and a fine, numerical model of the SRM and the optimization of the SRM taking into account the local material degradation. We observed that the space mapping technique, using a hybrid method and trust region strategies, can be used to perform an accurate and fast optimization. The results of the optimization show that the local material degradation has a significant influence on the dimensions of stator and rotor yoke. As a result, it is possible to optimize the design and production process of electrical machines. Keywords— optimization, space mapping, electromagnetic calculations, switched reluctance motor, material degradation

I. I NTRODUCTION

E

LECTRICAL machine laminations are shaped by cutting electrical steel plates. Frequently used cutting techniques are laser cutting, guillotine cutting and punching. However, the application of such techniques on electrical steel introduces significant local strains near the material border, which influence its magnetic properties drastically [1]. Producers take this effect into account by introducing an empirical factor in the machine design or by annealing the shaped components. Obviously, this has an impact on the optimality of the result. Currently used procedures to optimize electrical machines are either based on analytical models, which lead to a fast calculation but a low accuracy. On the other hand, when optimizing using numerical models, a high accuracy can be achieved at the cost of a high calculation time. These disadvantages are strengthened by the added complexity of the models when the local material degradation is taken into account. Therefore the need exists to develop an optimization procedure which guarantees a high accuracy but at a low calculation cost. This paper presents the use of the space mapping technique, which combines the fast computations of the coarse model with the accuracy of the fine model, in order to optimize the magnetic circuit of the SRM. The fine model uses finite element calculations and the coarse model solves the magnetic circuit in a semianalytical way. Both models take into account the local change of material properties. S. Uitterhaegen, L. Dupré and G. Crevecoeur are with the Department Electrical Energy, Systems and Automation, Ghent University (UGent)

II. S PACE M APPING The space mapping (SM) technique is based on the optimization using surrogate models by exploiting a coarse model to align with the computationally intensive fine model [2], [3]. We denote xc the geometrical parameters of the coarse model, c(xc ) the torque profile resulting from the evaluation of the coarse model. xf and f (xf ) are defined identically, applied on the fine model. The cost function to be minimized can strive to a maximized torque or can take the torque ripple into account. The SM technique starts by calculating the optimal value x∗c in the coarse model [2]. A parameter mapping p is then constructed, which approximates the fine model by a surrogate model based on the coarse model: f (xf ) ≈ c(p(xf )).

(1)

In this way, the optimum of the fine model is expected to be close to the optimum of c(p(xf )). The mapping p is determined for a certain value xf by the so-called parameter extraction: p(xf ) = arg min kc(p) − f (xf )k p∈Ωc

(2)

Ωc is the parameter space in the coarse model and thus limits the result to physically meaningful parameters. Every parameter extraction requires one optimization performed in the coarse model and one evaluation of the fine model. We implemented the aggressive space mapping (ASM) technique. ASM solves p(xf ) = x∗c for xf by using quasi-Newton iterations, using the classical Broyden formula. We extended the ASM algorithm with a hybrid method, based on [3], which guarantees convergence by switching to a direct optimization in the fine model when the ASM procedure does not converge. III. M ODELLING THE PROBLEM :

FORWARD PROBLEM

In order to apply SM, a coarse and a fine forward model were constructed which modelled the SRM including local material degradation. The fine and coarse model are presented in figure 1. The variable parameters xc , xf of the models are the width of the stator teeth bst and the rotor tooth brt , the internal diameter of the stator yoke Ds and the external diameter of the rotor yoke Dr . The material under consideration is high-quality, nonoriented electrical steel (JIS 35H230). We used a single sheet tester in order to determine the influence of plastic strains on the B-H-characteristics of the material and we fitted the measured

(a)

(a)

(b)

Fig. 1. Fine (a) and coarse (b) model

data with analytical formulas, fixing three parameters H0 , B0 and ν through the use of ν B B H = H0 + (3) B0 B0 Through the use of data, which characterize the change of the plastic strain as a function of the distance from the cutting edge, the width of the degraded borders and the B-H-characteristics in the subregions were selected. The coarse model is constructed in such a way that every geometrical part of the magnetic circuit is substituted by a certain reluctance. The number of windings w is dependent on the geometry of the SRM. The model can be solved in a fast, analytical way. Since the model is non-linear, we use a damped, iterative procedure. In every iteration, the permeability of every material sector corresponding to a certain reluctance in the magnetic circuit is calculated using the B-H-characteristics and the magnetic induction from the previous iteration. The magnetic flux is then determined in every reluctance by solving the magnetic circuit, using the permeabilities of the material in order to find the values of the reluctances, and the magnetic induction. Every step ∆B is damped by multiplication with α: π(∆B)a (4) α = tanh ∆B with (∆B)a the maximum allowed step to pass unchanged. This maximum step is decreased in the course of the algorithm. The fine model is implemented using finite element calculations. The degraded material borders are modeled as five discrete zones with different material properties. IV. O PTIMIZATION OF THE MAGNETIC CIRCUIT We have optimized the magnetic circuit of the SRM using the space mapping technique for enforcement currents I of 1 A, 2 A, 3 A and 4 A. We have performed this optimization with and without taking the local material degradation into account. The variable parameters are the stator tooth width bst , the rotor tooth width brt , the internal diameter of the stator yoke Ds and the external diameter of the rotor yoke Dr . The resulting geometrical parameters are presented in figure 2(a). We denote the cost of the evaluation of the torque with degradation and optimized with degradation as C (d) , the cost without degradation (and optimized as such) as C (nd) and the cost of the evaluation

(b) Fig. 2. Results of optimization. (a): Difference ∆x between the optimal parameter values with and without material degradation as a function of current I. (b): Torque profiles with optimal geometries at I = 1 A, 2 A, 3 A, 4 A. Blue: T (nd) , red: T (d) , purple: T (no)

with degradation but optimized without taking the degradation into account as C no . The same suffixes are valid for the torque profiles, presented in figure 2(b). We calculated the relative improvement α due to the adapted optimization for every current: α=

C (no) − C (d) C (no) − C (nd)

(5)

The values of α are ranging from 7% to 16%. The computational time needed when optimizing using SM is about 10-15 times faster when using traditional optimization techniques. V. C ONCLUSION We have developed a method that optimizes the geometrical parameters of electrical machines, using the space mapping technique. We also constructed a fine, numerical model and a coarse, semi-analytical model in order to solve the forward problem. In this way, we have succeeded to perform accurate and fast optimizations of the SRM under consideration. Both models are designed to implement local material degradation. The magnetic properties of the material have been determined experimentally. We have quantified the influence of the material degradation on the optimization of the electrical machine. We have observed a significant influence of the material degradation on the dimensions of the stator and rotor yoke. Using this method, it is possible to take the material degradation due to production processes into account for optimizing the machine. This feature can lead to a more efficient design or production process. R EFERENCES [1] V. Maurel, F. Ossart, and R. Billardon, “Residual stresses in punched laminations: Phenomenological analysis and influence on the magnetic behavior of electrical steels,” Journal of Applied Physics, vol. 93, pp. 7106–7108, 2005. [2] J.W. Bandler, Q.S. Cheng, S.A. Dakroury, A.S. Mohamed, M.H. Bakr, K Madsen, and J Søndergaard, “Space mapping: The state of the art,” IEEE Transactions on Microwave Theory and Techniques, vol. 52, pp. 337–361, 2004. [3] G. Crevecoeur, H. Hallez, P. Van Hese, Y. D’Asseler, L. Dupré, and R. Van De Walle, “EEG source analysis using space mapping techniques,” Journal of Computational and Applied Mathematics, 2006, To be published.

Inhoudsopgave Extended abstract

ii

1 Inleiding

1

2 Materiaaleigenschappen

4

2.1

Inleiding . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

2.2

Metingen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

2.2.1

Werkwijze . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

2.2.2

Mechanische vervorming . . . . . . . . . . . . . . . . . . . . . . . . . .

6

2.2.3

Elektromagnetisch meetcircuit . . . . . . . . . . . . . . . . . . . . . .

8

2.3

Fitten van de gemeten materiaalkarakteristieken . . . . . . . . . . . . . . . .

11

2.4

Ruimtelijke verdeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

13

2.4.1

Algemeen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

13

2.4.2

Verdeling van de materiaalkarakteristieken voor het fijne model . . . .

13

2.4.3

Verdeling van de materiaalkarakteristieken voor het ruwe model . . . .

15

Conclusie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

17

2.5

3 Ruw model

19

3.1

Inleiding . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

19

3.2

Magnetisch netwerk . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

19

3.2.1

Algemeen model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

19

3.2.2

Vervangreluctanties

. . . . . . . . . . . . . . . . . . . . . . . . . . . .

20

3.2.3

Bekrachtiging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

26

3.3

Fluxberekening . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

28

3.4

Koppelberekening . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

31

3.5

Resultaten . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

34

3.6

Conclusie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

34

v

INHOUDSOPGAVE

vi

4 Fijn model

35

4.1

Inleiding . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

35

4.2

Geometrie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

36

4.3

Meshing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

37

4.4

Fysische eigenschappen

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

40

4.5

Veldberekeningen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

42

4.5.1

Differentiaalvergelijkingen . . . . . . . . . . . . . . . . . . . . . . . . .

42

4.5.2

Variationele vorm

. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

44

4.5.3

Discretisatie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

45

4.5.4

Niet-lineaire solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

47

4.5.5

Lineair subprobleem . . . . . . . . . . . . . . . . . . . . . . . . . . . .

48

4.6

Koppelberekening . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

50

4.7

Resultaten en conclusie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

51

5 Directe optimalisatie

52

5.1

Inleiding . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

52

5.2

Implementatie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

53

5.2.1

Minimalisatie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

53

5.2.2

Kostenfunctie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

54

5.3

5.4

Algoritmes

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

55

5.3.1

Meerdere variabelen . . . . . . . . . . . . . . . . . . . . . . . . . . . .

55

5.3.2

Eén variabele . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

62

Conclusie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

65

6 Space Mapping

67

6.1

Inleiding . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

67

6.2

Space mapping: principe . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

68

6.3

Aggressive Space Mapping . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

69

6.4

Trust Region . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

72

6.5

Hybride Methode . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

73

6.6

Conclusie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

76

7 Optimaal ontwerp van machines

77

7.1

Inleiding . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

77

7.2

Geldigheid en belang van de resultaten . . . . . . . . . . . . . . . . . . . . . .

77

7.3

Ontwerp in ruw model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

79

7.3.1

Methode . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

79

7.3.2

Optimalisatie naar maximaal koppel . . . . . . . . . . . . . . . . . . .

79

7.3.3

Optimalisatie met inbegrip van koppelrimpel . . . . . . . . . . . . . .

81

Ontwerp met behulp van space mapping . . . . . . . . . . . . . . . . . . . . .

87

7.4

INHOUDSOPGAVE

7.5

vii

7.4.1

Methode . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

87

7.4.2

Optimalisatie naar maximaal koppel . . . . . . . . . . . . . . . . . . .

88

Conclusie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

92

8 Conclusie

93

Bibliografie

93

Symbolen en notaties Algemene notaties • Scalaire grootheden: R (weerstand), ε (rek),... • Eenheden: m (meter), nT (nanotesla),... • Operatoren en vaste functies: d/dx, sin(x) • Vectoren en matrices: A (vectorpotentiaal), J (jacobiaan), U0 ,... • Stapnummer, benadering: x(k) • Elementen van vector of matrix: x2 , Jk,l ,... • Schatting: x Symbolen A, A Magnetische vectorpotentieel, active set α, β Dempingsfactoren B, B Magnetische inductie B0 Parameter materiaalmodel b0 Breedte van het middelste deel van de meetsample bb Breedte van de zijdelen van de meetsample brt Breedte van de rotortanden bst Breedte van de statortanden d Diameter van de draden van de bekrachtigingswikkelingen d0 Dikte van het middelste deel van de meetsample dls Dikte van de luchtspleet Dls Diameter van de luchtspleet Dr Uitwendige diameter van het rotorjuk Ds Inwendige diameter van het statorjuk E Elasticiteitsmodulus viii

INHOUDSOPGAVE

ix

εe Elastische rek in het midden van de meetsample εp Plastische rek in het midden van de meetsample εt Totale rek in het midden van de meetsample εw Wenswaarde voor de rek in het midden van de meetsample F Magnetomotorische kracht φ Fysische magnetische flux, tweedimensionale basisfunctie van Lagrange Φ Gekoppelde magnetische flux g (x ) Functies met randvoorwaarden voor de optimalisatie Λ Magnetische permeantie H, H Magnetische veldsterkte, hessiaan H0 Parameter materiaalmodel I, i Bekrachtigingsstroom, stroomsterkte J, J Stroomdichtheid, jacobiaan L Totale lengte van de meetsample, Lagrange-functie L0 Lengte van het middelste deel van de meetsample Leq Equivalente lengte van de meetsample Lm Axiale lengte van de motor λ Lengte van het juk van de meetopstelling, Lagrange-vermenigvuldiger µ Absolute magnetische permeabiliteit µ0 Absolute magnetische permeabiliteit van het vacu¨ um µr Relatieve magnetische permeabiliteit ν Parameter materiaalmodel νp Modulus van Poisson ψ Eendimensionale basisfunctie Lagrange R Reluctantie; s: stator, spreiding, r: rotor, j: juk, t: tanden, ls: luchtspleet, 0: midden, 1: rand, k: kop S Oppervlakte Sd Oppervlakte van een draad in de bekrachtigingswikkeling Sw Oppervlakte beschikbaar voor bekrachtigingswikkeling σ Mechanische spanning in het middelste deel van de meetsample σys Opgegeven vloeigrens van het materiaal T Elektromagnetisch koppel

INHOUDSOPGAVE θ Rotorhoek U Oplossingenvector van de eindige-elementenberekening V Elektrische spanning w Aantal windingen Wem Elektromagnetische energie Wco Elektromagnetische co-energie wem Elektromagnetische energiedichtheid ξ Lokale parameter

x

Hoofdstuk 1

Inleiding Bij het produceren van een elektrische machine wordt het magnetisch circuit gevormd uit lamellen van elektrisch staal door lasersnijden, guillotinesnijden of ponsen. Het toepassen van deze processen op elektrisch staal kan echter een grote invloed hebben op de magnetische eigenschappen van dit materiaal [1]. Door het introduceren van plastische vervormingen aan de randen van de onderdelen wordt de magnetische karakteristiek immers lokaal gewijzigd [2, 3]. In het bijzonder bij lasersnijden van de onderdelen kan een grotere lokale materiaaldegradatie ontstaan dan bij mechanische processen [4]. Dit lasersnijden is echter een aantrekkelijke methode aangezien er geen fysisch contact vereist is, wat een grotere flexibiliteit toelaat. Traditioneel ontwerpt men het magnetisch circuit van een elektrische machine gebruik makend van de gemeten of gegeven materiaalkarakteristieken van het gebruikte elektrisch staal. Hierbij wordt rekening gehouden met de grootte van het resulterende of gewenste koppel en het koppelrimpel. Doordat bij het produceren echter een significante wijziging van de materiaaleigenschappen kan optreden, is de uiteindelijke machine niet noodzakelijk meer optimaal [5]. Op dit fenomeen wordt momenteel op twee manieren ingespeeld. In Europese bedrijven is het de gewoonte om een hoger dan vereiste kwaliteit van materiaal te gebruiken of de machine te overdimensioneren door middel van een empirische factor, zodat de machine ook na de productie nog aan de specificaties voldoet. In Amerika worden de onderdelen, indien dit mogelijk is, na de productie uitgegloeid om de spannings- en rekeffecten teniet te doen. Het spreekt voor zich dat beide methodes een zekere inefficiëntie in de productie introduceren. Daarom bestaat de noodzaak om een methode te ontwerpen die elektrische machines optimaliseert, rekening houdend met de lokale materiaaldegradatie. 1

HOOFDSTUK 1. INLEIDING

2

Bij het optimaliseren van het model van een elektrische machine moet een afweging gemaakt worden tussen de gewenste nauwkeurigheid en de vereiste rekentijd. Eén mogelijkheid is de machine voor te stellen als een magnetisch netwerk, wat tot een snelle berekening leidt maar geen grote nauwkeurigheid kan bieden, in het bijzonder als het model extra complex wordt door het in rekening nemen van de lokale materiaaldegradatie. Het alternatief is een numerieke berekening, bijvoorbeeld met de eindige-elementenmethode. Deze kan een heel precies resultaat geven, maar neemt ook een zeer grote rekentijd in beslag. Aangezien een optimalisatieprocedure meestal vele iteraties van het directe model vereist kan deze grote rekentijd het model zelfs onbruikbaar maken. Dit probleem wordt eveneens versterkt door de grotere complexiteit van een model waarbij de aangetaste materiaalranden in acht genomen worden. We zoeken dus een methode die een nauwkeurige optimalisatie van het magnetisch circuit kan uitvoeren in een aanvaardbare rekentijd. Deze optimalisatie moet uitgevoerd worden rekening houdend met de lokale materiaaldegradatie. Wij passen hiertoe de ‘space mapping’-methode toe. Dit is een recente methode die gebruik maakt van zowel een ruw als een fijn model. De vele evaluaties die eigen zijn aan het oplossen van een invers probleem, de optimalisatie van de machine, worden uitgevoerd in een ruw ‘surrogaatmodel’. Dit surrogaatmodel wordt echter voortdurend aangepast door middel van een afbeelding van de parameters tussen het ruwe en fijne model, zodat het ruwe model in een omgeving van het optimum het fijne model benadert. Opdat de geometrische parameters met een fysische werkelijkheid zouden overeenstemmen implementeren we ‘trust region’-strategieën. We bevorderen de convergentie door het gebruik van een hybride model, dat een overstap toelaat tussen het space-mappingalgoritme en een directe optimalisatie van het fijne model. De beide modellen berekenen oplossingen van het directe probleem: uitgaande van de machinevariabelen en de geometrische, elektrische en magnetische parameters wordt het koppelprofiel berekend voor opgegeven rotorhoeken en bekrachtigingsstromen. Het ruwe model is een semi-analytisch model waarbij de motor benaderd wordt door een magnetisch netwerk. In het ruwe model bepalen we de koppelkarakteristiek door het iteratief bepalen van de gekoppelde flux en vervolgens toepassen van de co-energiemethode. We werken met twee materiaalkarakteristieken, namelijk ‘gewijzigd’ materiaal over een zekere breedte aan de rand en ‘ongewijzigd’ materiaal in de rest van de machine. Het fijne model implementeren we als een eindige-elementenberekening. Hier wordt de wijziging van het materiaal aan de randen in fijnere stappen verdeeld. Beide modellen vatten we uiteraard parametrisch op, zodat een

HOOFDSTUK 1. INLEIDING

3

wijziging van een motorvariabele meteen tot een aangepaste oplossing leidt. De machine die we implementeren is de 6/4 reluctantiemotor. Bij deze machine is de invloed van de permeantie immers het meest uitgesproken, aangezien het koppel niet evenredig is met de permeantie maar met het verschil van de permeanties van de twee assen. We voeren alle berekeningen uit op een Intel Pentium 4 computer met een processorsnelheid van 3.2 Ghz en 3.11 GB intern geheugen. Het bepalen van de materiaalkarakteristieken gebeurt experimenteel. Door metingen op een trekbank onderzoeken we de invloed van homogene plastische vervormingen op de BH-karakteristiek van elektrisch staal. We stellen ook analytische benaderingen op van de opgemeten karakteristieken. Vervolgens gebruiken we de gegevens voor de plastische rek in het materiaal in functie van de afstand tot de rand om te bepalen welke magnetische karakteristieken op welke plaats moeten gebruikt worden in het ruwe en het fijne model. Deze scriptie stelt voor hoe we de magnetische eigenschappen van het elektrisch staal bepalen in de omgeving van de gedegradeerde randen. Vervolgens wordt de opbouw van het ruwe en het fijne model besproken. Verder wordt ingegaan op de verschillende optimalisatietechnieken die we gebruiken en op de implementatie van space mapping, inclusief de verschillende uitbreidingen ervan. Ten slotte wordt de invloed van de lokale materiaaldegradatie op het optimaal ontwerp van het magnetisch circuit van de machine geëvalueerd.

Hoofdstuk 2

Materiaaleigenschappen 2.1

Inleiding

De eerste stap die moet worden gezet om elektrische machines te optimaliseren met inbegrip van lokale materiaaldegradatie is het onderzoeken van de lokale materiaaldegradatie. Meer bepaald moeten de B-H-karakteristieken van het materiaal bepaald worden in functie van de afstand tot de rand. Deze gegevens kunnen dan gebruikt worden in het ruwe en het fijne model van de motor om het directe probleem op te lossen. Het gelamineerde magnetisch circuit van de machine onder beschouwing is opgebouwd uit niet-georiënteerd elektrisch staal (JIS 35H230) van hoge kwaliteit 1 . De materiaaleigenschappen zijn weergegeven in tabel 2.1. In dit hoofdstuk bespreken we enerzijds hoe metingen van de magnetische karakteristieken van het materiaal in functie van de plastische rek werden uitgevoerd en hoe deze gemeten karakteristieken analytisch benaderd werden. Anderzijds gaan we na hoe de gegevens voor de plastische rek in het materiaal in functie van de afstand tot de rand gebruikt werden om de uiteindelijke invloed van de productie van de machine op de magnetische eigenschappen te kwantificeren.

2.2

Metingen

2.2.1

Werkwijze

Het doel van de metingen is om de B-H-karakteristieken van het gegeven materiaal te bepalen bij verschillende plastische rekken. We voeren hiervoor metingen uit op een sample 1

Verliezen (50 Hz; 1.5 T) ≤ 2.3 W

4

HOOFDSTUK 2. MATERIAALEIGENSCHAPPEN

5

Tabel 2.1: Eigenschappen materiaal

material:

non-oriented electrical steel (Nippon Steel)

grade:

JIS 35H230 (corresp. with V230-35A)

Si-content:

≈ 3 wt%

thickness:

0.35 mm(±10%)

avg grain size:

≈ 150 µm

yield stress:

441 MPa

tensile strength:

539 MPa

elongation:

14 %

hardness:

213 HV(1)

density:

7600 kg/m3

electr. resistivity:

59 µΩ·cm

B25:

1.56 T

B50:

1.66 T

van het materiaal. De vorm en afmetingen van deze sample zijn weergegeven in figuur 2.1. Door middel van een trekbank leggen we een bepaalde mechanische rek op aan de sample. Vervolgens laten we het materiaal terug ontspannen zodat enkel plastische vervorming aanwezig blijft. Op deze plastisch vervormde sample meten we dan een reeks hysteresiscurves met toenemende amplitude. De eenwaardige B-H-karakteristiek van het materiaal wordt dan bekomen door de toppunten te nemen van de gemeten hysteresiscurves. Deze karakteristiek is dan uiteindelijk de maagdelijke curve van het materiaal.

Figuur 2.1: Afmetingen van de sample

Aangezien we een eenwaardige karakteristiek hanteren en we zodus de effecten van hysteresis verwaarlozen zullen we enigszins aan nauwkeurigheid inboeten. Het gebruik van zo een eenwaardige karakteristiek wordt echter opgelegd door de invoervereisten van onze eindigeelementenberekening.


2.2.2

6

Mechanische vervorming

De mechanische plastische rekken worden opgelegd door middel van een trekbank ontworpen voor cyclische mechanische belastingen [6], gedimensioneerd voor een maximale trekkracht van 17 kN en aangedreven door een stroomgestuurde brushless DC motor. Deze trekbank wordt hier echter gebruikt om een monotoon stijgende spanning aan te leggen. De verplaatsing van de trekbank wordt opgemeten door een zogenaamde Linear Variable Differential Transducer (LVDT) en de kracht kan opgemeten worden met behulp van een load cell. Het geheel is te zien in figuur 2.2.

Figuur 2.2: De trekbank

Aangezien we verwachten dat de grootste degradatie van de magnetische eigenschappen van het materiaal reeds voorkomt bij kleine plastische rekken, tussen 0 en 0.2 %, verrichten we veel metingen in dit gebied [2], zie ook figuur 2.5. We bepalen van tevoren bij welke plastische rekken we de metingen willen uitvoeren. Aan de hand van deze plastische rekken berekenen we dan de totale verplaatsing (∆L)t die we aan de sample moeten opleggen om tot een welbepaalde plastische vervorming te komen. Dit wordt verduidelijkt in figuur 2.3. Hierbij verwaarlozen we het kleine verschil in elastische rek tussen de huidige en de vorige meting, dat ontstaat door de iets hogere benodigde spanning in de huidige meting. We kunnen de totale verplaatsing ontleden als plastische en elastische verplaatsing (of rek): (∆L)t = (∆L)p + (∆L)e

(2.1)

De gewenste plastische verplaatsing kan bepaald worden als de gewenste plastische rek εp,w vermenigvuldigd met de oorspronkelijke lengte L0 van het plastisch vervormbare deel van de sample (het deel met vernauwde sectie). De elastische verplaatsing kan benaderd worden


7

door de werkelijke elastische verplaatsing in het vorige meetpunt, gedefinieerd als het verschil tussen de totale en plastische verplaatsing. Samengevat kunnen we de totale verplaatsing in het n-de meetpunt schrijven als:

(n)

(∆L)t

(n−1)

(n) = εp,w · L0 + ((∆L)t

− (∆L)(n−1) ) p

(2.2)

Hierbij staat de index w voor ‘wenswaarde’. Het komt er visueel op neer (zie figuur 2.3) dat we de totale op te leggen verplaatsing (rood), berekenen als de som van de gewenste plastische rek (groen) en de elastische rek in het vorige meetpunt (twee identieke zwarte lijntjes)

Figuur 2.3: Spanning-rekdiagramma

Voor het bepalen van de op te leggen verplaatsing om de gewenste plastisch rek te verkrijgen in het eerste meetpunt moeten we de elasticiteitsconstante E = σ/ε berekenen aan de hand van metingen in het elastische gebied. Plastische vervorming treedt uiteindelijk enkel op in het centrale, smallere gebied van de sample, aangezien daar een hogere spanning actief is. De elastische rek geldt echter over de volledige lengte en is evenredig met de plaatselijke spanning. Daarom definiëren we een equivalente lengte Leq , die rekening houdt met de grotere breedte aan de basis.

Leq = L0 + (L − L0 ) ·

b0 bb

(2.3)

In het elastische gebied is de rek in het midden van de sample εe = ∆L/Leq . De eerste op te leggen uitrekking is dan:

(1)

(∆L)t

= ε(1) w,p L0 +

σys Leq E

(2.4)


8

Hierbij is σys de opgegeven vloeigrens van het materiaal. Een zekere onnauwkeurigheid wordt ge¨ıntroduceerd, ten eerste door de onnauwkeurigheid van de opgegeven vloeigrens, verder door het feit dat bij de vloeigrens in feite al ongeveer 0.2 % plastische vervorming is opgetreden. Ten slotte hebben we ook geen rekening gehouden met de afrondingen aan het einde van het centrale deel. De fout op deze waarde leidt echter niet tot een systematische meetfout, aangezien de precieze verplaatsing gemeten wordt en aan de basis ligt van de volgende stap. Na het opleggen van deze verplaatsing ontspannen we het materiaal terug tot de spanning tot (n)

(n)

nul herleid is. Dan meten we de nieuwe breedte b0 en dikte d0 in het midden van de sample met behulp van respectievelijk een elektronische schuifmaat en een micrometer, waarna we de magnetische metingen uitvoeren. De dikte en breedte zijn nodig bij het bepalen van de magnetische eigenschappen, zie paragraaf 2.2.3. Enige onnauwkeurigheid wordt hierdoor opnieuw ge¨ıntroduceerd. Voor de dikte van de sample meten we immers constante, zelfs zeer licht stijgende waarden, in tegenstelling tot wat de theorie voorspelt:

∆d0 ∆b0 ∆L0 = = −νp d0 b0 L0

(2.5)

Hierbij is νp ≈ 0.3 de modulus van Poisson. Een test met een ongecoate sample wijst uit dat deze afwijking te wijten is aan het loskomen van de coating. Bovenstaande cyclus herhalen we voor alle gewenste plastische rekken tot er breuk optreedt bij εp = 17.9 %. We bekomen de volgende waarden: E = 191 MPa en Leq = 179 mm. De resultaten zijn samen met de resultaten van de magnetische metingen te vinden in tabel 2.2. We werken consequent met de technische rek en de technische spanning.

2.2.3

Elektromagnetisch meetcircuit

Bij elke plastische rek willen we een eenwaardige magnetische karakteristiek van het materiaal bekomen. Dit bereiken we door telkens een reeks van een dertigtal quasistatische metingen uit te voeren. Bij elke quasistatische meting wordt een laagfrequente mmk over een deel van de sample aangelegd terwijl het resulterende fluxverloop gemeten wordt. Op deze manier bekomen we voor elk meetpunt een reeks hysteresislussen waarvan de toppunten de gevraagde eenwaardige karakteristiek genereren. De kern van het meetcircuit is gebaseerd op de single sheet tester (SST). Essentieel is de single sheet tester een set van twee gekoppelde spoelen die rond de sample aangebracht worden, samen met een sluitjuk. Door de excitatiespoel wordt een bepaalde stroom gestuurd, die een mmk opwekt over het deel van de sample dat binnen het juk ligt. Afhankelijk van de


9

reluctantie wordt hierdoor een flux opgewekt, en de fluxverandering dφ/dt vertaalt zich in een spanning in de meetspoel. We bekijken eerst het circuit van de excitatiespoel, we gebruiken index 1. Voor een schematische weergave van de meetspoelen, zie figuur 2.4. De magnetische veldsterkte wordt gegeven door:

H=

I1 · w1 λ

(2.6)

Met I1 de stroom, w1 het aantal windingen en λ de magnetische weglengte in de sample (de binnenafstand tussen de poolschoenen) [7]. Hierbij wordt verondersteld dat er geen magnetomotorische-krachtval optreedt over de sluitjukken (deze hebben een veel grotere sectie dan de te bemeten sample en zijn gebouwd uit hoogpermeabel materiaal), zodat de aangelegde magnetomotorische kracht (I1 w1 ) bij benadering gelijk kan gesteld worden aan de magnetomotorische-krachtval over de sample, H ·λ. De stroom I1 wordt als volgt gestuurd. In het stuurprogramma, geprogrammeerd in LabView, wordt een bepaalde gewenste stroomvorm ingegeven. Deze wordt, onder de vorm van een spanningssignaal, via de data-acquisitiekaart naar een lineaire vermogensversterker gestuurd, die een stroom produceert volgens de ratio I1 = 1.2 Vin . Deze stroom wordt dan uiteindelijk naar de excitatiespoel gestuurd. We leggen een laagfrequent signaal aan, in casu 0.5 Hz, zodat zo weinig mogelijk dynamische effecten (wervelstromen) optreden. Voor elke reeks metingen moet het materiaal gedemagnetiseerd worden, dit is vooral belangrijk voor de eerste metingen met kleine amplitude. Hiertoe wordt een wisselstroom met exponentieel afnemende amplitude door de excitatiewikkeling gestuurd. Bekijken we nu het circuit van de meetspoel, we gebruiken index 2. De opgewekte spanning in de spoel wordt gegeven door:

V2 = w2 ·

dφ dt

(2.7)

met w2 het aantal windingen in de meetspoel en φ de fysische flux. Het spanningssignaal V2 vormt de ingang van een analoge integrator, die voor elke meting gereset wordt. Dit is ingebouwd in het programma. De opamp van de analoge integrator is voorzien van de mogelijkheid om met een potentiometer de offset te compenseren en zo de drift te beperken. De uitgang van de integrator is een spanningssignaal evenredig met de flux, dat dan via een data-acquisitiekaart in de computer wordt ingelezen.


10

Figuur 2.4: Schematische weergave van meetspoelen

Om uiteindelijk de inductie B = kBk te bepalen moeten we nog corrigeren met de doorsnede van de meetspoel. Bij grote veldsterktes is de flux in de lucht φl , tussen het materiaal van de sample en de meetspoel, immers niet meer verwaarloosbaar ten opzichte van de flux in het materiaal φmat . De relatieve permeabiliteit van het materiaal neemt immers af door verzadiging. Meer bepaald geldt:

φmat = φtot − φl

(2.8)

φtot Smat

(2.9)

dus met:

Bmeting = geldt:

Bmat =

Bmeting · Smat − µ0 H · Sl Smat

(2.10)

of:

Bmat = Bmeting − µ0 H ·

Sspoel −1 Smat

(2.11)

Met Sl de oppervlakte van de lucht tussen meetspoel en materiaal, Sspoel de oppervlakte van de spoel en Smat de oppervlakte van het materiaal van de sample. Op deze manier kunnen we voor elk meetpunt de veldsterkte en de inductie bepalen. Figuur 2.5 geeft de B-H-karakteristieken weer voor verschillende plastische rekken. Hierin zien we duidelijk dat


11

een kleine plastische vervorming (0.23 %) een grote invloed heeft op de B-H-karakteristiek. De zwarte karakteristiek heeft betrekking op de sample zonder plastische vervorming.

Figuur 2.5: B-H karakteristieken materiaal

2.3

Fitten van de gemeten materiaalkarakteristieken

Om de bekomen materiaaleigenschappen te gebruiken in onze modellen moeten ze eerst aangepast worden aan de gewenste invoervorm. Dit betekent dat we de discrete meetwaarden van de eenwaardige karakteristiek omzetten naar analytische formules. We gebruiken hiervoor de algemene formule:

H = H0 ·

B + B0

B B0

ν (2.12)

We moeten dus voor elke meting, of voor elke plastische rek, de parameters B0 , H0 en ν bepalen. Dit gebeurt door middel van het programma Microcal Origin 6.0. Door formule (2.12) te gebruiken houden we geen rekening met de vorm van de curve in het Rayleighgebied, in de buurt van de oorsprong, zie figuur 2.6. In het Rayleigh-gebied wordt de B-Hkarakteristiek benaderd door B = aH +bH 2 . Het resultaat van het bepalen van de parameters H0 ,B0 en ν als functie van de plastische rek wordt weergegeven in tabel 2.2. De relatieve permeabiliteit µr uitgaande van formule (2.12) wordt dus bepaald door:


B0 · µr = µ0 H0

1+

B B0

12

ν−1 !−1 (2.13)

Figuur 2.6: Karakteristieken in het Rayleigh-gebied

Tabel 2.2: Materiaalparameters in functie van plastische rek

n

εp (%)

H0

B0

ν

00

0

57.70

1.156

10.54

01

0.23

484.52

1.351

8.167

02

0.38

535.66

1.337

7.518

03

0.68

596.84

1.318

6.843

04

0.98

613.60

1.302

6.467

05

2.0

666.82

1.270

5.776

06

3.0

678.93

1.235

5.353

07

4.0

699.64

1.216

5.167

08

5.0

690.57

1.193

4.969

09

7.0

679.33

1.160

4.732

10

9.0

666.03

1.120

4.573

11

12.0

664.15

1.077

4.312

12

14.0

689.94

1.056

4.200

13

16.0

661.78

1.019

4.005


2.4 2.4.1

13

Ruimtelijke verdeling Algemeen

In deze paragraaf passen we de resultaten uit vorige paragrafen toe om de invloed van de plastische vervorming van het materiaal op de magnetische eigenschappen in de buurt van de randen te bepalen. In het onaangetaste materiaal, ver van de rand, gebruiken we uiteraard de meting van de materiaaleigenschappen voor het oorspronkelijke materiaal, namelijk meting n = 00 (zie tabel 2.2). Als basis voor de berekeningen gebruiken we de vergelijking voor de plastische vervorming in functie van de afstand tot de rand die vooropgesteld wordt in [8]: εp (%) =

y d

1 +b +a

(2.14)

Hierbij is d de dikte van de lamel. Deze wordt als parameter behouden om de flexibiliteit van de modellen niet onnodig te beperken. De andere constanten krijgen de volgende waarden: a = 0.0969

(2.15)

b = −0.3229

(2.16)

Vergelijking (2.14) resulteert in figuur 2.7, met de relatieve afstand tot de rand als abscis. Meer informatie over het bepalen van de materiaaleigenschappen aan de rand kan gevonden worden in [4]. De ruimtelijke verdeling van de materiaaleigenschappen zullen we op drie manieren moeten implementeren. In het fijne model werken we met vijf randgebieden, die allemaal verschillende eigenschappen kunnen krijgen. In het ruwe model wordt een bepaalde randdikte als ‘aangetast’ gedefinieerd en de rest van het materiaal als ‘onaangetast’. Als inputparameter is keuze tussen twee mogelijkheden. Ofwel worden de werkelijke materiaalkarakteristieken gebruikt, ofwel worden voor deze beide gebieden bilineaire karakteristieken gebruikt. In deze paragraaf zullen we voor elk geval de correcte parameters afleiden.

2.4.2

Verdeling van de materiaalkarakteristieken voor het fijne model

Het fijne model maakt gebruik van de eindige-elementenmethode, waarin een geometrische verdeling van de parameters gekarakteriseerd moet worden door gebruik te maken van subdomeinen met constante materiaaleigenschappen. Ruimtelijk verdelen we de magnetische


14

Figuur 2.7: Plastische vervorming in % in functie van y/d

materiaalparameters in zes verschillende subdomeinen (vijf randen en het ongewijzigde materiaal in het midden), waarbij in ieder subdomein de materiaalparameters constant zijn. We moeten dus enerzijds de diktes van de vijf randdomeinen bepalen en anderzijds de waarde van de magnetische materiaalparameters in ieder domein. In de praktijk moeten we de curve die de plastische rek weergeeft in functie van de afstand tot de rand dus discretiseren. Hiervoor lossen we deze functie voor alle meetpunten op naar de afstand tot de rand y, zodat elk meetpunt meteen gerelateerd wordt aan een zekere afstand tot de rand. We zullen ervoor kiezen om niet te interpoleren tussen verschillende meetpunten, om met andere woorden enkel plastische rekken te aanvaarden waarvoor magnetische metingen uitgevoerd zijn, en geen tussenliggende punten. Dit omdat het gedrag van de magnetische functie bij interpolatie van de parameters niet gekend is. Het algoritme dat we zullen volgen is als volgt: we overlopen de afstand- rekcurve (figuur 2.7) van rechts naar links, van kleine rekken naar grote rekken. We beginnen met het eerste punt. Hierbij berekenen we tot welke abscis, welke afstand tot de rand, de benaderende rechthoek zich moet uitstrekken zodat de oppervlaktes onder de rechthoek en onder de curve gelijk zijn. Vanaf deze abscis gaan we dan verder en bekijken we de punten die links ervan liggen. Voor elk van de eerstvolgende punten berekenen we tot welke abscis de corresponderende rechthoek zich zou uitstrekken, en dan kan een keuze gemaakt worden welke waarde gebruikt wordt. We zoeken met andere woorden telkens de y-waarde waarvoor:


Z

y (n−1) x d

y

15

1 + b − ε(n) dx = 0 +a

(2.17)

Hierbij is y de gezochte linkergrens van het huidige interval waarin de plastische rek door ε(n) benaderd wordt. Bij het laatste punt, dus bij hoge rek of dicht bij de rand van het materiaal, ligt de breedte van de rechthoek vast, namelijk tot de rand van het materiaal. Hier wordt opgelost naar de gewenste ordinaat of rek. Het resultaat wordt benaderd door een nabije waarde waarvoor magnetische metingen voorhanden zijn. Voor de eerste waarde, voor lage plastische rek, is van dit algoritme afgeweken. Dit omdat de plastische rek over een vrij grote afstand lager blijft dan de eerste rek waarvoor meetwaarden voorhanden zijn. Dit zou resulteren in een eerste interval dat een groot deel van het volledige domein beslaat. Toch kan, met figuur 2.5 in gedachten, aangenomen worden dat de materiaaleigenschappen bij zeer kleine mechanische vervormingen al snel naar de waarden van de eerste meting neigen. Daarom werd hier een arbitraire startwaarde gekozen voor de berekeningen. De resultaten zijn te zien in figuur 2.8 en tabel 2.3. In de tabel staat n voor het nummer van de gebruikte magnetische meting, zie tabel 2.2. Ook de relatieve afstand tot de rand die overeenstemt met de rek waarop het interval gebaseerd is en waarvoor metingen gebeurd zijn, dit is eigenlijk een soort centrum van het interval, wordt weergegeven, net als de breedte van elk interval en de rechtergrens ervan.

Tabel 2.3: Discretisering van de rek- afstandscurve

2.4.3

n

ε

Centrum

Breedte

Rechtergrens

01

0.23

1.71

1.6

3

03

0.68

0.90

0.88

1.4

05

1.95

0.33

0.33

0.52

08

5.0

0.09

0.15

0.19

10

9.0

0.01

0.04

0.04

Verdeling van de materiaalkarakteristieken voor het ruwe model

In het ruwe model kunnen twee verschillende sets materiaaleigenschappen ge¨ımplementeerd worden. We zullen eerst de randdikte en de te gebruiken materiaalkarakteristieken bepalen. Voor het model waarin de werkelijke materiaalkarakteristiek gebruikt wordt zijn de parameters bepaald aan de hand van tabel 2.2. Voor het bilineaire model echter zullen zowel de


16

Figuur 2.8: Discretisering van de rek- afstandscurve

randkarakteristiek als de karakteristiek van het onaangetaste materiaal nog gebilineariseerd worden. Voor de randdikte is het, rekening houdend met figuur 2.5, logisch om de volledige rand te nemen waar enige plastische vervorming voorkomt. De materiaaleigenschappen veranderen reeds zeer drastisch bij zeer kleine mechanische vervormingen. Volgens vergelijking (2.14) is de randdikte dus drie keer de lameldikte. Om de te gebruiken materiaalkarakteristiek te bepalen zal een procedure toegepast worden die gelijkaardig is aan wat gedaan werd bij het bepalen van de laatste karakteristiek voor het fijne model. We zoeken de meetwaarde waarbij een constante plastische rek over de hele rand dezelfde oppervlakte beslaat als de werkelijke curve, of de εp waarvoor geldt: Z 0

3 x d

1 + b − εp dx = 0 +a

(2.18)

Het resultaat is: εp = 0.83 ≈ 0.98

(2.19)

Aangezien we niet wensen te interpoleren kiezen we de dichtstbijzijnde meetwaarde. Dit is de meetwaarde met n = 04, die overeenstemt met een plastische rek van 0.98 %. De bilineaire benaderingen van de gekozen, aangetaste karakteristiek en de onaangetaste karakteristiek (n = 00) zullen we bepalen met behulp van de functie van Microsoft Excel


17

voor lineaire interpolatie. We werken op basis van de originele meetwaarden, en niet op basis van de gefitte analytische benaderingen. We nemen van beide sets de eerste 10 punten als basis voor een interpolatie, waarbij vereist wordt dat de resulterende rechte door de oorsprong gaat (zie figuur 2.9). Het is belangrijk om voldoende punten te nemen en de rechte door de oorsprong te definiëren, aangezien de vorm van de curve in het Rayleigh-gebied anders voor een aanzienlijke afwijking zou kunnen zorgen. Het resultaat van deze functie is de richtingscoëfficiënt van het eerste deel van de bilineaire karakteristieken.

Figuur 2.9: De bilineaire materiaaleigenschappen

Om het tweede deel van de karakteristiek te bepalen gebruiken we de laatste 6 punten. Deze keuze werd gemaakt op basis van de grafiek, zie figuur 2.9. We kiezen er dus voor om de karakteristieken te baseren op het begin en het einde van de curves. Hierdoor zal een vrij grote afwijking ontstaan in het midden van de curve, ter hoogte van de ‘knie’, doch dit wordt verkozen boven fundamentele afwijkingen bij alle waarden. Het resultaat van een lineaire interpolatiefunctie op deze verzameling punten is de richtingscoëfficiënt van het tweede deel van de karakteristieken en het snijpunt met de ordinaatas. De resultaten voor implementatie in het ruwe model zijn weergegeven in tabel 2.4.

2.5

Conclusie

Door middel van metingen op een trekbank hebben we de invloed onderzocht van de plastische rek van het gebruikte materiaal op de magnetische eigenschappen ervan (zie figuur 2.5). Vervolgens hebben we alle bekomen B-H-karakteristieken gefit met analytische formules, waardoor elke karakteristiek bepaald wordt door drie parameters zoals gedefinieerd in


18

Tabel 2.4: Materiaalparameters voor het ruwmodel

Centrum

Rand

n

00

04

ε(%)

0

0.98

Optie 1: werkelijke mat. kar. H0

57.7

613.6

B0

1.156

1.302

ν

10.54

6.467

Optie 2: bilineaire mat. kar. m0

0.011

0.0030

m1

5.6e-5

1.05e-4

b1

1.49

1.25

vergelijking (2.12). Uiteindelijk hebben we het verband tussen de afstand tot de snijrand van het materiaal en de mechanische eigenschappen ervan gezocht. Deze gegevens hebben we gebruikt om de materiaalkarakteristieken te selecteren die in het ruwe en fijne model gebruikt worden en om de diktes van de gedegradeerde randen in deze modellen te bepalen. Zie tabellen 2.3 en 2.4. De fysische eigenschappen die in de modellen van het directe probleem gebruikt worden zijn hiermee bepaald. Hierdoor kan ook het inverse probleem opgelost worden en kan dus ook de invloed van de productie van de machine op het optimaal ontwerp ervan bepaald worden.

Hoofdstuk 3

Ruw model 3.1

Inleiding

Indien we enkel zouden optimaliseren gebruik makende van een numeriek fijn model, dan zou dat zeer veel rekentijd vergen. Met het oog op space mapping, waarbij in de optimalisatie een minder accuraat maar vlugger werkend model dan het fijne model wordt ge¨ıncorporeerd, ontwikkelen we zo een ruw model. Toch moet dit model een vrij grote nauwkeurigheid bezitten, aangezien de afwijking tussen het ruwe en fijne model niet te groot mag zijn. Het ruwe model, dat in dit hoofdstuk besproken wordt, werkt met een benadering van de motor als magnetisch netwerk. Ons model zal als ingang enerzijds de vast genomen parameters, rotorhoeken en bekrachtigingsstromen hebben en anderzijds de variërende geometrische parameters. De uitgang van het model zal het berekende koppel zijn voor de verschillende parameters. Om de efficiëntie te verhogen is het model matrix- en vectorgewijs ge¨ımplementeerd: de berekening wordt tegelijk uitgevoerd voor elke stroom en rotorhoek waarvoor we een resultaat wensen.

3.2 3.2.1

Magnetisch netwerk Algemeen model

De motor die we behandelen wordt analytisch benaderd door een magnetisch netwerk. Hierbij wordt elk fysisch deel van de motor gesubstitueerd door een theoretische vervangreluctantie. De stroomvoerende wikkelingen van het bekrachtigde tandenpaar zijn vervangen door een theoretische magnetomotorische-krachtbron. De geometrie van de motor is te zien in figuur 3.1. Het volledige magnetisch netwerk wordt weergegeven in figuur 3.2. 19

HOOFDSTUK 3. RUW MODEL

20

Figuur 3.1: De geometrie van de motor

Kenmerkend voor het ruwe model in het kader van dit project is dat de degradatie van de materiaaleigenschappen aan de rand weergegeven wordt door een ‘binair’ model: er wordt een kunstmatige randdikte gedefinieerd waarin de materiaaleigenschappen uniform gedegradeerd zijn; in de rest van het materiaal wordt met de ongewijzigde materiaaleigenschappen gewerkt. Deze aangetaste randen worden als aparte reluctanties gemodelleerd. De gebruikte materiaalparameters werden bepaald in het hoofdstuk dat de materiaaleigenschappen behandelt. Daarin hebben we vermeld dat we werken met twee verschillende sets van materiaalkarakteristieken, namelijk met het bilineaire materiaalmodel of met de werkelijke B-H-karakteristieken. Aan de hand van een invoerparameter in het model is het mogelijk om tussen beide te kiezen.

3.2.2

Vervangreluctanties

Juk- en tandreluctanties In figuur 3.2 zien we de jukreluctanties van stator en rotor met indices ‘sj’ respectievelijk ‘rj’. De tanden zijn respectievelijk met ‘st’ en ‘rt’ aangeduid. Het laatste deel van de index heeft de volgende betekenis: ‘0’ staat voor het ongewijzigde materiaal in het centrum van een metaaldeel. ‘1’ wijst op het gedegradeerde materiaal van de randen en ‘k’ op de gedegradeerde koppen van de tanden. De index ‘s’ duidt de spreidingsreluctantie aan. We beschouwen in dit model enkel de twee bekrachtigde tanden van de stator, de bijhorende tanden van de rotor en de jukken van stator en rotor. De niet-essentiële tanden worden met andere woorden verwaarloosd.


21

Figuur 3.2: Het magnetische netwerk van het ruwe model met ls de verschillende delen van de luchtspleet en hierbij de indices overeenstemmend met de delen aangegeven op figuren 3.4 en 3.5. rj: rotorjuk, sj: statorjuk, s: spreiding, rt: rotortand, st: statortand en hierbij 1: aangetaste rand, 0: niet aangetast, k: aangetaste kop

Het statorjuk wordt behandeld als twee parallelle balken met als lengte het gemiddelde van de bekomen lengtes op basis van de binnen- en buitendiameter en aan beide zijden eindigend na één vierde van de tandbreedte, aannemend dat de flux over de breedte van de tanden uniform wegvloeit uit het juk, zie figuur 3.3. Elke balk wordt gemodelleerd als een parallelschakeling van twee reluctanties die de aangetaste randen vertegenwoordigen en een centrale reluctantie. Met andere woorden geldt voor de reluctantie van het volledige statorjuk:

Rsj

1 = · 2

2 1 + Rsj1 Rsj0

−1 (3.1)


22

Met hierbij:

lsj µsj1 · Ssj1 lsj = µsj0 · Ssj0

Rsj1 =

(3.2)

Rsj0

(3.3)

Hierbij is l de lengte van het gebied dat overeenstemt met een reluctantie. Binnen één reluctantie wordt de permeabiliteit van het materiaal dus als constant beschouwd. De permeabiliteit wordt voor elke reluctantie bepaald aan de hand van de materiaaleigenschappen (de B-H-karakteristiek) en de magnetische flux in die reluctantie. Dit is uiteraard een iteratief proces, dat in paragraaf 3.3 meer in detail besproken wordt. Het rotorjuk wordt volledig analoog behandeld.

Figuur 3.3: De flux door statortand en -juk

De statortanden zijn iets gecompliceerder aangezien ook het uiteinde ervan, dat uit gedegradeerd materiaal bestaat, in rekening gebracht moet worden. We kunnen de statortanden modelleren op twee manieren. In het ene scenario zien we het lichaam en de kop van de tand als een serieschakeling. Hierbij nemen we in feite aan dat de flux volledig uniform door de tandkop stroomt. In het andere geval maken we een parallelschakeling van het centrale deel en de twee randen, telkens inclusief het kopdeel. We nemen dan eigenlijk aan dat de fluxen in de respectievelijke delen van kop en lichaam gelijk blijven. In werkelijkheid ‘verspreidt’ de flux zich opnieuw doorheen de tandkop, zoals te zien in figuur 3.4. De werkelijkheid ligt dus tussen beide modellen in. Wij gebruiken het tweede model.

Luchtspleetreluctantie We modelleren de luchtspleetreluctantie alsof de tanden rechthoekig zijn en een lineaire, horizontale translatie ten opzichte van elkaar ondergaan. Deze translatie is uiteraard gekoppeld


23

Figuur 3.4: De kop van de statortand

met de positie van de rotor volgens:

∆x = θ ·

Dls 2

(3.4)

Hierbij is x de translatie, θ de rotatie van de rotor ten opzichte van de gealigneerde toestand van stator- en rotortand en Dls de diameter van de luchtspleet. Hierbij is impliciet een coördinatentransformatie gebeurd. De fysische, horizontale translatie (∆x)lin is immers evenredig met sin(θ) in plaats van met θ. De horizontale coördinaat x, en ook ∆x, beschouwen we echter als de booglengte, en niet als een lineaire afmeting. Dus is ∆x = (∆x)lin / sin(θ). De totale luchtspleetreluctantie bestaat uit verschillende onderdelen, afhankelijk van de positie van de rotor. De verschillende onderdelen zijn weergegeven voor een rotorhoek van 40◦ in figuur 3.5 en voor 0◦ in figuur 3.4.

Figuur 3.5: De luchtspleetflux

Het eerste deel van de luchtspleetflux is de flux in de zone waar beide tanden overlappen. De reluctantie die met deze flux overeenstemt is:


24

Rls1 =

dls µ0 · Sls1

(3.5)

Hierbij is Sls1 het product van de overlapping en de motordiepte (Lm ) en is dls de dikte van de luchtspleet. Het model is ook bruikbaar als de rotortand breder is dan de statortand. Het tweede deel van de flux vloeit tussen beide tanden, links van de overlapping. We nemen hierbij aan dat de fluxlijnen de vorm van cirkelbogen van 90◦ hebben. Zolang er overlap is tussen beide tanden (zie figuur 3.4) begint deze flux aan de linkeronderhoek van de statortand. Ook hier is het model geldig voor een bredere statortand dan rotortand, we bespreken verder echter het normale geval. Door de cirkelvorm is er in het model in dit geval altijd een breedte ter grootte van de luchtspleetdikte waar geen flux naar de rotortand vloeit. De flux stopt aan het einde van de rotortand. Wanneer de tanden niet meer overlappen wordt de breedte van de flux volledig bepaald door het begin en het einde van de rotortand. De permeantie Λls2 die hiermee overeenstemt wordt berekend door integratie: Z µ0 ·

Λls2 = Sls2

dS l

(3.6)

met, als y de integratieparameter is over de breedte van de rotortand:

dS = dy · Lm l=

2π ·y 4

(3.7) (3.8)

Hierbij is l de booglengte van de fluxlijnen. We bekomen dus, met Lm de axiale lengte van de motor:

Z

y2

Λls2 = y1

2µ0 · Lm 2µ0 Lm y2 dy = · ln( ) π·y π y1

(3.9)

Het derde deel van de flux is analoog, zie figuur 3.5. Dit deel wordt, behalve voor zeer kleine rotorhoeken, niet beperkt door de breedte van de rotortand, maar door die van de statortand. We moeten echter rekening houden met het feit dat de rotortanden meestal een vrij beperkte hoogte hebben. Op het moment dat de cirkelbogen van de flux het rotorjuk raken (zie figuur 3.5) stopt dit boogmodel en ontstaat een vierde soort flux, namelijk deze tussen de statortand en het rotorjuk. Deze stroomt gewoon rechtlijnig, alsof er in die zone een grote luchtspleet is tussen rotorjuk en statortand.


25

De eerste drie reluctanties worden in het model in parallel geschakeld. Dit betekent dat, in het model, deze delen van de flux door de volledige stator- en rotortand stromen. Dit is uiteraard een afwijking van de werkelijkheid. Het tweede en derde deel van de flux bypassen immers een deel van een tand, en missen ook de tandkop van één der tanden. In de plaats van een tandkop te doorstromen, doorstromen ze de eveneens aangetaste zijkant van de tand. Daarom nemen we aan dat de afwijking die ontstaat aanvaardbaar is in dit ruwe model. De vierde reluctantie wordt parallel geschakeld met de rest van de luchtspleet én de rotortand. Het gedegradeerde metaal dat de flux doorstroomt bij het binnentreden in het juk wordt niet gemodelleerd. Een nauwkeurige modellering van de intrede door de gedegradeerde rand zou zeer gecompliceerd zijn, aangezien zowel bij het juk als bij de zijkant van de tanden een gezamenlijk effect in rekening gebracht zou moeten worden van een flux die binnentreedt en een flux die in de loodrechte richting stroomt. Aan de hand hiervan zou de verzadigingstoestand van het materiaal berekend kunnen worden. De beperkte nauwkeurigheid van ons model voor de fluxlijnen in de lucht rechtvaardigt echter dat we dit niet gedetailleerd modelleren.

Spreidingsreluctantie De spreidingsreluctantie van de bekrachtigingswikkelingen staat in parallel met het hele verdere netwerk. We gebruiken het fijne model, dat in FEMLab ge¨ımplementeerd is (zie hoofdstuk 4) om deze reluctantie te berekenen als het quotiënt van de fysische flux door de wikkelingen en de bekrachtiging. We beschouwen de spreidingsflux als het deel van de flux dat de rotor niet binnentreedt, zie figuur 3.6. Daarom isoleren we dit deel van de flux door de rand van de rotor als een magnetische isolatie te definiëren. We leggen met andere woorden aan de rand van de rotor een Dirichlet-randvoorwaarde op. Dit levert het fluxpatroon van figuur 3.7 op.

Figuur 3.6: De spreidingsflux, samen met de hoofdflux


26

Figuur 3.7: De spreidingsflux ge¨ısoleerd

We bemerken dat de sterkte van de bekrachtigingsstroom geen invloed heeft op de spreidingsreluctantie. Dit is in feite onnauwkeurig: door het niet-lineaire systeem kan de spreidingsflux niet zomaar op de hoofdflux gesuperponeerd worden, maar hangt deze ook af van de verzadiging die door de hoofdflux veroorzaakt wordt. De spreidingsreluctantie blijkt wel afhankelijk te zijn van de rotorhoek. We hebben een berekening uitgevoerd voor de hoeken van 0 tot 45◦ met telkens een interval van 5◦ . De resultaten zijn te zien in tabel 3.1. In het ruwe model bepalen we door lineaire interpolatie van deze waarden voor elke opgegeven rotorhoek de spreidingsreluctantie.

Tabel 3.1: Spreidingsreluctantie

3.2.3

α (◦ )

Aw R 106 Wb

0

8.540

5

8.250

10

8.152

15

8.189

20

8.312

25

8.584

30

9.151

35

9.602

40

9.950

45

10.069

Bekrachtiging

De stroomvoerende spoelen introduceren in het magnetisch circuit van de werkelijke motor een magnetomotorische kracht. Volgens de wet van Ampère geldt voor de magnetomotorische kracht opgewekt door een spoel het volgende:


27

F =w·I

(3.10)

I is de stroom door de spoel en w is het aantal windingen van de spoel. Deze uitdrukking wordt bekomen door te integreren over een pad rond de windingen. We vatten de bekrachtiging op als twee spoelen in serie, rond de twee overstaande statortanden. De bekrachtigingswikkelingen van de andere tanden voeren geen stroom en worden hier dus ook niet beschouwd. We nemen aan dat de volledige ruimte tussen het statorjuk, de statortanden en de luchtspleet beschikbaar is voor windingen. Uiteraard wordt de ruimte voor de wikkelingen van de andere tanden wel vrijgelaten. De wikkelingsoppervlakte Sw is dus afhankelijk van de geometrische parameters van de machine. We geven de draaddiameter d als parameter aan de procedure waardoor de oppervlakte Sd van één draad berekend kan worden. Het aantal windingen kan dan berekend worden als het quotiënt van de wikkelingsoppervlakte en de draadoppervlakte. Hierbij nemen we een correctiefactor in acht, die eveneens als parameter kan opgegeven worden. Meestal is niet de volledige oppervlakte immers opgevuld met stroomvoerend materiaal: bij ronde draden blijft vrij veel lucht vrij tussen de draden. Twee uiterste gevallen zijn getekend in figuur 3.8. Figuur 3.8a illustreert draden die vierkant gestapeld zijn. De oppervlakte van één draad wordt gegeven door:

Sd = d 2 ·

π 4

(3.11)

De totale oppervlakte die door een draad in beslag genomen wordt is echter Sd = d2 , zodat we de volgende correctiefactor voor het aantal windingen krijgen:

Ca =

π = 0.785 4

(3.12)

Het aantal windingen is dus:

w = Ca ·

Sw Sd

(3.13)

In figuur 3.8b zijn de windingen compact gestapeld. Hier geldt een wikkelfactor van:

Cb =

−1 4 o · sin(60 ) = 0.907 π

(3.14)


28

Bij het berekenen van het aantal windingen houden we tevens rekening met het feit dat slechts de ene kant van de spoel unieke windingen vertegenwoordigt en dat de spoel bestaat uit de wikkelingen van beide tanden in serie.

Figuur 3.8: Stapelingen van de bekrachtigingswindingen

3.3

Fluxberekening

Om tot een koppelprofiel te komen moeten we eerst de gekoppelde flux door de bekrachtigingswikkelingen berekenen. Deze berekening kan echter niet rechtstreeks uitgevoerd worden, aangezien veel elementen van het magnetisch netwerk niet-lineair zijn, en dus afhankelijk zijn van de flux. Een eenvoudige iteratieve oplossing kan hier het best toegepast worden [9]. De procedure is voorgesteld in figuur 3.9. We berekenen eerst, aan de hand van opgegeven startwaarden voor de magnetische inducties B = kBk in elke niet-lineaire reluctantie en de opgegeven B-H-karakteristieken van het materiaal van de reluctanties, alle nieuwe µ-waarden. Als parameter geven we op of we met werkelijke materiaalkarakteristieken werken of met bilineaire benaderingen, zoals besproken in paragraaf 2.4.3. Als we met de werkelijke materiaalkarakteristieken werken zijn ν, B0 en H0 als parameters opgegeven voor het aangetaste en onaangetaste materiaal. Deze parameters bepalen de reeds besproken B-H-karakteristiek (met H = kH k):

H = H0 · Aangezien we de waarden voor µ

B + B0

B B0

ν (3.15)


29

Figuur 3.9: Procedure om flux te berekenen

µ= moeten bepalen, berekenen we:

B H

(3.16)


30

B0 · µ= H0

1+

B B0

ν−1 !−1 (3.17)

Indien we opteren om te werken met de bilineaire karakteristiek, zijn m1 , m2 en b de op te geven parameters. m1 is de richtingscoëfficiënt van het eerste lineaire deel van de karakteristiek, m2 is de richtingscoëfficiënt van het tweede lineaire deel en b2 is het snijpunt met de ordinaatas van het tweede lineaire deel van de karakteristiek. µ(B) is dus een stuksgewijs lineaire functie. Het snijpunt wordt gevonden door in de vergelijkingen:

Bs = m1 · Hs

(3.18)

B s = m2 · H s + b 2

(3.19)

H te elimineren: Bs =

b2 m1 m1 − m2

(3.20)

Voor het eerste deel van de karakteristiek is µ = m1 . Voor het tweede deel hebben we:

µ=

b2 b2 m2 B = m2 + = m2 + H H B − b2

(3.21)

Een andere deelprocedure bepaalt dan uit de resulterende permeabiliteiten de totale reluctantie waardoor de totale gekoppelde flux door de spoel berekend kan worden:

Φ = w · φtot =

w·F Rtot

(3.22)

Deze totale flux, weer omgerekend naar de fysische flux, wordt samen met de waarden van de verschillende reluctanties van het magnetische netwerk gebruikt om de flux door elke reluctantie te bepalen. Deze fluxen worden omgerekend naar de bijhorende magnetische inducties, met behulp van de doorsneden van de onderdelen die door de reluctanties vertegenwoordigd worden. Met deze vernieuwde waarden voor de inducties in de onderdelen van het netwerk worden dan opnieuw de permeabiliteiten van het materiaal berekend met behulp van de eerste deelprocedure. Deze iteraties gaan verder tot de oplossing convergeert. We hebben echter demping moeten voorzien opdat de convergentie bevorderd zou worden. Deze is ge¨ımplementeerd in de deelprocedure die de fluxen en inducties berekent. De demping berekent telkens de stap ∆B


31

van de inductie in elke reluctantie, tussen de nieuwe, berekende waarde en de vorige waarde. Deze stap wordt vermenigvuldigd met een dempingsfactor α, die opgesteld is als volgt: π · (∆B)t α = tanh ∆B

(3.23)

(∆B)t is de grens voor de stapgrootte die ongewijzigd toegelaten wordt. De initiële (∆B)t wordt als parameter opgegeven aan de functie. Wanneer elke stap, dus voor elke stroom, hoek en reluctantie, onder een bepaalde grens komt te liggen wordt (∆B)t verkleind met een op te geven factor. Ook als een op te geven aantal iteraties wordt overschreden bij dezelfde (∆B)t wordt deze factor aangepast. We besluiten dat convergentie bereikt is als de absolute waarde van de stap (∆B), relatief ten opzichte van de inducties, of meer bepaald het maximum van de stappen voor elke stroom, rotorhoek en reluctantie, kleiner is dan een bepaalde waarde. Ook deze tolerantie is op te geven als invoerparameter. Dit systeem zou echter tot een valse convergentie kunnen leiden, aangezien de stap na een bepaald aantal iteraties verkleind wordt, ook als bij die stapgrootte nog niet voldoende convergentie opgetreden is. Om dit te vermijden stellen we voorop dat niet de werkelijke, gedempte stap kleiner moet worden dan de tolerantie, maar wel de stap die zou worden gezet als er geen demping aanwezig was, dus de stap die met de dempingsfactor vermenigvuldigd wordt. Verder eisen we dat deze stap gedurende een klein aantal opeenvolgende iteraties binnen de gestelde limieten moet vallen.

3.4

Koppelberekening

De berekende gekoppelde flux Φ dient nu de ingang te zijn voor de procedure die het koppel berekent. Het koppel wordt berekend met de co-energiemethode. Deze methode kunnen we als volgt inzien. Bij een vaste rotorhoek wordt de magnetische energie in het systeem opgebouwd door toename van de bekrachtigingsstroom:

dWe = i · dΦ = dWm

(3.24)

De magnetische energie wordt dus gegeven door (bij een vaste rotorhoek):

Z Wm (θ0 , Φ0 ) =

Φ0

i (θ0 , Φ) dΦ = We 0

(3.25)


32

Wanneer de rotor in beweging komt moet van de toename van de magnetische energie, te wijten aan de veranderende bekrachtigingsstroom, de mechanische arbeid die de machine levert worden afgetrokken.

dWm = i · dΦ − T · dθ

(3.26)

De co-energie wordt gedefinieerd als [10] (zie ook figuur 3.10):

Wco (θ0 , Φ0 ) = I0 · Φ0 − Wm (θ0 , Φ0 )

(3.27)

dWco = Φ · di + T · dθ

(3.28)

Figuur 3.10: Definitie van de co-energie

Het koppel wordt dan gegeven door: T (θ0 , I0 ) =

∂ Wco (θ0 , I0 ) ∂θ

(3.29)

Om de co-energie te berekenen, voor elke stroom, moeten we dus Φdi integreren tot die stroomwaarde. Dit uiteraard voor elke rotorhoek. Die integratie vatten we aan met de regel van Simpson. Zie figuur 3.11 en [11]. Hierbij wordt elke set van drie punten, of van twee intervallen, benaderd door een tweedegraadsfunctie. De benadering luidt:

Z a

b

b−a f (x)dx ≈ 6

f (a) + 4 · f

a+b 2

+ f (b)

(3.30)

We passen deze integratie toe voor elke stroomwaarde. Indien het aantal punten waarvoor Φ gekend is binnen het integratie-interval oneven is benaderen we het eerste interval door


33

Figuur 3.11: Regel van Simpson

een driehoek. Het laten wegvallen van tussenpunten is immers niet toegelaten aangezien de abscissen equidistant moeten zijn. De ingegeven stroomwaarden waarbij het ruwe model moet uitgevoerd worden, worden hier ook op getest bij invoer. Er ontstaat wel een probleem als de berekening bijvoorbeeld maar bij één of twee (grote) stroomwaarden uitgevoerd wordt. Dan wordt deze integratie immers onnauwkeurig aangezien geen tussenpunten meer aanwezig zijn. In dit geval, als er minder dan 4 bekrachtigingsstromen ingegeven worden, worden aan het begin van de procedure een aantal tussenpunten ingevoerd waarbij de berekening ook uitgevoerd wordt. Deze worden niet als uitgang gegeven, maar wel gebruikt om een nauwkeurige integratie te bekomen. Om het koppel te berekenen moet de co-energie afgeleid worden naar de rotorhoek. Ook deze differentiatie moet discreet uitgevoerd worden:

(k+1) (k−1) ∂ Wco − Wco Wco (θ(k) , I) ≈ . ∂θ 2∆θ

(3.31)

Ook de rotorhoeken waarbij de berekening dient te gebeuren moeten dus equidistant zijn. Bij 0◦ en 45◦ is het koppel wegens symmetrie gelijk aan nul. Bij een andere begin- of eindwaarde van de hoekenvector waarvoor het koppelprofiel moet bepaald worden kan eenzijdig worden afgeleid:

(2) (1) Wco − Wco ∂ Wco (θ(1) , I) ≈ . ∂θ ∆θ (n) (n−1) ∂ Wco − Wco (n) Wco (θ , I) ≈ . ∂θ ∆θ

(3.32) (3.33)


3.5

34

Resultaten

Op een performante PC is de rekentijd van het ruwe model in de orde van enkele seconden, ondanks het feit dat in sterk niet-lineaire gevallen vaak meer dan honderd iteraties nodig zijn om convergentie te bereiken. Het ruwe model blijkt betrouwbaar te functioneren als fysisch verantwoorde parameters worden opgegeven. Ook de nauwkeurigheid ervan is zeer aanvaardbaar, in het bijzonder bij kleinere bekrachtigingsstromen. Bij grotere bekrachtigingsstromen is een onderschatting van de flux te zien bij kleine rotorhoeken, wat tot een te laag koppel leidt bij kleine rotorhoeken en een te hoog koppel bij grote rotorhoeken. Figuur 3.12 toont een vergelijking de koppelprofielen van het ruwe model (rood) met de exact veronderstelde koppelprofielen van het fijne model (blauw).

Figuur 3.12: Nauwkeurigheid van koppelprofiel ruw model

3.6

Conclusie

We hebben een ruw model opgesteld voor het oplossen van het directe probleem. Het model is volledig parametrisch, zoals vereist wordt om een optimalisatie te kunnen uitvoeren met behulp van space mapping of een andere optimalisatieprocedure. Het model voldoet tevens aan de vereisten wat betreft snelheid en nauwkeurigheid.

Hoofdstuk 4

Fijn model 4.1

Inleiding

Opdat we de fluxen en koppels van de geschakelde reluctantiemotor op een accurate manier zouden kunnen berekenen, bouwen we een numeriek model. De numerieke methode die de veldberekeningen uitvoert is de eindige-elementenmethode. In het kader van de optimalisatie met space mapping is dit model het zogenaamde fijne model. Het is een nauwkeuriger model dan het ruwe model, maar heeft meer tijd nodig om de veldberekeningen uit te voeren. We hebben het model net als het ruwe model ge¨ımplementeerd in Matlab. Hierbij worden procedures van het programma FEMLab 3.1 aangeroepen die de eindige-elementenmethode toepassen. Ook dit model hebben we zeer parametrisch opgevat. We hebben de procedure die het fijne model implementeert eenzelfde interface gegeven als die van het ruwe model. We hebben ervoor gekozen om één algemene gegevensstructuur op te stellen die zowel als input kan dienen voor het ruwe als het fijne model. Zowel alle geometrische kenmerken als de beoogde nauwkeurigheid, fysische eigenschappen en andere voorkeuren kunnen als parameter opgegeven worden. Ook de rotorhoeken en bekrachtigingsstromen kunnen opgegeven worden. De uitgangen van de twee modellen hebben eveneens dezelfde vorm zodanig dat space mapping op een eenvoudige manier ge¨ımplementeerd kan worden. In dit hoofdstuk bekijken we alle aspecten van het fijne model zoals ze doorlopen worden door de routines van ons opgebouwd programma. Zo bekijken we eerst het opstellen van de geometrie. Daarna nemen we het principe en de implementatie van het meshen onder de loep, gevolgd door de fysische eigenschappen van de verschillende domeinen. Vervolgens beschouwen we het principe van het feitelijke oplossen van de niet-lineaire partiële differentiaalvergelijkingen. Uiteindelijk hebben we het over de berekening van het mechanische koppel. 35

HOOFDSTUK 4. FIJN MODEL

36

Ten gepaste tijde gaan we iets dieper in op de theoretische achtergrond van een bewerking, anderzijds bekijken we ook de praktische uitvoering.

4.2

Geometrie

Het opbouwen van de geometrie gebeurt bij het uitvoeren van het programma slechts één keer per opgegeven rotorhoek. Elke opgebouwde geometrie gebruiken we dus om bij die hoek voor elke opgegeven bekrachtigingsstroom de berekening uit te voeren. De hele geometrie is tweedimensionaal opgevat. De volledige geometrie van de motor is te zien in figuur 4.1.

Figuur 4.1: De volledige geometrie

In FEMLab kunnen drie types geometrische elementen worden gedefinieerd: punten, curves en oppervlakken. Wij hebben de volledige geometrie opgebouwd gebruik makend van oppervlakken of unies en doorsnedes van oppervlakken. De basiselementen zijn eenvoudige cirkels en rechthoeken. In FEMLab zelf wordt elk geometrisch element echter opgeslagen als een gegevensstructuur met de punten die de verschillende delen van de randen scheiden, de randen zelf, gedefinieerd als Bézier-curves van eerste, tweede of derde orde en eigenschappen die subdomeinen van het bewuste element definiëren. Achteraf worden alle elementen samen geanalyseerd en worden de eigenschappen van de volledige geometrie bepaald. We hebben ervoor gekozen om elk oppervlak te laten samenvallen met een fysisch subdomein. Dit is een samenhangend gebied met constante materiaaleigenschappen. Een andere mogelijkheid was geweest, gewoon overlappende oppervlakken te definiëren. Deze zouden


37

dan door het programma in subdomeinen worden ingedeeld. De reden voor onze werkwijze wordt verduidelijkt in paragraaf 4.4. Metalen onderdelen hebben we ook hier telkens afgezoomd door een aangetaste rand. We vatten dit op als een opeenvolging van vijf trappen met op te geven diktes, meer informatie daarover is te vinden in paragraaf 2.4.2. Verder is het concept van dit model analoog aan wat aan de basis ligt van het ruwe model. Het gaat eveneens om een 6/4 reluctantiemotor waarvan alle geometrische kenmerken parametrisch te bepalen zijn. We nemen aan dat de hele vrije ruimte tussen de statortanden beschikbaar is voor bekrachtigingswikkelingen. Uiteraard modelleren we effectief enkel de wikkeling van de bekrachtigde tanden, zoals weergegeven in figuur 4.1.

4.3

Meshing

Het meshen van een tweedimensionale geometrie betekent essentieel het verdelen van deze geometrie in driehoekige elementen. De elementen mogen enkel een volledige zijde of een hoekpunt gemeenschappelijk hebben, er mogen met andere woorden geen meshknopen in een zijde van een meshelement liggen. De elementen moeten de volledige geometrie bedekken zonder te overlappen en één element moet volledig uit één fysisch medium bestaan. Het spreekt voor zich dat voor elke verschillende geometrie en dus voor elke verschillende rotorhoek een mesh opgesteld wordt. Met één mesh kunnen wel voor alle bekrachtigingsstromen de berekeningen uitgevoerd worden. De volledige mesh is te zien in figuur 4.2, een detail in figuur 4.3 (voor θ = 0). Bij het opstellen van de mesh worden eerst de bekende punten van de geometrie als meshknopen gedefinieerd. Vervolgens worden de grenzen tussen twee subdomeinen gediscretiseerd door enkele knopen op deze krommes te leggen. Door deze knopen te verbinden ontstaat een initiële, zeer ruwe mesh. Hierin worden dan inwendig knopen toegevoegd om de uiteindelijke mesh te bekomen. Om de fijnheid en omvang van de mesh te bepalen kunnen verschillende parameters worden opgegeven. We hebben ervoor gekozen in ons fijne model slechts drie parameters als ‘input’ te aanvaarden. De eerste parameter, een natuurlijk getal van 1 tot 5, definieert een keuze voor een voorgedefinieerde set waarden voor een aantal interne parameters, in overeenstemming met de wijze waarop deze gegevens in FEMLab zelf kunnen ingevoerd worden. Een grotere waarde voor deze parameter leidt tot een ruwere mesh. De (interne) meshparameters die hierdoor bepaald worden zijn de volgende.


38

Figuur 4.2: De volledige mesh

Figuur 4.3: Een detail van de mesh

De maximum-elementgrootteschaalfactor bepaalt de maximum toegelaten elementgrootte, relatief ten opzichte van de standaard maximum elementgrootte. Deze standaard maximum elementgrootte is 1/15 van de grootste afstand in de geometrie parallel aan de assen. Deze parameter varieert bij onze voorgedefinieerde waarden van 1 tot 5. Door de intrinsieke fijnheid van de mesh, aangezien er zeer fijne gebieden in de geometrie voorkomen, zal deze instelling meestal echter slechts invloed hebben in de omgeving van de motor, dus buiten de motor. De elementgroeiratio bepaalt hoe snel aanliggende elementen in grootte kunnen toenemen, bijvoorbeeld als een zeer fijne sector aanligt bij een grote sector waar een fijne mesh niet noodzakelijk is. Deze situatie doet zich veelvuldig voor in onze geometrie gezien de zeer fijne materiaalranden. De waarden die wij gebruiken variëren van 1.3 tot 2, dit betekent dat een


39

element 30% tot 100% groter kan zijn dan een aanliggend element. De meshkromtefactor bepaalt de grootte van de lijnstukken waarin de geometrische randen onderverdeeld worden (en dus de driehoeken die erop gebaseerd zijn), relatief ten opzichte van de kromtestraal van de randen. Wij hebben waarden van 0.3 tot 1, maar aangezien de kromtestralen in onze geometrie steeds vrij groot zijn zal dit geen bepalende factor zijn. De meshkromte-afsnijwaarde vermijdt het voorkomen van zeer veel elementen bij sterk gekromde delen van de geometrie. Wanneer de kromtestraal kleiner is dan deze waarde, vermenigvuldigd met de grootste afstand in de geometrie parallel aan de assen, wordt voor het meshen gebruik gemaakt van dit product in plaats van de kromtestraal. Ook deze optie zal over het algemeen slechts een beperkte invloed hebben gezien de relatief grote kromtestraal van de randen in onze geometrie. De waarde varieert van 0.001 tot 0.05. Deze vier meshparameters worden dus gezamenlijk aangepast door één invoerparameter. De tweede invoerparameter is een geavanceerde mesheigenschap die in ons geval wel een grote invloed heeft. Het gaat om de resolutie in nauwe gebieden, dit is het aantal lagen van meshelementen die gevormd worden in nauwe gebieden. In onze geometrie heeft dit betrekking op alle gedegradeerde randen van de materiaalgebieden. Door de beperkte groeisnelheid van de mesh, zoals vermeld bij de bespreking van de elementgroeiratio, heeft de resolutie in die nauwe gebieden invloed op grote delen van de mesh. Voor de invoerparameter die de resolutie in nauwe gebieden bepaalt is elk scalair getal toegelaten. In ons geval gaat het echter typisch om waarden tussen 0 en 1. Door de zeer fijne geometrie met zeer grote schaalverschillen leidt dit in het algemeen toch nog tot een zeer fijne mesh, enkel de nauwste rand wordt gemesht met elementen die anisotroop, dus zeer langwerpig, gevormd zijn. Een analyse van de invloed van de meshfijnheid op de uiteindelijk bekomen gekoppelde flux toont echter aan dat in ons geval de invloed zeer sterk beperkt blijft. Een randresolutie van 0.5 in plaats van 1 levert typisch slechts een relatieve afwijking op van de orde 10−4 . Ook de algemene meshfijnheid, zoals bepaald door middel van de eerste invoerparameter, heeft slechts een effect van deze orde. Op de snelheid van de berekeningen heeft de grootte van de mesh echter wel een grote invloed. De derde invoerparameter voor de mesh bepaalt de maximaal toegelaten grootte van de mesh. Als maatstaf werd gekozen voor het aantal elementen. Op de computer waar de berekeningen werden uitgevoerd, die nochtans over een uitgebreid geheugen beschikt, was een mesh van 250000 elementen het maximum toelaatbare om de berekeningen te kunnen uitvoeren. Onze subroutine is voorzien van een lus waarbij een te grote mesh verworpen wordt. In dat


40

geval wordt de parameter ‘resolutie in nauwe gebieden’ met 0.1 verminderd, waarna opnieuw gemesht wordt. Deze procedure herhaalt zich tot een geschikte mesh bekomen wordt. Deze beveiliging is nodig omdat de meshgrootte niet alleen van de meshparameters afhangt, maar ook zeer sterk varieert met de geometrische parameters. Onder meer de rotorhoek heeft een sterke invloed, en de diktes van de aangetaste randen kunnen de mesh een orde 5 laten variëren.

4.4

Fysische eigenschappen

Eenmaal de geometrie gegenereerd is, moeten aan de verschillende delen ervan fysische eigenschappen toegekend worden. Hier ontstaat een praktisch probleem. Fysische kenmerken worden immers toegekend aan subdomeinen, en niet aan geometrische gebieden. Het nummeren van de subdomeinen gebeurt echter door FEMLab op een wijze die voor een complex model als het voorliggende niet zeer voorspelbaar is. Door de parametrische opvatting verandert de geometrie immers vrij, waardoor ook de nummering van de subdomeinen aangepast wordt. We kunnen dus niet zomaar een reeks materiaaleigenschappen vooropstellen, aangezien deze na wijziging van de geometrie aan de verkeerde subdomeinen zouden toegekend worden. FEMLab bevat wel een procedure om een geometrie te analyseren en de nummers van de subdomeinen te vergelijken met een voordien geanalyseerde geometrie. Zo hebben we een standaardgeometrie opgesteld waarvoor de volgorde van de subdomeinen gekend is. De verschillen tussen de voorliggende geometrie en de gestelde standaard worden dan geanalyseerd. De analyse-procedure blijkt bij dit complexe model echter geen foutvrij resultaat te garanderen. Er bestaat wel een mogelijkheid om de relaties tussen de geometrische gebieden en de fysische subdomeinen bekend te maken bij de verwerking van de geometrie. Hierbij wordt door middel van een matrix aangeduid welke subdomeinen deel uitmaken van een geometrisch gebied. Wij maken gebruik van deze methode, waarbij we ervoor zorgen dat de geometrische gebieden meteen overeenstemmen met een subdomein en dus niet overlappen. Op deze manier krijgen we een matrix met slechts één positief element op elke rij en elke kolom, zodat unieke relaties bestaan tussen de nummers van de geometrische gebieden, zoals ingegeven, en de nummers van de subdomeinen. FEMLab is voorzien van applicatiemodes, die een kader bieden om op een eenvoudige manier fysische betekenis te geven aan een model. In dit geval werken we met de mode ‘Magneto-


41

statics’ binnen de module ‘Electromagnetics’. We hebben meer bepaald opgegeven dat het om externe stromen loodrecht op het vlak gaat, wat meteen ook bepaalt dat alle inducties en veldsterktes in het tweedimensionale vlak van de geometrie liggen. Wij moeten enkel de gepaste materiaaleigenschappen van het materiaal invullen, in casu de relatieve permeabiliteit µr , de elektrische conductiviteit σ en de externe stroomdichtheid Jz . Het materiaal is isotroop. Alle velden en eigenschappen worden constant verondersteld in de tijd. De wetten van Maxwell en de constitutieve wetten worden opgelost:

∇·B =0

(4.1)

∇×H =J

(4.2)

B =µ0 µr H

(4.3)

Uit de eerste vergelijking weten we dat er een magnetische vectorpotentiaal bestaat zodat

B =∇×A

(4.4)

Uit deze vergelijkingen kunnen we de Poisson-vergelijking afleiden:

∇ × µ−1 ∇ × A = J

(4.5)

Het systeem wordt opgelost naar de vectorpotentiaal. Aangezien we ervan uitgaan dat alle stromen loodrecht op het vlak staan, en B en H dus in het vlak liggen, weten we dat de vectorpotentiaal A loodrecht op het vlak staat (en dus enkel een z-component heeft) en dus als een scalair kan weergegeven worden, A = Az = kAk, net als de stroomdichtheid J = Jz = kJ k.

∂ ∂x

−1

µ

∂ ∂ −1 ∂ A(x, y) + µ A(x, y) = −J(x, y) ∂x ∂y ∂y ∇ · µ−1 ∇A = −J

(4.6) (4.7)

De magnetische eigenschappen van het aangetaste materiaal in de vijf randen en deze van het onaangetaste materiaal werden besproken in paragraaf 2.4.2. Uit de formule voor de B-H-karakteristiek volgt voor de relatieve permeabiliteit, met B = kBk:

B0 µr = · µ0 H0

1+

B B0

ν−1 !−1 (4.8)


42

Dit is opnieuw een functie van de magnetische inductie B. De formule kan in deze vorm in FEMLab ingegeven worden. De elektrische conductiviteit van het materiaal wordt uit de materiaaleigenschappen in tabel 2.1 gehaald. De elektrische stroomdichtheid in de gebieden die de bekrachtigingswindingen vertegenwoordigen wordt als volgt berekend. Het aantal windingen van de bekrachtigingsspoel wordt berekend zoals in het ruwe model, paragraaf 3.2.3. De stroomdichtheid J is dan, met Sw de wikkelingsoppervlakte:

J=

w·I Sw

(4.9)

Het is duidelijk dat het toekennen van de fysische eigenschappen voor elke berekening opnieuw moet gebeuren, dus voor elke geometrie, elke hoek en elke stroom. Praktisch worden enkel verschillende soorten materiaal gedefinieerd. Dan moet via een vector van indices elk subdomein aan een bepaald type materiaal gekoppeld worden. Deze vector wordt bij het uitvoeren van het programma opgesteld aan de hand van de materiaaleigenschappen die aan de geometrische gebieden moeten gekoppeld worden, de vaste volgorde van de geometrische gebieden en de unieke relaties tussen de geometrische gebieden en de fysische subdomeinen. Verder moeten ook de randvoorwaarden worden ingesteld. We hebben een ‘frame’ voorzien dat als omgeving van het magnetische circuit dient. Enkel de buitenranden van dit frame moeten dus behandeld worden. De rand werd gedefinieerd als een magnetische isolatie, wat neerkomt op een Dirichlet-randvoorwaarde. Dit betekent dat A over de rand constant is, of dat het magnetische veld raakt aan de rand. Meer bepaald wordt A gelijk gesteld aan nul, aangezien de vectorpotentiaal op een constante na willekeurig kan gekozen worden. Om een natuurlijke overeenkomst met onze ronde motor te krijgen kenden we het frame een cirkelvorm toe. Deze randvoorwaarde is echter tevens de standaardwaarde in FEMLab, dus moest het toekennen ervan niet expliciet gebeuren en hoefde hier geen rekening gehouden te worden met het verband tussen geometrische en fysische randen. Figuur 4.4 geeft de magnetische vectorpotentiaal A weer. Figuur 4.5 is een ietwat artistieke afbeelding van B, de norm van de magnetische inductie.

4.5 4.5.1

Veldberekeningen Differentiaalvergelijkingen

Zoals vermeld in paragraaf 4.4 moet het volgende stelsel van partiële differentiaalvergelijkingen worden opgelost:


43

Figuur 4.4: De magnetische vectorpotentiaal

Figuur 4.5: De norm van de inductie

∇ · µ−1 ∇A = −J

(4.10)

We kiezen echter voor de ‘algemene’ formulering van de vergelijkingen. Dit betekent dat als basis aan het algoritme het volgende stelsel van vergelijkingen ligt:


44

∇ · Γ = −J

(4.11)

A = 0,

δΩ

(4.12)

Hierbij is:

Γ = µ−1 ∇A

(4.13)

en δΩ is de rand van het domein. De andere mogelijkheid was de ‘coëfficiënt’-vorm. Hierbij wordt A expliciet als variabele gebruikt in plaats van Γ . Deze vorm is echter minder geschikt voor sterk niet-lineaire problemen [12].

4.5.2

Variationele vorm

Vermenigvuldigen we beide leden van vergelijking (4.11) met een willekeurige testfunctie v en integreren we over het volledige domein van het model Ω , dan krijgen we:

Z (∇ · Γ + J) v dS = 0

(4.14)

Ω

Partiële differentiatie (met behulp van de functie van Green) geeft:

Z

Z (Γ · ∇v − Jv) dS −

Ω

(n · Γ v) ds = 0

(4.15)

δΩ

Hierbij is n de uitwendige normaalvector op de buitenranden van het domein. We stellen de testfunctie v echter gelijk aan 0 op δΩ . Hierdoor valt deze tweede term weg. Er geldt nog steeds dat A = 0 op de buitenranden. Vermenigvuldigen we ook deze vergelijking met een testfunctie w en integreren we de bekomen uitdrukking over de rand, dan krijgen we de volgende uitdrukking.

Z Aw ds = 0 δΩ

Deze vergelijkingen worden de zwakke of variationele vorm van het systeem genoemd.

(4.16)


4.5.3

45

Discretisatie

De essentie van de eindige-elementenmethode is dat de variationele vergelijkingen geprojecteerd worden op een eindig-dimensionale functieruimte. De functies waarover het gaat zijn gewoonlijk (ruimtelijk één- of tweedimensionale) stuksgewijze veeltermfuncties, men spreekt over Lagrange-elementen. In ons geval wordt de werkelijke functie op elk mesh-element benaderd door een kwadratische functie. Een kwadratische functie wordt, in één dimensie, bepaald door drie punten. Daarom wordt op elke rand van een driehoekig mesh-element, dus tussen twee meshknopen in, een nieuwe Lagrange-knoop ingevoerd. Elke knoop betekent een vrijheidsgraad Ui in ons systeem, aan elke knoop wordt ook een basisfunctie gekoppeld. Basisfuncties verbonden aan een meshknoop zijn enkel verschillend van nul in deze meshknoop en in de meshelementen die deze meshknoop als hoekpunt hebben. Basisfuncties verbonden aan een knoop in het midden van een zijde zijn enkel verschillend van nul in deze knoop en in de twee meshelementen die deze zijde delen. Dit betekent dat in elk meshelement de interpolatiefunctie ontstaat als de lineaire interpolatie van zes basisfuncties. Bekijken we deze basisfuncties van iets dichter, zie figuur 4.6. Elke basisfunctie is dus gelijk aan 1 in de Lagrange-knoop waaraan ze verbonden is en gelijk aan nul in elke andere knoop. We definiëren lokale parameters ξ1 , ξ2 en ξ3 in een driehoekig meshelement, waarbij ξi = 1 op meshknoop i, ξi = 0 op de zijde recht tegenover meshknoop i en ξi lineair varieert tussen deze beide waarden. Dan is de basisfunctie ϕ die met meshknoop 1 verbonden is gelijk aan:

ϕ1 = ξ1 (2ξ1 − 1)

(4.17)

De basisfunctie die met Lagrangeknoop 4 verbonden is, heeft als vergelijking:

ϕ4 = 4ξ1 ξ2

(4.18)

Het feit dat we met kwadratische lagrange-elementen werken in plaats van met de traditionele lineaire elementen heeft twee gevolgen. Enerzijds stijgt het aantal vrijheidsgraden uiteraard, ongeveer met een factor 2, wat tot extra geheugenvereisten en rekentijd leidt. Anderzijds is de nauwkeurigheid ook groter. Men kan aantonen dat de benaderingsfout van hogere orde wordt naarmate de graad van de stuksgewijze veeltermfuncties stijgt [13]. Op deze manier wordt A benaderd door:


46

Figuur 4.6: De twee soorten basisfuncties

A≈

X

Ui ϕi

(4.19)

i

Op de rand van het domein definiëren we aparte, eendimensionale basisfuncties ψn die zullen dienen om de randvoorwaarden te evalueren. Deze zijn op dezelfde manier gedefinieerd als de gewone basisfuncties, er wordt ook voor deze basisfuncties gekozen voor kwadratische Lagrange-elementen. De lokale parameters zijn hier uiteraard eendimensionaal opgevat, per element zijn er ook slechts twee lokale parameters gedefinieerd, namelijk de lineaire afstand tot de beide eindpunten van het element, en drie kwadratische Lagrange-elementen. Op de rand van het domein kunnen we dus stellen:

A=

P

Un ψn , δΩ

(4.20)

n

Als testfuncties vk kiezen we de basisfuncties ϕk . De testfuncties op de rand wm kiezen we gelijk aan de basisfuncties ψm . We krijgen als nieuwe vergelijkingen, met τ de integratie over alle Gauss-punten (dit zijn de punten waarin Un gedefinieerd is, dus de Lagrange-knopen) en xj de j-de ruimtelijke coördinaat:

Z X τ

Γj

j

∂φk ∂xj

− Jφk dS = 0

(4.21)

Z Aψm ds = 0

(4.22)

δτ

∀ m, k Deze tweede vergelijking behandelen we vormelijk verder; gezien de aard van dit probleem komt deze er echter gewoon op neer dat de variabele A gelijk aan nul is in elke lagrangeknoop


47

op de rand. Samengevat moet het volgende niet-lineaire systeem van partiële differentiaalvergelijkingen worden opgelost:

L(U ) = 0

(4.23)

M (U ) = 0

(4.24)

Met L en M functies van de vector U met alle vrijheidsgraden Ui . Hierbij houden we vergelijkingen (4.13), (4.19) en (4.20) in gedachten. We krijgen specifiek:

Lk =

Z X τ

j

Z Mm = δτ

∂φk Γj ∂xj

X

− Jφk dS

Un ψn ψm ds

(4.25) (4.26)

n

en met de definitie van Γ zoals in vergelijking (4.13). De integralen in de definities van L en M worden berekend met behulp van kwadratuurformules (zie ook de formule van Simpson in paragraaf 3.4). Deze methode komt erop neer dat een integraal over een meshelement benaderd wordt door de gewogen som van de waarden in een eindig aantal punten in dit meshelement.

4.5.4

Niet-lineaire solver

Het gegeven stelsel van niet-lineaire stationaire partiële differentiaalvergelijkingen wordt opgelost door een gedempte Newton-methode. Bij elke iteratie k wordt, uitgaand van de huidige waarde U (k) , een correctie (∆U )(k) berekend door het oplossen van het gelineariseerde systeem:

K (U (k) ) (∆U )(k) = L(U (k) )

(4.27)

N (U (k) ) (∆U )(k) = M (U (k) )

(4.28)

Hierbij zijn K en N gelijk aan de jacobianen van L en M in U (k) :

∂L(U (k) ) ∂U ∂M (U (k) ) N (U (k) ) = ∂U K (U (k) ) =

(4.29) (4.30)


48

Deze jacobianen worden analytisch berekend, aan de hand van de waarden van onder meer µ(B(U (k ) )). Er is ook een mogelijkheid om de jacobianen numeriek te berekenen, indien er bijvoorbeeld geen analytische uitdrukkingen voorhanden zijn. Dit kost echter veel rekentijd. Verder is het mogelijk om bepaalde herberekeningen van de jacobiaan te vervangen door Broyden-updates, zie ook paragraaf 6.3. De stap die bij elke iteratie wordt gezet is:

U (k+1) = U (k) + α(k) (∆U )(k)

(4.31)

Hierbij is α(k) een dempingsfactor die per stap wordt bepaald. De iteratie gaat door tot een voldoende nauwkeurige oplossing bekomen wordt. De vereiste nauwkeurigheid kan opgegeven worden aan het algoritme. Hierbij kan best ook rekening gehouden worden met de opgegeven fijnheid van de mesh: het heeft geen zin een zeer nauwkeurige berekening uit te voeren als door een ruwe mesh de oplossing toch een relatief grote afwijking vertoont. Een nauwkeurigheid van 10−3 blijkt toch een nodige voorwaarde te zijn om een bruikbare koppelkarakteristiek te bekomen. Als α = 1 is deze procedure de gewone Newton-methode. Deze convergeert kwadratisch als de initiële waarde voldoende dicht bij de oplossing ligt. Als dit niet het geval is wordt door een gepaste keuze van de dempingsfactoren toch convergentie bereikt. Als startwaarde van het algoritme worden in de procedure twee mogelijkheden voorzien. Als de berekening voor dezelfde geometrie al eens uitgevoerd is, bijvoorbeeld bij een andere bekrachtigingsstroom, wordt de bekomen waarde als startwaarde opgegeven. Om de integratie van de flux mogelijk te maken in het ruwe model en de compatibiliteit te bewaren, zie paragraaf 3.4, moeten de stromen in stijgende volgorde opgegeven worden. Dit levert ook hier voordelen op aangezien de startwaarde hierdoor telkens de oplossing vrij dicht benadert. Als de berekening voor een bepaalde geometrie voor de eerste keer wordt uitgevoerd, wordt als startwaarde een functie opgegeven met niet-nul afgeleide, bijvoorbeeld x2 + y 2 . Als we dit niet zouden voorzien zou de startwaarde gewoon constant zijn, wat, aangezien de coëfficiënten van de differentiaalvergelijkingen afhankelijk zijn van de gradiënt van deze waarde, tot een singuliere jacobiaan zou leiden en uiteindelijk tot een foutmelding.

4.5.5

Lineair subprobleem

Voor elke iteratie van de niet-lineaire solver moet dus het volgende lineair probleem worden opgelost naar ∆U .


49

K (U0 ) ∆U = L(U0 )

(4.32)

N (U0 ) ∆U = M (U0 )

(4.33)

Het oplossen van het lineaire subprobleem gebeurt met behulp van de lineaire solver. Hiervoor moet een algoritme geselecteerd worden. We hebben ervoor geopteerd deze keuze in de procedure vast te leggen en dus niet uitwendig bereikbaar te maken. De complexiteit van het probleem maakt het vinden van een oplossing met beperkte geheugenvereisten immers niet triviaal. Als algoritme voor het oplossen van het lineaire subprobleem gebruiken we de Cholesky-decompositie. We bespreken deze hier kort. Stel dat het volgende lineaire probleem moet worden opgelost, we gebruiken hier algemene benamingen:

Ax = b

(4.34)

A = LLT

(4.35)

Eerst wordt de matrix A ontbonden als:

Hierbij is L een matrix met enkel elementen verschillend van nul op en onder de diagonaal. Vervolgens worden de volgende twee stelsels van vergelijkingen opgelost, naar y en uiteindelijk x:

Ly = b

(4.36)

LT x = y

(4.37)

Deze twee problemen zijn triviaal op te lossen aangezien de matrices driehoekig zijn (en de elementen dus sequentieel bepaald kunnen worden). Er zijn enkele verschillende doch gelijkaardige algoritmes om de Cholesky-decompositie uit te voeren. We lichten kort het Cholesky-Banachiewicz-algoritme toe. Als we de vergelijking A = LLT uitschrijven krijgen we de volgende uitdrukkingen voor de elementen van L:


Li,j =

50

1 Lj,j

Ai,j −

j−1 P

Li,k Lj,k

, i>j

(4.38)

k=1

s Li,i =

!

Ai,i −

i−1 P k=1

L2i,k

(4.39)

Een element van de matrix L kan dus bepaald worden als alle elementen links ervan en erboven gekend zijn. Het Cholesky-Banachiewicz-algoritme begint met het eerste element L1,1 en berekent alle elementen vervolgens rij per rij. Voor dit gebeurt worden de rijen wel op een andere manier geordend om de efficiëntie te verhogen, het gaat immers om ‘sparse’ matrices. Om de Cholesky-decompositie uit te voeren moet de matrix symmetrisch en positief-definiet zijn. In ons geval gaat het in het bijzonder om de matrix K (we stappen opnieuw over op onze opgebouwde nomenclatuur). De vergelijking met matrix N is immers triviaal, zoals we al aangestipt hebben. De vereiste dat de matrix positief-definiet moet zijn ligt aan de basis van de startwaardes, zoals besproken in paragraaf 4.5.4. Symmetrie van de systeemmatrix is in normale omstandigheden voldaan bij het oplossen van elektromagnetische vergelijkingen.

4.6

Koppelberekening

Bij het berekenen van het koppel kunnen we dezelfde methode toepassen als in het ruwe model. Dit betekent dat we de gekoppelde flux integreren om zo de co-energie te bekomen en dat we deze dan afleiden om het koppel te verkrijgen. Het nauwkeurig integreren van de co-energie is echter niet haalbaar, aangezien wanneer een klein aantal bekrachtigingsstromen opgegeven wordt, de berekening bij extra tussenliggende stroomwaarden moet worden uitgevoerd om een goede integratie te bekomen. Dit is gezien de grote kost van een berekening niet verdedigbaar. De co-energie kan echter wel rechtstreeks worden berekend (zie paragraaf 3.4):

Z

Z JA − wm dV =

Wco = V

(JA − wm )Lm dS

(4.40)

S

De integraal van J · A berekenen we door integratie over de subdomeinen van de bekrachtigingswikkelingen. Aangezien J constant is integreren we enkel A. Dit heeft als voordeel dat deze berekening meteen ook aanzien kan worden als een berekening van de fysische flux, na deling door de wikkelingsoppervlakte en rekening houdend met het feit dat de integratie beide wikkelingen omvat. We willen de gekoppelde flux immers eveneens beschikbaar als uitgang.


51

De magnetische energiedichtheid wordt berekend als:

Z

B

H · dB

wm =

(4.41)

0

De aldus bekomen co-energie wordt zoals in paragraaf 3.4 gedifferentieerd om het koppel te bekomen.

4.7

Resultaten en conclusie

We hebben een nauwkeurig eindige-elementenmodel ontworpen voor het oplossen van het directe probleem van een elektrische machine. Dit model zal verder gebruikt worden door de space-mappingprocedure om tot een optimalisatie te komen. De rekentijd van het fijne model varieert (voor één bekrachtigingsstroom) op een snelle computer van een half uur tot enkele uren, afhankelijk van de vereiste nauwkeurigheid, van de meshfijnheid en van de niet-lineariteit of de grootte der bekrachtigingsstroom. Het model is betrekkelijk betrouwbaar, maar aangezien enkele procedures inwendig in het programma FEMLab geprogrammeerd zijn is het niet mogelijk om volledige betrouwbaarheid te controleren. In het bijzonder bij het meshen werd een geval gemeld dat het programma er niet in slaagde een oplossing te vinden. Verder wordt, met aanvaardbare parameters, wel steeds naar een oplossing geconvergeerd. Gezien de complexiteit van het model is wel een performante computer nodig om het model uit te voeren.

Hoofdstuk 5

Directe optimalisatie 5.1

Inleiding

We hebben een procedure ontworpen die de geometrische parameters van de motor optimaliseert aan de hand van één model. Deze wordt gebruikt om een optimalisatie uit te voeren gebruik makende van evaluaties in het ruwe model. De optimalisatie is direct, dit betekent dat de evaluaties in het model geforceerd worden om te convergeren naar een vooropgestelde waarde, al kan dit een onbereikbare waarde zijn als bijvoorbeeld een maximalisatie van het koppel uitgevoerd wordt. Deze procedure is van belang bij het uitvoeren van de zogenaamde parameter-extractie in de space mapping, evenals bij het bepalen van het absolute optimum in het ruwe model. In theorie kan deze optimalisatieprocedure ook uitgevoerd worden wanneer we willen optimaliseren in het fijne model. Dat zou dan meteen tot een nauwkeurig resultaat leiden. Aangezien echter vaak honderden evaluaties van het gebruikte model nodig zijn, in het bijzonder als naar meerdere parameters tegelijk geoptimaliseerd wordt, zou de toepassing van deze procedure op het fijne model tot een ontoelaatbare (minstens 1 week) rekentijd leiden. De vijf parameters waarnaar we optimaliseren zijn de diameter van de luchtspleet Dls , de breedte van stator- en rotortand bst en brt , de inwendige diameter van het statorjuk Ds en de uitwendige diameter van het rotorjuk Dr . Als invoer leggen we bepaalde beperkingen op. Deze worden gebruikt om alle variabelen tussen bepaalde grenzen te laten variëren, maar ook om één of meer parameters a priori te definiëren. Op deze manier kunnen we er voor kiezen om slechts één variabele vrij te laten. Dit speciale geval wordt met behulp van een bijzonder algoritme aangepakt. In dit hoofdstuk bekijken we eerst hoe we de procedure ge¨ımplementeerd hebben. Hierbij 52

HOOFDSTUK 5. DIRECTE OPTIMALISATIE

53

komen het geval met één variabele en met meerdere variabelen tegelijk aan bod. Vervolgens gaan we in op de algoritmes die gebruikt worden voor de beide gevallen.

5.2 5.2.1

Implementatie Minimalisatie

Het optimaliseren van een functie komt neer op het minimaliseren van een kostenfunctie, die de afwijking van deze functie ten opzichte van een bepaalde streefwaarde weergeeft. Deze minimalisatie moet binnen bepaalde grenzen gebeuren, aangezien het mogelijk moet zijn om de modellen te evalueren en de aangezien resultaten een fysische betekenis moeten hebben. In onze implementatie hebben de begrenzingen een plaats gevonden in de algemene invoerstructuur die alle procedures gemeenschappelijk hebben. Er zijn verschillende invoerparameters die verwacht worden. Ten eerste zijn er de numerieke onder- en bovengrenzen van de parameters, lo ≤ x ≤ lb . Het is aan de hand van deze grenzen dat achterhaald wordt of er naar meer dan één variabele geoptimaliseerd wordt. Om dus naar één variabele te optimaliseren moeten vier elementen van de onder- en bovengrens identiek zijn. Als blijkt dat slechts één variabele voorkomt definiëren deze begrenzingen meteen de hele parameterruimte. Als we echter naar meer dan één variabele optimaliseren, verwachten we nog andere begrenzingen. De matrix A en de vector B leggen de volgende ongelijkheid op:

A·x ≤B

(5.1)

met x de vector der parameters. Hierdoor kunnen we relaties tussen de verschillende afmetingen opleggen. We kunnen bijvoorbeeld vastleggen dat de rotor altijd kleiner moet zijn dan de stator, of dat dit verschil minstens een bepaalde, eindige waarde moet bedragen. De standaardverhoudingen zijn uiteraard reeds voorgedefinieerd in het initialisatiescript voor de invoerstructuur. Verder bepalen de matrix Aeq en de vector Beq de eventuele gelijkheden. De vorm is:

Aeq · x = Beq

(5.2)

Deze vergelijkingen worden standaard niet gebruikt, er worden enkel nulmatrices gedefinieerd. Ten slotte is er ook nog de mogelijkheid om niet-lineaire begrenzingen op te leggen,


54

hiertoe moet een procedure worden opgegeven die de begrenzingen geeft als functie van de parameters. Wanneer meer dan één variabele voorkomt, verwacht de procedure ook een beginwaarde voor de variabelen. De functie voorziet twee mogelijkheden. Enerzijds kan de startwaarde worden opgegeven als argument, om het mogelijk te maken verschillende opeenvolgende optimalisaties snel uit te voeren. Als echter geen expliciete beginwaarde opgegeven wordt zal in de invoerstructuur voor de parameters gezocht worden naar een startwaarde. Deze kan tevens standaard ge¨ınitialiseerd worden. Het is ook mogelijk om, als we naar meerdere variabelen optimaliseren, geen startwaarde op te geven maar de procedure zelf startwaarden te laten genereren. Hierbij moeten we enkel opgeven worden hoeveel waarden gegenereerd moeten worden. Deze worden dan willekeurig verdeeld tussen de boven- en ondergrenzen en eventueel aangepast om ook aan de relatieve begrenzingen te voldoen. Deze aanpassing gebeurt met de procedure die in paragraaf 6.4 uit de doeken gedaan wordt. Vervolgens wordt voor elke startwaarde een optimalisatie uitgevoerd, zodat een aantal sets ‘optimale’ parameters bekomen worden. Er bestaat immers een kans dat de procedure convergeert naar een lokaal minimum. Daarom wordt van de bekomen sets parameters deze gekozen die de kleinste kostenfunctie oplevert. Daarmee kan meestal aangenomen worden dat een absoluut minimum bereikt is. Het minimaliseren zelf gebeurt met behulp van interne functies van Matlab. Wanneer meerdere variabelen voorkomen gaat het om de functie ‘fmincon’, bij één variabele wordt naar de functie ‘fminbnd’ gegrepen. De precieze werking van deze procedures wordt in paragraaf 5.3 in detail besproken. Verder kan als optie ingegeven worden of na elke iteratie een output naar het scherm moet verstuurd worden. Ook de minimumwaarde van de evaluatie van de kostenfunctie en de minimumwaarde voor de stap van de variabelen kunnen als optie worden ingesteld.

5.2.2

Kostenfunctie

De kostenfunctie is de functie die effectief geminimaliseerd wordt. Het is een functie die de afwijking aangeeft van de evaluatie van het model met de huidige waarden voor de parameters ten opzichte van een zekere streefwaarde. We maken gebruik van het kwadraat van de norm van deze afwijking. De waarde waarmee vergeleken wordt, waarnaar dus geoptimaliseerd wordt, kan opgegeven worden. Dit kan gebruikt worden om te optimaliseren naar een bepaald koppelprofiel, of om


55

het koppel te maximaliseren door een grote waarde in te geven. Er is tevens de mogelijkheid om rekening te houden met de koppelrimpel. Dan is de afwijking voor een bepaalde bekrachtigingsstroom gedefinieerd als de gewogen som van enerzijds de afwijking van het gemiddeld koppel tot een streefwaarde en anderzijds, bij elke rotorhoek, de absolute waarde van het verschil tussen het gemiddelde koppel en het werkelijke koppel. Als slechts één bekrachtigingsstroom opgegeven wordt, wordt slechts één koppelprofiel geoptimaliseerd. Als meerdere stromen opgegeven worden, dan zullen alle waarden meegeteld worden in de norm om de kost te berekenen.

5.3 5.3.1

Algoritmes Meerdere variabelen

Algemeen Het minimaliseren van een niet-lineaire en niet-kwadratische scalaire functie f (x ), met begrenzingen

gi (x ) = 0, i = 1 . . . me

(5.3)

gi (x ) ≤ 0, i = me + 1 . . . m

(5.4)

werd traditioneel aangepakt door sequentieel onbegrensde problemen op te lossen, waarbij er kostenfuncties voorzien werden voor overtredingen van de begrenzingen. Op deze manier werd geconvergeerd naar een oplossing die binnen de grenzen lag. Deze methodes blijken echter relatief inefficiënt te zijn. Daarom wordt heden vooral gefocust op het oplossen van de zogenaamde Kuhn-Tucker-vergelijkingen die het oorspronkelijke probleem anders uitdrukken. We zoeken x ∗ waarvoor:

∇f (x ∗ ) +

m X

λ∗i ∇gi (x ∗ ) = 0

(5.5)

i=1

λ∗i · gi (x ∗ ) = 0,

i = 1 . . . me

(5.6)

i = me + 1 . . . m

(5.7)

gi (x ) = 0,

i = 1 . . . me

(5.8)

gi (x ) ≤ 0,

i = me + 1 . . . m

(5.9)

λ∗i ≥ 0,


56

Als het probleem convex is, als de tweede afgeleiden dus positief zijn, zijn de Kuhn-Tuckervergelijkingen een nodige en voldoende voorwaarde voor een globale oplossing. In de eerste vergelijking is gezorgd voor een neutralisering van de gradiënten van de te minimaliseren functie en de actieve begrenzingen. Enkel de actieve begrenzingen mogen hierin voorkomen, de Lagrange-vermenigvuldigers λi van de niet-actieve begrenzingen moeten dus nul zijn. Dit wordt impliciet gesteld in de tweede, derde en laatste vergelijking. De twee laatste vergelijkingen zijn de bekende begrenzingen.

Sequential Quadratic Programming De meeste algoritmes proberen uit de Kuhn-Tucker-vergelijkingen direct de Lagrange vermenigvuldigers te berekenen. Begrensde quasi-Newton methodes-leiden tot een meer dan lineaire convergentie door tweede-orde-informatie van de vergelijkingen in acht te nemen. Deze methodes worden gewoonlijk Sequential Quadratic Programming (SQP) methodes genoemd, aangezien herhaaldelijk een kwadratisch subprobleem opgelost wordt [14]. Dit is gebaseerd op een kwadratische benadering van de Lagrange-functie L(x , λ):

L(x , λ) = f (x ) +

m X

λi · gi (x )

(5.10)

i=1

Het k-de subprobleem is de minimalisatie naar d (een vector van dezelfde dimensie als de vector der variabelen) van de volgende functie:

T 1 min d T H (k) d + ∇f (x (k) ) d x 2 T ∇gi (x (k) ) d + gi (x (k) ) = 0, T ∇gi (x (k) ) d + gi (x (k) ) ≤ 0,

met :

(5.11)

i = 1 . . . me

(5.12)

i = me + 1 . . . m

(5.13)

Dit is de kwadratische benadering van de te minimaliseren functie, rond het punt x (k) . De matrix H (k) is een positief definiete benadering van de hessiaan van de Lagrange-functie. De hessiaan van een (scalaire) functie is de matrix der tweede afgeleiden naar de variabelen. De opstelling van deze matrix bespreken we in de volgende paragraaf. De nieuwe waarde van x na elke oplossing van het subprobleem wordt gegeven door:

x (k+1) = x (k) + α(k) d (k)

(5.14)


57

Hierbij is α(k) een dempingsfactor. Hierover verder meer. De procedure stopt als de stap ∆x kleiner wordt dan een opgegeven tolerantiewaarde, de stap van de objectieffunctie ∆f (dus de kostenfuntie) klein wordt of een op te geven maximaal aantal iteraties overschreden wordt. Het volledige verloop van de procedure is weergegeven in figuur 5.1. Meer informatie over deze optimalisatiemethode is te vinden in [15] en [16].

Oplossing subprobleem Het oplossen van het kwadratische subprobleem gebeurt in drie stappen. Eerst wordt de hessiaan H bijgewerkt. Vervolgens wordt het feitelijke ‘Quadratic Programming’-subprobleem opgelost om de richting van de stap te bepalen en ten slotte wordt een interne kostenfunctie opgesteld die toelaat de staplengte te bepalen. De hessiaan wordt bij elke iteratie bijgewerkt volgens de BFGS-methode (Broyden-FletcherGoldfarb-Shanno) (zie [17]):

H

(k+1)

=H

(k)

T T q (k) q (k) H (k) s (k) s (k) H (k) + − T T q (k) s (k) s (k) H (k) s (k)

(5.15)

met:

s (k) = x (k+1) − x (k)

(5.16)

q (k) = ∇L(x (k+1) , λ(k+1) ) − ∇L(x (k) , λ(k) )

(5.17)

Hierbij is λ(k) de k-de schatting van de vector der Langrange-vermenigvuldigers. De hessiaan wordt in dit algoritme echter positief definiet gehouden. Hiertoe moet de startT waarde positief definiet zijn en moet bij elke update q (k) s (k) positief zijn. Wanneer dit laatste niet het geval is wordt q (k) element per element bijgewerkt totdat aan de voorwaarde voldaan is. Hierbij wordt als eerste stap het meest negatieve element van het elementgewijze T T product q (k) . ∗ s (k) gehalveerd. Dit wordt herhaald tot q (k) s (k) kleiner wordt dan T 10−5 . Als deze methode niet tot een positieve waarde q (k) s (k) leidt wordt een tweede methode gebruikt. In dat geval wordt een vector v , vermenigvuldigd met een scalair w opgeteld bij q :

(k)

qmod = q (k) + w(k) v (k) De elementen van de vector v worden gedefinieerd als (met i = 1 . . . m):

(5.18)


Figuur 5.1: Het optimalisatie-algoritme

58


(k)

vi

59

= ∇gi (x (k+1) ) · gi (x (k+1) ) − ∇gi (x (k) ) · gi (x (k) )

(5.19)

wanneer voldaan is aan:

q (k) · w < 0, i (k) · s (k) < 0 q i

i

(5.20) (5.21)

In elk ander geval is vi = 0 Dan wordt w stelselmatig verhoogd tot q (k)

T

s (k) positief wordt. Als gevraagd wordt om

de iteraties van de optimalisatie op het scherm weer te geven wordt ook weergegeven of de hessiaan bijgewerkt is om positief-definiet te worden en of dat enkel met de eerste of ook met de tweede methode gebeurd is. Als deze wijzigingen noodzakelijk zijn is dat gewoonlijk een indicatie dat het probleem zeer niet-lineair is en dat de convergentie dus misschien trager dan verwacht zal verlopen. Vervolgens wordt, bij elke iteratie, het feitelijke subprobleem opgelost. Noemen we in vergelijking (5.11):

c = ∇f (x (k) )

(5.22)

1 q(d ) = d T H (k) d + c T d 2

(5.23)

A = ∇g (x (k) )

(5.24)

b = g (x (k) )

(5.25)

d is de kleine afwijking ten opzichte van x (k) in het gekwadratiseerde probleem. Ai is de i-de rij van de matrix A, bi is het i-de element van de vector b. q(d ) is dus de te minimaliseren functie, de oplossing is d (k) . De begrenzingen worden uitgedrukt als:

Ai d = bi

i = 1 . . . me

(5.26)

Ai d ≤ bi

i = me + 1 . . . m

(5.27)

Het oplossen van het subprobleem gebeurt met een zogenaamde active set methode, zoals beschreven in [18]. Dit is een iteratief proces. Als we hier spreken over een iteratie zal dus


60

een iteratie binnen de procedure die het subprobleem oplost bedoeld worden. In die context wordt de index n gebruikt. De matrix A(n) is een schatting van de ‘active set’. Dit is een schatting van de begrenzingen die actief zijn, deze waarbij de huidige oplossing effectief op de grens ligt. De matrix wordt bij elke iteratie bijgewerkt en wordt gebruikt om een basis te vormen voor de zoekrichting dˆ (n) . We gebruiken deze notatie om een verschil te maken met d (k) , de oplossing van deze algemene iteratie van de SQP methode. De zoekrichting dˆ (n) wordt berekend om de functie te minimaliseren, maar hierbij wordt op de actieve begrenzingen gebleven. Ze worden dus niet overschreden maar ook niet verlaten. Hiertoe wordt dˆ (n) opgebouwd in een ruimte met basis Z (n) . De kolommen van Z (n) staan orthogonaal op de schatting van de active set A(n) . Elke zoekrichting die opgebouwd wordt als een lineaire combinatie van de kolommen van Z (n) zal dus op de actieve grenzen blijven. De matrix Z (n) wordt opgebouwd als de laatste m−l kolommen van de QR-decompositiematrix van A(n) , met l het aantal actieve begrenzingen. Dus:



Q T A(n)

T

=

R 0

 

Z (n) = Q:,l+1...m

(5.28) (5.29)

Vervolgens moet de zoekrichting dˆ (n) gezocht worden die q(d ) minimaliseert. Deze moet bestaan uit een lineaire combinatie van de kolommen van Z (n) . We kunnen dus schrijven:

dˆ (n) = Z (n) p (n)

(5.30)

Hierbij is p (n) een vector die dˆ (n) definieert in de ruimte Z (n) . We kunnen q nu in functie van p schrijven, zie vergelijking (5.23):

T 1 q(p) = p T Z (n) HZ (n) p + c T Z (n) p 2

(5.31)

Deze functie afleiden naar p geeft:

T T ∇q(p) = Z (n) HZ (n) p + Z (n) c

(5.32)

∇q(p) noemt met de geprojecteerde gradiënt van de kwadratische functie. Het is inderdaad de gradiënt van de functie, geprojecteerd in de subruimte met Z (n) als basis. Aangezien de


61

hessiaan H hier positief definiet is, wordt het minimum van de functie q(p) in de subruimte gedefinieerd door Z (n) gevonden als ∇q(p) = 0. De vector p (n) wordt gevonden door een stelsel lineaire vergelijkingen op te lossen. Met p (n) is meteen ook dˆ (n) = Z (n) p (n) bekend. Wanneer d ˆ(n) bekend is wordt (binnen het kwadratisch probleem) een stap gezet, waarbij y (n) de n-de benadering is van d (k) , de k-de stap in het SQP-algoritme.:

y (n+1) = y (n) + β (n) dˆ (n)

(5.33)

Als een stap langs dˆ (n) kan gezet worden met β = 1 is het minimum bereikt. Het kan echter zijn dat een nieuwe grens ontmoet wordt, die nog geen deel uitmaakt van de active set. De relatieve afstand tot de grenzen van het toelaatbare gebied in een bepaalde richting dˆ (n) is: ( β (n) = min i

(n)

−(Ai y (n) − bi ) (n) A dˆ (n)

) (5.34)

i

Dit quotiënt is gedefinieerd voor begrenzingen die niet in de active set zitten en waarvoor dˆ (n) in de richting van de grens is. Nu kan het zijn dat bepaalde grenzen moeten verlaten worden. Als l onafhankelijke begrenzingen in de active set zitten wordt een vector van Lagrange vermenigvuldigers λ(n) berekend zodat aan het volgende stelsel van lineaire vergelijkingen voldaan is:

A(n)

T

λ(n) = c

(5.35)

Als alle elementen van λ(n) positief zijn en β = 1 is y (n) de optimale oplossing van dit kwadratische probleem, d k = y (n) Als één element echter negatief is, en op de corresponderende component is geen gelijkheid van toepassing, wordt dit element uit de active set verwijderd en wordt een nieuwe iteratie gestart. Dit algoritme gebruikt als startwaarde gewoonlijk de vorige waarde van de SQP-procedure. Als deze echter niet aan de voorwaarden voldoet wordt een punt berekend dat aan de gelijkheids- en ongelijkheidsvoorwaarden voldoet.

Staplengte We keren terug naar het SQP-niveau. De term ‘iteratie’ wordt dus opnieuw toegepast op de SQP-procedure. De oplossing van het kwadratisch probleem is de vector d (k) . Hiermee wordt een stap gezet:


62

x (k+1) = x (k) + α(k) d (k)

(5.36)

De factor α wordt bepaald aan de hand van een interne kostenfunctie die per stap voldoende moet verminderd worden. Er wordt een startwaarde voor α genomen, waarna de kostenfunctie in het bekomen punt berekend wordt. Als de afname van de kostenfunctie ten opzichte van de waarde in α = 0 zeer klein is, is α te klein genomen. Is de kostenfunctie toegenomen, dan is α te groot genomen. In beide gevallen wordt α aangepast tot een geschikte waarde gevonden wordt. De kostenfunctie is:

Ψ(x ) = f (x ) +

me X

ri · gi (x ) +

i=1

m X

ri · max{0, gi (x )}

(5.37)

i=me +1

Met als strafparameter ri :

1 (r(k+1) )i = max λi , ((r(k) )i + λi ) , 2

i = 1...m

(5.38)

De startwaarde voor de strafparameter is:

ri =

k∇f (x )k k∇gi (x )k

(5.39)

Dit is een toepassing van de line search methode, zie [17].

5.3.2

E´ en variabele

Het algoritme dat gebruikt wordt om een functie te minimaliseren met slechts één vrije variabele wordt de methode van Brent genoemd. Meer informatie over het algoritme kan worden gevonden in [19] en [20]. Essentieel begint het algoritme telkens met het insluiten van het minimum in een zogenaamd ‘tripelpunt’. Dit is een set van drie punten a < b < c waarbij f (a) > f (b) en f (c) > f (b). Inderdaad zal, als een dergelijk tripelpunt bestaat en de functie continu is, zeker een minimum bestaan tussen a en c. Het punt b is te beschouwen als de huidige benadering voor het minimum. Het is de bedoeling bij elke iteratie deze punten dichter naar het minimum te brengen. Dit kan echter op twee manieren gebeuren. Een langzaam convergerend, maar stabiel algoritme met behulp van gulden secties wordt gecombineerd met een snel convergerend maar instabiel algoritme, de inverse parabolische interpolatie of parabolische benadering. In elke stap wordt er gecontroleerd of het mogelijk


63

is de inverse parabolische interpolatie toe te passen. Als dit niet mogelijk is wordt een guldensnede-iteratie toegepast. We bespreken eerst de parabolische benadering, inclusief de voorwaarden die gesteld worden om dit algoritme toe te passen. Vervolgens gaan we in op de guldensnedemethode.

Inverse parabolische interpolatie Met inverse parabolische interpolatie bedoelen we dat het minimum van een functie benaderd wordt door het minimum van een parabool die als interpolatie van de drie gekende punten (het tripelpunt) van de functie geldt. De interpolatie wordt invers genoemd omdat naar een abscis gezocht wordt in plaats van een ordinaat. De methode wordt weergegeven in figuur 5.2 [19]. Er wordt gestart met het tripelpunt P1 , P2 en P3 . Dit levert een parabool op. De functie wordt geëvalueerd in het minimum van de parabool, wat leidt tot punt P4 . De punten P1 , P4 en P2 definiëren een nieuwe parabool, die aan de basis ligt van de volgende stap. Aangezien f (P4 ) < f (P2 ) kan immers een tripelpunt gevormd worden zonder punt P3 . Vervolgens wordt opnieuw het minimum van de nieuwe parabool berekend, dit levert punt P5 op. Het nieuwe tripelpunt bestaat uit de punten P1 ,P5 en P4 .

Figuur 5.2: Parabolische benadering

Het minimum van de parabool die door de punten P1 (a, f (a)), P2 (b, f (b)) en P3 (c, f (c)) gaat is:

x=b−

1 (b − a)2 [f (b) − f (c)] − (b − c)2 [f (b) − f (a)] 2 (b − a) [f (b) − f (c)] − (b − c) [f (b) − f (a)]

(5.40)

Dit zien we in als volgt. Een parabool die zeker door de punten P1 , P2 en P3 gaat is de volgende:


f (x) = f (a)

(x − b)(x − c) (x − a)(x − c) (x − a)(x − b) + f (b) + f (c) (a − b)(a − c) (b − a)(b − c) (c − a)(c − b)

64

(5.41)

Aangezien een parabool uniek gedefinieerd wordt door drie punten is dit de parabool die als benadering gebruikt wordt. Om het minimum te zoeken berekenen we de afgeleide:

f 0 (x) = f (a)

(x − b) + (x − c) (x − a) + (x − c) (x − a) + (x − b) + f (b) + f (c) (a − b)(a − c) (b − a)(b − c) (c − a)(c − b)

(5.42)

Gelijkstellen aan nul en oplossen naar x geeft vergelijking (5.40). Om aanvaardbaar te zijn, en dus als stap uitgevoerd te worden, moet het resultaat van een parabolische iteratie in de eerste plaats binnen het initiële interval [a, b] liggen. De tweede voorwaarde is dat de stapgrootte kleiner moet zijn dan de helft van de voorlaatste stap. Deze vereiste stelt eenvoudig dat de procedure moet convergeren. Dat de stap vergeleken wordt met de voorlaatste stap in plaats van met de laatste is een heuristische stelregel. Het is ook niet toegelaten de functie te evalueren in een kleine omgeving van een punt dat reeds geëvalueerd is. Dit zou immers geen nieuwe significante informatie opleveren. Als aan één van deze voorwaarden niet voldaan is wordt een gulden-snedestap uitgevoerd.

Gulden snede De gulden-snedemethode is intu¨ıtief vergelijkbaar met de bissectiemethode voor het vinden van een nulpunt van een functie [21]. Ook hier gaan we echter uit van een tripelpunt met abscissen (a, b, c). Bij elke iteratie wordt een nieuw punt x gekozen, ofwel tussen a en b ofwel tussen b en c. Stel dat het nieuwe punt x tussen b en c ligt. Dan evalueren we f (x). Als f (b) < f (x) is het nieuwe tripelpunt (a, b, x), is f (b) > f (x) dan is het nieuwe tripelpunt (b, x, c). Het middelste punt van het tripelpunt is altijd het voorlopig beste minimum. We moeten dus beslissen waar het nieuwe punt x genomen wordt. Stel dat b op een fractie w van de afstand tussen a en c ligt.

b−a c−a c−b 1−w = c−a w=

(5.43) (5.44)

Stel nu dat het volgende punt x een fractie z (groter of kleiner dan 0) verder dan b ligt:


z=

x−b c−a

65

(5.45)

Het volgende tripelpunt is, zoals besproken, als x rechts van b ligt ofwel (a, b, x), ofwel (b, x, c). De lengte van het buitenste interval is dus w + z of 1 − w. Als x links van b ligt (z < 0) wordt het tripelpunt (a, x, b) of (x, b, c), met als lengtes w of 1 − w − z. Om de slechtste mogelijkheid te minimaliseren nemen we beide mogelijke lengtes gelijk. Dan is z = 1 − 2w. Op deze manier komt x vanzelf in het grootste van beide intervallen te liggen. We zien dat |b − a| = |x − c|, x ligt met andere woorden symmetrisch ten opzichte van b in het interval van het tripelpunt. We gaan er echter van uit dat ook in de vorige stap dezelfde strategie gevolgd is, dus om b te bepalen. Er geldt dus dat x dezelfde fractie van de afstand van b naar c is (we nemen even aan dat z > 0), als b van de afstand van a tot c. Dus:

w=

z 1−w

(5.46)

Met z = 1 − 2w geeft dit de kwadratische vergelijking w2 − 3w + 1 = 0, of voor w: √ 3− 5 w= ≈ 0.38197 2

(5.47)

In elk tripelpunt (a, b, c) zal b dus een relatieve afstand 0.38197 tot de ene grens hebben en 0.61803 tot de andere grens. Dit zijn de fracties van de gulden snede. Bij elke iteratie zal, indien het vorige interval eveneens aan deze verhouding voldoet, het tripelpuntinterval dus afnemen tot 61.803% van zijn vorige waarde. Als gestart wordt met een tripelpunt dat niet aan deze verhouding voldoet (zoals wanneer deze iteratie volgt op een iteratie met paraboolbenadering) zal door opeenvolgende punten te nemen volgens de vooropgestelde regel vanzelf geconvergeerd worden naar de guldensnedeverhouding.

5.4

Conclusie

We gebruiken een procedure die in staat is het absolute optimum van de geometrie van een motormodel te bepalen. We hebben deze routine theoretisch onderbouwd. Deze kan als subprocedure van space mapping gebruikt worden, maar kan ook gebruikt worden om enkel het ruwe model te optimaliseren. Het aantal vereiste iteraties van het model is zeer groot, in


66

het bijzonder als meerdere parameters tegelijk geoptimaliseerd worden. Een optimalisatie in het fijne model zal dus niet haalbaar zijn. De gebruiker kan kiezen voor een optimalisatie naar een bepaald koppelprofiel of voor een maximalisering van het koppelprofiel. Ook met de koppelrimpel kan rekening gehouden worden. Als naar meerdere variabelen tegelijk geoptimaliseerd wordt, kan de optimalisatie meerdere malen uitgevoerd worden bij willekeurige startwaarden die tussen de opgegeven grenzen liggen, om het vinden van lokale minima te vermijden.

Hoofdstuk 6

Space Mapping 6.1

Inleiding

De space-mappingtechniek, ge¨ıntroduceerd in 1994 door Bandler [22], is een optimalisatietechniek die oorspronkelijk gebruikt werd voor het ontwerp van radiofrequentie- en microgolfcomponenten. Directe optimalisatietechnieken, zoals besproken in hoofdstuk 5, maken gebruik van vele evaluaties van het model. Indien we een rekenintensief model gebruiken, kan de rekentijd aanzienlijk groot worden bij de optimalisatie. Dan is het idee ontstaan om in de optimalisatieprocedure een minder rekenintensief model te implementeren. De space mapping zorgt voor de manier waarop dit ruwe model ge¨ımplementeerd wordt in de optimalisatieprocedure van het fijne model zodanig dat er weinig evaluaties nodig zijn in het fijne model. Bij space mapping wordt immers het merendeel van het rekenwerk uitgevoerd met behulp van een ruw model zoals ontwikkeld in hoofdstuk 3. De nauwkeurigheid van het resultaat van de optimalisatie met behulp van space mapping is groter dan de nauwkeurigheid van het resultaat van de optimalisatie met het ruwe model, aangezien de optimalisatie ‘gevalideerd’ wordt door een fijn model. In dit hoofdstuk bespreken we eerst de algemene principes van space mapping. Vervolgens gaan we in op het principe en onze implementatie van het gebruikte ‘Aggressive Space Mapping’-algoritme. Het gebruik van dit algoritme garandeert echter nog geen convergentie. Daarom maken we gebruik van een hybride methode die we bespreken en implementeren.

67

HOOFDSTUK 6. SPACE MAPPING

6.2

68

Space mapping: principe

Intu¨ıtief is het principe van space mapping als volgt in te zien. Noemen we xf de vector der parameterwaarden van het fijne model, in de parameterruimte Ωf ⊂
xf∗ = arg min K(f (xf )) xf ∈Ωf

(6.1)

Het vinden van xf∗ is het einddoel van het algoritme. Het snelle, ruwe model noemen we r (xr ) : Ωr → <m met xr ∈ Ωr ⊂
xr∗ = arg min K(r (xr )) xr ∈Ωr

(6.2)

De onnauwkeurigheid van het ruwe model wordt opgevangen door niet met het werkelijke ruwe model te werken, maar wel met surrogaatmodellen, modellen die het fijne model benaderen maar op basis van het ruwe model werken. Deze surrogaatmodellen worden bekomen door een afbeelding of mapping p : Ωf → Ωr te zoeken tussen de parameterruimtes van de beide modellen. Het surrogaatmodel is dan gegeven door:

r (p(xf )) ≈ f (xf )

(6.3)

Het bepalen van de afbeeldingsfunctie p(xf ) noemen we de parameterextractie. Dit is de kern van het space-mappingalgoritme. De parameterextractie gebeurt voor een bepaalde waarde van xf door een snelle minimalisatie uit te voeren in het surrogaatmodel:

p(xf ) = arg min kr (p) − f (xf )k p∈Ωr

(6.4)

Hierbij is f (xf ) bepaald door middel van één evaluatie van het fijne model bij deze waarde xf . Als gevolg van het gebruik van het surrogaatmodel wordt de xf∗ bekomen in vergelijking (6.1) benaderd door door x f :


69

x f = arg min K(r (p(xf )))

(6.5)

xf ∈Ωf

Deze x f wordt bepaald met behulp van space mapping door p(xf ) − xr∗ = 0 op te lossen (k)

(k)

naar xf . Dit is een iteratief proces waarbij we in de k-de iteratie werken met xf . f (xf ) moet dan geëvalueerd worden in het fijne model, waarna door middel van vergelijking (6.4) (k) de waarde p xf bepaald kan worden. Vergelijking (6.4) wordt berekend met behulp van de optimalisatieprocedure die in hoofdstuk 5 ontwikkeld werd. Vervolgens moet een nieuwe (k+1)

xf

bepaald worden. Deze procedure gaat door tot p (xf ) naar xr∗ convergeert of een (k+1)

maximum aantal iteraties bereikt wordt. Het bepalen van de volgende iteratie xf

gebeurt

in ons geval met behulp van de ‘Aggressive Space Mapping’-methode. Dit is het onderwerp van de volgende paragraaf. Het is uiteraard de bedoeling dat beide modellen niet te sterk verschillen zodat formule (6.3) een voldoende nauwkeurige benadering is. In het bijzonder is het best dat de nietlineariteit van het fijne model ook voor een deel ge¨ımplementeerd is in het ruwe model, zodat de afbeelding niet te zeer niet-lineair moet zijn [23]. Een meer uitgebreide uiteenzetting over de principes van space mapping is te vinden in [24].

6.3

Aggressive Space Mapping (k+1)

In deze paragraaf bespreken we hoe de stap bepaald wordt naar de waarde xf

. De

‘Aggressive Space Mapping’-methode gebruikt hiervoor quasi-Newton-iteraties [25]. De quasiNewton-methode is een eenvoudige methode om niet-lineaire stelsels van vergelijkingen op te lossen. We houden hier nog geen rekening met de begrenzingen. Deze komen aan bod in paragraaf 6.4. De vergelijking die moet worden opgelost naar xf drukken we uit als:

e(xf ) = p(xf ) − xr∗ = 0

(6.6) (1)

e : Ωf → Ωr noemen we de foutfunctie. Als startwaarde stellen we xf

= xr∗ . Voor een k-de

iteratie nemen we een stap h (k) :

(k+1)

xf

(k)

= xf

+ h (k)

De stap h k wordt bepaald volgens de methode van Newton:

(6.7)


70

(k)

B (k) h (k) = −e(xf ) −1 (k) h (k) = − B (k) e(xf )

(6.8) (6.9) (k)

De matrix B (k) ∈ <m×m is de benadering van de jacobiaan Jp (xf ) van de afbeeldingsfunctie: (k)

(k) Jp (xf )

=

∂p(xf ) ∂xf

(6.10)

Met Jf de jacobiaan van het fijne model en Jr de jacobiaan van het ruwe model kunnen we ook zeggen:

Jf = Jr Jp

(6.11)

We berekenen deze matrix B (k) met de recursieve eersterangsformule van Broyden (e (k) = (k)

e(xf )):

B (k+1) = B (k) +

e (k+1) − e (k) − B (k) h (k) (k) T h T h (k) h (k)

(6.12)

Als startwaarde voor B nemen we de eenheidsmatrix. Meer informatie over de Broydenmethode is te vinden in [19] en [26]. Deze methode om de matrix B te bepalen vereist geen enkele bijkomende evaluatie van de modellen. Ze is opgesteld om te voldoen aan twee eigenschappen. We willen ten eerste bereiken dat:

(k+2)

e(xf (k+2)

We willen immers dat e(xf

(k+1)

) − e(xf

) ≈ B (k+1) h (k+1)

(6.13)

) → 0, wat samen met vergelijking (6.8) tot vergelijking (6.13)

leidt. Dit doel kan slechts benaderend bereikt worden aangezien het om een niet-lineair systeem gaat. Hiertoe bepalen we B (k+1) zo dat:

(k)

(k+1)

e(xf ) = e(xf

) − B (k+1) h (k)

(6.14)

We zorgen er met andere woorden voor dat B (k+1) , de benadering van de jacobiaan voor de volgende stap, in de richting van de vorige stap (een stap wordt immers gezet in een vijfdimensionale variabelenruimte) gelijk is aan de werkelijke helling van de vorige stap. We


71

zien dat vergelijking (6.12) hieraan voldoet. Dit principe is, voor één dimensie, weergegeven in figuur 6.1.

Figuur 6.1: De Broyden-update

De tweede eigenschap van deze methode is dat de benadering van de jacobiaan, B (k+1) , enkel bijgewerkt wordt in de richting van de vorige stap. Dus geldt voor elke willekeurige v ⊥h (k) :

B (k+1) v = B (k) v

(6.15)

Ook hieraan wordt voldaan door vergelijking (6.12). Het alternatief voor deze methode is om bij elke iteratie de jacobiaan numeriek te berekenen door eindige differentie naar de n verschillende elementen van xf . Dit vereist n evaluaties in het fijne model, bovenop de evaluatie die uitgevoerd werd om f (xf ) te berekenen, dat aan de basis lag van de berekening van p(xf ). Deze methode neemt een veel grotere rekentijd in beslag. Rekening houdend met vergelijking (6.8) is vergelijking (6.12) te vereenvoudigen tot:

B (k+1) = B (k) +

e (k+1) h (k)T h (k)T h (k)

(6.16)

Als ke (k) k voldoende klein wordt, stopt het algoritme. Dit gebeurt ook als de staplengte zeer klein wordt of als een maximaal aantal iteraties overschreden wordt. De waarden van de verschillende criteria kunnen opgegeven worden.


6.4

72

Trust Region

We gebruiken de uitdrukking ‘Trust Region’ voor het opleggen van begrenzingen aan de waarden van de variabele parameters. Deze waarden moeten immers fysisch verantwoord zijn. Het is niet mogelijk om, zoals vroeger toegepast werd, een onbegrensde optimalisatie uit te voeren waarbij de kostenfunctie rekening houdt met de overtredingen van de begrenzingen. Een waarde die niet met een fysische werkelijkheid overeenstemt zal immers niet uitvoerbaar zijn in de modellen die het directe probleem oplossen. Deze toepassing is niet identiek aan deze die in [27] besproken wordt, aangezien daar bij elke iteratie trust region toegepast wordt als dempingsmethode. (k)

De trust-regionprocedure wordt telkens toegepast als een nieuwe xf -waarde bepaald is. Deze waarde ligt niet meteen aan de basis van de volgende evaluatie van het fijne model, maar moet eerst aangepast worden om aan alle fysische begrenzingen te voldoen. Deze begrenzingen kunnen vrij worden opgegeven, maar er wordt een standaardset ingesteld bij de initialisatie. De hier besproken trust-regionmethode wordt ook gebruikt in de optimalisatieprocedure met één model met meerdere startwaarden, namelijk om de willekeurige startwaarden in overeenstemming te brengen met de begrenzingen. Zie paragraaf 5.2. Deze procedure, met meerdere startwaarden, wordt in het space mapping algoritme ook gebruikt om bij het begin van het algoritme het absolute optimum xr∗ te bepalen in het ruwe model evenals bij de parameterextractie. We noemen x (k) de waarde die moet worden getest. Het algoritme begint door te controleren of de waarde x (k ) voldoet aan de begrenzingen. Zoals eerder vermeld bestaan deze begrenzingen uit de volgende vergelijkingen:

Ax (k ) ≤ B

(6.17)

Aeq x (k ) = Beq

(6.18)

lo ≤ x (k ) ≤ lb

(6.19)

Als aan alle begrenzingen voldaan is, is de procedure ten einde. Als echter aan één of meer begrenzingen niet voldaan is wordt eerst weergegeven welke parameterwaarden welke begrenzingen overtreden. Deze output kan uiteraard worden uitgeschakeld. Vervolgens worden de parameters aangepast om binnen de grenzen te komen. Hiervoor zijn twee methodes ge¨ımplementeerd, waartussen door middel van een invoerparameter kan gekozen worden. We

HOOFDSTUK 6. SPACE MAPPING (k)

noemen xn

73

de nieuwe, aangepaste waarde voor de parameters die bepaald wordt in deze

procedure. De eerste methode wordt standaard gebruikt. Hierbij wordt de loodrechte projectie van x (k) op de rand van het domein bepaald. Dit gebeurt door een klassieke minimalisatie met het algoritme van paragraaf 5.3.1, uiteraard rekening houdend met de begrenzingen.

2 De kostenfunctie die geminimaliseerd wordt is y − x (k) , met y de tijdelijke variabele waarnaar geminimaliseerd wordt. Er geldt dus:

xn(k) = arg min ky − x k2

(6.20)

y

(k)

Na het uitvoeren van deze minimalisatie wordt gecontroleerd of de nieuwe waarde xn voldoet aan de begrenzingen. Als dit niet het geval is wordt de procedure herhaald. Deze situatie kan normaal echter niet voorkomen. De tweede methode verwacht als invoer niet alleen de huidige, te testen waarde x (k) , maar ook een waarde x (k−1) waarvan we weten dat ze voldoet aan de begrenzingen. Deze methode laat de variabele y een rechte tussen beide waarden overlopen, van x (k) naar x (k−1) , waarbij telkens een kleine stap gezet wordt en na elke stap gecontroleerd wordt of de nieuwe y -waarde aan de begrenzingen voldoet. Wanneer een y -waarde gevonden wordt die voldoet stopt het (k)

algoritme en is de nieuwe, aangepaste waarde xn

= y . Als we voor de gekende waarde

effectief x (k−1) gebruiken beperkt deze methode gewoon de lengte van de stap h (k−1) , terwijl de richting wel behouden blijft. Deze methode stemt overeen met de Levenberg-Marquadtmethode.

6.5

Hybride Methode

Het voorgestelde ‘Aggressive Space Mapping’-algoritme heeft te kampen met een aantal beperkingen [24]. Ten eerste kan convergentie niet gegarandeerd worden, in het bijzonder in punten waar een relatief groot verschil is tussen het ruwe en het fijne model en wanneer het algoritme de oplossing zeer dicht benadert. Ten tweede bestaat er ook geen garantie dat het bereikte optimum exact gelijk is aan het optimum van het fijne model. Daarom stellen we een hybride methode voor. Deze gebruikt het space-mappingalgoritme zolang dit convergeert. Indien wordt vastgesteld dat geen convergentie meer optreedt wordt echter een procedure aangeroepen die een directe optimalisatie in het fijne model gebruikt. Het algoritme is weergegeven in figuur 6.2.


74

Figuur 6.2: De hybride methode

De criteria die gecontroleerd worden om te bepalen of wordt overgestapt naar een directe optimalisatie zijn de volgende [28]. Als eerste criterium stellen we dat de kostenfunctie toe-

HOOFDSTUK 6. SPACE MAPPING (k)

(k−1)

neemt: K(f (xf )) > K(f (xf

75

)). Bovendien moet aan één van volgende twee voorwaarden (k)

(k)

voldaan zijn. De te minimaliseren functie e(xf ) = p(xf ) − xr∗ kan kleiner worden dan een opgegeven minimumwaarde. Hiermee wordt dus de nabijheid bij de oplossing gedetecteerd. Ofwel kan het voorkomen dat de werkelijke afname van de foutfunctie kleiner is dan 1% van de verwachte afname. Dit wordt gecontroleerd als volgt:

(k−1) (k) (k−1) (k−1) (k−1) (k−1) ) − e(xf ) < 0.01 e(xf ) − e(xf )+B h

e(xf

(6.21)

Indien aan het eerste criterium en aan één van de beide laatste criteria beantwoord wordt, roepen we de procedure aan die een directe optimalisatie uitvoert. Deze procedure is gebaseerd op de methode die in paragraaf 5.3.1 besproken werd. Ze begint met als startpunt de (k)

iteratie k te kiezen waar de kostenfunctie K(f (xf )) het kleinst is. Vanaf dit punt worden de (m)

jacobianen van de kostenfuncties JK

(m)

= ∇K(f (xf

)) numeriek berekend. Hierbij overloopt

de waarde m alle iteraties vóór iteratie k waarvoor de berekening van de jacobiaan nog niet is uitgevoerd in een eerdere aanroeping van de directe optimalisatie. Aan de hand van al deze waarden wordt de hessiaan van iteratie k berekend met de formule van Broyden, Fletcher, Goldfarb en Shanno, vergelijking 5.15. Hiertoe wordt een startwaarde H (0) , de eenheidsmatrix, of een waarde uit de vorige aanroeping van de directe optimalisatie, voor elke iteratie k −m ge¨ updatet om zo uiteindelijk H (k) te vinden. Aangezien de benodigde jacobianen reeds voorhanden zijn gebeurt dit zonder evaluaties van het fijne model. Het numeriek berekenen van de jacobianen gebeurt met behulp van eindige differenties. Hierbij wordt een zeer kleine, eindige stap εej gezet in de richting van elk van de vrije variabelen, waarbij ej gelijk is aan de eenheidsvector in de richting van de j-de variabele. In deze punten wordt een evaluatie uitgevoerd in het fijne model en wordt de kostenfunctie van deze evaluatie (m) berekend. De bijhorende waarde JK in de jacobiaan is dan gelijk aan: j

(m)

(m) JK

j

=

K(f (xf

(m)

)) − K(f (xf εej

+ εej ))

(6.22)

Voor variabelen die door middel van de begrenzingen vastgelegd zijn wordt deze berekening niet uitgevoerd, de bijhorende waarde in de jacobiaan stellen we gewoon gelijk aan 0. De jacobiaan is hier een vector, aangezien de kostenfunctie een scalaire functie is. De evaluatie (k)

in het punt xf

en de berekening van de kostenfunctie van deze evaluatie zijn reeds eerder

gebeurd, er zijn om één jacobiaan te berekenen dus even veel evaluaties in het fijne model nodig als er vrije variabelen zijn. De relatieve stapgrootte ε kan opgegeven worden, maar


76

wordt wel aangepast als de opgegeven waarde kleiner is dan de vierkantswortel van de machinenauwkeurigheid, van de orde 10−16 . Het product van ε met de elementen van de vector (k)

mag immers niet kleiner zijn dan deze nauwkeurigheid. Daarom wordt een element van

(k)

vervangen door ε als de absolute waarde ervan kleiner is dan ε.

xf xf

Vervolgens berekent de procedure de stap die gezet moet worden door de kwadratische benadering van de kostenfunctie te minimaliseren. Dit geeft als resultaat:

(k+1)

xf (k+1)

Op xf

(k)

= xf

−1 (k) ∇K(f (xf )) − H (k)

(6.23)

wordt trust region toegepast zodat de waarde met een fysische werkelijkheid over-

eenstemt, waarna het fijne model geëvalueerd wordt in de nieuwe waarde. Hierna wordt de kostenfunctie berekend. Vervolgens wordt, om aansluiting te behouden met het space-mappingalgoritme, een parame(k +1 )

terextractie toegepast om p(xf

(k +1 )

) te bepalen en wordt de foutfunctie e(xf

) berekend.

Uiteindelijk worden de voorwaarden om de space-mappingprocedure te stoppen overlopen en herbegint het algoritme, waarbij zoals eerder besproken beslist wordt of voor een spacemappingiteratie gekozen wordt of voor een iteratie met directe optimalisatie.

6.6

Conclusie

We hebben een algoritme opgesteld om een invers probleem op te lossen, gebruik makend van een ruw en een fijn model. We hebben hierbij gebruik gemaakt van de ‘Aggressive Space Mapping’-principes, uitgebreid tot een hybride methode om convergentie te garanderen en voorzien van trust-regionstrategieën om geen fysisch onmogelijke situaties te simuleren. De procedure is ontworpen in het kader van onze zoektocht naar een methode om een geschakelde reluctantiemotor te optimaliseren. Ze kan echter in elke context worden toegepast waar een ruw en een fijn model van eenzelfde systeem voorhanden zijn. Deze procedure brengt onze beide directe modellen samen, waardoor het mogelijk wordt om een snelle en nauwkeurige optimalisatie van de elektrische machine uit te voeren.

Hoofdstuk 7

Optimaal ontwerp van machines 7.1

Inleiding

We onderzoeken de invloed van de lokale materiaaldegradatie van het elektrisch staal op het optimale ontwerp van de machine. We analyseren wat de verschillen zijn in de optimale parameterwaarden indien we al dan niet werken met aangepaste randen. Hierbij bekijken we eerst het ontwerp enkel werkend met het ruwe model. Verder gaan we na wat het optimale ontwerp is indien we werken met het numerieke model. Voor deze optimalisatie maken we gebruik van space mapping.

7.2

Geldigheid en belang van de resultaten

In paragraaf 2.3 hebben we analytische benaderingen bepaald van de B-H-karakteristieken van het materiaal. Het zijn deze analytische benaderingen die in de voorwaartse modellen gebruikt worden, we werken dus niet met het biineaire materiaalmodel. De vergelijkingen die we opgesteld hebben zijn echter slechts geldig binnen een bepaald gebied voor de magnetische inductie, namelijk het gebied waar de karakteristieken direct gebaseerd zijn op werkelijke meetwaarden. Het verloop van de analytische benaderingen van de materiaalkarakteristieken die in het ruwe model worden gebruikt is weergegeven in figuur 7.1. We stellen vast dat het fysisch niet zinvol is deze karakteristieken te gebruiken buiten het geldige werkingsgebied, aangezien deze een snijpunt vertonen. Dit betekent dat, voor zeer hoge veldsterktes, het aangetaste materiaal van hogere kwaliteit zou blijken te zijn dan het originele materiaal. De werkelijkheid is uiteraard dat beide karakteristieken elkaar benaderen voor hoge veldsterktes. 77

HOOFDSTUK 7. OPTIMAAL ONTWERP VAN MACHINES

78

Figuur 7.1: Analytische benaderingen van gebruikte materiaalkarakteristieken in ruw model, blauw: onaangetast, rood: aangetast

Onderzoek in het ruwe model heeft aangetoond dat de optredende magnetische inducties binnen het geldige gebied blijven bij bekrachtigingsstromen tot 5 A. We zullen het ontwerp van de machine dan ook beperken tot deze stromen. Dit is geen grote inhoudelijke beperking: voor grotere stromen zal in kritische zones toch geen wezenlijk verschil zijn tussen de gevallen met of zonder lokale materiaaldegradatie, aangezien de eigenschappen van beide materialen elkaar benaderen bij verzadiging. Uit de resultaten zal ook blijken dat de invloed van de materiaaldegradatie op het ontwerp afneemt naarmate de bekrachtigingsstroom groter wordt. Om de invloed van de lokale materiaaldegradatie op het geleverde koppel aan te tonen, voeren we een koppelberekening uit bij een bepaald ontwerp van de motor, dus met dezelfde geometrische parameterwaarden voor elke bekrachtigingsstroom en al dan niet met inbegrip van de lokale materiaaldegradatie. In figuur 7.2 zien we voor de gehele bekrachtigingsstromen van 1 A tot 5 A het koppelprofiel met en zonder aangetaste randen. We stellen vast dat de invloed van de lokale materiaaldegradatie significant is. Bij kleine rotorhoeken loopt de afwijking op tot meer dan 20%. Het is dus zinvol om aparte optimalisaties uit te voeren voor al dan niet aangetaste randen.


79

Figuur 7.2: Koppelprofiel bij één geometrie, bij de bekrachtigingsstromen I = 1 A, 2 A, 3 A, 4 A, 5 A, blauw: zonder lokale materiaaldegradatie, rood: met lokale materiaaldegradatie

7.3 7.3.1

Ontwerp in ruw model Methode

De benodigde rekentijd voor een optimalisatie met behulp van het ruwe model, met vijf startpunten en één bekrachtigingsstroom, varieert op de gebruikte computer met stijgende bekrachtigingsstromen tussen 150 s en 500 s. Voor één optimalisatie, bij één startwaarde en naar vier parameters worden tussen 15 en 35 iteraties uitgevoerd, hierbij wordt het ruwe model tussen 100 en 320 keer geëvalueerd. Bij de resultaten die we in deze paragraaf beschrijven is geen enkele begrenzing actief. Elk bereikt optimum bevindt zich binnen de opgelegde grenzen. Dit is voordelig omdat een actieve begrenzing van één parameter invloed heeft op de optimale waarden van alle parameters, en zodoende een negatieve impact heeft op de consequentie van de resultaten.

7.3.2

Optimalisatie naar maximaal koppel

We optimaliseren naar maximaal koppel (zie paragraaf 5.2.2) voor elke gehele bekrachtigingsstroom van 1 A tot 5 A en zowel met inbegrip van de materiaaldegradatie als zonder materiaaldegradatie. De optimalisatie gebeurt naar vier parameters: de inwendige diameter van het statorjuk Ds , de uitwendige diameter van het rotorjuk Dr , de breedte van de rotortand brt en de breedte van de statortand bst . De diameter van de luchtspleet houden


80

we constant: Dls = 60 mm. De resultaten zijn weergegeven in figuur 7.3. Het verschil ∆x tussen de parameterwaarden met en zonder inbegrip van de lokale materiaaldegradatie is weergegeven in figuur 7.4. De maximale afwijking tussen het geval met en het geval zonder materiaaldegradatie is 8 %.

Figuur 7.3: Optimale parameterwaarden naar maximaal koppel, rood: met materiaaldegradatie, blauw: zonder materiaaldegradatie. Ds : inwendige diameter van statorjuk, Dr : uitwendige diameter van rotorjuk, brt : breedte van rotortand, bst : breedte van statortand

We stellen inderdaad vast dat de invloed van de materiaaldegradatie afneemt met toenemende bekrachtigingsstromen. Zoals fysisch te verwachten is, zien we dat de uitwendige diameter van het rotorjuk Dr , en dus de breedte van het rotorjuk, toeneemt wanneer de materiaaldegradatie in acht wordt genomen. De inwendige diameter van het statorjuk Ds neemt af bij inbegrip van de materiaaldegradatie, wat tot een verbreding van het statorjuk leidt. We kunnen de materiaaldegradatie immers beschouwen als een ‘versmalling’ van het magnetische pad, aangezien de permeabiliteit afneemt. Om hetzelfe effect te krijgen moet de werkelijke breedte dus toenemen. We stellen ook vast dat een toename van de bekrachtigingsstroom dezelfde invloed heeft op Ds en Dr als de materiaaldegradatie. Verzadiging van het materiaal heeft immers hetzelfde effect als lokale materiaaldegradatie: de reluctantie van een materiaaldeel neemt toe. De invloed op de breedtes van stator- en rotortanden, bst en brt , is vrij beperkt. De rotortand is breder dan de meeste andere materiaaldelen, de reluctantie ervan zal dus relatief klein zijn. De optimale breedte van de rotortand wordt dan ook eerder bepaald door de invloed ervan


81

Figuur 7.4: Verschil van optimale parameterwaarden naar maximaal koppel met materiaaldegradatie en zonder materiaaldegradatie. Ds : inwendige diameter van statorjuk, Dr : uitwendige diameter van rotorjuk, brt : breedte van rotortand, bst : breedte van statortand

op de luchtspleetreluctantie dan door de materiaaleigenschappen. De optimale statortandbreedte wordt eveneens deels vastgelegd door de luchtspleetreluctantie. Hier zijn verder twee tegengestelde effecten actief. Enerzijds zal de versmalling van het magnetische pad door de materiaaldegradatie aanleiding geven tot een toename van de tandbreedte. Het aantal bekrachtigingswindingen zou hierdoor echter worden beperkt, wat een negatieve invloed heeft op het gemiddelde koppel. Dit laatste effect is ook geldig voor Ds , wat meteen verklaart dat de variatie van Ds beperkter is dan deze van Dr . In figuur reffig:koppelmax zijn de koppelprofielen weergegeven bij de optimale parameterwaarden, zowel met (rood) als zonder (blauw) inbegrip van materiaaldegradatie en bij bekrachtigingsstroman van 1 A tot 5 A.

7.3.3

Optimalisatie met inbegrip van koppelrimpel

De koppelprofielen in figuur 7.5 zijn opvallend smal. Dit volgt uit het feit dat we naar het maximale (kwadratisch) gemiddelde koppel geoptimaliseerd hebben. Het spreekt voor zich dat een dergelijk koppelprofiel in werkelijke toepassingen over het algemeen niet wenselijk is. Daarom onderzoeken we het effect van een optimalisatie waarbij rekening gehouden wordt met de koppelrimpel. Dit gebeurt door middel van een gewichtsfactor w in de kostenfunc-


82

Figuur 7.5: Koppelprofielen bij optimale parameterwaarden naar maximaal koppel, bij de bekrachtigingsstromen I = 1 A, 2 A, 3 A, 4 A, 5 A, blauw: zonder lokale materiaaldegradatie, rood: met lokale materiaaldegradatie

tie, zie paragraaf 5.2. De gewichtsfactor geeft het belang van een groot gemiddeld koppel in de optimalisatie ten opzichte van het belang van een kleine koppelrimpel. De optimalisatie is uitgevoerd voor een gewichtsfactor die de waarden 1, 0.9, 0.7, 0.5 en 0.3 aanneemt. De opgelegde bekrachtigingsstroom is 4 A. De resulterende optimale parameterwaarden zijn weergegeven in figuur 7.6, de bijhorende koppelprofielen in figuur 7.7. We zien dat, wanneer een meer ‘normaal’ koppelprofiel gevraagd wordt, dus voor kleinere w, enkele parameterwaarden substantieel gewijzigd worden. De breedtes van stator- en rotorjuk nemen sterk toe (dit vertaalt zich in de wijzigingen van Ds en Dr ). Hierdoor neemt de verzadiging bij kleine rotorhoeken (en dus bij lage totale reluctantie) af, waardoor de flux in dit gebied stijgt. Dit zorgt dan weer voor een verbreding van het koppelprofiel aan de linkerkant, en dus een lagere koppelrimpel. Verder neemt de breedte van de statortand fel toe. Dit zorgt inderdaad voor een sterke verbreding van het koppelprofiel, niet alleen door de verlaging van de reluctantie in de tand, maar vooral omdat de variatie van de luchtspleetreluctantie over een groter gebied van de rotorhoek gespreid is. We selecteren w = 0.7 als gewichtsfactor voor onze verdere optimalisaties. We kunnen nu de invloed van de materiaaldegradatie en de bekrachtigingsstroom op de optimale parameterwaarden opnieuw onderzoeken. De optimale waarden van de geometrische parameters zijn weergegeven in figuur 7.8. Het verschil tussen de optimale parameterwaarden met en zonder rekening te houden met de aangetaste randen is weergegeven in figuur 7.9. In figuur 7.10 is het


Figuur 7.6:

83

Optimale parameterwaarden in functie van gewichtsfactor w voor koppelrimpel bij I = 4 A. Ds : inwendige diameter van statorjuk, Dr : uitwendige diameter van rotorjuk, brt : breedte van rotortand, bst : breedte van statortand

Figuur 7.7: Koppelprofielen bij I = 4 A en bij optimale parameterwaarden, in functie van gewichtsfactor w waarbij de optimalisatie uitgevoerd werd

koppelprofiel weergegeven bij optimale geometrie, dus met een verschillende geometrie voor elke bekrachtigingsstroom en met of zonder materiaaldegradatie. Het koppel is weergegeven voor 1 A, 2 A, 3 A, 4 A en 5 A.


Figuur 7.8:

84

Optimale parameterwaarden met w = 0.7, rood: met materiaaldegradatie, blauw: zonder materiaaldegradatie. Ds : inwendige diameter van statorjuk, Dr : uitwendige diameter van rotorjuk, brt : breedte van rotortand, bst : breedte van statortand

Figuur 7.9: Verschil van optimale parameterwaarden met w = 0.7 met materiaaldegradatie en zonder materiaaldegradatie. Ds : inwendige diameter van statorjuk, Dr : uitwendige diameter van rotorjuk, brt : breedte van rotortand, bst : breedte van statortand

Opvallend is hier het verloop van Ds . Deze neemt voor bepaalde stroomwaarden toe wanneer de materiaaldegradatie in acht genomen wordt, voor andere stroomwaarden is een negatieve invloed merkbaar. Dit kan verklaard worden als een gevolg van het samenspel tussen de


85

Figuur 7.10: Koppelprofielen bij optimale parameterwaarden met w = 0.7, bij de bekrachtigingsstromen I = 1 A, 2 A, 3 A, 4 A, 5 A, blauw: zonder lokale materiaaldegradatie, rood: met lokale materiaaldegradatie

twee tegenstrijdige effecten waarvan sprake: enerzijds wordt de toename van de reluctantie opgevangen door meer bekrachtigingswindingen toe te laten, en dus een smaller statorjuk te voorzien, anderzijds wordt de reluctantie beperkt door de breedte van het statorjuk te vergroten. Aangezien de materiaaleigenschappen sterk veranderen tussen de situaties met 1 A en 5 A bekrachtigingsstroom, is het aanvaardbaar dat de invloed van de materiaaldegradatie op de optimale parameterwaarden niet bij elke bekrachtigingsstroom dezelfde richting heeft. De andere parameterwaarden verlopen ongeveer gelijk als bij de optimalisatie naar maximaal koppel. In het bijzonder zien we opnieuw een sterke toename van Dr wanneer de materiaaldegradatie in acht genomen wordt. In figuur 7.11 zien we, naast de koppelprofielen zonder materiaaldegradatie en geoptimaliseerd zonder aangetaste randen (blauw) en deze met materiaaldegradatie en geoptimaliseerd met aangetaste randen (rood), ook de koppelprofielen die bekomen worden bij simulatie met materiaaldegradatie, maar met de parameters uit een optimalisatie die geen rekening houdt met deze materiaaldegradatie (paars). We observeren dat de sterk negatieve invloed van de materiaaldegradatie op het koppelprofiel voor een deel gecompenseerd wordt door deze materiaaldegradatie in acht te nemen bij de optimalisatie, rekening houdend met de afweging tussen koppelrimpel en maximaal koppel. Deze compensatie varieert wel met verschillende bekrachtigingsstromen. We hebben eveneens een optimalisatie uitgevoerd voor alle gehele bekrachtigingsstromen


86

Figuur 7.11: Koppelprofielen bij optimale parameterwaarden met w = 0.7, bij de bekrachtigingsstromen I = 1 A, 2 A, 3 A, 4 A, 5 A, blauw: zonder lokale materiaaldegradatie, geoptimaliseerd zonder materiaaldegradatie, rood: met materiaaldegradatie, geoptimaliseerd met materiaaldegradatie, paars: met lokale materiaaldegradatie, geoptimaliseerd zonder lokale materiaaldegradatie

van 1 A tot 4 A tegelijk. De resultaten zijn weergegeven in tabel 7.1. Deze resultaten zijn consequent met de resultaten die eerder bekomen werden; de grootte van de wijzigingen aan de parameterwaarden is echter vrij beperkt.

Tabel 7.1:

Resultaat van optimalisatie naar I=1 A, 2 A, 3 A en 4 A tegelijk bij w = 0.7. Ds : inwendige diameter van statorjuk, Dr : uitwendige diameter van rotorjuk, brt : breedte van rotortand, bst : breedte van statortand

bst (mm)

brt (mm)

Ds (mm)

Dr (mm)

Met mat. deg.

16.45

21.89

118.84

44.02

Zonder mat.

17.04

21.79

118.63

43.70

∆x

-0.59

0.10

0.21

0.32


7.4 7.4.1

87

Ontwerp met behulp van space mapping Methode

We hebben de optimalisaties uitgevoerd met behulp van de space-mappingprocedure die we ontworpen hebben zoals besproken in hoofdstuk 6. Een optimalisatie met behulp van space mapping duurt, voor toenemende bekrachtigingsstroom, tussen 4 en 13 uren. Hierbij wordt het fijne model ongeveer 10 keer geëvalueerd. Van het ruwe model worden enkele duizenden evaluaties uitgevoerd. Het is duidelijk dat dit een sterke verbetering is ten opzichte van de directe optimalisatie in het fijne model, waarbij minstens een week nodig is om tot een resultaat te komen. We hebben bij deze berekeningen geen gebruik gemaakt van de hybride methode. Deze zou tot meer nauwkeurige resultaten leiden, maar de benodigde rekentijd zou eveneens met ongeveer een factor 5 toenemen. We beschouwen bij wijze van voorbeeld van het space-mappingalgoritme het verloop van (k)

de kostenfunctie K(f (xf )) bij de optimalisatie voor een bekrachtigingsstroom van 3 A, zie figuur 7.12.

Figuur 7.12: Verloop van de kostenfunctie bij space-mappingoptimalisatie

(k)

(k) ∗

De norm van de foutfunctie e = p(xf ) − xr is weergegeven in figuur 7.13. Beide worden verwacht te dalen, de kostenfunctie moet naar een minimum streven en de foutfunctie naar nul. We zien dat dit ook het geval is. We bemerken één uitschietende


88

Figuur 7.13: Verloop van de foutfunctie bij space-mappingoptimalisatie

waarde. Deze is echter te wijten aan het feit dat de gewenste stap in de betreffende iteratie niet kon worden gezet omwille van een fysische begrenzing.

7.4.2

Optimalisatie naar maximaal koppel

We voeren optimalisaties uit naar maximaal koppel (zie paragraaf 5.2.2) voor de bekrachtigingsstromen 1 A, 2 A, 3 A en 4 A en zowel met als zonder inbegrip van de lokale materiaaldegradatie aan de randen. De optimalisatie gebeurt naar dezelfde vier parameters als deze die we gebruikten bij de optimalisatie in het ruwe model: de inwendige diameter van het statorjuk Ds , de uitwendige diameter van het rotorjuk Dr , de breedte van de rotortand brt en de breedte van de statortand bst . De diameter van de luchtspleet houden we constant: Dls = 60 mm. De bekomen optimale parameters zijn weergegeven in figuur 7.14. Het verschil tussen de parameters bekomen door optimalisatie met inbegrip van materiaaldegradatie en deze bekomen door optimalisatie zonder inbegrip van materiaaldegradatie is weergegeven in figuur 7.15. We zien dat de inwendige statordiameter Ds kleiner wordt, dus het statorjuk breder wordt als de materiaaldegradatie in rekening genomen wordt. Hetzelfde gebeurt met de breedte van de rotortanden. Deze wijzigingen stemmen overeen met het effect dat in paragraaf 7.3.2 besproken werd. In rekening brengen van lokale materiaaldegradatie aan de rand komt erop neer dat de reluctantie van het magnetische pad toeneemt. Dit wordt gecorrigeerd door de breedte van het betreffende onderdeel te laten toenemen.


89

Figuur 7.14: Optimale parameterwaarden, rood: met materiaaldegradatie, blauw: zonder materiaaldegradatie. Ds : inwendige diameter van statorjuk, Dr : uitwendige diameter van rotorjuk, brt : breedte van rotortand, bst : breedte van statortand

Figuur 7.15: Verschil van optimale parameterwaarden met materiaaldegradatie en zonder materiaaldegradatie. Ds : inwendige diameter van statorjuk, Dr : uitwendige diameter van rotorjuk, brt : breedte van rotortand, bst : breedte van statortand

Het effect op de breedte van de rotortand brt is weliswaar vrij beperkt omdat deze relatief breed is in vergelijking met de andere onderdelen van het magnetisch circuit. Dit volgt uit


90

de invloed die deze breedte heeft op de flux in de luchtspleet. De inwendige statordiameter ondervindt dan weer invloed van twee tegenstrijdige effecten. Enerzijds moet de stijging van de reluctantie, te wijten aan de materiaaldegradatie, worden gecompenseerd door een verbreding van het juk. Anderzijds zou dit ten koste gaan van het aantal windingen van de bekrachtigingswikkeling. We zien dat, door de kleine breedte van het statorjuk, het eerste effect de bovenhand haalt. De optimale breedte van de statortand bst ondervindt weinig invloed van de materiaaldegradatie. Aangezien de breedte van de statortand groter is dan de breedte van het statorjuk, houden beide effecten elkaar hier ongeveer in evenwicht. De uitwendige diameter van het rotorjuk Dr neemt voor kleine bekrachtigingsstromen sterk toe als de materiaaldegradatie in rekening genomen wordt. Dit beantwoordt aan de fysische verwachtingen, volgens het principe dat een afname van de reluctantie gecompenseerd wordt door een verbreding van het onderdeel. We zien echter dat de invloed van de materiaaldegradatie afneemt met toenemende bekrachtigingsstroom. Dit is deels te verklaren door het feit dat de materiaalkarakteristieken van het aangetaste materiaal en deze van het onaangetaste materiaal elkaar benaderen in het verzadigde gebied. Voor I = 4 A wordt Dr echter kleiner bij inachtname van de materiaaldegradatie. Een mogelijke verklaring hiervoor is te vinden in de invloed van de flux φtj die rechtstreeks tussen de statortand en het rotorjuk stroomt. Deze invloed manifesteert zich bij grote rotorhoeken. De flux φtj zorgt ervoor dat de totale flux minder afneemt bij toenemende rotorhoeken dan het geval zou zijn als φtj niet aanwezig zou zijn. Hierdoor verlaagt het koppel voor grotere rotorhoeken. Bij grote bekrachtigingsstromen is het magnetisch materiaal meer verzadigd, en speelt deze flux een grotere rol. Ook als de lokale materiaaldegradatie in rekening genomen wordt neemt de invloed van φtj toe, aangezien het magnetisch materiaal dan een lagere permeabiliteit vertoont. Dit effect kan gecompenseerd worden door een afname van de uitwendige diameter van het statorjuk, aangezien de reluctantie tussen statortand en rotorjuk dan groter wordt en φtj dus afneemt. We merken tevens op dat de resulterende variabelen, zoals weergegeven in figuur 7.15 een vrij grote overeenkomst vertonen met de variabelen die we bekomen bij directe optimalisatie in het ruwe model, zie figuur 7.4. In figuur 7.16 zijn drie sets koppelprofielen weergegeven voor de bekrachtigingsstromen 1 A, 2 A, 3 A en 4 A. Ten eerste zijn er de koppelprofielen die bekomen worden door evaluatie van het fijne model, zonder lokale materiaaldegradatie en met de geometrische parameters die bekomen worden uit een optimalisatie waarin eveneens geen materiaaldegradatie in rekening genomen wordt (blauw). Dan zijn er de koppelprofielen die bekomen worden door een


91

evaluatie met aangetaste randen, maar met de geometrische parameters uit een optimalisatie waarbij deze materiaaldegradatie niet in rekening genomen wordt (paars). Ten slotte zijn er de koppelprofielen die bekomen worden uit het fijne model, met inbegrip van de materiaaldegradatie en met de geometrische parameters waarbij ook de optimalisatie uitgevoerd werd rekening houdend met de lokale materiaaldegradatie (rood).

Figuur 7.16: Koppelprofielen bij optimale parameterwaarden bij de bekrachtigingsstromen I = 1 A, 2 A, 3 A, 4 A, blauw: zonder lokale materiaaldegradatie, geoptimaliseerd zonder materiaaldegradatie, rood: met materiaaldegradatie, geoptimaliseerd met materiaaldegradatie, paars: met materiaaldegradatie, geoptimaliseerd zonder materiaaldegradatie

Om de positieve invloed van de optimalisatie met inbegrip van de lokale materiaaldegradatie te illustreren, zijn de waarden van de kostenfunctie voor de koppelprofielen uit figuur 7.16 weergegeven in tabel 7.2. We noemen K (a) de kost bij een aangetaste rand, geoptimaliseerd met aangetaste karakteristieken. K (na) is de kost zonder materiaaldegradatie, waarbij de optimalisatie eveneens zonder materiaaldegradatie gebeurde. K (no) is de kost met materiaaldegradatie maar waarbij de optimalisatie zonder materiaaldegradatie uitgevoerd werd. Er wordt een identieke term van elke kost afgetrokken en een factor 1000 toegepast om overzichtelijke getalwaarden te bekomen. We geven voor elke bekrachtigingsstroom ook de fractie α, gedefinieerd als:

α=

K (no) − K (a) K (no) − K (na)

(7.1)

Dit is de fractie van de toename van de kostenfunctie door de materiaaldegradatie die gecompenseerd wordt door bij de optimalisatie eveneens de materiaaldegradatie in acht te nemen.


92

Tabel 7.2: Kostenfuncties van koppelprofielen met en zonder aangetaste randen en al dan niet met inbegrip van aangetaste randen in de optimalisatie

1A

2A

3A

4A

K (na)

9.2503

7.4055

4.9813

2.4884

K (no)

9.3438

7.6071

5.2826

2.8050

K (a)

9.3292

7.6057

5.2597

2.7797

α (%)

15.63

0.69

7.61

8.01

We merken op dat α een zeer kleine waarde aanneemt voor I = 2 A. Dit is te wijten aan een onvolledige convergentie, wat een gevolg is van het feit dat deze berekeningen zijn uitgevoerd zonder gebruik te maken van het hybride model.

7.5

Conclusie

We hebben de invloed onderzocht van de lokale materiaaldegradatie op de optimale geometrische parameters van de machine. Deze invloed blijkt vooral sterk uitgesproken te zijn voor de afmetingen van het stator- en rotorjuk. Bovendien kunnen we besluiten dat de procedures die we opgesteld hebben in staat zijn deze minima te bepalen.

Hoofdstuk 8

Conclusie We hebben een methode ontworpen om de geometrische parameters van een elektrische machine te optimaliseren. Hierbij hebben we de ‘space mapping’-techniek toegepast. Op deze manier zijn we erin geslaagd nauwkeurige optimalisaties uit te voeren in een beperkte rekentijd. We hebben twee modellen opgebouwd die door de space-mappingprocedure gebruikt worden om het directe probleem op te lossen. Het eerste is een ruw model dat de motor benadert door een magnetisch netwerk. Dit leidt tot een snelle evaluatie maar een beperkte nauwkeurigheid. Het tweede model maakt gebruik van eindige-elementenberekeningen. In onze modellen is rekening gehouden met de lokale materiaaldegradatie die ontstaat ten gevolge van het productieproces van de machine. De precieze magnetische eigenschappen die gebruikt worden in de modellen hebben we experimenteel bepaald. We hebben ook de invloed van de lokale materiaaldegradatie op het ontwerp van de elektrische machines gekwantificeerd. Uit onze berekeningen blijkt dat in het bijzonder de optimale afmetingen van stator- en rotorjuk een relatief grote invloed ondergaan van de materiaaldegradatie. Door middel van deze methode is het mogelijk om bij het ontwerp van de machine meteen op een verantwoorde manier rekening te houden met de lokale materiaaldegradatie. Deze werkwijze kan de bijkomende thermische bewerkingen of de overdimensioneringen die heden gebruikt worden overbodig maken en dus tot een efficiënter productieproces leiden.

93

Bibliografie [1] A.J. Moses, N. Derebasi, G. Loisos en A. Schoppa. Aspects of the cut-edge effect stress on the power loss and flux density distribution in electrical steel sheets. Journal of Magnetism and Magnetic Materials, vol. 215-216 (2000), pp. 690–692. [2] V. Maurel, F. Ossart en R. Billardon. Residual stresses in punched laminations: phenomenological analysis and influence on the magnetic behavior of electrical steels. Journal of Applied Physics, vol. 93 (2005), pp. 7106–7108. [3] V. Permiakov. 1D and 2D magnetization in electrical steels under uniaxial stress. Doctoraatsthesis, Universiteit Gent (2005). [4] G. Crevecoeur, L Dupré en L. Vandenbossche. Reconstruction of local magnetic properties of steel sheets by needle probe methods using space mapping techniques. Journal of Applied Physics, (2006). To be published. [5] M. De Wulf, E. Hoferlin en L. Dupré. Finite element modelling of induction machines taking into account manufacturing processes. ICEM (2004). [6] L. Vandenbossche, L. Dupré en J. Melkebeek. Preisach based magnetic evaluation of fatigue damage progression. Journal of Magnetism and Magnetic Materials, vol. 290291 (2005), pp. 486–489. [7] IEC-norm 404-3. Methods of measurement of the magnetic properties of magnetic sheet & strip by means of a single sheet tester. International Electrotechnical Commission (IEC) (1992). [8] M. De Wulf, L. Dupré, D. Makaveev en J. Melkebeek. Needle-probe techniques for local magnetic flux measurements. Journal of Applied Physics, vol. 93 (2003), pp. 8271–8273. [9] P. Sergeant. Laagfrequente magnetische afscherming van elektrische installaties. Doctoraatsthesis, Universiteit Gent (2006). 94

BIBLIOGRAFIE

95

[10] L. Dupré. Elektromagnetische energieomzetting. Syllabus voor studenten eerste proef burgerlijk werktuigkundig-elektrotechnisch ingenieur (2003). Universiteit Gent, Faculteit Ingenieurswetenschappen. [11] W. Cheney en D. Kincaid. Numerical mathematics and computing. Science and Math, vijfde editie (2004). [12] G. Strang en G.J. Fix. An analysis of the finite element method. Prentice-Hall (1973). [13] Y. Stry, M. Hainke en Th. Jung. Comparison of linear and quadratic shape functions for a hybrid control-volume finite element method. International Journal of Numerical Methods for Heat & Fluid Flow, vol. 12 (2002), pp. 1009–1031. [14] E. Haber, U.M. Ascher en D. Oldenburg. On optimization techniques for solving nonlinear inverse problems. Inverse problems, vol. 16 (2000), pp. 1263–1280. [15] M.C. Biggs. Constrained minimization using recursive quadratic programming. In Towards global optimization, pp. 341–349. North Holland (1970). [16] K. Madsen, H.B. Nielsen en O. Tingleff. Optimization with constraints. Technical University of Denmark, Informatics and mathematical modelling, tweede editie (2004). [17] K. Madsen, H.B. Nielsen en O. Tingleff. Methods for non-linear least squares problems. Technical University of Denmark, Informatics and mathematical modelling, tweede editie (2004). [18] P.E. Gill, W Murray en M.H. Wright. Numerical linear algebra and optimization, vol. 1. Addison Wesley (1991). [19] W.H. Press, B.P. Flannery, S.A. Teukolsky en W.T. Vetterling. Numerical recipes in C: the art of scientific computing. Cambridge University Press (1993). [20] R.P. Brent. Algorithms for minimization without derivatives. Prentice-Hall (1973). [21] G.K. Smyth. Optimization. In Encyclopedia of environmetrics, vol. 3, pp. 1481–1487. John Wiley & Sons (2002). [22] J.W. Bandler, R.M. Biernacki, S.H. Chen, P.A. Grobelny en R.H. Hemmers. Space mapping techniques for electromagnetic optimization. IEEE Transactions on Microwave Theory and Techniques, vol. 42 (1994), pp. 2536–2544.

BIBLIOGRAFIE

96

[23] G. Crevecoeur, H. Hallez, P. Van Hese, Y. D’Asseler, L. Dupré en R. Van De Walle. EEG source analysis using space mapping techniques. Journal of Computational and Applied Mathematics, (2006). To be published. [24] J. Søndergaard. Optimization using surrogate models - by the space mapping technique. Doctoraatsthesis, Technical University of Denmark, Informatics and mathematical modelling (2003). [25] E. Haber. Quasi-Newton methods for large-scale electromagnetic inverse problems. Inverse problems, vol. 21 (2005), pp. 305–323. [26] C.G. Broyden. A class of methods for solving nonlinear simultaneous equations. Mathematics of Computation, vol. 19 (1965), pp. 577–593. [27] J.W. Bandler, Q.S. Cheng, S.A. Dakroury, A.S. Mohamed, M.H. Bakr, K Madsen en J Søndergaard. Space Mapping: the state of the art. IEEE Transactions on Microwave Theory and Techniques, vol. 52 (2004), pp. 337–361. [28] M.H. Bakr, J.W. Bandler, M.B. Radoslaw, S. Chen en K. Madsen. A trust region aggressive space mapping algorithm for EM optimization. IEEE Transactions on Microwave Theory and Techniques, vol. 46 (1998), pp. 2412–2425.

Universiteit Gent Faculteit Ingenieurswetenschappen. Promotor: Prof. dr. ir. L. DUPRE Begeleider: ir. G. CREVECOEUR

Recommend Documents