Benoît, Grégoire, Manon Judith Simon

Damien Laurent

Doorian

Adeline Morgan

Arnould Nicola Pascal

Marie In

Fatigue induced by deterministic and random loading

Keywords: Isogeometric analysis, Finite element method, Random vibration, Numerical and experimental tests

Nowadays, mechanical systems are composed of numerous components, leading to specific requirements for numerical analysis methods in terms of accuracy and time computation. The finite element analysis (FEA) is essential in simulation technology, as it can help engineers predict results to avoid experimental cost testing. However, the conventional FEA has some mesh generation and geometric representation weaknesses. The time spent in the mesh generation process is often much longer than the analysis time. It is mandatory to communicate with the original geometry during each mesh refinement process, which takes time and can lead to complex problems. Furthermore, it is challenging to accurately represent complex geometric models based on the Lagrangian basis function in classical FEA. Thus, the approximate models would result in inaccurate analysis results.

Failures in mechanical structures are mostly due to fatigue. Fatigue life prediction is essential to know the actual life before a part is used in operation. These operating parts experience random fatigue loads and the most appropriate approach is to use a probabilistic method to develop the fatigue analysis. This analysis can be performed in the time or frequency domain for Gaussian random processes. However, the time domain analysis is time consuming compared to using power spectral density (PSD) in the frequency domain.

Thus, this study considered the random vibration fatigue analysis in the frequency domain using the isogeometric analysis (IGA) method. Geometric models can be described by non-uniform rational B-spline curves (NURBS), and the corresponding mesh generation and refinement processes are extremely fast compared to FEA.

The main content of the thesis is divided into four parts. In the first part, the IGA was developed on a simple plate model to be compared with FEA based on the analytical solution. In the second part, fatigue analyses by Dirlik method were performed on a plate model. FEA was used to check the results and experimental tests were used to give a relative conclusion on the methods. An L-shaped plate and a wind turbine tower model were also studied in random fatigue in the third part. The results obtained in terms of damage were validated by FEA and a Matlab program that we developed.

In the fourth part, isogeometric optimization was developed on the tower model by combining Matlab programming and LS DYNA.

In conclusion, the IGA is suitable for numerical simulations of mechanical structures since the results obtained highlight:

• IGA is efficient in terms of mesh generation compared to FEA.

• IGA can provide similar results to FEA with an acceptable relative error given that for damage calculation, IGA requires fewer integration points and mesh elements.

• Optimization with design variables as control points can provide consistent models.

Résumé

De nos jours, les systèmes mécaniques sont composés de nombreux éléments, ce qui entraîne des exigences spécifiques pour les méthodes d'analyse numérique en termes de précision et de temps de calcul. La méthode d'analyse par éléments finis (MEF) est essentielle dans la technologie de simulation, car elle peut aider les ingénieurs à prévoir les résultats afin d'éviter les essais expérimentaux coûteux. Cependant, la MEF classique présente certaines faiblesses en matière de génération de maillage et de représentation géométrique. Le temps consacré au processus de génération de maillage est souvent beaucoup plus long que le temps d'analyse. Il est nécessaire de communiquer avec la géométrie d'origine pendant chaque processus de raffinement du maillage, ce qui prend du temps et peut rendre le problème complexe. En outre, il est difficile de représenter avec précision des modèles géométriques complexes basés sur la fonction de base Lagrangienne dans la MEF classique. Ainsi, des modèles approximatifs donnent des résultats d'analyse inexacts.

Les défaillances des structures mécaniques sont principalement dues à la fatigue. La prédiction de la durée de vie en fatigue est essentielle pour connaître la durée de vie effective avant qu'une pièce ne soit utilisée en fonctionnement.

Ces pièces subissent des charges de fatigue aléatoires, et l'approche la plus appropriée est d'utiliser une méthode probabiliste pour développer l'analyse de fatigue. L'analyse de la fatigue peut être effectuée dans le domaine temporel ou fréquentiel pour les processus aléatoires Gaussiens. Cependant, l'analyse dans le domaine temporel prend beaucoup de temps par rapport à l'utilisation de la densité spectrale de puissance (DSP) dans le domaine fréquentiel.

Ainsi, cette étude est basée sur l'analyse de la fatigue par vibration aléatoire dans le domaine fréquentiel en utilisant la méthode d'analyse isogéométrique (IGA). Les modèles géométriques peuvent être décrits par des courbes B-spline rationnelles non uniformes (NURBS) et les processus de génération et de raffinement du maillage correspondant sont extrêmement rapides par rapport à la MEF.

Le contenu principal de la thèse est divisé en quatre parties. Dans la première, l'IGA a été développée sur un modèle de plaque simple afin d'être comparée aux résultats de la MEF en se basant sur la solution analytique. Dans la deuxième partie, des analyses de fatigue par la méthode de Dirlik ont été réalisées sur un modèle de plaque. La MEF a été utilisée pour vérifier les résultats obtenus par fatigue et des tests expérimentaux ont été utilisés pour donner une relative conclusion sur les méthodes. Dans la troisième partie, une plaque en forme de L et un modèle de tour d'éolienne ont également été étudiés en fatigue aléatoire dans la troisième partie. Les résultats obtenus en termes de dommages ont été validés par la MEF et un programme Matlab que nous avons développé. Dans la quatrième partie, une optimisation isogéométrique a été développées sur le modèle de tour en combinant la programmation Matlab et LS DYNA.

Pour conclure, l'IGA est adaptée aux simulations numériques de structures mécaniques étant donné que les résultats obtenus mettent en évidence :

• l'IGA est efficace en termes de génération de maillage par rapport à la MEF.

• l'IGA peut fournir des résultats similaires à la MEF avec une erreur relative acceptable étant donné que pour le calcul des dommages, l'IGA nécessite moins de points d'intégration et d'éléments de maillage. ii • L'optimisation avec des variables de conception comme points de contrôle peut fournir des modèles cohérents.

Mots clés

Analyse isogéométrique ; Méthode des éléments finis ; Vibration aléatoire ; Essais numériques et expérimentaux iii

General Introduction

Based on industrial data, 80 to 95% of mechanical structures fail to function due to fatigue [START_REF] Ringeval | Random vibration fatigue analysis with ls dyna[END_REF], in which most loadings are random. It's essential to predict the structural fatigue strength before its application in a working environment.

Under most circumstances, the fatigue analysis relies on the FEA, which has some weaknesses. Firstly, the time spent in the mesh generation process is often much longer than analysis time [START_REF] Owen | An immersive topology environment for meshing[END_REF]. It is assumed that 80% of overall analysis time is spent on the mesh generation in automotive, aerospace, and shipbuilding industries [START_REF] Hughes | Isogeometric analysis: Cad, finite elements, nurbs, exact geometry and mesh refinement[END_REF]. Furthermore, it is necessary to communicate with the original geometry model during each mesh refinement process. Secondly, the approximate finite element mesh model often results in inaccurate analysis results.

On the other hand, the geometry and analysis models have different representations in the FEA field. The shape optimization with design variables as element nodes would result in zigzag and unrealistic shapes. Moreover, an extra post-processing step is required to exchange data with the CAD system to obtain the final design results, making the structural optimization more complicated.

At present, mechanical structures often consist of numerous components. It is vital to develop structural fatigue analysis and design optimization with a compact, time-efficient analysis process, precise and ideal analysis results to improve the effectiveness and efficiency. Pitoiset et cl. [START_REF] Pitoiset | Tools for a multiaxial fatigue analysis of structures submitted to random vibrations[END_REF] proposed two kinds of frequency domain fatigue analysis methods, the equivalent von Mises stress and multiaxial rainflow methods, to evaluate the high cycle fatigue damage of a supported rectangular plate subjected to a band-limited white noise random pressure. The results show that the two approaches can provide similar PSD functions and damage ratios, allowing a designer to localize the most critical elements. After, Pitoiset et cl. [START_REF] Pitoiset | Spectral methods for multiaxial random fatigue analysis of metallic structures[END_REF] also applied frequency domain fatigue analysis methods on an L-shaped plate model subjected to multiaxial random loading to estimate the high-cycle fatigue life and compared the analysis results with the critical plane time-domain method. It has been found that compared with the time-domain approach, the frequency domain methods are extremely time-saving (17 seconds for the equivalent von Mises stress, 1.6 minutes for the multiaxial rainflow, and one week for the critical plane) and have a good agreement in localizing critical areas in the structure. E.Dowling [START_REF] Dowling | Mean stress effects in stress-life and strain-life fatigue[END_REF] investigated the different fatigue life correction methods, including Goodman, Morrow, and Smith Watson Topper (SWT) and Walker as their ability to correlate stress and strain life on several sheets of steel and nonferrous metals, and obtained the following conclusions from stress-life correlation: 1) the modified Goodman equation with ultimate tensile strength was found to be highly inaccurate; 2 ) the Morrow equation using the true fracture strength worked well on certain materials; viii 3) the Morrow expression with fatigue strength coefficient led to non-conservative values; 4) the SWT method was a rational choice that can avoid above shortcomings; 5) Walker method with an adjustable exponent can fit very well with the test date, allowing higher accuracy. Seong-in Moon et cl. [START_REF] Moon | Fatigue life evaluation of mechanical components using vibration fatigue analysis technique[END_REF] proposed a methodology to decide the optimum vibration fatigue test, which gives an equivalent failure mode with driving test condition, through a series of vibration fatigue analyses by changing acceleration directions and magnitudes. Arshad et cl. [START_REF] Sarkar | Random vibration analysis and fatigue life evaluation of auxiliary heater bracket[END_REF] developed frequency domain random vibration fatigue analysis on an auxiliary heater bracket in Abaqus commercial software. By modifying the geometry model, an infinite fatigue life in 1σ and 2σ level of confidence was achieved. It has been concluded that Abaqus is efficient in modeling random vibration fatigue. Adam et cl. [START_REF] Niesłony | Mean stress value in spectral method for the determination of fatigue life[END_REF] firstly proposed to take into account the effects of mean stress on the fatigue life calculation using the spectral method, in which the PSD has been determined based on the transformed stress considering the mean value, and the detailed steps have been presented.Essam et cl. [START_REF] Essam | Nvh and random vibration fatigue analysis of a landing gear ' s leg for an unmanned aerial vehicle using ls-dyna ®[END_REF] developed random vibration fatigue analysis on a landing gear based on spectral methods with the use of different ways, including Dirlik, Steinberg, Hancock, Wirsching, Tunna, Chaudhury, and Dover in Ls Dyna commercial software.

Through comparison of analysis results, it was observed that there were many discrepancies among obtained fatigue life based on these methods. Yusuf et cl. [START_REF] Eldoǧan | Vibration fatigue analysis of a cantilever beam using different fatigue theories[END_REF] developed finite element random vibration fatigue analysis on a cantilever beam using different damage models and compared fatigue life of the beam obtained from the frequency domain, time domain, and experimental test. On the other hand, they also studied the effects of the damping ratio on damage results. It has been observed that the damping ratio is vital to determine correct fatigue damage, and compared with other damage models, the Dirlik method can provide more rational results. Giovanni de Morai Teixeira [START_REF] Teixeira | Random vibration fatigue analysis of a notched aluminum beam[END_REF] developed random vibration fatigue analysis on a notched plate in the software Fe-safe. The analysis results were validated by the time-domain method based on an equivalent approach, predicting that the fatigue life from frequency and time domain fatigue analysis at critical locations differs by 20%, in the acceptable range. It was concluded that the random fatigue analysis in fe-safe is robust and can provide accurate results. Yun Huang et cl. [START_REF] Huang | Recent updates in fatigue analysis with ls-dyna[END_REF] developed frequency and time domain fatigue analysis on different cases in commercial software Ls Dyna, in which in the frequency domain, random and steady-state vibration fatigue analysis were introduced, and in the time domain, stress and strain-based fatigue analysis were applied. In the frequency domain fatigue analysis, the mean stress correction methods, such as Goodman, Soderberg, Gerber, were developed on a bracket model. From analysis results, it can be shown that the adoption of mean stress correction can lead to obtaining shorter fatigue life. Demirel [START_REF] Demirel | Implementation of dirlik's damage model for the vibration fatigue analysis[END_REF] developed random vibration fatigue analysis on a rectangular cross-section notched beam based on Dirlik's model, studied the effects of modal damping ratio on fatigue analysis results, and validated the fatigue analysis results with an experimental test. However, some discrepancies in fatigue analysis and experimental test results have been found, which have been explained by the different ultimate tensile strength values, damping ratio, and the difficulties of observing the crack initiation time in a real experimental test. Böhm et cl. [START_REF] Böhm | Fatigue life estimation with mean stress effect compensation for lightweight structures-the case of glare 2 composite[END_REF] developed frequency and time domain random vibration fatigue analysis with the application of mean stress compensation in hybrid materials consisting of three layers of aluminum and two layers of Glass Fiber Reinforced Polymer, in which the fatigue life was calculated under three different fiber orientations. It has been observed that in a high-cycle fatigue regime, the fatigue life calculated from Rainflow and spectral methods show good agreement; in a low-cycle fatigue regime, compared with the Goodman model, the fatigue life based on the ix Gerber model can provide good compatibility between Rainflow and spectral calculations.

Based on the papers above, it can be concluded that the fatigue analysis mainly focused on the numerical simulation and experimental test. In terms of numerical analysis: some mean stress compensation methods, including Goodman and Gerber, were proposed; different fatigue analysis methods like Dirlik, Steinberg, Tuna were applied to analyze fatigue damage in the frequency domain, and the results were validated by the Rainflow counting method in the time domain or experimental tests. It can also be found that almost all of the fatigue analysis simulations were based on the FEA, in which the approximate mesh models may lead to inaccurate results.

Therefore, this work addresses to develop a random vibration fatigue analysis and design optimization on a complex mechanical structure based on a new numerical analysis method, IGA. In the beginning, the IGA was carried on plate models to validate the static, convergence, and fatigue analysis results by an analytical solution, experimental test, FEA, and Matlab programming. Then, the isogeometric random vibration fatigue analysis and design optimization were developed on a wind turbine model and validated by FEA and Matlab programming. This thesis consists of 6 chapters, organized as follows. Chapter 1 focuses on the theories related to the fatigue damage calculation induced by deterministic and random loading, followed by the chapter conclusion. The isogeometric modeling is described in Chapter 2. It starts from the description of fundamental theories, concepts, and research state of IGA. At the same time, the IGA is considered on a simple plate model, based on the validation by an analytical solution.

Chapter 3 presents the IGA on a specimen and an L-shaped plate model. The isogeometric static, convergence, and random vibration fatigue analysis are respectively taken into account in which the fatigue analysis results are validated by FEA, Matlab programming, and experimental tests. The isogeometric random vibration fatigue analysis, size, and shape optimization are carried on a wind turbine tower model in Chapter 4 with the design variables of tower segment thicknesses and control point position, respectively. Chapter 5 summarizes this work. IGA formulations, wind loads calculation, the IGA procedures in Ls Dyna, and this work's general analysis procedures are given in the appendices.

x

B i,p (ξ), B j,q (η), B k,z (ζ): ith, jth, kth Bernstein basis functions of order p, q, z along ξ, η, ζ parametric direc- tions N i,p (ξ), M j,q (η), R k,z (ζ): ith, jth, kth B-spline basis functions of order p, q, z along ξ, η, ζ parametric direc- tions R i,p (ξ), R p,q i,j (ξ, η), R p,q,z i,j,k (ξ, η, ζ): Univariate,

Introduction

In the mid-nineteenth century, with the invention of steam locomotives, it was discovered that wheel structures failed under the limit stresses far less than their static strength, promoting the occurrence of general fatigue problems. Then, it was observed that under the effects of structural resonance frequency, a vehicle would also experience failure, even though the load was much smaller than the general fatigue load, and this contributed to the occurrence of vibration fatigue problems (or dynamic fatigue problems) [START_REF] Yao | Vibration fatigue in engineering structures[END_REF]. Fatigue is defined as the progressive and localized structural damage caused by repeated loads [START_REF] Moon | Fatigue life evaluation of mechanical components using vibration fatigue analysis technique[END_REF][START_REF] Eldoǧan | Vibration fatigue analysis of a cantilever beam using different fatigue theories[END_REF]. In the beginning, only static fatigue analysis that doesn't consider the inertial effects has attracted much attention. In 1963 S.H.Crandall and W.D.Mark firstly proposed that vibration fatigue is irreversible and is a cumulative damage process caused by dynamic loads [START_REF] Crandall | 1 -characterization of random vibration[END_REF]. However, it's not enough for a vibration fatigue analysis. In recent decades, Yao Qihang, Yao Jun, Zhang Azhou et cl. [START_REF] Yao | Vibration fatigue in engineering structures[END_REF][START_REF] Azhou | Vibration Environment Engineering[END_REF][START_REF] Yao | The characteristics and analysis methods of structural vibration fatigue problems[END_REF][START_REF] Yao | Simplified analysis of structural random vibration response[END_REF][START_REF] Yao | Structural vibration fatigue problems and airplane vibration assessment[END_REF] proposed that "vibration fatigue" refers to the fatigue damage caused by applied dynamic loads including vibration, noise, harshness approximate the natural frequency of a system.

The fatigue failures are featured by: 1) under the repetitive action of alternating load, fatigue may also occur, even though the obtained component stress is much smaller than the material's strength limit (sometimes elastic limit);

2) regardless of whether a material is brittle or plastic, macroscopically, the sudden fatigue fracture is featured by no obvious plastic deformation and low-stress brittle fracture, which make fatigue damage difficult to predict, and dangerous; 3) there are two main zones close to the fatigue crack areas, the smooth and the rough zone, also termed 1 as fatigue crack growth zone and fast fracture zone, and these are important criteria in determining fatigue failure; 4)

fatigue failure often has a local feature and does not involve all materials of the entire structure. Most of the time, as long as the local detailed design or process measures are changed, the fatigue life can be increased significantly.

Therefore, the ability of a structure or component to resist fatigue damage not only depends on the materials used but also sensitively depends on the structural shape, size, components connection and matching form, performance shape, and environmental conditions, etc; 5) fatigue failure is a cumulative damage process, which usually involves crack formation, cracks propagation, and rapid fracture when the crack propagates to a critical size, requiring a certain time history, even a long time course.

The fatigue performance of materials or components under alternating loads can be measured by fatigue strength.

The intensity of fatigue strength represented by fatigue limit or endurance limit refers to the maximum stress at which a material or component can withstand unlimited cycles under stress cycle characteristic R. As the fatigue limit of a material varies with the loading method and stress ratio, the fatigue limit under symmetrical stress cycles is used as the basic fatigue limit of the corresponding material. Usually, fatigue life is given as the cycle numbers of stress or strain experienced during a failure, generally denoted by N . The fatigue life of a test piece depends on the mechanical properties of the material and the level of applied stress. Generally speaking, the higher the strength limit of the material and the lower the applied stress level, the longer the fatigue life of a specimen; otherwise, the shorter the fatigue life.

Based on the highness of cyclic stress or the condition of occurrence in plastic strain, the fatigue can be classified into stress fatigue or strain fatigue. When the obtained maximum cyclic stress is much less than the corresponding material yield one, and no plastic strain occurs, the fatigue is called stress fatigue. In this situation, the number of stress cycles is relatively high, generally higher than 10 4 times. Therefore, the stress fatigue is also called high cycle fatigue (HCF). On the other hand, if the maximum cyclic stress is higher than the material yield limit, and there is an occurrence of plastic strain, it is more appropriate to use strain as the fatigue life estimation parameter. And due to the lower number in stress cycles, generally less than 10 4 times, the strain fatigue is also called low cycle fatigue.

This chapter provides a general introduction to the theoretical backgrounds of fatigue damage induced by deterministic and random loading. As shown in Fig. 1.1, it begins with the introduction in Section 1.1 in which the historical backgrounds in fatigue analysis, and the concepts in fatigue strength, fatigue limit, fatigue life, stress and strain fatigue are respectively described. Section 1.2 presents the deterministic fatigue analysis theories in S-N curves, Miner's rule, Rainflow counting, etc. Subsequently, Section 1.3 proposes the random vibration fatigue analysis theories such as the definition, classification of random processes, the power spectral density (PSD) and probability density function (PDF), and the fatigue analysis formulation in the frequency domain, etc. The chapter conclusions are given in Section 1.

Deterministic fatigue analysis

Fatigue is defined as a progressive change in the material properties following the application of loading cycles, the repetition of which can lead to the fracture. Schijve [START_REF] Schijve | Four lectures on fatigue crack growth[END_REF] defines fatigue life as the addition of the initiation period and the crack growth period. The author found that more than 90% of the fatigue life is usually spent before detectable cracks under HCF conditions. A cumulative damage calculation is used to quantify the total damage characterized by the fraction of life consumed by the crack initiation.

This study focused on HCF thus the material behaviour remains elastic during the fatigue loading until crack initiation. Therefore, an elastic computation was performed to obtain stress states used for the damage calculations.

For deterministic excitation, the linear model of Palmgren-Miner damage accumulation [START_REF] Miner | Cumulative damage in fatigue[END_REF], [START_REF] Palmgren | Die lebensdauer von kugellagern (life length of roller bearings or durability of ball bearings)[END_REF] knowing as Miner's rule, is one of the most widely used methods [START_REF] Krasnowski | Application of damage tolerance to increase safety of helicopters in service[END_REF][START_REF] Strauss | Life-cycle and sustainability of civil infrastructure systems[END_REF][START_REF] Risitano | Cumulative damage by miner's rule and by energetic analysis[END_REF]. Let us consider m stress amplitudes denoted σ a,i (i = 1, . . . , m), which characterize the time history response calculated at the maximum stress location of a structure (by FEA or IGA).

The authors assume that each significant stress amplitude extracted from the stress-time history (by using Rainflow counting techniques, see [START_REF] Lee | Cycle counting techniques[END_REF], [START_REF] Milne | Cyclic loading and fatigue[END_REF]) produces an individual damage and the total damage D is obtained by the following summation:

D = m i=1 1 N f,i N i (1.1)
where N f,i is the number of cycles allowable at particular stress before a material fails by fatigue (the number of cycles to failure at constant stress amplitude σ a ) defined from the S-N curve (representing the stress amplitudes versus the number of cycles to failure). N i the ith applied stress cycle number. For steel structures, if the cumulative fatigue damage reaches a critical value of 1, thus the structure is considered damaged.

The relationship between the fatigue limit σ RN (a particular value of stress amplitude, expressed as σ a , for the simplification in the following) and cycle number N is called the fatigue curve or S-N curve in which the stress cycle characteristic R is defined as the proportion between the minimum and maximum stresses. As shown in Fig. 1.2, when the cycle number N < 10 3 -10 4 , the fatigue limit almost approximates the ultimate strength σ u and generally can be calculated based on the standard of static strength. When the cycle number 10 3 -10 4 < N < N b , the fatigue limit would decrease with the increase of cycle number N . For most materials, the relationship between the fatigue limit σ a and cycle number N f follows the exponent Eq. (1.2). When the fatigue cycle N > N b , or the fatigue limit σ < σ b , it is traditionally assumed that the fatigue would not happen.

σ β a N f = C (1.2)
which can also be written as:

σ a = σ f N b f (1.3) with σ f = C 1 β and b = -1 β . lgσ a lgN f 10 3 (10 4 ) σ b A B N b σ u

Low cycle

High cycle finite life region infinite life region endurance limit

0 σ β a N f = C Figure 1.2: Schematic illustration of a material S-N curve
Therefore, with the consideration of Eq. (1.2), Eq. (1.1) becomes:

D = m i=1 C -1 σ β a,i N i (1.4)
For a deterministic loading, the fatigue analysis is performed in the time domain while for Gaussian random excitation, it can be developed in the time domain or in the frequency domain.

Cumulative damage for variable amplitudes of stresses

The cumulative damage is easy to calculate when a load signal behaves in a sine function. However, when the load amplitudes are variable, the rainflow cycle counting method, firstly proposed by Matsuishi and Endo [START_REF] Matsuichi | Fatigue of metals subjected to varying stress[END_REF], is used to count the number of fatigue cycles from a given stress or strain time history. There are several rainflow counting methods, such as hysteresis filtering, Peak-valley Filtering, four-point counting method, etc. The Hysteresis Filtering removes tiny oscillations from the load-time history that contribute negligible damage through the definition of a gate of a specific amplitude. Any cycle with amplitude smaller than the gate would be removed from the load-time history.

The purpose of peak-valley filtering is to save data points that are reversals in direction. Any intermediate data points between the maximum and minimum values of a given cycle can be deleted. The four-point counting method roughly requires the following steps: firstly, define inner and outer stress; secondly, based on the inner and outer stress range, the stress cycle is counted; thirdly, store the obtained cycles in a rainflow matrix, which is an n × n matrix of data, and n equals the number of bins. Then the combination of rainflow counting, material S-N curve, and Miner's rule allows calculating the fatigue damage of a component subjected to complex loads.

Factors affecting fatigue strength

Under the effects of alternating loads, the fatigue strength of a mechanical component is affected by several factors, such as stress concentration, component dimension, surface state, environmental medium, load sequence, load characteristics, etc., in which the first three ones are most significant. Thus, the component fatigue curve is different from the corresponding material counterpart.

• I) Component stress concentration coefficient K σ

With the demands of structures, a component usually has abrupt changes in holes, fillets, grooves, etc., leading to a phenomenon of stress concentration, under a load application. The component stress concentration coefficient K S is applied to represent the reduced certainties of the material fatigue limit and defined by the Eq. (1.7).

K σ = σ -1 (σ -1 ) k (1.5)
where σ -1 and (σ -1 ) k are respectively the material and component fatigue limit under a symmetric cycle.

• II) Component dimension coefficient S

Under the same conditions, the larger the component dimension, the higher the chance of fatigue defects. This is because when a component size is large, the material particles are coarse, the probability of defects is high, and the fatigue cracks are easy to form. Component dimension coefficient S is defined as:

σ = (σ -1 ) d (σ -1 ) d0 (1.6)
where (σ -1 ) d and (σ -1 ) d0 are respectively the fatigue limit with the specimen diameters of d and d 0 (normally

d 0 = 6 -10 mm).
• III) Component surface state coefficient β σ

The component surface state represent the component surface roughness and surface processing state. It can produce significant effects on fatigue failure in the aspects of crack initiation, stress concentration, and anti-fatigue ability. The Component surface state coefficient β σ is defined as:

β σ = (σ -1 ) βσ (σ -1 ) β0 (1.7)
where (σ -1 ) βσ is the fatigue limit under certain surface state, (σ -1 ) β0 is the fatigue limit under fine polishing.

• IV) The allowable stress Under considerations of these factors mentioned above and safety coefficient s d , the allowable stress of the symmetric and pulsating cycle can be defined by the Eqs. (1.8) and (1.9), respectively.

[σ -1 ] = σ β σ σ -1 k σ s d (1.8) [σ 0 ] = σ β σ σ 0 k σ s d (1.9) • V) Component σ a -σ m
It has been found that the component stress concentration, dimensions, and surface state coefficients have effects mainly only on stress amplitude σ a . Thus, through the combination of these three factors, the comprehensive impact coefficient (K σ ) D is decided by Eq. (1.10).

(K σ ) D = k σ σ β σ (1.10)
The component fatigue limit diagram is obtained after the consideration of comprehensive impact coefficient (K σ ) D and fatigue life coefficient K N . The comprehensive impact coefficient (K σ ) D has effects on the component stress amplitude. The fatigue life coefficient K N has effects both on the component stress amplitude and mean stress.

Summary

As shown in Fig. 1.3, a deterministic fatigue analysis usually would require Miner's rule, stress cycle counting algorithm, and component S-N curve to calculate the total fatigue damage. In the first step, obtain the stress time history of a structure from numerical simulation methods in which for a multiaxial stress state, the equivalent von Mises stress can be calculated, based on the Eq. (1.41); Then, according to the stress cycle counting method like the rainflow counting approach, compute the stress cycle numbers under different mean and alternating stresses in which for a material S-N curve obtained under the condition of zero mean stress, the component S-N curve with non-zero mean stress can be computed considering the comprehensive impact factor, as shown in Eq. (1.10), and mean stress correlation algorithm like Goodman [START_REF] Smith | The effect of range of stress on the fatigue strength of metals[END_REF], as shown in Eq. (1.11). In the final step, with the utilization of obtained component S-N curves under different mean stresses, Miner's rule, and stress cycle numbers from the algorithm of stress cycle, calculate the fatigue damage.

σ a σ ar = 1 - σ m σ u (1.11)
where σ ar is the obtained stress amplitude, σ u is the ultimate strength.

Random vibration fatigue analysis

For a deterministic fatigue loading, the fatigue analysis is ideally suited to the time domain. However, the fatigue analysis in the time domain would often require extensive time records to obtain an accurate analysis for a random loading, so the analysis process is always time-consuming. In reality, most real fatigue loads are unexpected. And for the case of a random loading signal, there are several advantages for a fatigue analysis in the frequency domain rather than the time domain. The main advantage is that the fatigue damage calculation is time-efficient [START_REF] Al-Bahkali | Fatigue life estimate of landing gear's leg using modal analysis[END_REF].

Gaussian and stationary random process

When random loading is considered, the damage and fatigue life are random variables. The random process σ(t)

corresponds to a random variation in time of a stress component or effective stress. At each instant t i , i = 1, 2, . . . , n, the process is composed of a number of random variables σ(t i ) that follow a Gaussian distribution. If the time t is continuous, the process is termed as a continuous-time random process. Similarly, if t is discrete, the process is called a discrete-time random process. Over a duration of observation T , the probability that a random process σ(t) does not exceed, in absolute value, a given value denoted by r can be given by the cumulative distribution function (CDF) of the amplitudes (or maxima for zero-mean stress) denoted by F σ (r):

F σ (r) = r -∞ p σ (r)dr = prob(|σ(t)| ≤ r), t ∈ [0, T ] (1.12)
where p σ (r) is the probability density function (PDF).

The mth order statistical moments of the continuous random variable σ(t i ) is obtained from the expected value operator:

E[σ(t) m ] = +∞ -∞
r m p σ (r)dr (1.13) m = 1 corresponds to the expected value and m = 2 gives the second moment or variance for the zero mean.

The random variables are completely characterized by their statistical properties and more precisely by the CDF of each variable such as:

F σ(t1),••• ,σ(tn) (r 1 , • • • , r n ) = prob(σ(t 1 ) ≤ r 1 , σ(t 2 ) ≤ r 2 , • • • , σ(t n ) ≤ r n ) (1.14)
where

0 < t 1 < t 2 < • • • < t n .
The autocorrelation function measures the similarity between two different observation times t 1 and t 2 at the random process σ(t), defined as:

R σσ (t 1 , t 2 ) = E[σ(t 1 )σ(t 2 )] (1.15)
A stationary random process can be categorized into a strict or wide sense counterpart based on its statistical properties. If its CDF, given in Eq. (1.14) is invariant to a time shift τ , the random process can be defined as a strict sense stationary (SSS) random process:

F σ(t1),••• ,σ(tn) (r 1 , r 2 , • • • , r n ) = F σ(t1-τ ),••• ,σ(tn-τ ) (r 1 , r 2 , • • • , r n ) (1.16)
The random process is a wide sense stationary (WSS) random process if its moments and autocorrelation functions do not depend on absolute time, such as :

E[σ(t 1 ) m ] = E[σ(t 2 ) m ] = • • • = E[σ(t n ) m ] = E[σ(t) m ] (1.17)
and,

R σσ (t, t -τ ) = R σσ (τ ) (1.18)
Note, the autocorrelation function reaches his maximum for τ = 0. In this case R σσ (0

) = E[σ(t) 2 ].
In addition, a random process is ergodic if the time average of the sample function and the statistical average are equals:

E[σ(t)] = lim T →∞ 1 2T +T -T σ(t)dt (1.19)
and the time autocorrelation is identical to the statistical autocorrelation:

R σσ (τ ) = lim T →∞ 1 2T +T -T σ(t)σ(t -τ )dt (1.20)

Fatigue analysis induced by random acceleration

The study of the dynamic response exploiting vibration phenomena aimed to determine the dynamic properties which are directly connected to the geometrical and mechanical characteristics of a system. Hence, some concepts of structural dynamic response and how to obtain the dynamic parameters will be summarized in this section. The excitation with base motion (e.g., car suspension, earthquake ground motion or electrodynamic shaker, etc.) does not provide the force applied information. Therefore the formalism related to the base motion is different from the classical force-excitation theory.

For random vibration loading, the global matrices M, C and K (and mode shapes) given in Eq.(A.6) remain deterministic. In time domain, the force and displacement field are written in terms of expected values denoted, respectively, E[F(t)] and E[u(t)]. In spectral domain, these two quantities are connected by the following relation:

S uu (f ) = H uF (f ) 2 S FF (f ) (1.21)
S uu (f ) and S FF (f ) are the PSDs of control point displacement and external forces. |H uF (f ) corresponds to the frequency response magnitude of the structure obtained from the ratio between the displacement and the force given in frequency domain. In time domain, the stress expected value of the structure can be, therefore, deduced from the force expected value by using the following expression:

E[σ(t)] = C B K -1 E[F(t)] (1.22) K -1 = H uF (0)
represents the static part. The transformation of Eq.(1.22) into the spectral domain leads to a relationship between the stress and displacement PSDs:

S σσ (f ) = C B S uu (f ) B T C T (1.23)

Cumulative fatigue and frequency formulation for random excitation

To work in the frequency domain, it is necessary to make certain assumptions on a random excitation to ensure its stability. Firstly, a random signal should be stationary, meaning the statistical properties, including expected value, variance do not change with time. Secondly, it should be ergodic to ensure that an arbitrary sample of the excitation can be applied to substitute for the whole sample. Thus, the Fourier transform applied on the autocorrelation function R σσ (τ ) to obtain the power spectral density (PSD) defined by:

S σσ (f ) = +∞ -∞ R σσ (τ )exp(-j2πf τ )dτ (1.24)
When a Gaussian stationary random process is considered, these parameters are mainly the kth order spectral moments given by:

m k = +∞ 0 |2πf | k S σσ (f )df (1.25)
S σσ (f ) is the stress PSD function. When the mechanical structure is submitted to a zero-mean Gaussian stationary random excitation, the expected damage per unit time can be obtained from the cumulative expression given in Eq.

(1.1): it corresponds to βth order statistical moment.

E[D] = C -1 E |σ(t)| β E[N ] (1.26
Knowing that for period T , E[N ] = T × E[N + 0 ], so Eq. (1.26) becomes:

E [D] = C -1 T E[N + 0 ] +∞ 0 r β p σ (r)dr (1.28) E[N + 0 ] = m2 m0
is the expected number of zero-crossings with positive slope per second. m 0 and m 2 are respectively the zeroth-order and second-order spectral moments of the PSD.

√ m 0 is the root mean square (RMS) and represents the square root of the area under the PSD-frequency graph. The expected fatigue life denoted E[T f ] can be then deduced from Eq. (1.28):

E[T f ] = T E[D] = 1 C -1 E[N + 0 ] +∞ 0 r β p σ (r)dr (1.29)
Before the damage evaluation, it is necessary to define a probability density function p σ (r) of the maxima. The different spectral methods for fatigue damage assessment were first characterized by the cycle counting procedure. Since the Rainflow method has been recognized to give the best predictions, the spectral methods were mainly focused on the Rainflow cycle distribution. There is no analytical formula to establish the cycle distribution, therefore the existing spectral methods evaluate these distributions approximately or empirically. Some of the most commonly used methods include the narrowband approach based on Rayleigh approximation, Dirlik's amplitude distribution [START_REF] Dirlik | Application of computers in fatigue analysis[END_REF]. In this work, Dirlik's distribution is used to evaluate the expected value of damage from Eqs. (1.28) and (1.29). ). The irregularity factor varies from 0 to 1. As it approaches 1, the random stress signal is close to NB. It comes 0, the stress signal approaches WB.

γ = E[N + 0 ] E[P ] = m 2 2 m 0 • m 4 (1.30)
In 1964, Bendat [START_REF] Bendat | Probability functions for random responses: Prediction of peaks, fatigue damage, and catastrophic failures[END_REF] firstly provided the frequency domain fatigue damage calculation through the definition of the PDF function, as shown in Eq. (1.31). But the damage values from Bendat's equation give conservative results for Wide-Band application. Therefore, Bandat's solution is also termed a narrow-band solution.

p σB (r) = r 4m 0 exp -r 2 8m 0 (1.31) E[D] NB = E[N + 0 ]T C +∞ 0 r β r 4m 0 exp -r 2 8m 0 dr (1.32)
To solve this problem, the conservative results for wide-band application, Wirsching [START_REF] Wirsching | Fatigue under wide band random stresses[END_REF] firstly proposed the expression for wide-band solution based on Bendat's, as shown in Eq. (1.33).

E[D] W = E[D] NB (a + (1 -a)(1 -ε) c ) (1.33)
where E[D] NB is the expected damage determined by the narrow band solution and a = 0.926 -0.033b, c =

1.587b -2.323, ε = 1 -γ 2 .
In the field of narrow-band fatigue analysis, Tuna [START_REF] Tunna | Fatigue life prediction for gaussian random loads at the design stage[END_REF] and Dirlik [START_REF] Dirlik | Application of computers in fatigue analysis[END_REF] respectively proposed the new PDF. Tuna's PDF is defined by:

p σT (r) = r 4γm 0 exp -r 2 8γm 0 (1.34)
Dirlik's PDF for a normalized variable Z = σa √ m0 is:

p σD (Z) = D1 Q exp -Z Q + D2Z R 2 exp -Z 2 2R 2 + D 3 Q S exp -Z 2 2 2 √ m 0 (1.35)
with m 4 the fourth-order spectral moment of the PSD. x m is defined by Dirlik as the mean frequency and is expressed such as:

x m = m 1 m 0 m 2 m 4 (1.36)
and the other parameters are obtained from:

D 1 = 2(x m -γ 2 ) 1 + γ 2 D 2 = 1 -γ -D 1 + D 2 1 1 -R D 3 = 1 -D 1 + D 2 Q S = 1.25(γ -D 3 -D 2 ) D 1 R = γ -x m -D 2 1 1 -γ -D 1 + D 2 1
The expected damage ratio can also be expressed in the form of equivalent stress, given by Eq. (1.37). Hancock [START_REF] Kam | Fast fatigue assessment procedure for offshore structures under random stress history[END_REF] firstly gave the equivalent stress expression by the Eq. (1.38). Later Chaudhuri and Dover (C&D) [START_REF] Chaudhury | Fatigue analysis of offshore platforms subject to sea wave loadings[END_REF] and Steinberg [START_REF] Steinberg | Vibration analysis for electronic equipment[END_REF] respectively gave the expression also, as shown in Eqs. (1.39) and (1.40). The time-series generation method from a given spectrum can be classified into deterministic or random amplitude scheme, in which the former is commonly applied [START_REF] Merigaud | Free-surface time-series generation for wave energy applications[END_REF]. Generally, the deterministic amplitude scheme has several procedures for generating random signals from the given PSD function. That are respectively 1) the selection of the frequencies f i and random phase angles ϕ i from the PSD function; 2) evaluation of the deterministic amplitudes:

E[D] = E[N + 0 ]T C σ eq (1.37) (σ eq ) Hancock = (2 √ 2m 0 ) γΓ β 2 + 1 1 β (1.38) (σ eq ) C and D = (2 √ 2m 0 ) β+2 2 √ pi Γ β 2 + 1 + γ 2 Γ β 2 + 2 + erf (γ) γ 2 Γ β 2 + 2 1 β (1.39) (σ eq ) Steinberg = 0.683(2 √ m 0 ) β + 0.271(4 √ m 0 ) β + 0.043(6 √ m 0 ) β 1 β (1.
A i = √ 2G i × ∆f i
, where ∆f i is the frequency bandwidth; 3) summation for each time t: Y (t) = n i A i sin(2πf i + ϕ i ).

Equivalent stress for random loading

In triaxial stress state, the von Mises stress σ c is given as:

σ c = 1 2 [(σ 1 -σ 2 ) 2 + (σ 2 -σ 3 ) 2 + (σ 3 -σ 1 ) 2 + 6(σ 2 12 + σ 2 23 + σ 2 31 )] (1.41)
where σ 1 , σ 2 , σ 3 and σ 12 , σ 23 , σ 31 are respectively normal and shear stresses. Defining the stress vector as σ = (σ 1 , σ 2 , σ 3 , σ 12 , σ 23 , σ 31 ) T , Eq. (1.41) can be written as:

σ 2 c = σ T Qσ = T raceQ[σσ T ] (1.42)
where:

Q =               1 -1 2 -1 2 0 0 0 -1 2 1 -1 2 0 0 0 -1 2 -1 2 1 0 0 0 0 0 0 3 0 0 0 0 0 0 3 0 0 0 0 0 0 3               (1.43)
In recent decades, multiaxial fatigue damage models were defined in the stress, strain, or energy time histories [START_REF] You | A critical review on multiaxial fatigue assessments of metals[END_REF].

However, these methods are often time-consuming for a finely discretized structure, as they often require complete stress-strain time histories. 1996, Pitoiset et cl. [START_REF] Pitoiset | Tools for a multiaxial fatigue analysis of structures submitted to random vibrations[END_REF] proposed a frequency domain fatigue analysis method for multiaxial stress states based on the definition of equivalent von Mises stress (an equivalent uniaxial stress PSD).

Take the expectation on Eq. (1.42):

E[σ 2 c ] = Trace{QE[σσ T ]} (1.44)
where E[σσ T ] is the covariance matrix of the stress vector σ, which is related to the PSD Φ σσ (ω) of the stress vector:

E[σσ T ] = +∞ -∞ Φ σσ (ω)dω (1.45)
where Φ σσ characterizes a zero mean Gaussian random stress field.

Thus, combine the Eqs. (1.44) and (1.45), obtain:

E[σ 2 c ] = +∞ -∞ Trace{QΦ σσ (ω)}dω (1.46) E[σ 2 c ] also equals to: E[σ 2 c ] = +∞ -∞ Φ c (ω)dω (1.47)
Combining the Eqs. (1.46) and (1.47), obtain the PSD Φ c (ω) of equivalent von Mises stress as:

Φ c (ω) = Trace{QΦ σσ (ω)} (1.48)

Conclusion

In this chapter, the fatigue analysis backgrounds were described in the introduction. Then deterministic fatigue analysis theories, such as the explanations in Miner's rule, S-N curves were presented in Section 1.2. The dynamic fatigue analysis theories, including the definition of the random process and formulations in the frequency domain damage calculation, were given in Section 1.3, and based on these theories, the random vibration fatigue analysis will be developed on different models in chapter 3. In the chapter 2, I will present the theoretical backgrounds of the IGA approach and perform IGA on a simple plate model.

Chapter 2

Isogeometric Analysis Approach

Introduction

Isogeometric Analysis (IGA) proposed by T.J.R Huges et cl. [START_REF] Hughes | Isogeometric analysis: Cad, finite elements, nurbs, exact geometry and mesh refinement[END_REF] is a relatively new numerical analysis method, which integrates Computer-Aided Design (CAD) and Computer-Aided Engineering (CAE) to save the time spent in mesh generation, and also to obtain more precise analysis results, through exact geometric representation and direct utilization of CAD-based NURBS described geometric information.

Numerical modelling of a complex structure is a time-consuming aspect. Nowadays, mechanical systems are an assembly of many components, leading to specific requirements on numerical analysis methods in terms of accuracy and speed of analysis. Conventional FEA has some weaknesses. The time spent in the mesh generation process is often much longer than analysis time [START_REF] Owen | An immersive topology environment for meshing[END_REF]. It is estimated that 80% of analysis time is devoted to the mesh generation in some fields, such as automotive or shipbuilding industries [START_REF] Hughes | Isogeometric analysis: Cad, finite elements, nurbs, exact geometry and mesh refinement[END_REF]. On the other hand, it is often necessary to communicate with the original geometry during each mesh refinement, and this process is time-consuming. The second disadvantage lies in geometric approximation. Indeed, it is challenging to accurately represent complex geometric models based on the Lagrangian basis function in classical FEA. Thus, the approximate models would result in inaccurate analysis results. A possible alternative to finite elements is IGA, often based on non-uniform rational B-splines (NURBS) basis functions.

This chapter mainly introduces the IGA approach regarding theoretical backgrounds and applications on a simple plate model. As shown in Fig. 2

C(ξ, η) (ξ ∈ [0, 1], η ∈ [0, 1]) and solid C(ξ, η, ζ) (ξ ∈ [0, 1], η ∈ [0, 1], ζ ∈ [0, 1]
) are respectively defined by the bivariate and trivariate Bernstein polynomials and control points.

1) Bernstein basis function

A pth degree Bernstein polynomial B i,p (ξ) is defined as:

B i,p (ξ) = C i p ξ i (1 -ξ) p-i ξ ∈ [0, 1] (2.1)
where p is the polynomial order in ξ parametric direction, C i p is the ith binomial coefficient and expressed by:

C i p = p! i!(p -i)! (2.2)
Considering polynomial order p = 4, the formulation process of the Bernstein basis function would be presented following. Beginning with i = 0, obtain:

B 0,4 (ξ) = C 0 4 ξ 0 (1 -ξ) 4 = (1 -ξ) 4 (2.3) i = 1, obtain B 1,4 (ξ) = C 1 4 ξ 1 (1 -ξ) 3 = 4ξ(1 -ξ) 3 (2.4) i = 2, obtain B 2,4 (ξ) = C 2 4 ξ 2 (1 -ξ) 2 = 6ξ 2 (1 -ξ) 2 (2.5) i = 3, obtain B 3,4 (ξ) = C 3 4 ξ 3 (1 -ξ) 1 = 4ξ 3 (1 -ξ) (2.6) i = 4, obtain B 4,4 (ξ) = C 4 4 ξ 4 (1 -ξ) 0 = ξ 4 (2.7)
Based on obtained Eqs. from 2.3 to 2.7, the corresponding basis functions can be constructed in Matlab, as shown in Fig. 2.2.

Figure 2.2: Bernstein basis functions

The important properties of Bernstein basis functions can be concluded as:

(1) partition of unity ξ ∈ [0, 1],

p+1 i=1 B i,p (ξ) = 1 (2) non negative ξ ∈ [0, 1], B i,p (ξ) ≥ 0,
(3) kronecker's delta B i,p (ξ j ) = δ ij which is satisfied at element boundary.

(4) linear independence

p+1 i=1 α i B i,p = 0, if α 1 = α 2 = • • • = α n = 0
(5) satisfy variation diminishing property, which leads to Bézier curves behaving monotonous with the increase of polynomial orders of basis functions.

(6) C 0 inter element continuity.

2) Bézier curve

The general form of a p degree Bézier curve C(ξ), which requires p + 1 number of control point is expressed as:

C(ξ) = p+1 i=1 B i,p (ξ)P i , ξ ∈ [0, 1],
(2.8)

where P i = (x i , y i , z i ) T (i = 1, . . . , p + 1) is the ith control point vector. As shown in Fig. 2.3, the polynomial order of 4 Bézier curve is defined by 5 control points, which begins at control point P 1 and ends at P 5 . B i,p (ξ)B j,q (η)P i,j (2.9)

A simplified form of Eq. 2.9 can be written as:

C(ξ, η) = p+1 i=1 q+1 j=1 B p,q i,j P i,j (2.10)
A Bézier solid can be defined as:

C(ξ, η, ζ) = p+1 i=1 q+1 j=1 z+1 k=1 B i,p (ξ)B j,q (η)B k,z (ζ)P i,j,k (2.11)
Similarly, a simplified form of Eq. (2.11) can be written as:

C(ξ, η, ζ) = p+1 i=1 q+1 j=1 z+1 k=1
B p,q,z i,j,k P i,j,k (2.12)

where ξ, η, ζ are the parameters defined in the range [0, 1], p, q, z represent the polynomial orders in ξ, η, ζ parametric directions, P i,j and P i,j,k are control point net and control point volume arrays, respectively.

According to Eq. (2.9), a Bézier surface can be constructed in Matlab, as shown in Fig. 2.4, in which the polynomial order in x and y directions are respectively 1 and 2. The main innovation of the B-spline basis functions comparing to Bernstein polynomial is in the knot vector, a set of coordinates defined in parametric space. In univariate Bernstein polynomial, we use the parameter ξ generally in the range [0, 1] to determine the basis function. In the B-spline basis function, however, we use a knot vector.

1) Knot vector

A knot vector in one dimension is defined as a series of non-decreasing coordinates in parametric space, denoted by Ξ = {ξ 1 , ξ 2 , ..., ξ n+p+1 }, where ξ i ∈ R is the ith knot (or coordinate), and i is the knot index from 1, 2, . . . , n + p + 1, in which n is the number of B-spline basis function along ξ parametric direction, and p is the polynomial order of B-spline basis function. In constructing B-spline surface and solid, it is necessary to use 2 and 3-knot vectors respectively directed along ξ and η directions. Each knot or coordinate of the knot vector is used to divide the parametric space of a geometrical model to obtain elements. Thus all of the mesh elements can be selected by knot values of knot vectors.

In terms of the space between different knots, a knot vector can be a uniform or non-uniform knot vector. In a uniform knot vector, the knots are equally spaced in the parametric space, such as

Ξ = {0.1, 0.2, 0.3, • • • , ξ n+p+1 }.
Similarly, in a non-uniform knot vector, the knots are unequally spaced in the parametric space, such as Ξ = {0.1, 0.15, 0.25, 0.3, . . . , ξ n+p+1 }. In a knot vector, there can be repeated knots, and a knot vector is said to be open if its first and last knots repetition are equal to the p + 1, in which p is the polynomial order of the basis function.

In one dimension, the basis functions constructed by an open knot vector interpolate the ends of parametric space. For example, in the following case, the knot vector and polynomial order are respectively Ξ = {0, 0, 0, 0.5, 1, 1, 1} and 2.

The polynomial order p + 1 is 3, which equals the number of first and last repeated knots, leading to the B-spline curve to interpolate the first and last control points, as shown in Fig. 2.8.

2) B-spline basis function

A B-spline basis function is defined by following Eqs. (2.13) and (2.14).

For polynomial order p = 0:

N i,0 (ξ) =    1 if ξ i ξ < ξ i+1 0 otherwise (2.13)
for p > 1:

N i,p (ξ) = ξ -ξ i ξ i+p -ξ i N i,p-1 (ξ) + ξ i+p+1 -ξ ξ i+p+1 -ξ i+1 N i+1,p-1 (ξ) (2.14)
The building processes of B-spline basis functions corresponding to the knot vector Ξ = {0, 0, 0, 0.5, 1, 1, 1} would be presented to explain the effects of the knot vector.

Starting with i = 1, polynomial order p = 0, obtain

N 1,0 (ξ) =    1 if ξ 1 ξ < ξ 2 0 otherwise (2.15)
as ξ 1 = ξ 2 = ξ 3 = 0, it can be observed that there is no value of ξ, such that 0 ξ and ξ < 0, and therefore N 1,0 (ξ) = 0.

Using the same logic to the remaining indices, the following piecewise constant functions can be obtained Now, consider polynomial order p = 1. Beginning with i = 1, the following linear functions can be obtained

N 2,0 (ξ) = 0 (2.16) N 3,0 (ξ) =    1 if 0 ξ < 0.5 0 otherwise (2.17) N 4,0 (ξ) =    1 if 0.5 ξ < 1 0 otherwise (2.
N 1,1 (ξ) = ξ -0 0 -0 N 1,0 (ξ) + 0 -ξ 0 -0 N 2,0 (ξ) (2.21) as N 1,0 (ξ) = N 2,0 (ξ) = 0, obtain N 1,1 (ξ) = 0. when i=2, obtain N 2,1 (ξ) = ξ -ξ 2 ξ 3 -ξ 2 N 2,0 (ξ) + ξ 4 -ξ ξ 4 -ξ 3 N 3,0 (ξ) =    (1 -2ξ) if 0 ξ < 0.5 0 otherwise (2.22)
when i=3, obtain

N 3,1 (ξ) = ξ -ξ 3 ξ 4 -ξ 3 N 3,0 (ξ) + ξ 5 -ξ ξ 5 -ξ 4 N 4,0 (ξ) = 2ξ    1 if 0 ξ < 0.5 0 otherwise + 2(1 -ξ)    1 if 0.5 ξ < 1 0 otherwise =          2ξ if 0 ξ < 0.5 2(1 -ξ) if 0.5 ξ < 1 0 otherwise (2.23) when i = 4, obtain N 4,1 (ξ) = ξ -ξ 4 ξ 5 -ξ 4 N 4,0 (ξ) + ξ 6 -ξ ξ 6 -ξ 5 N 5,0 (ξ) =    (2ξ -1) if 0.5 ξ < 1 0 otherwise (2.24)
when i = 5, obtain

N 5,1 (ξ) = ξ -ξ 5 ξ 6 -ξ 5 N 5,0 (ξ) + ξ 7 -ξ ξ 7 -ξ 6 N 6,0 (ξ) (2.25)
similarly, as N 5,0 (ξ) = N 6,0 (ξ) = 0, obtain N 5,1 (ξ) = 0.

Based on the obtained Eqs. Now, consider polynomial order p = 2. Beginning with i = 1, the following quadratic functions can be obtained

N 1,2 (ξ) = ξ -ξ 1 ξ 2 -ξ 1 N 1,1 (ξ) + ξ 4 -ξ ξ 4 -ξ 2 N 2,1 (ξ) =    ξ4-ξ ξ4-ξ2 ξ4-ξ ξ4-ξ3 if 0 ξ < 0.5 0 otherwise =    4( 1 2 -ξ) 2 if 0 ξ < 0.5 0 otherwise (2.26) when i = 2, obtain N 2,2 (ξ) = ξ -ξ 2 ξ 4 -ξ 2 N 2,1 (ξ) + ξ 5 -ξ ξ 5 -ξ 3 N 3,1 (ξ) =          ξ-ξ2 ξ4-ξ2 ξ4-ξ ξ4-ξ3 + ξ5-ξ ξ5-ξ3 ξ-ξ3 ξ4-ξ3 if 0 ξ < 0.5 ξ5-ξ ξ5-ξ3 ξ5-ξ ξ5-ξ4 if 0.5 ξ < 1 0 otherwise =          -6ξ 2 + 4ξ if 0 ξ < 0.5 2(1 -ξ) 2 if 0.5 ξ < 1 0 otherwise (2.27) when i = 3, obtain N 3,2 (ξ) = ξ -ξ 3 ξ 5 -ξ 3 N 3,1 (ξ) + ξ 6 -ξ ξ 6 -ξ 4 N 4,1 (ξ) =          ξ-ξ3 ξ5-ξ3 ξ-ξ3 ξ4-ξ3 if 0 ξ < 0.5 ξ-ξ3 ξ5-ξ3 ξ5-ξ ξ5-ξ4 + ξ6-ξ ξ6-ξ4 ξ-ξ4 ξ5-ξ4 if 0.5 ξ < 1 0 otherwise =          2ξ 2 if 0 ξ < 0.5 -6ξ 2 + 8ξ -2 if 0.5 ξ < 1 0 otherwise (2.28) when i = 4, obtain N 4,2 (ξ) = ξ -ξ 4 ξ 6 -ξ 4 N 4,1 (ξ) + ξ 7 -ξ ξ 7 -ξ 5 N 5,1 (ξ) =    ξ-ξ4 ξ6-ξ4 ξ-ξ4 ξ5-ξ4 if 0.5 ξ < 1 0 otherwise    4ξ 2 -4ξ + 1 if 0.5 ξ < 1 0 otherwise (2.29)
Based on the obtained Eqs. Important properties of B-spline basis functions can be concluded as:

(1) partition of unity ∀ξ,

n i=1 N i,p (ξ) = 1 (2) non negative ∀ξ, 0 ≤ N i,p (ξ)
(3) kronecker's delta N i,p (ξ j ) = δ ij , which is satisfied at patch boundary.

(4) linear independence

n i=1 α i N i,p = 0, if α 1 = α 2 = • • • = α n = 0
(5) variation diminishing property, leading to B-spline curves behaving monotonous with the increase of polynomial orders of basis functions.

(6) basis function continuity. pth order basis functions have p -1 continuous derivative C p-1 across the knots, element boundaries, under the condition of no repeated knots. With repeated knots, the pth order basis functions would constitute p -k -1 continuous derivatives C p-k-1 in which k is the multiplicity number of knots in a knot vector. For an open knot vector in which the multiplicity of the first and last knot is p + 1, the basis functions would possess C 0 continuity at the boundary of a patch.

3) B-spline curve

A B-spline curve is defined by the linear combination of B-spline basis function N i,p (ξ) and the corresponding control point P i , i = 1, 2, . . . , n, as shown in Eq. (2.30).

C(ξ) = n i=1 N i,p (ξ)P i (2.30)
Based on Eq. (2.30), a B-spline curve can be constructed in Matlab, as shown in Fig. 2.8, in which the polynomial order and number of control points are respectively 2 and 4.

C(ξ, η) = n i=1 m j=1 N i,p (ξ)M j,q (η)P i,j (2.31)
where N i,p (ξ) and M j,q (η) are the pth and qth order B-spline basis functions, n and m are respectively number of basis functions in ξ and η parametric directions. A compact form of Eq. (2.31) can be written as:

C(ξ, η) = n i=1 m j=1 N p,q i,j (ξ, η)P i,j (2.32)
For a control point volume array P i,j,k , i = 1, 2, ..., n, j = 1, 2, ..., m, k = 1, 2, ..., l, a B-spline solid is defined as:

C(ξ, η, ζ) = n i=1 m j=1 l k=1 N i,p (ξ)M j,q (η)R k,z (ζ)P i,j,k (2.33)
where R k,r (ζ) is rth order B-spline basis function, and l is the number of basis functions in ζ parametric direction. The simplified form of Eq. (2.33) can be written as:

C(ξ, η, ζ) = n i=1 m j=1 l k=1 N p,q,z i,j,k (ξ, η, ζ)P i,j,k (2.34)
Based on Eq. (2.31), a B-spline surface can be constructed in Matlab, as shown in Fig. 2.9. • III) NURBS curve, surface, and volume

1) NURBS basis function

Univariate NURBS basis function is described by the rationale of weighted B-spline basis functions as:

R i,p (ξ) = ω i N i,p (ξ) W (ξ) = ω i N i,p (ξ) n i=1 ω i N i,p (ξ) (2.35)
where ω i denotes the weight value of the ith control point P i , and W (ξ) is the weighted linear combination of B-spline basis functions.

The building procedure of the NURBS basis function with knot vector Ξ = {0, 0, 0, 0.5, 1, 1, 1}, and polynomial order p = 0 would be presented.

As the related equations of B-spline basis functions have been reviewed before, we only focus on the NURBS basis function.

R 1,2 (ξ) = ω 1 N 1,2 (ξ) W (ξ) = ω 1 N 1,2 (ξ) n4 i=1 ω i N i,2 = ω 1 N 1,2 (ξ) ω 1 N 1,2 (ξ) + ω 2 N 2,2 (ξ) + ω 3 N 3,2 (ξ) + ω 4 N 4,2 (ξ) =    4ξ 2 -4ξ + 1 if 0 ξ < 0.5 0 otherwise (2.36) R 2,2 (ξ) = ω 2 N 2,2 (ξ) W (ξ) = ω 2 N 2,2 (ξ) n4 i=1 ω i N i,2 = ω 2 N 2,2 (ξ) ω 1 N 1,2 (ξ) + ω 2 N 2,2 (ξ) + ω 3 N 3,2 (ξ) + ω 4 N 4,2 (ξ) =          -6ξ 2 +4ξ 4ξ 2 -4ξ+2 if 0 ξ < 0.5 2(1 -ξ) 2 if 0.5 ξ < 1 0 otherwise (2.37) R 3,2 (ξ) = ω 3 N 3,2 (ξ) W (ξ) = ω 3 N 3,2 (ξ) n4 i=1 ω i N i,2 = ω 1 N 1,2 (ξ) ω 1 N 1,2 (ξ) + ω 2 N 2,2 (ξ) + ω 3 N 3,2 (ξ) + ω 4 N 4,2 (ξ) =          2ξ 2 4ξ 2 -4ξ+2 if 0 ξ < 0.5 (-6ξ 2 + 8ξ -2) if 0.5 ξ < 1 0 otherwise (2.38) R 4,2 (ξ) = ω 4 N 4,2 (ξ) W (ξ) = ω 4 N 4,2 (ξ) n4 i=1 ω i N i,2 = ω 1 N 1,2 (ξ) ω 1 N 1,2 (ξ) + ω 2 N 2,2 (ξ) + ω 3 N 3,2 (ξ) + ω 4 N 4,2 (ξ) =    4ξ 2 -4ξ + 1 if 0.5 ξ < 1 0 otherwise (2.
C(ξ, η) = n i=1 m j=1 R p,q i,j (ξ, η)P i,j (2.41)
where the bivariate NURBS basis function R p,q i,j (ξ, η) is given as: R p,q i,j (ξ, η) = N i,p (ξ)M j,q (η)w i,j n i=1 m j=1 N i,p (ξ)M j,q (η)w i,j

(2.42)

where N i,p (ξ) and M j,q (ξ) are respectively the pth and qth order B-spline basis function, which are defined by the Eqs.

(2.21) and (2.27).

Similarly, NURBS solids are defined as

C(ξ, η, ζ) = n i=1 m j=1 l k=1 R p,q,z i,j,k (ξ, η, ζ)P i,j,k (2.43)
where the trivariate NURBS basis function R p,q,z i,j,k (ξ, η, ζ) is given as

R p,q,z i,j,k (ξ, η, ζ) = N i,p (ξ)M j,q (η)R k,z (ζ)w i,j,k n i=1 m j=1 l k=1 N i,p (ξ)M j,q (η)R k,z (ζ)w i,j,k (2.44)
Based on Eq. (2.41), a NURBS surface can be constructed in Matlab, as shown in Fig. 2.12.

Basic concepts related to IGA

In this part, some definitions of fundamental concepts related to IGA will be presented based on figures created in Ls Dyna and the paper of Agrawal [START_REF] Gautam | Iga: A simplified introduction and implementation details for finite element users[END_REF].

• I) Definitions and transformation of different spaces

Index space in two dimensions is equally divided no matter with knot values of knot vectors. In this case, the knot vectors are respectively Ξ 1 = {0, 0, 0, 0.5, 1, 1, 1} and Ξ 2 = {0, 0, 0, 0.5, 1, 1, 1}, in which the index space ranges from [0, 1] ( see Fig. 2.13 (a)). Parameter space in two dimensions is [0, 1] ⊗ [0, 1] domain where NURBS basis functions are defined ( see Fig. 2.13 (b)). Control points, physical mesh, and control mesh are defined in physical space ( see Fig.

(c))

.

ξ = 1 2 (ξ i+1 -ξ i )ξ + (ξ i+1 + ξ i ) η = 1 2 [(η i+1 -η i )η + (η i+1 + η i )] (2.45)
where ξ and η are the known integration points in the parent domain, which can be obtained based on the Gauss-Legendre quadrature rule.

In order to compute numerical integration, in the first step the mapping from the physical space to parameter space would be developed, based on Eq. (2.46) [START_REF] Gautam | Iga: A simplified introduction and implementation details for finite element users[END_REF]. where n e cp is the number of control points in an element Ω e .

J 1 =   ∂x
In the second step, the mapping from parameter to master spaces can be obtained, based on Eq. (2.49) [START_REF] Gautam | Iga: A simplified introduction and implementation details for finite element users[END_REF].

J 2 = ∂ξ ∂ξ ∂η ∂η (2.49)
The overall mapping processes can be explained by the following example where a random physical quantity f (x, y) is computed over the physical space, and the Gauss integration is performed over the parent domain [START_REF] Gautam | Iga: A simplified introduction and implementation details for finite element users[END_REF].

Ω f (x, y)dΩ = nel e=1 Ω e f (x, y)dΩ = nel e=1 Ωe f (ξ, η)|J 1 |dξdη, ( φ : Ω e → Ωe ) = nel e=1 Ω e f (ξ, η)|J 1 ||J 2 |dξdη, (φ : Ωe → Ω e ) = nel e=1 1 -1 1 -1 f (ξ, η)|J 1 ||J 2 |dξdη = nel e=1   n e gp i=1 f (ξ i , η i )ωp i |J 1 ||J 2 |   (2.50)
where n e gp and ωp i are respectively the number of Gauss points and their weights in element Ω e of the parent domain. In physical space of Fig. 2.16, there are two kinds of mesh models that are respectively control mesh and physical mesh. The control mesh is described by control points, and it controls an actual geometry model. The physical mesh is a real representation of an actual geometry model.

State of the art in IGA

At present, the researches in IGA are mainly in the aspect of isogeometric basis functions algorithm [START_REF] Sederberg | T-splines and t-nurccs[END_REF][START_REF] Bazilevs | Isogeometric analysis using t-splines[END_REF][START_REF] Buffa | Isogeometric analysis in electromagnetics: B-splines approximation[END_REF][START_REF] Dörfel | Adaptive isogeometric analysis by local h-refinement with t-splines[END_REF][START_REF] Wang | Converting an unstructured quadrilateral mesh to a standard t-spline surface[END_REF][START_REF] Nguyen-Thanh | Isogeometric analysis using polynomial splines over hierarchical t-meshes for two-dimensional elastic solids[END_REF][START_REF] Campen | Similarity maps and field-guided t-splines: a perfect couple[END_REF],

contact mechanics [START_REF] El-Abbasi | On the modelling of smooth contact surfaces using cubic splines[END_REF][START_REF] Jianli | Advances in isogeometric analysis[END_REF][START_REF] Lu | Isogeometric contact analysis: Geometric basis and formulation for frictionless contact[END_REF][START_REF] İlker Temizer | Contact treatment in isogeometric analysis with nurbs[END_REF][START_REF] İlker Temizer | Three-dimensional mortar-based frictional contact treatment in isogeometric analysis with nurbs[END_REF], fluid mechanics [START_REF] Bazilevs | Nurbs-based isogeometric analysis for the computation of flows about rotating components[END_REF][START_REF] Bazilevs | Isogeometric fluid-structure interaction analysis with applications to arterial blood flow[END_REF][START_REF] Bazilevs | Isogeometric fluid-structure interaction: theory, algorithms, and computations[END_REF], structural optimization [START_REF] Hassani | Application of isogeometric analysis in structural shape optimization[END_REF][START_REF] Shojaee | Isogeometric structural shape optimization using particle swarm algorithm[END_REF][START_REF] Wall | Isogeometric structural shape optimization[END_REF][START_REF] Cho | Isogeometric shape design optimization: exact geometry and enhanced sensitivity[END_REF][START_REF]Isogeometric shape design optimization using nurbs basis functions[END_REF][START_REF] Ha | Numerical method for shape optimization using t-spline based isogeometric method[END_REF][START_REF] Nagy | Isogeometric sizing and shape optimisation of beam structures[END_REF][START_REF] Nagy | Isogeometric design of elastic arches for maximum fundamental frequency[END_REF][START_REF] Seo | Shape optimization and its extension to topological design based on isogeometric analysis[END_REF][START_REF] Seo | Isogeometric topology optimization using trimmed spline surfaces[END_REF][START_REF] Manh | Isogeometric shape optimization of vibrating membranes[END_REF][START_REF] Qian | Full analytical sensitivities in nurbs based isogeometric shape optimization[END_REF], shell analysis [START_REF] Hartmann | About isogeometric analysis and the new nurbs-based finite elements in ls-dyna[END_REF][START_REF] Kiendl | Isogeometric shell analysis with kirchhoff-love elements[END_REF][START_REF] Benson | Isogeometric shell analysis: The reissner-mindlin shell[END_REF][START_REF] Benson | A large deformation, rotation-free, isogeometric shell[END_REF][START_REF] Uhm | T-spline finite element method for the analysis of shell structures[END_REF][START_REF] Kiendl | The bending strip method for isogeometric analysis of kirchhoff-love shell structures comprised of multiple patches[END_REF], damage and fracture mechanics [START_REF] Moon | Fatigue life evaluation of mechanical components using vibration fatigue analysis technique[END_REF][START_REF] Eldoǧan | Vibration fatigue analysis of a cantilever beam using different fatigue theories[END_REF][START_REF] Peng | Isogeometric boundary element methods for three dimensional static fracture and fatigue crack growth[END_REF], structural vibration analysis [START_REF] Cottrell | Isogeometric analysis of structural vibrations[END_REF][START_REF] Zhang | Efficient isogeometric formulation for vibration analysis of complex spatial beam structures[END_REF][START_REF] Wang | Novel higher order mass matrices for isogeometric structural vibration analysis[END_REF][START_REF] Shojaee | Free vibration analysis of thin plates by using a nurbs-based isogeometric approach[END_REF][START_REF] Gondegaon | Static structural and modal analysis using isogeometric analysis[END_REF][START_REF] Hien | Stochastic isogeometric analysis of free vibration of functionally graded plates considering material randomness[END_REF]. In the following, I will mainly present the IGA in the algorithm, structural optimization, shell, and vibration analysis.

• I) Isogeometric algorithm research NURBS is widely used in CAD systems due to its distinct advantages, in which the most promising is in representing the free-form shapes, such as cylinders, spheres, ellipsoids, and conical sections [START_REF] Bazilevs | Isogeometric analysis using t-splines[END_REF]. However, with the wide application, the NURBS basis function has been exposed to have some shortcomings. For example, NURBS using tensor product form makes local mesh refinement inefficient. Some structures like a cylinder surface cannot be represented by a simple sealed NURBS face, as there are gaps or overlaps at the junction of NURBS faces. These disadvantages must be improved so it can get a broader range of applications. Sederberg et cl. Through comparison with NURBS, T-splines have several advantages. However, the complexity of knot insertion for local refinement is far from trivial, particularly in 3D problems [START_REF] Dörfel | Adaptive isogeometric analysis by local h-refinement with t-splines[END_REF][START_REF] Nguyen-Thanh | Isogeometric analysis using polynomial splines over hierarchical t-meshes for two-dimensional elastic solids[END_REF]. Moreover, the linear independence of the T-spline basis functions is not available for generic T-meshes. Nguyen-Thanh et cl. [START_REF] Nguyen-Thanh | Isogeometric analysis using polynomial splines over hierarchical t-meshes for two-dimensional elastic solids[END_REF] proposed a new basis function for isogeometric analysis, that is, PHT-spline, which inherits almost all of the advantages of NURBS basis functions, including the non-negative, partition of unity, local support, etc., and developed isogeometric analysis based on PHT-spline basis functions on four models, and obtain that PHT-splines are superior to finite element cubic quadrilateral elements (FEA-Q16) and NURBS in terms of convergence rates and total errors, also PHT-splines support simpler local mesh refinement.

On the other hand, for obtaining a smooth surface from an input surface, usually, there are two processes, parametrization and construction of spline-based surface, in which the major challenge is in the parametrization process in designing an efficient algorithm to obtain the same topological structure as the input field. In 2017, Campen et cl. [START_REF] Campen | Similarity maps and field-guided t-splines: a perfect couple[END_REF] present a surface construction technique, T-splines with half-edge knots, which is a straightforward generalization of classical T-spline, and a class of parametrizations, a seamless similarity maps, from which a smooth piecewise rational surface with precisely the same input structure can be achieved for any given designs.

• II) Isogeometric design optimization

Structural design optimization is the eternal pursuit of the engineering community. It was initially achieved using classical analytical methods, such as variational or differential methods [START_REF] Jianli | Advances in isogeometric analysis[END_REF]. In 1960 Schmit [START_REF]Structural design by systematic synthesis[END_REF] firstly combined finite element analysis with nonlinear mathematical planning methods for structural design. With the widespread application of computer technology in structural analysis, structural design optimization has evolved from low-level size and shape optimization to high-level topological, structural, and multidisciplinary optimization stages [START_REF] Jianli | Advances in isogeometric analysis[END_REF]. In most processes of structural optimization, the finite element method is often used to compute structural response and sensitivity.

However, the approximate geometry description in classical FEA can cause accuracy issues. At present, B-spline and NURBS have been widely adopted in shape parameterization for structural optimization, since the work [START_REF] Braibant | Shape optimal design using b-splines[END_REF], in which the B-spline curves were applied to define design elements. The main reasons are because NURBS has a powerful capability in representing complex free form shapes and has a widespread application in CAD systems [START_REF] Jamshid | A survey of shape parameterization techniques[END_REF]. On the other hand, for finite element shape and topology optimization, the final design results depend on the mesh elements and extra post-processing steps to exchange data with the CAD system. The geometry and analysis models have different representations in the FEA field, making structural optimization more complicated. This urgently requires a close combination of CAD and CAE, and isogeometric analysis provides a promising way to combine design and analysis models. In 2005, Huges et cl. [START_REF] Hughes | Isogeometric analysis: Cad, finite elements, nurbs, exact geometry and mesh refinement[END_REF] firstly proposed the isogeometric shape optimization method and mentioned that it has a distinct advantage when compared with classical finite element counterparts. Wall et cl. [START_REF] Wall | Isogeometric structural shape optimization[END_REF] presented a framework of isogeometric shape optimization and applied it to the optimization problem of 3 classical cases, cantilever beam, plate with a hole, and a fully open end spanner structure, in which for the case of the cantilever beam, through the validation by an analytical solution, it has been obtained that the proposed method can provide a good result. Cho S et cl. [START_REF] Cho | Isogeometric shape design optimization: exact geometry and enhanced sensitivity[END_REF] proposed a continuum-based adjoint sensitivity method for isogeometric shape optimization and demonstrated the effectiveness and applicability of the isogeometric shape optimization method through the development of several cases.

Ha S H [START_REF]Isogeometric shape design optimization using nurbs basis functions[END_REF][START_REF] Ha | Numerical method for shape optimization using t-spline based isogeometric method[END_REF] introduced T-splines in isogeometric shape optimization to obtain optimal solution more efficiently and developed the T-spline and NURBS based isogeometric shape optimization on a bracket model. Nagy et cl. [START_REF] Nagy | Isogeometric sizing and shape optimisation of beam structures[END_REF] proposed isogeometric structural sizing and shape optimization approach on curved beam structures and also introduced a kind of multilevel approach, the combination of sizing and shape optimization, to solve the difficulties induced by structural sizing and shape optimization method. After Nagy et cl. extended the isogeometric design framework to optimize elastic arches for fundamental frequency maximization [START_REF] Nagy | Isogeometric design of elastic arches for maximum fundamental frequency[END_REF]. Yu-Deok Seo et cl. [START_REF] Seo | Shape optimization and its extension to topological design based on isogeometric analysis[END_REF] developed isogeometric shape optimization on three cases, a fillet, a plate with a hole, and a cantilever beam, in which for the cantilever beam problem, T-splines were adopted to solve the geometric and numerical instabilities, and then extended this approach to topology optimization, showing the capability to combine shape and topology optimization based splines. However, they found problems in terms of the computation time and systematic strategy for topological changes. To address these challenges, Yu-Deok Seo et cl. [START_REF] Seo | Isogeometric topology optimization using trimmed spline surfaces[END_REF] also proposed a new spline-based topology optimization based on isogeometric analysis. Through numerical analysis, it has been observed that the proposed approach is beneficial in expanding flexible design space, treating design-dependent load problems, and integrating design optimization framework. Manh N D et cl. [START_REF] Manh | Isogeometric shape optimization of vibrating membranes[END_REF] developed isogeometric shape optimization on vibrating membranes with the practical extension of a B-spline parametrization from the boundary of a domain onto its interior. Through the experiment with two numerical methods, respectively based on constructing a quasi-conformal mapping and a spring-based mesh model, it has been concluded that isogeometric analysis fits very well with shape optimization.

• III) Isogeometric shell analysis J. Kiendl et cl. [START_REF] Kiendl | The bending strip method for isogeometric analysis of kirchhoff-love shell structures comprised of multiple patches[END_REF] developed isogeometric bending analysis in Ls DYNA software on a thin shell structure comprised of multiple patches. As Kirchhoff-Love shell theory depends on higher-order basis functions and NURBS elements are usually c 1 or higher-order continuities, J. Kiendl et cl. [START_REF] Kiendl | The bending strip method for isogeometric analysis of kirchhoff-love shell structures comprised of multiple patches[END_REF] proposed Kirchhoff-Love shell theory based on isogeometric analysis, in which the bending strip method was introduced to treat c 0 patch boundaries, through adding strips of fictitious material with uni-directional bending stiffness and zero membrane stiffness at patch boundaries. The proposed framework demonstrated good analysis results through the application of a series of examples.

2011, Stefan et cl. [START_REF] Hartmann | About isogeometric analysis and the new nurbs-based finite elements in ls-dyna[END_REF] developed isogeometric shell analysis on an underbody cross member and compared the isogeometric analysis with finite element analysis under the condition of different polynomial orders and element size.

Through the comparison of draw-in length, contact forces between upper die and blank ( as shown in Fig. 2.19 (a))and CPU time, it was realized that IGA with higher polynomial order sometimes would not lead to better results, if the spacing of the control points will not be changed significantly. And under the condition of comparable discretization with standard linear finite elements, the NURBS shell elements can produce as good results as FEA with less CPU time.

For example, for the comparison of contact forces and CPU time, as shown in Fig. 2.19 (b), a higher polynomial order will not produce better results because of the constant control point space, and from Fig. 2.19 (c), it can be observed that the computation time with quadratic NURBS elements with a mesh size of 4 mm (P2-4 mm) is about 30% faster than the one with linear standard shell elements with a mesh size of 2 mm (Std-2 mm). It is also necessary to note that these analysis results were obtained without optimized isogeometric code conditions. Through comparison with FEA, it can be concluded that under adaptive mesh refinements, isogeometric analysis can produce as good results as a comparable discretization with standard linear finite elements. • IV) Isogeometric vibration analysis Cottrell et cl. [START_REF] Cottrell | Isogeometric analysis of structural vibrations[END_REF] firstly developed isogeometric vibration analysis on several basic structural modes such as rods, beams, membranes, plates, and three-dimensional solids to verify the effectiveness of IGA on vibration problems; from the comparison point of view, with the finite element vibration analysis, it can be observed that the isogeometric analysis results are better than the finite element analysis results. On the other hand, Cottrell et cl. [START_REF] Cottrell | Isogeometric analysis of structural vibrations[END_REF] also developed isogeometric vibration analysis on a NASA aluminum tested cylinder, made a comparison between numerical analysis and experimental results, and obtained a good agreement. Sangmesh et cl. [START_REF] Gondegaon | Static structural and modal analysis using isogeometric analysis[END_REF] presented complete details of calculation, explanation, and examples in isogeometric static and vibration analysis with the use of Matlab code and verified the results by FEA.

Xiangkui Zhang et cl. [START_REF] Zhang | Efficient isogeometric formulation for vibration analysis of complex spatial beam structures[END_REF] presented isogeometric modeling with global Timoshenko beam theory, then developed modal analysis on a complex beam structure with random curvature and torsion, in which a multi-patch non-uniform B-spline represents the beam. The isogeometric analysis results were verified by finite element analysis. The results confirmed the efficiency and precision of this formulation. From analysis results, it can be observed that the IGA produces results as same as the FEA with the use of fewer elements; the locking problem can be solved by the increase of the basis function orders; the number of the patch doesn't have effects on natural frequencies, which furtherly prove the effectiveness of this isogeometric formulations.

The locking problem is a notorious phenomenon in numerical simulation of beams and plates when elements used to analyze deep beams or thick plates are utilized to analyze slender or thin plates. This is featured that the obtained displacement is unreasonably small when calculating beams with large length-radius ratio [START_REF] Zhang | Efficient isogeometric formulation for vibration analysis of complex spatial beam structures[END_REF]. At present, several methods were introduced to solve the locking problem in static isogeometric analysis, including the new B method for volume and shear locking problems, RI method, order reduction method, the mixed formulation method, the DSG method, the single-variable method, assumed natural strain method, high-order non-uniform rational B-splines interpolation method, and couple polynomial field method.

2013, Dongdong et cl. [START_REF] Wang | Novel higher order mass matrices for isogeometric structural vibration analysis[END_REF] presented a method of higher-order mass matrices to improve the accuracy of isogeometric vibration analysis based on NURBS results, in which the higher-order mass matrix is obtained through a mixed mass matrix that was formulated through a linear combination of reduced bandwidth mass matrix and consistent mass matrix, and verified this method on several models. The analysis results showed that the higher-order matrix method is efficient and can produce high order of accuracy. 2017, Ta Duy et cl. [START_REF] Hien | Stochastic isogeometric analysis of free vibration of functionally graded plates considering material randomness[END_REF] applied isogeometric vibration analysis on a stochastic structure with random material properties in which Young's modulus and mass density were modeled as homogeneous Gaussian random fields. In this study, the IGA vibration formulation is presented in the appendices A.

• V) Summary on present studies in IGA From the papers above, obtain the following conclusions: I) in the isogeometric algorithm, the current research mainly focus on the improvement of basis functions properties to obtain a gap-free model and also develop a simple local mesh refinement method; II) in isogeometric shape optimization, 1) compared with finite element optimization, the isogeometric optimization has advantages in terms of the geometric representation and combination of CAD and CAE models; 2) most isogeometric design optimization is based on NURBS, but in some cases, T-spline based isogeometric optimization has also developed, as T-spline basis functions support local refinement and also seamless geometry, which are difficult to be accomplished using NURBS; 3) isogeometric topology optimization is available for specimen cases. 4) In IGA, both the control point location and weights have effects on a geometry representation, which can be explained by Fig. 2.20 which describes a NURBS surface consisting of the 4 × 3 control net. As shown in Figs. 2.20

(b) and (c), when control point position and weigh are altered from Q a to Q b and ω a to ω b , the related knot span and NURBS surface has also changed. However, until recently, for isogeometric shape optimization [START_REF] Cho | Isogeometric shape design optimization: exact geometry and enhanced sensitivity[END_REF][START_REF]Isogeometric shape design optimization using nurbs basis functions[END_REF], the analytical sensitivities are mainly calculated based on the design variables of control point location. This is mainly due to the complexity in the derivation of analytical sensitivities for control point weights [START_REF] Qian | Full analytical sensitivities in nurbs based isogeometric shape optimization[END_REF]. An analytical formulation of calculating sensitivities for both control point position and weight is necessary, this is because [START_REF] Qian | Full analytical sensitivities in nurbs based isogeometric shape optimization[END_REF]: a) to obtain more precise and effective computation of derivative information essential in gradient-based optimization. b) to achieve more flexible and compact shape representation, and increase design space, which can lead to obtaining a class of optimal design, including conic curves, circles, ellipses, and cylinders, spheres, ellipsoids; 5) Some scholars present formulation for the utilization of control point position and weights as design variables. It can be found that when compared with formulation just based on control point weight as a design variable, the proposed approach can obtain more robust optimization results.

However, it can also be found that isogeometric design optimization is mainly focused on some basic models. The main contribution of this study is that isogeometric design optimization is developed on an actual industrial model, and the elaborate process of the shape optimization is presented in chapter 4, which would provide easy access for a widespread application of the isogeometric design optimization. Ls DYNA software was chosen for this study. It is a high industrial software that has implemented IGA, and I was in touch with Ls DYNA developing team during my thesis preparation.

Comparison with FEA

The differences between IGA and FEA are the mesh elements formulation, choice of basis functions, and calculation of field variables in the governing equation. As shown in Figure 2.21, the mesh model obtained in FEA consists of mesh elements represented in physical space or parent domain [START_REF] Zienkiewicz | The finite element method: Its basis and fundamentals[END_REF]. In physical space, the mesh elements are defined by their nodal coordinates, and finite element basis functions, i.e. Lagrangian functions [START_REF] Gautam | Iga: A simplified introduction and implementation details for finite element users[END_REF], that interpolate the coordinate between the nodes. Each element in the physical space can be mapped to the parent domain through coordinate transformation, in which the width and length of the parent element are respectively 1, and the Gauss integration is performed on the parent element.

In IGA besides the usual physical mesh and mesh in the parent domain, the control mesh is introduced. The control mesh is defined by control points. It defines geometry, however, it does not have to coincide with the real geometry of a studied object. In two dimensions the control mesh is a bilinear quadrilateral element [START_REF] Hughes | Isogeometric analysis: Cad, finite elements, nurbs, exact geometry and mesh refinement[END_REF]. The physical mesh is a representation of actual geometry. It is obtained by the projection of control points with NURBS basis functions, and the discretization is governed by knot vectors discretization in the parent domain. A physical mesh model can consist of several patches, e.g. subdomains in which the section and material properties are the same. Next, IGA usually adopts the NURBS as its basis functions. Compared to Bernestain or B-spline basis functions, the NURBS basis functions possess more flexible properties and can exactly represent conic, circular and sphere models due to the effects of B-spline basis function and weight points associations [START_REF] Gautam | Iga: A simplified introduction and implementation details for finite element users[END_REF]. For FEA, the Lagrangian basis functions shown in Eq. (2.51) are used to approximate solutions at element nodes [START_REF] Zienkiewicz | The finite element method: Its basis and fundamentals[END_REF]:

L i,p (ξ) = p+1 k=1,k =i ξ -ξ k ξ i -ξ k , 1 ≤ k ≤ p + 1 and -1 ≤ ξ ≤ 1 (2.51)
Due to the different nature of basis functions, the description model also gives different results. The critical properties are in variation diminishing characteristics and inter-element boundary continuity. The variation diminishing property (the number of sign changes) is used to characterize the smoothness of a curve. For FEA, the Lagrangian basis functions can have any sign (-1 ≤ ξ ≤ 1, -1 ≤ L i,p (ξ) ≤ 1), so the oscillation of the fitting curve would be increased with the increase of the polynomial order, resulting in a non-smooth representation of the fitting curve (can not satisfy the variation diminishing property) and leading to contact problems between different description models. In IGA, the NURBS basis functions can satisfy non-negativity (-1 ≤ ξ ≤ 1, 0 ≤ R i,p (ξ)). The obtained results, same as the NURBS curves, are less sensitive to the polynomial orders and can present a smooth representation of the geometry (possess the variation diminishing property) and the contact surfaces.

Additionally, the NURBS basis functions present c p-1-k continuity, in which k is the number of repeating knots in a knot vector. However, the finite element basis functions are restricted to only c 0 continuity, leading to a non-smooth representation of the physical derivative quantities like stresses or strains.

Thirdly, in IGA, the combination of the control point P i and NURBS basis function are used to define NURBS-based elements, and field variables, such as displacement u in Eq. (A.4), are performed based on control points. For FEA, the Langrangian based finite elements are defined by the combination of elements node variable X and Langrangian basis function. The field variable U in finite element weak formulation is performed based on element nodes [START_REF] Gautam | Iga: A simplified introduction and implementation details for finite element users[END_REF][START_REF] Zienkiewicz | The finite element method: Its basis and fundamentals[END_REF].

Application on a simple example

The IGA is developed on a plate model, in which the FEA and analytical solution verify the isogeometric analysis results. The geometric model of the plate is created in Ls Dyna, as shown in Fig. 2.22. Due to the symmetric effects in geometry and applied load, only a quarter of the model is analyzed. During convergence analysis, different densities of mesh models are constructed to investigate the CPU time and maximum stress in y-direction to demonstrate the differences between IGA and FEA. The isogeometric and finite element mesh models are presented in Fig. 2.23. In Ls Dyna, a load can be applied on either the NURBS element, chosen by using the keyword, LOAD NURBS SHELL, or on control points. In this section, the load is carried on control points to simulate the traction effects on the plate model. In terms of element formulation, the isogeometric shell element with NURBS is adopted. For FEA, considering analysis accuracy, quadrilateral four nodes elements are used to create mesh models. Because of the applied load condition, it is not necessary to consider the load effects in the z-direction. On the other hand, the aspect ratio S w = W1 e of this plate is higher than 20, so plane stress shell element formulation, in which the normal stress σ z and shear stress τ zx and τ zy are assumed to be 0 is selected in static analysis.

During numerical simulation, the boundary conditions are the same for IGA and FEA. The translational constraints in the x-direction and y-direction are respectively imposed on the bottom and right edges of the model, which are marked by black bracket, as shown in Fig. 2.23.

Then 1 Mpa traction load is applied on the control point (nodes in FEA) on the left edge in the y direction to develop convergence analysis, which is marked by red bracket. Here it is necessary to note that in Ls Dyna, the keyword, LOAD NODE SET, is used to develop force load for each control point, which means that the pressure load should be transformed to force load. Then the obtained force load is divided by the number of control points (for FEA, divided by the number of elements nodes) to calculate the forces for each control point. (2.53). The CPU time for IGA and FEA are respectively 74 and 720 seconds. It has been found that IGA is more time-efficient compared with FEA, with a time difference of 89%. This is mainly because of mesh density, as IGA can provide similar analysis results with FEA using fewer NURBS elements.

All the numerical calculations presented are run on the computer Intel(R) Core(TM) i5-6440HQ CPU 2.60GHz with RAM 8Gb.

Relative error =

IGAresult -F EAresult F EAresult (2.53) The following well-known application [START_REF] Gautam | Iga: A simplified introduction and implementation details for finite element users[END_REF][START_REF] Milić | Isogeometric structural analysis based on nurbs shape functions[END_REF] is used to validate the isogeometric and finite element models developed in LS-DYNA software based on Kirsch's solution. A rectangular plate with a hole of radius a is loaded in the plane by a one-direction tension denoted by σ ∞ . The plane stress condition is considered. The plate parameters are depicted in Fig. 2.26 and the related analytical solution of the stress around the hole is given by [START_REF] Kirsch | Die theorie der elastizitat und die bedurfnisse der festigkeitslehre[END_REF]:

σ rr = σ ∞ 2 1 - a r 2 + σ ∞ 2 1 -4 a r 2 + 3 a r 4 cos 2θ σ θθ = σ ∞ 2 1 + a r 2 - σ ∞ 2 1 + 3 a r 4 cos 2θ τ rθ = - σ ∞ 2 1 + 2 a r 2 -3 a r
σ rr = 0 σ θθ = σ ∞ (1 -2 cos 2θ) τ rθ = 0 (2.55)
in addition, when θ=±90°, the proportion between the obtained maximum stress σ θθ and the applied stress σ ∞ is 3. This analytical solution is used in the following section to compare FEA and IGA.

Conclusion

In this chapter, the definition and theoretical backgrounds of IGA are reviewed. Then, IGA is performed on a simple plate model, and FEA and the analytical solution verify the analysis results. By comparing the results, such as the proportion between the obtained maximum stress and applied stress, it is observed that the IGA and FEA results have a good agreement with the analytical solution. Through the comparison of the CPU time, we can find that IGA is time-efficient, compared to FEA, with a time difference of 89%. The plate model (as shown in Fig. 3.2), loaded in bending in the static analysis, is firmly clamped on the bottom side and kept free on the other. The reduced section is used to localize the stress away from the clamp. A random acceleration is applied to the clamping fixture during fatigue analysis.

1) Convergence analysis

For convergence analysis, the von Mises stresses are computed to decide the convergence points. From Fig. 3.4, it can be observed that with the increase of control points and element nodes number, the von Mises stresses tend to reach stable value, in which fitting curves are used to approximate the obtained stresses. The computed correlation coefficients r between stress values and fitting curves are respectively 0.92 and 0.99.

For IGA, the maximum stress values start to converge from the mesh density of 52 (the number of control points);

for FEA, from the mesh density of 133 (the number of element nodes). After the isogeometric and finite element convergence points are decided, the analysis results of these points are obtained to compare the differences between IGA and FEA. comparison, it can be observed that the obtained isogeometric and finite element natural frequencies and vibration modes match very well. The relative error for the first isogeometric and finite element natural frequency is 0.9%. We can also observe that the relative error for the second and fourth natural frequencies is higher. In our opinion, the difference is due to the low correlation observed in Fig. 3.4 (a). Therefore, the number of control points chosen for the IGA leads to a small deviation of the von Mises stress in static (see Figure 3.5) and to a more consequent deviation in dynamic especially for the pair modes (e.g. torsion). However, only the first mode is useful for this study because it causes the most damage. For this mode, the deviation of the displacements given in Figs 3.6 and 3.7 is small. The frequency response corresponding to the ratio of the displacement over the input acceleration were used to follow the resonant frequency of the system. If a crack occurs, a shift and a decrease of magnitude of the frequency response is observed. Therefore, it is possible to detect a change in the response of the specimen by the evolution of the resonant frequency deviation defined by:

RF D() = 100 × f 0 r -f i r f 0 r (3.1)
where f 0 r corresponds to the initial resonant frequency of undamaged specimen and f i r is the ith measure of the resonant frequency during the fatigue test.

The strain signal in time recorded by the HBM device via the strain gauge is represented in Fig. 3.15 for the specimen S 1 (with a sampling rate 4800 per second). The related spectrogram given in Fig. 3.16 shows the evolution of resonant frequency during the exposure time. The deviation of this resonant frequency is clearly visible. pa respectively, leading to the relative error of -0.10%. From Fig. 3.26, it can be seen that the obtained isogeometric and finite element cumulative damage ratios are respectively 7.571e-3 and 7.341e-3, leading to the relative error of 3.13%.

The maximum damage ratios were located on similar elements. According to the Eq. (1.29), the expected isogeometric and finite element fatigue life E[T f ] are 1.5850e+03, and 1.6347e+03 seconds respectively. Based on own developed Matlab programming, the isogeometric and finite element damage ratios are 7.6e-03 and 7.4e-03, which have a good agreement to the damage ratios obtained from Ls Dyna.

Summary

In this part, the isogeometric and finite element convergence and random vibration fatigue analyses were developed on the L-shaped plate model, in which the obtained frequencies were compared with the ones from Pitoiset [START_REF] Pitoiset | Spectral methods for multiaxial random fatigue analysis of metallic structures[END_REF], and fatigue analysis results were validated by the Matlab programming.

In convergence analysis, the plate is clamped on the two edges (see Fig. The wind turbine tower model is created by a series of different thickness cylinders and conical shell sections, in which the geometry parameters such as the height, thickness, etc. are respectively displayed in the form of mm in the The obtained maximum, minimum, and original (no consideration of optimization) frequencies are displayed in the Tab. 4.4. When compared with the original state in terms of natural frequency, it can be found that through optimization, the natural frequency can be improved with a relative error of 1.76%. The obtained maximum, minimum, and original (no consideration of optimization) effective stress RMS are presented in Tab. 4.5. When compared with the original state in terms of natural frequency and stress values, it can be found that through optimization, the RMS can be reduced with a relative error of 6.07%, meaning that under the same applied load condition, based on material S-N curve, the optimized model have a longer fatigue life.

Conclusion

In this chapter, the isogeometric random vibration fatigue analysis, size and shape optimizations were developed on the wind turbine tower model. The design objectives were the maximum first natural frequency and minimum effective stress RMS. The design variables were respectively the tower segment thicknesses and control point position of the NURBS curves.

For the fatigue analysis, the tower model was clamped on the base flange, and random force PSDs were applied on the concentrated mass element. From modal analysis, it can be found that the obtained first five natural frequencies and vibration modes from IGA and FEA have a good agreement. Fatigue analyses show that the obtained isogeometric and finite element maximum effective stress RMS are 70.62 and 72.39 pa with a relative error of -2.45%, and cumulative damage ratios are 4.178e-5 and 4.112e-5 with a relative error of 1.61%. Based on the Matlab program, the isogeometric and finite element damage are respectively 4.2e-5. Thus, the damage comparison between IGA and FEA, numerical simulation, and Matlab programming shows that the obtained values agree very well.

From the sizing and shape optimization results, it can be observed that the natural frequency can be improved to 0.6996 (isogeometric size optimization), 0.6881 (finite element size optimization), and 0.692 Hz (isogeometric first shape optimization) from the original value of 0.680 Hz. It can be seen that the design variables have fewer effects on the first natural frequency, however have significant effects the effective stress RMS, determining a structure's damage ratio. From the isogeometric shape optimization results, it can be found that under consideration of mass constraints, the effective stress RMS value can be reduced with relative errors of 6.07% and 1.94%. Furthermore, based on proper fixed control points, we can develop shape optimization considering different door shapes.

Chapter 5

Conclusions

In this studying, isogeometric static, convergence, fatigue, size and shape optimization are respectively developed, in which the static, convergence and fatigue analysis results are validated by FEA, and Matlab programming.

Firstly, from the aspect of mesh generation, it is obvious that IGA is time-efficient compared with FEA. Normally, in the FEA field, the more complex a geometry, the more time spent in mesh generation. For example, for the L-shaped and wind turbine tower model, in FEA, the mesh generation time occupies respectively 50% and 80% of whole analysis time; in IGA, almost 10%.

In terms of convergence analyses, compared with FEA, the IGA normally predict the results with the use of less elements and CPU time. For example, for the convergence analysis on the simple plate model, the CPU time are respectively 231 seconds (IGA) and 462 seconds (FEA).

From the fatigue analysis on plate and tower models, it can be found that IGA can provide similar fatigue life with FEA, normally using fewer control points (compared with the nodes of FEA) and integration points in the thickness direction of the models. For example, for the plate model, during IGA and FEA, the number of integration points through-thickness direction are respectively 2 and 4.

From the isogeometric optimization, it can be observed that firstly, thorough the change of control point position, an ideal and realistic door shape can be easily obtained, which is normally difficult with the design variables of element nodes; secondly, because of the combination of CAE and CAD models, after calculation, there is no need to post-processing steps to exchange date with the CAD systems.

In the future, it is necessary to consider the IGA in the following aspects: 1) complex geometrical modeling. As the tower flange is simplified in this work, next step, we will consider the IGA on a more realistic tower model; 2) During isogeometric static analysis, a load can be applied directly on control points, and also on NURBS elements.

The two cases would be explained in the following respectively. The setting process in modal analysis and fatigue analysis will be presented.

1) Modal analysis

In the first step, the number of eigenvalues to extract can be set, based on the keyword, IMPLICIT_EIGENVALUE (

Figure 1 :

 1 Figure 1: Schematic illustration of the thesis outline

2 . 5

 25 Constant basis functions corresponding to Ξ = {0, 0, 0, 0.5, 1, 1, 1} . . . . . . . . . . . . . . . . . . 2.6 Linear basis functions corresponding to Ξ = {0, 0, 0, 0.5, 1, 1, 1} . . . . . . . . . . . . . . . . . . . 2.7 Quadratic basis functions corresponding to Ξ = {0, 0, 0, 0.5, 1, 1, 1} . . . . . . . . . . . . . . . . . 2.8 B-spline curve corresponding to the knot vector Ξ = {0, 0, 0, 0.5, 1, 1, 1} . . . . . . . . . . . . . . . 2.9 B-spline surface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) R 1,2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) R 2,2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.10 NURBS basis functions corresponding to Ξ = {0, 0, 0, 0.5, 1, 1, 1} . . . . . . . . . . . . . . . . . . (c) R3,2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (d) R 4,2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.11 NURBS curve . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.12 NURBS surface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.13 Schematic illustration of different spaces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) Index space . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) Parameter space . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (c) Physical space . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.14 Schematic illustration of the space transformations . . . . . . . . . . . . . . . . . . . . . . . . . . 2.15 Schematic illustration of patch elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) Model 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) Model 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.16 Schematic illustration of mesh elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.17T-spline repairment (from [16]) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) B-spline surfaces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) T-splines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.18 NURBS, T-spline control mesh and surfaces (from [17]) . . . . . . . . . . . . . . . . . . . . . . . . (a) NURBS and T-spline control mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) NURBS and T-spline surfaces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xii (a) Underbody cross member . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.19 Isogeometric and finite element analysis comparison on a underbody cross member (from [18]) . . . (b) Comparison of contact forces between upper die and blank . . . . . . . . . . . . . . . . . . . . (c) Comparison of necessary CPU-time . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.20 NURBS surfaces (from [19]) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) Initial NURBS surface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) Modified NURBS surface after the position of control point Q is changed from Q a to Q b . . . . (c) Modified NURBS surface after the weight of control point ω is changed from ω a to ω b . . . . . 2.21 Meshing process in isogeometric and finite element analysis. . . . . . . . . . . . . . . . . . . . . . 2.22 The schematic diagram of of the plate model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.23 Mesh models and boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) IGA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) FEA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.24 Convergence analysis results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) IGA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) FEA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.25 Maximum stresses in y directions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) IGA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) FEA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.26 Schematic diagram of a plate model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Schematic illustration of the third Chapter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 The plate model with reduced section . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Mesh models and boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) IGA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) FEA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Convergence analysis results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) IGA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) FEA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 The effective stress . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) IGA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) FEA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 Isogeometric first five vibration modes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) First . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) Second . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (c) Third . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (d) Fourth . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (e) Fifth . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) First . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) Second . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (c) Third . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.7 Finite element first five vibration modes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (d) Fourth . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (e) Fifth . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.8 Applied acceleration PSD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.9 Effective stress PSD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.10 The effective stress RMS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) IGA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) FEA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1 β

 1 xvi (a) Wind force PSD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) Normalized wind force PSD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C.5 Eurocode 1 force PSDs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) Longitudinal force PSDs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) Lateral force PSDs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (c) Vertical force PSDs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (d) Vertical force PSDs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . D.1 Schematic illustration of the overall analysis procedures . . . . . . . . . . . . . . . . . . . . . . . . D.2 Schematic illustration of the fatigue analysis procedure . . . . . . . . . . . . . . . . . . . . . . . . xvii List of Tables2.1 Material properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Material properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 The first five natural frequencies (Hz) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 The isogeometric and finite element cumulative damage ratios . . . . . . . . . . . . . . . . . . . . 3.4 The fatigue life E[T f ] results related to the 5% variation of σ f variation (in the units of hours: minutes: seconds) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Results of the fatigue life E[T f ] obtained from the resonant frequency deviation . . . . . . . . . . . 3.6 Material properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.7 The first five natural frequencies from [5] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.8 The first five naural frequencies (Hz) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .4.1 Main parameters of the wind turbine tower . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Material properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 The first five natural frequencies (Hz) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 The first natural frequencies (Hz) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 The effective stress RMS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.6 The effective stress RMS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xviii Nomenclature FEA: Finite Element Analysis IGA: Isogeometric Analysis NURBS: Non-Uniform Rational B-spline PSD: Power Spectrum Density FRF: Frequency Response Function CDF: Cumulative Distribution Function PDF: Probability Density Function SSS: Strict Sense Stationary random process WSS: Wide Sense Stationary random process R: Fatigue cycle characteristic N : Fatigue cycle number S rN : Fatigue limit under the cycle characteristic r and cycle number N C, β, σ f = C , b = -1 β : Material constants K N : Fatigue life coefficient N b : Base cycle number σ b : Base fatigue limit (endurance limit) σ u : Ultimate strength N : Applied cycle number N f : Number of cycles allowable at a particular stress amplitude [σ -1 ]: Allowable stress under a symmetric cycle [σ 0 ]: Allowable stress under a pulsating cycle σ -1 , (σ -1 ) k : Material and component fatigue limit under symmetric cycles (σ -1 ) d , (σ -1 ) d0 : Symmetric cycle fatigue limit under the specimen diameters of d and d 0 (σ -1 ) β S , (σ -1 ) β0 : Symmetric cycle fatigue limit under certain and fine polishing surface states K σ : Component stress concentration coefficient σ : Component dimension coefficient β σ : Component surface state coefficient (K σ ) D : Comprehensive impact coefficient σ a : Stress amplitude σ ar : Obtained stress amplitude from the Goodman equation σ m : Mean stress s d : Safety coefficient σ c : von Mises stress σ = (σ 1 , σ 2 , σ 3 , σ 12 , σ 23 , σ 31 ) T : Stress vector σ 1 , σ 2 , σ 2 : Normal stresses σ 12 , σ 23 , σ 13 : Shear stresses Q: Matrix used in definition of von Mises stress σ(t)(t ∈ T ): A random process σ i , i = 1, 2, . . . : Random variable t ∈ T : Time parameter τ : Time shift xix E[σσ T ]: Covariance matrix of the stress vector σ Φ σσ (ω): PSD of the stress vector Φ c (ω): PSD of the equivalent von Mises stress D: Total damage for deterministic loading E[N + 0 ]: Expected number of zero-crossing per unit time E[T f ]: Fatigue life for random loading (sec.) F σ (r): Cumulative density function of stress E[p]: The number of peak per second E[D] N B : Cumulative damage ratio defined by Benda's PDF E[D] W : Wirching's cumulative damage ratio expression (σ eq ) Hancock , (σ eq ) C and D , (σ eq ) Steinberg : Equivalent stress expression given by Hancock, Chaudhuery and Dover, and Steinberg p σB (r), p σT (r), p σD (Z): Benda, Tuna, Dirlik's PDF γ: Irregularity factor e, cp: Element and control point f : Frequency (Hz) √ m 0 : Root mean square value of stress m k : kth order spectral moment n e : Number of element n e cp : Number of control points in an element e p, q, z: Polynomial order along ξ, η and ζ parametric directions

  bivaraite, trivariate NURBS basis functions W (ξ): Weighted linear combination of B-spline basis functions n, m, l: Number of basis functions along ξ, η, ζ parametric directions P i , P i,j , P i,j,k : Control point vector, net and volume arrays n e cp : Number of control points in element Ω e n e : Number of nodes in element Ω e Ω e , Ωe , Ω e : An element in the physical space, parameter space, and parent domain n e gp : Number of Gauss points K: Stiffness matrix M: Mass matrix F: Force vector e = (e 1 , e 2 , e 3 ) T : Orthonormal basis f: Deterministic volume force n: Outward normal vector t: Deterministic surface force u, u e : Global and element displacement vector u, ü: Velocity and acceleration vectors x = (x, y, z) T : Global Cartesian reference system B: Strain displacement matrix C: Damping matrix C: Fourth order material constitutive tensor F, F e : Global and element external force vector |H(f ) : Magnitude of frequency response matrix K, K e : Global and element stiffness matrix P i : ith control point vector xx R: NURBS basis functions matrix S FF (f ): External force PSD matrix S uu (f ): Displacement PSD matrix S σσ (f ): Stress PSD matrix ε: Strain tensor η: Parametric coordinates given by ξ in one dimension, (ξ,λ) in 2 dimensions and (ξ, λ, ζ) in 3 dimensions σ: Cauchy stress tensor Ξ, Λ: Knot vectors E[.]: Expected value operator δ : Use for virtual quantities ∇:

Figure 1 . 1 :

 11 Figure 1.1: Schematic illustration of the first chapter

Figure 1 . 3 :

 13 Figure 1.3: Deterministic fatigue analysis procedures

Figure 1 . 4 :

 14 Figure 1.4: Schematic illustration of a random process

  ) where C and β are deterministic. To determine the expected value of |σ(t)| β , we can use the expectation definition of random variable in HCF domain such as: E |σ(t)| β = +∞ 0 r β p σ (r)dr (1.27)

1. 3 . 4

 34 Fatigue analysis under different probability density functionsA signal can be classified into the narrow-band (NB) or wide-band (WB), based on the bandwidth value (the difference between the upper and lower frequencies). The NB refers to the signal, whose maximum bandwidth is 3400 Hz, as illustrated by the solid line in Fig.1.5. The WB is the combination of NB and upper-band (UP) and has the bandwidth of 50-7000 Hz, defined by the international telecommunication union (ITU)[START_REF] Ding | Wideband audio over narrowband low-resolution media[END_REF].

Figure 1 . 5 :

 15 Figure 1.5: Bandwidth definition

Figure 1 . 6 :

 16 Figure 1.6: Vibration fatigue analysis procedures in the frequency domain

Figure 2 . 3 :

 23 Figure 2.3: Bézier curve

Figure 2 . 4 :

 24 Figure 2.4: Bézier surface

18 )

 18 Fig. 2.5, based on Eqs. (2.15) to (2.20), the constant basis functions can be plotted in corresponding index space, in which it can be observed that function N 3,0 (ξ) and N 4,0 (ξ) are 1 in index space 0 ξ < 0.5 and 0.5 ξ < 1, respectively.

Figure 2 . 5 :

 25 Figure 2.5: Constant basis functions corresponding to Ξ = {0, 0, 0, 0.5, 1, 1, 1}

  (2.21) to (2.25), the corresponding basis functions can be constructed in Matlab, as shown in Fig. 2.6.

Figure 2 . 6 :

 26 Figure 2.6: Linear basis functions corresponding to Ξ = {0, 0, 0, 0.5, 1, 1, 1}

( 2 .

 2 [START_REF] Schijve | Four lectures on fatigue crack growth[END_REF] to(2.29), the corresponding basis functions can be constructed in Matlab, as shown in Fig.2.7.

Figure 2 . 7 :

 27 Figure 2.7: Quadratic basis functions corresponding to Ξ = {0, 0, 0, 0.5, 1, 1, 1}

Figure 2 . 8 :

 28 Figure 2.8: B-spline curve corresponding to the knot vector Ξ = {0, 0, 0, 0.5, 1, 1, 1}

Figure 2 . 9 :

 29 Figure 2.9: B-spline surface

39 ) 2 Figure 2 . 11 :

 392211 Figure 2.11: NURBS curve

Figure 2 .

 2 Figure 2.12: NURBS surface

Figure 2 . 13 :

 213 Figure 2.13: Schematic illustration of different spaces

Figure 2 . 14 :

 214 Figure 2.14: Schematic illustration of the space transformations

Figure 2 . 15 :Figure 2 . 16 :

 215216 Figure 2.15: Schematic illustration of patch elements

Figure 2 . 17 :

 217 Fig.2.17 (b), which can save much time spent in repairing gaps and so improve the precision in analysis results.

  (a) NURBS and T-spline control mesh (b) NURBS and T-spline surfaces

Figure 2 . 18 :

 218 Figure 2.18: NURBS, T-spline control mesh and surfaces (from [17])

  (a) Underbody cross member (b) Comparison of contact forces between upper die and blank (c) Comparison of necessary CPU-time

Figure 2 . 19 :

 219 Figure 2.19: Isogeometric and finite element analysis comparison on a underbody cross member (from [18])

Figure 2 . 20 :

 220 Figure 2.20: NURBS surfaces (from [19])

Figure 2 . 21 :

 221 Figure 2.21: Meshing process in isogeometric and finite element analysis.

2. 5 . 1

 51 Analysis preparation• I) Geometric modelAs shown in Fig.2.22, the first analyzed model is a rectangular plate with the thickness e, length L 1 , width w 1 , and radius R r of 0.001, 10, 2, and 0.25 (m).The material of this model is steel DC01, and properties are given in Tab. 2.1.

Figure 2 . 22 :

 222 Figure 2.22: The schematic diagram of of the plate model

Figure 2 . 24 :

 224 Figure 2.24: Convergence analysis results

Figure 2 . 25 :

 225 Figure 2.25: Maximum stresses in y directions

Figure 2 . 26 :

 226 Figure 2.26: Schematic diagram of a plate model

Figure 3 . 1 :

 31 Figure 3.1: Schematic illustration of the third Chapter

Figure 3 . 2 : 3 3)

 323 Figure 3.2: The plate model with reduced section

Figure 3 . 3 :

 33 Figure 3.3: Mesh models and boundary conditions

Figure 3 . 4 :

 34 Figure 3.4: Convergence analysis results

Figure 3 . 5 :

 35 Figure 3.5: The effective stress

Figure 3 . 6 :

 36 Figure 3.6: Isogeometric first five vibration modes

Figure 3 . 7 :

 37 Figure 3.7: Finite element first five vibration modes

Figure 3 . 8 :

 38 Figure 3.8: Applied acceleration PSD

Figure 3 . 9 :

 39 Figure 3.9: Effective stress PSD

Figure 3 . 10 : 60 (

 31060 Figure 3.10: The effective stress RMS

Figure 3 . 11 :

 311 Figure 3.11: Isogeometric cumulative damage ratios

3. 2 . 2 Figure 3 . 13 :

 22313 Figure 3.13: Schematic representation of the system

Figure 3 . 14 :

 314 Figure 3.14: Frequency responses during fatigue test (4 measures)

  (a) Complete signal during fatigue test Table 3.5: Results of the fatigue life E[T f ] obtained from the resonant frequency deviation Id of specimen Fatigue life (in hour) time (b) Versus number of cycle

Figure 3 . 17 :

 317 Figure 3.17: Resonant frequency deviation during fatigue test

Figure 3 . 21 : 71 (

 32171 Figure 3.21: Isogeometric first five vibration modes

Figure 3 . 22 :

 322 Figure 3.22: Finite element first five vibration modes

Figure 3 . 23 :

 323 Figure 3.23: Applied load PSD

Figure 3 . 24 :Figure 3 . 25 :Figure 3 . 26 :

 324325326 Figure 3.24: Effective stress PSD

Fig. 4 . 3 . 3 Figure 4 . 3 :

 43343 Fig. 4.3. The tower model is assembled by five flange connections, whose base, second, third, fourth, and top flange thicknesses are 330, 200, 200, 200, and 110 mm. At the base of the tower, a door opening is taken into account. The detailed geometrical information can be found in the reference [97].

Figure 4 . 5 :

 45 Figure 4.5: Isogeometric first five vibration modes

Figure 4 . 6 : 2 )

 462 Figure 4.6: Finite element first five vibration modes

Figure 4 . 7 :

 47 Figure 4.7: Effective stress PSD

Figure 4 . 8 :

 48 Figure 4.8: Effective stress RMS

  (a) Sizing optimization result (b) Sizing optimization result with consideration of mass constraint

Figure 4 . 12 :Fig. 4 .

 4124 Figure 4.12: Finite element sizing optimization results

Table 4 . 4 :Fig. 4 .

 444 Fig. 4.15 (a) present the shape optimization results, in which z label shows the obtained effective stress RMS under different door shapes. It can be seen that under the mass constraint (see Fig. 4.16 (b)), the minimum and maximum effective stresses are situated on the points (0.008, -0.002, 6.633e+1) and (-0.046, 0.058, 1.505e+2), marked by red diamonds in the Fig. 4.15 (c).

Table 4 . 5 :

 45 The effective stress RMS

Figure 4 . 16 :

 416 Figure 4.16: Isogeometric shape optimization results

Figure 4 . 17 :Fig. 4 .

 4174 Figure 4.17: Illustration of fixed control points

Figure 4 . 19 :

 419 Figure 4.19: Isogeometric shape optimization results

Figure B. 3 :

 3 Figure B.3: Element shell NURBS patch

1) case 1 :(Figure B. 4 :

 14 Fig. B.5 (b)). In the final step, the defined load can be applied on the selected control points though the keyword, LOAD ( Fig. B.5 (c)), in which the load can be applied in different directions through the selection of Degree of Freedom (DOF).2) case 2: load applied on NURBS elements

Figure B. 6 :

 6 Figure B.6: load applied on NURBS elements

Fig. B. 7 (

 7 Fig. B.7 (a)), which would be used to activate eigenvalue analysis to compute eigen modes, which would be saved in a binary database "d3eigv"; the analysis method can be selected, using the keyword IMPLICIT_GENERAL( Fig. B.7 (b)), in which the time step for implicit analysis can also be defined; the analysis solution for implicit analysis can be chosen by the keyword, IMPLICIT_SOLUTION, ( Fig. B.7 (c)).

Figure B. 7 :

 7 Figure B.7: Modal analysis precedures

Figure B. 8 :

 8 Figure B.8: Fatigue analysis procedures

1

  Schematic illustration of the thesis outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi 1.1 Schematic illustration of the first chapter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Schematic illustration of a material S-N curve . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Deterministic fatigue analysis procedures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Schematic illustration of a random process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1.5 Bandwidth definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.6 Vibration fatigue analysis procedures in the frequency domain . . . . . . . . . . . . . . . . . . . . 2.1 Schematic illustration of the second section . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Bernstein basis functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Bézier curve . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Bézier surface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  xiii (a) Steinberg . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) Dirlik . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (c) Narrow band . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (d) Wirshing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (e) Chaudhury and Dover . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (f) Tunna . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.11 Isogeometric cumulative damage ratios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (g) Hancock . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.12 Finite element cumulative damage ratios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) Steinberg . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) Dirlik . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (c) Narrow band . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (d) Wirshing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (e) Chaudhury and Dover . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (f) Tunna . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (g) Hancock . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.13 Schematic representation of the system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .14 Frequency responses during fatigue test (4 measures) . . . . . . . . . . . . . . . . . . . . . . . . . .16 Spectrogram obtained from strain measure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.17Resonant frequency deviation during fatigue test . . . . . . . . . . . . . . . . . . . . . . . . . . . . Fourth . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.22 Finite element first five vibration modes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (e) Fifth . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.23 Applied load PSD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.24 Effective stress PSD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv 3.25 Effective stress RMS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) IGA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) FEA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.26 Cumulative damage ratio . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) IGA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Mesh models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) IGA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) FEA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (c) Flange model for IGA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (d) Flange model for FEA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) First . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) Second . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (c) Third . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 Isogeometric first five vibration modes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (d) Fourth . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (e) Fifth . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) First . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) Second . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (c) Third . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.6 Finite element first five vibration modes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Cumulative damage ratio . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) IGA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) FEA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.10 Sizing optimization procedures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.11 Isogeometric sizing optimization results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) Sizing optimization result . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) Sizing optimization result with consideration of mass constraint . . . . . . . . . . . . . . . . . 4.12 Finite element sizing optimization results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Fixed control points . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) Without constraint . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.14 Isogeometric shape optimization results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) With mass constraint . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (c) With mass constraint . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Without constraint . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.15 Isogeometric shape optimization results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) With mass constraint . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (c) With mass constraint . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.16 Isogeometric shape optimization results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . With mass constraint . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.18 Isogeometric shape optimization results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (c) With mass constraint . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.19 Isogeometric shape optimization results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) With optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) With optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) Plate construction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) Circle construction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (c) Geometry trimming . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.1 The geometry construction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (d) Created model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (e) Final model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.2 The mesh model construction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) Rough mesh model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Boundary condition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.6 load applied on NURBS elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.7 Modal analysis precedures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Wind velocity against the height in Normandy region . . . . . . . . . . . . . . . . . . . . . . . . . C.2 The wind speed PSDs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C.3 The normalized wind speed PSDs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	(a) Versus time . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) Versus number of cycle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.18 Dimensions of the plate (in mm) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.19 Mesh models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) IGA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) FEA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (c) Splitted geometry model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.20 Convergence analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) IGA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) FEA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) First . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) Second . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (c) Third . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (d) Fourth . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.21 Isogeometric first five vibration modes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (e) Fifth . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) First . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) Second . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (c) Third . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1 Schematic illustration of the fourth chapter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Layout of a 45 KW wind turbine tower with a height of 38 m . . . . . . . . . . . . . . . . . . . . . 4.3 Geometric model (in mm) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 (d) Fourth . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (e) Fifth . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.7 Effective stress PSD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.8 Effective stress RMS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) IGA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) FEA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.9 (a) Sizing optimization result . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) Sizing optimization result with consideration of mass constraint . . . . . . . . . . . . . . . . . 4.13 NURBS curves for the door . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) NURBS curves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) Same control points . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) (a) Original state . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) With optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (c) With optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.17 Illustration of fixed control points . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) Without constraint . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) (b) Mesh refinement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (c) Final mesh model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.3 Element shell NURBS patch . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.4 Material, section properties and boundary condition setting . . . . . . . . . . . . . . . . . . . . . . (a) Material properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) Section properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (c) Boundary condition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.5 The setting procedures of load application on control points . . . . . . . . . . . . . . . . . . . . . . (a) Material properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) Section properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (c) (a) Material properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) Section properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (c) Boundary condition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.8 Fatigue analysis procedures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) First . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) Second . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (c) Third . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (d) (b) FEA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (c) xv C.1

(a) The vibration-based bending fatigue bench . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) The schematic details . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3(a) Complete signal during fatigue test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.15 Strain measure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) A focus of signal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3C.4 The force PSD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  The IGA and FEA are developed on a simple plate model in Section 2.5, in which the convergence analysis is carried on with the criterion of maximum stress, and the results are validated by the analytical solution.

		Chapter 2:
		IGA approach
	Section 2.2	Section 2.3	Section 2.4	Section 2.5
	Theoretical	State of the	Comparison	Application on a
	backgrounds	art in IGA	with FEA	simple example
	The definition and properties of Bézier,B-spline, NURBS basis function, curve, surface and solids; Basic concepts related to IGA	Isogeometric algo-rithm; isogeometric vibration analy-sis; isogeometric shape optimiza-tion; Isogeometric shell analysis	The differences between IGA and FEA in terms of mesh elements, basis functions, and calculation of field variables	IGA and FEA developed on a plate model: convergence analysis with the criterion of stress, comparison with analytical solution
		Figure 2.1: Schematic illustration of the second section
	2.2 Theoretical background

.1, in Section 2.2, the isogeometric theoretical backgrounds are presented. The present research state in IGA is reviewed in Section 2.3 through the discussion of the isogeometric algorithm, shape optimization, shell and vibration analysis. In Section 2.4, the differences between IGA and FEA are described.

2.2.1 Bézier, B-spline, NURBS

• I) Bézier curve, surface, and solid A Bézier curve C(ξ) (ξ ∈ [0, 1]) has been widely employed in practical engineering including mechanical design and computer graphics

[START_REF] Choi | Path planning based on bézier curve for autonomous ground vehicles[END_REF][START_REF] Elhoseny | Bezier curve based path planning in a dynamic field using modified genetic algorithm[END_REF]

. It is a parametric curve defined by Bernstein polynomials and control points. A Bézier surface

Table 2 . 1 :

 21 Material properties

	Mass density Young's modulus poisson's ratio
	7800 Kg/m 3	2.05e+11 Pa	0.3
	• II) Mesh models and boundary conditions

Table 3 .

 3

		2: The first five natural frequencies (Hz)
		1	2	3	4	5
	IGA	211 1282 1510 1623 4253
	FEA	209 1141 1506 1579 4266
	relative error(%) 0.9 12.4 0.3	2.8	-0.3
			56

Table 3 . 3 :

 33 The isogeometric and finite element cumulative damage ratios Fatigue methods Steinberg Dirlik Narrow band Wirsching Chaudhurry and Dover Tunna Hancock

	IGA	1.180	1.449	34.25	22.14	32.61	1.424	32.73
	FEA	1.197	1.467	35.10	21.49	32.26	1.440	32.24
	relative errors (%)	1.44	1.24	-2.42	3.02	1.08	1.12	1.52

Table 3 . 4 :

 34 The fatigue life E[T f ] results related to the 5% variation of σ f variation (in the units of hours: minutes: seconds)

	σ f (Mpa) 581.685	612.3	642.915
	IGA	2:19:40 3:27:02 4:48:43
	FEA	2:17:06 3:24:30 4:46:48

Table 3 . 8 :

 38 The first five naural frequencies (Hz)

	Method CPU	1	2	3	4	5
	IGA	1 second 28.96 112.28 132.93 298.10 304.88
	FEA	1 second 28.63 114.42 131.30 299.63 303.78

Table 4 . 1 :

 41 3.3), and based on obtained first natural frequency, the convergence points are decided. The results show that the IGA and FEA convergence analyses predict similar modal analysis results. The isogeometric and finite element first natural frequency agree very well with the relative error of 1.15%. On the other hand, the obtained first natural frequencies from IGA and FEA are close to the one from Pitoiset, and the values are respectively 28.96, 28.63, and 27.7 Hz.In fatigue analysis, a random vibration acceleration in the vertical direction to the plate surface is applied on the clamping fixture to simulate the base acceleration. The simulation results show that the obtained maximum RMS and cumulative damage ratios from IGA and FEA have a good agreement, with the relative errors of -1.03%, 3.1% respectively. Based on the Matlab program, the predicted isogeometric and finite element damage ratios are 7.6e-03 and 7.4e-03, which are close to the damage values calculated from Ls Dyna, 7.571e-3 and 7.314e-3. Main parameters of the wind turbine tower

	Item	Value
	Rated power(KW)	450
	Number of blades	3
	Tower height(m)	38
	Tower top diameter(m)	1.82
	Tower bottom diameter(m)	2.8
	Rotor and nacelle assembly mass (g) 4.023e+7
	• I) Analysis preparation
	1) Geometric model, and material properties

Table 4 . 6 :

 46 The effective stress RMS

	Maximum minimum original
	162.7	68.84	70.62

sin 2θ(2.54)

(b) With mass constraint (c) With mass constraintFigure 4.14: Isogeometric shape optimization results

(b) With mass constraint (c) With mass constraintFigure 4.15: Isogeometric shape optimization results

(b) With mass constraint (c) With mass constraintFigure 4.18: Isogeometric shape optimization results

Acknowledgement

From the construction of NURBS basis functions, the following essential properties can be obtained:

(1) a partition of unity ∀ξ, n i=1 R i,p (ξ) = 1 (2) non-negative ∀ξ, 0 ≤ R i,p (ξ)

(3) kronecker's delta R i,p (ξ j ) = δ ij which is satisfied at patch boundary.

(4) linear independence

(5) variation diminishing property, leading to NURBS curves behaving monotonous with increased polynomial orders of basis functions. [START_REF] Dowling | Mean stress effects in stress-life and strain-life fatigue[END_REF] basis function continuity. It can be seen that NURBS basis functions are the rationale of B-spline basis functions, which result in sharing some same properties with B-spline basis functions like basis functions continuity. At the same time, because of the combination of the weight values, the NURBS basis functions are more flexible than the B-spline basis functions. [START_REF] Hughes | Isogeometric analysis: Cad, finite elements, nurbs, exact geometry and mesh refinement[END_REF][START_REF] Gautam | Iga: A simplified introduction and implementation details for finite element users[END_REF].

2) NURBS curve NURBS curve is defined by the linear combination of univariate NURBS basis function R i,p (ξ) and control point P i by the following expression [START_REF] Hughes | Isogeometric analysis: Cad, finite elements, nurbs, exact geometry and mesh refinement[END_REF]:

where the NURBS basis function R i,p (ξ) is defined by the Eq. (2.35), P i ={P x , P y , P z } is the ith control point vector.

Based on the Eq. (2.40), a NURBS curve can be constructed in Matlab, as shown in Fig. 2.11, in which the polynomial order, number of control points are respectively 2 and 4, and the corresponding knot vector and weight values are ξ = {0, 0, 0, 0.5, 1, 1, 1} and ω = {1, 1, 1, 1}.

• I) Convergence analysis

The maximum stress in the y-direction obtained from the static analysis is chosen to decide the convergence rate.

From Figure 2.24, it can be seen that with the increase of control points number of IGA and element nodes number of FEA, the maximum stresses tend to reach a stable value. In order to show the isogeometric convergence analysis intuitively, a fitting curve is used to approximate the obtained maximum stresses from different analyses. Furthermore, the minimum correlation coefficient r, which is defined by Eq. (2.52) between the maximum stress values and the values obtained from the fitting curve, is 0.92.

where Y i and Y f it,i are the ith values obtained from the numerical analysis and fitting curves, Y and Y f it are corresponding mean values, n is the extracted number sets and in this case n is 7. In this work, the convergence point is decided from the onset of the constant fitting curve. Thus, for IGA, the maximum stress values start to converge from the mesh density of 933 (the number of control points); for FEA, from the mesh density of 7701 (the number of element nodes). After the analysis results of these convergence points are obtained to compare the differences between IGA and FEA. On the other hand, we can observe some decreased stresses from isogeometric convergence analysis such as the stresses under mesh density of 2000 and 2200 (the final point from Figure 2.24 (a)) This is mainly because of the decreased applied stress on control points (in Ls DYNA, the applied load needs to be divided by the number of control points to obtain the load for each control point in this case). Under the mesh density of 2200, the control point load is smaller than the one under the mesh density of 2000, leading to more minor effects on the corner of the edge curve of the geometry and thus obtaining decreased stress. But, when the mesh density is large or enough, like in FEA, the load effects from element nodes on the corner may be minor. That's why in FEA, we can obtain a more smooth curve.

Chapter 3

Random Vibration Fatigue Analysis on Mechanical Structures

Introduction

This chapter develops the IGA and FEA on a specimen with a reduced section and an L-shaped plate model. The own-developed Matlab programming validates the fatigue analysis results.

The static convergence analysis is fulfilled for the specimen to compare the IGA and FEA differences. Then, a random acceleration load is applied to the clamping fixture to assess the random vibration fatigue life, based on the fatigue analysis methods of Steinberg's three bands, Dirlik, Narrowband, Wirsching, Chaudhury and Dover, Tunna,

Hancock. And the same excitation is considered in experimental tests to check the numerical results.

For the L-shaped plate model, the modal analysis is developed to compare the analysis results from the reference [START_REF] Pitoiset | Spectral methods for multiaxial random fatigue analysis of metallic structures[END_REF].

A random acceleration in the vertical direction to the plate surface is applied to the clamping fixtures to develop random vibration fatigue analysis. The Matlab program verifies the results.

Summary

In this section, isogeometric and finite element convergence and random vibration fatigue analyses were developed on the plate with notches model, in which the fatigue analysis results were validated by Matlab programming and experimental test.

In convergence analysis, the plate was clamped from the bottom side (see Fig. 3.3), and a 100 pa load was applied on top row elements to simulate the static bending problem. The results show that the IGA and FEA convergence analyses predict similar maximum stress values with the relative error of 0.75%, which are situated in similar locations; the CPU time of the IGA and FEA is respectively 231 and 462 seconds. It can be seen that IGA is more time-efficient compared with FEA.

The modal analysis results show that the obtained natural frequencies and vibration mode from IGA and FEA match well, in which the relative error of the first natural frequency from IGA and FEA is 0.9%.

In fatigue analysis, a random vibration acceleration in the vertical direction to the plate surface was applied to the clamping fixture to simulate the base acceleration. There were three kinds of comparisons in damage and fatigue life.

Firstly the simulation results show that the obtained maximum effective stress RMS and cumulative damage ratio from IGA and FEA have a good agreement, with the relative errors of -0.72%, 1.2% respectively. Secondly, based on Matlab programming, the predicted isogeometric and finite element damage ratios are 1.435 and 1.481, close to the damage values calculated from Ls Dyna, 1.449 and 1.467. Thirdly, the average fatigue life obtained from the experimental test is 3:30, which agrees very well with the isogeometric and finite element fatigue life calculated from Ls Dyna under σ f = 612.3 MPa, 3:27:02 and 3:24:30 (hours: minutes: seconds).

L-shaped specimen 3.3.1 The analysis preparation

Isogeometric and finite element convergence and random vibration fatigue analyses are developed on an L-shaped plate model [START_REF] Pitoiset | Spectral methods for multiaxial random fatigue analysis of metallic structures[END_REF] to compare the IGA and FEA differences in terms of mesh generation time, modal analysis, and fatigue analysis results.

• I) Geometric model

The model is featured by a hole, two notches as shown in Fig. 3.18, in which the thickness is 0.5 mm. The first five natural frequencies obtained from the X.Pitoiset et al. [START_REF] Pitoiset | Spectral methods for multiaxial random fatigue analysis of metallic structures[END_REF] are presented in Tab. 3.7. For IGA, the first natural frequency starts to converge from the mesh density of 738 (the number of control points);

for FEA, from the mesh density of 704 (the number of element nodes). Then the natural frequencies and vibration modes deciding from convergence points are obtained to compare the differences between IGA and FEA.

Conclusion

This chapter considered the static, convergence, and random vibration fatigue analysis on the two different models based on IGA and FEA. For the specimen with reduced section model, an experimental test and Matlab program are developed to check the fatigue analysis results. For the L-shaped plate, the isogeometric and finite element damage results are validated by the Matlab program.

For the plate model with reduced section, the model is clamped from one side during simulation, and random acceleration in a vertical direction to the plate surface is applied to the clamping fixture. Numerical simulations are verified by convergence analysis on the static bending problem. It is shown that the IGA and FEA convergence analyses lead to similar maximum stress values with a relative error of 0.75%, situating at similar locations of the reduced section.

The CPU time of the analysis, respectively 231 and 462 seconds, leading to the time difference of 50%, show that IGA is more time-efficient than FEA. The obtained first five natural frequencies and vibration modes have a good agreement.

Fatigue analyses show that the obtained isogeometric and finite element maximum effective stress RMS are 5.24e+7 and 5.278e+7 pa with a relative error of -0.72%, and the cumulative damage ratios are 1.449 and 1.467 with a relative error of 1.2%. Based on the Matlab program, the isogeometric and finite element damage are respectively 1.435 and 1.481. Thus, the damage comparison between IGA and FEA, numerical simulation, and Matlab programming shows that the obtained values agree very well. On the other hand, the fatigue life from the experimental test, 3:30:00 (hours: minutes: seconds), has a good agreement with the isogeometric and finite element fatigue life from LS Dyna, 3:27:02 and 3:24:30 under S f = 612.3 Mpa.

In the case of complex structures, the L-shaped plate model is clamped on two edges, and random acceleration in a vertical direction to the plate surface is applied to the clamping fixture. Based on obtained first natural frequencies, the convergence points are decided, and obtained frequencies are compared with those from Pitoiset. It can be observed that the IGA and FEA convergence analyses predict similar results, the first natural frequencies from IGA and FEA are respectively 28.96 and 28.63 Hz, leading to a relative error of 1.1%, which agree well with the first natural frequency from Pitoiset of 27.7 Hz. Fatigue analyses show that the obtained isogeometric and finite element maximum effective stress RMS are 1.919e+8 and 1.921e+8 pa with a relative error of -1.03%, and cumulative damage ratios are 7.571e-3 and 7.314e-3 with a relative error of 3.1%. Based on the Matlab program, the isogeometric and finite element damage are 7.6e-3 and 7.4e-3, leading to the relative error of 2.7%. Thus, the damage comparison between IGA and FEA, numerical simulation, and Matlab programming shows that the obtained values agree very well.

On the other hand, from the convergence analysis of these models, the numerical model development highlighted that for IGA during each mesh refinement step, it is not necessary to create mesh elements on the original geometry model. It is sufficient to develop mesh elements on the previous mesh model so that the mesh refinement time can be essentially saved. However, for the FEA, it is mandatory to communicate with the original geometric model for mesh refinement, and so this process was more time-consuming in LS Dyna software.

Chapter 4

Wind Turbine Tower optimization

Introduction

The aim of this chapter is to demonstrate the possibility of the developed modelling to increase the fatigue life of the realistic industrial structures through the shape optimization. The isogeometric random vibration fatigue analysis and design optimization are developed on a wind turbine tower model. The analysis results are verified by the FEA and Matlab programming. As shown in Fig. 4.1, this chapter mainly consists of three sections. In Section 4.2, a random force load is applied on the top concentrated mass node of the tower model to simulate the random vibration fatigue analysis. In Section 4.3, the isogeometric size optimization is considered with the design variables of the tower segment thicknesses, and the design objective is the direct maximization of the first natural frequency. In Section 4.4, the isogeometric shape optimization is fulfilled with the design variables of the control point position, in which the design objectives are respectively the direct maximization of the first natural frequency and direct minimization of the effective stress RMS. To obtain a regular door shape, in this case, only two control points are fixed. On the other hand, to furtherly investigate the different door shapes on the effects of effective stress, the shape optimization is also developed under the consideration of four fixed control points. During analysis, the tower was clamped on the base flange. From modal analysis results, it can be found that the isogeometric and finite element first natural frequencies are the same.

In fatigue analysis, longitudinal, lateral, and vertical random vibration force PSDs in the x, y, z directions to the tower model were applied on the concentrated mass node. The results show that firstly, the obtained maximum effective stress RMS and cumulative damage ratios have a good agreement, with the relative errors of -2.45%, and 1.61% respectively;

secondly, based on the Matlab program, the predicted isogeometric and finite element damage ratios are respectively 4.2e-5, which are close to the damage values calculated from Ls Dyna, 4.178e-5 and 4.112e-5. The damage is located near the door, thus providing the directions for the following shape optimization.

Sizing optimization

The mathematical formulation of a structural design optimization can be stated as [START_REF] Fiacco | Nonlinear programming: Sequential unconstrained minimization techniques[END_REF]:

subject to g(x) = 0 and h(X) ≥ 0

where x = (x 1 , x 2 , . . . , x n ) is the design variable with a set of components, f (x) is the objective function, g(x) and h(x) are the constraint functions.

Analysis preparation • I) Basic assumptions

The tower model assembled by a series of the cylinder and conical shell segments have different thicknesses. To reduce the dimensionality in the optimization problem, the following variables will be considered fixed:

1) tubular tower configuration 2) the height and diameter of each segment 3) material properties, concentrated mass, etc.

4) some analysis parameters like shell element formulation, integration points, boundary condition, etc.

• II) Design objectives

Several tower design objectives existed, including lightweight design, high stiffness, high ( stiffness/ mass )-ratio, etc.

In this studying, we choose the direct maximization of the system's first natural frequency as the tower design objective, which is beneficial for the long fatigue life, high stability, low noise, etc.

• III) Design variables

The design variables are selected to be the thickness of the segments. For the tower model, there are four kinds of segment thicknesses (t i , i = 1, 2, 3, 4), that are respectively 16, 14, 12, 10 mm, in which we only considered the effects of thicknesses 16 and 14 mm on the first natural frequencies.

• IV) Design constraints

The design constraints are in terms of mass and thickness.

1) Mass limitations

where M u and M l are respectively allowable upper and lower limiting value of the tower mass. For the size optimization, the values are respectively 6.8605E + 07 ≤ M ≤ 6.9475E + 07 (g)

2) Side constraints

where t iu and t il are respectively allowable upper and lower limiting value of the ith thickness.In this studying, the values are respectively 15 ≤ t 1 ≤ 17 and 13 ≤ t 2 ≤ 15 mm.

• V) Optimization procedures

There are three main steps for wind turbine size optimization, as shown in Fig.

The first shape optimization 4.4.1 Analysis preparation

The purpose is to find the optimal shape of the door, which will reduce the RMS stress and thus increase the service life of the wind tower. In IGA, the wind turbine door shape is controlled by the 10 NURBS curves (from 5 to 14), as shown in the Fig. 4.13 (a). To carry out the shape optimization, it is necessary to ensure that the control points connecting different NURBS curves should be the same. For example, for the NURBS curve 11 and 12, the curves are connected We can observe some asymmetric and large area stress distribution from the analysis results on the plate model (like Fig. 3.5), which is mainly because of the less mesh density in critical areas. We will try to use the new basis functions such as T spline to support the local mesh refinement; 3) design optimization on the whole structure. In this studying, the shape optimization on the wind turbine tower model is developed only considering the control point position of the door shape. Next step, we will perform isogeometric shape optimization on the whole tower model under consideration of multi-objective and advanced optimization algorithms; 4) advance fatigue criteria.

For the fatigue analysis, the Sine's or Crossland's fatigue criterion will be considered; 5) new material. In the future, we will consider IGA on new materials, especially on 3d printing materials.

Appendix A

IGA Formulation in vibration

Let's consider a body and the global Cartesian reference system x = (x, y, z) T and its associated orthonormal basis e = (e 1 , e 2 , e 3 ) T . The body is subjected to the volume forces f = f i e i and traction vector acting on the body surface t = t i e i . The domain of the body Ω is bounded by prescribed displacement Γ u and traction Γ t boundaries. We can introduce the displacement field as u = u i e i and the stress σ = σ ij e i ⊗ e j (i, j = 1...3) satisfying the equation in strong form:

where ∇ represents the divergence operator, ρ is the density assumed constant and n is the outward normal vector.

∂t 2 = ü is the acceleration vector of the body. This term stands for the inertial effect, which is considered negligible for a static problem.

The displacement field is approximated using NURBS basis functions to construct the CAD geometry of the structure.

The discretization of the domain Ω into a number of sub-domains

] is carried out by using the B-spline and NURBS formulations described in Section 2.1, which transform the parametric coordinates to physical coordinates. A mapping is introduced to perform the analysis on each control point cp of any element e:

x e (η) =

R e cp (η)P e cp (A.3) where η contains the parametric coordinates given by ξ in one dimension and (ξ, λ) in two dimensions. n e cp is the number of control points over element e. Using the Galerkin method, the displacement and virtual displacement fields can be deduced as follows:

R e cp (η)u e cp (A.4)

R e cp (η)δu e cp (A.5) where u cp and δu cp correspond to the values of the displacement and virtual displacement fields at the control point P cp . The momentum equation of Eq.(A.1) can be used in a weak form based on the virtual work principle which can be built for each element e to obtain the governing dynamic equilibrium of motion for a structure:

u e is the control point displacement and F e is external force vector that represents the contribution of the total body forces. M e , C e and K e are, respectively, the isogeometric element mass, damping using the damping property of an element e denoted by κ, and stiffness matrices given by:

)

)

B eT C B e dΩ e (A.9)

where C = C ijkl e i ⊗ e j ⊗ e k ⊗ e l and ε = ε ij e i ⊗ e j are, respectively, the fourth order material elastic and the strain tensors. The components of strain tensor can also be formulated from displacement as follows:

R is the matrix of NURBS basis function and B is the strain displacement matrix. It can be obtained by computing the derivatives of the basis functions R cp (η) for any element e. In bidimensional problem, the matrix has the following form:

The global matrices can be assembled by their corresponding elemental matrices as follows:

IGA in Ls Dyna

In this part, I will present the isogeometric static and dynamic analysis processes in Ls Dyna. Appendix C

Wind load calculation

The wind load applied in this work is calculated based on the reference [START_REF] Nguyen | Réponse des lignes de transport au vent turbulent par la méthode spectrale stochastique[END_REF]. A wind turbine is subjected to various loads from different sources such as aerodynamic loads. These loads depend on the environment of wind turbine systems. Indeed, there are many different climatic regions depending on the geography (mountain/plains, land/sea, type of vegetation) leading to a large variability of climatic conditions. The cyclic part of the load caused by wind shear is important to consider because it can lead to fatigue damage. The wind velocity depends on the time and the position of the calculation point. In the Cartesian space with parallel to the wind direction, the wind velocity V ent(t, P ) at the time t, point P decomposed by the longitudinal component -→ U (t), the lateral component -→ v (t) and the vertical component -→ w (t) is given by the equation following:

Where U is the longitudinal mean velocity which can be calculated, using the Eq. (C.2) or (C.3) for a height z denoted by U (z), and the corresponding wind-induced force F (z) can be calculated with the method of static equivalent, as shown in Eq. (C.4).

Mean velocity U (z) is extrapolated (in the hypothesis of atmosphere neutral stability) from measurements at lower heights to a height z according to following logarithmic law:

where z 0 is the roughness parameter (here z 0 = 0.05 m), U (z r ) is the average longitudinal wind speed at the reference height above the ground, normally given as z r = 10 m. The following power law can also be used after identification of α:

Where α is the power exponent depending on the surface roughness.

Based on the Eq. (C.2), the wind velocity at different height can be obtained, as shown in The wind-induced average load F (z) is given by:

Where ρ is the air density, C d is the coefficient of the wind resistance, S is the cross-sectional area of a structure perpendicular to the wind direction. u(t), v(t) and w(t) are considered as zero-mean Gaussian stationary processes. These processes are characterized by their standard deviations σ u , σ v and σ w , as shown in the equations following. And the effects of wind turbulence {u(t), v(t), w(t)} T are described by the corresponding wind speed PSD S v (f ) and wind force PSD S F (f ).

Where k 1 is the turbulence coefficient defined by the Eq. (C.9), k r is the field factor, determined by the Eq. (C.8).

K r = 0.19( z 0 z 0,11

) 0.07 (C.8)

The turbulent intensities are given by For Eurocode 1, for heights z less than 200 m, the expression is shown as:

where β = 0.67 + 0.05ln(z 0 ). Approximation of these parameters are often used such as L y u = 0.3L x u , L z u = 0.5L x u . Under the assumption of isotropic turbulence, the following Eq. can be given:

112

The longitudinal component of the wind is often decomposed into 3 zones. The low frequencies related to a large anisotropic vortex and leading to the energy generation. The mean frequencies leading to the energy transferred to small vortices. The high frequencies where the turbulent kinetic energy gets converted into heat by viscous dissipation due to viscous shear stresses.

The power spectral density (PSD) represents the energy distribution of the turbulence in different frequency ranges.

Numerous studies adopt the Kaimal spectrum or Von Karman spectrum recommended in various standards. The normalized wind PSDs f Su(n) σ 2 u are expressed by the normalized frequency which can simulate wind with or without turbulence. The wind applied in this work is considered with turbulence, the normalized frequency given by

include the turbulence scales.

For 0.007 Hz ≤ f ≤ 1 Hz, the unilateral von Karman's PSD of the longitudinal velocity is given by:

(C.17)

For the mean and high frequencies, the following Kaimal PSD is used:

With L x u,k = 2.329L x u . For the Davenport:

Where L x u,d = 1200. For the Eurocode 1: U (zr) 2 can be respectively given by: For the Von Karmen:

For the Kaimal:

Based on the equations above, the normalized wind velocity PSD in x direction can be obtained, as shown in the where C d is the dynamic coefficient.

The PSD of the wind-induced turbulent force at a point P can be described by: S F (P, f ) = (ρC d S) 2 U (z) 2 S u (P, f ) (C.26)

Where S u (P, f ) is the wind speed PSD at the point P .

Meanwhile, the normalized wind force PSD f S F (P,f ) U (zr) 2 can be given as: