Trade off Between Accuracy and Efficiency

CHAPTER SIX

Trade off Between Accuracy and Efficiency 6.1 DIGITAL REPRESENTATION OF MICROSTRUCTURE The concept of the Digital Material Representation (DMR) was proposed recently and is constantly evolving. The definition according to [325] states that DMR is a material description based on measurable quantities that provides the necessary link between simulation and experiment. The DMR is expected to create a possibility for analyzing material behavior in the conditions, which are difficult or even impossible to be monitored experimentally at the present state of equipment. The main objective of the DMR is creation of the digital representation of microstructure with its features represented explicitly to match real microstructure morphology. The DMR is used for calculations related to a unit cell or Representative Volume Element (RVE), depending on what kind of information is required either local or global. The RVE is a model of the material to be used to determine the corresponding effective properties for the homogenized macroscopic model. The RVE should be large enough to contain sufficient information about the microstructure in order to be representative, however it should be much smaller than the macroscopic body [107,133]. The unit cell is a part of RVE that enables obtaining results for particular part of the material. Thus the unit cell is not representative for the whole numerical model. Application of the unit cell idea, enables analyzing material behavior in particular location, for example, crack initiation along the inclusion, without focusing on the rest of the material. As a result the unit cell provides data only for local analysis, not for whole material behavior (unlike the RVE). Usually several unit cells can be considered as the RVE. It can be summarized that the DMR approach makes it possible to obtain a completely virtual description of the deformation process and provides better quality of numerical simulations to material scientists. Additionally this leads to a reduction in research costs and simplifies design of new materials. Computational Materials Engineering DOI: http://dx.doi.org/10.1016/B978-0-12-416707-0.00006-5

© 2015 Elsevier Inc. All rights reserved.

255

256

Computational Materials Engineering

The DMR concept of polycrystalline metallic materials based on the most commonly used approaches including the earlier described cellular automata method is presented in detail in the following chapters.

6.1.1 Voronoi The Voronoi tessellation is commonly used method dedicated to generation of initial digital microstructures, as well as to interpretation of some metallographic phenomena. The idea of this method is based on the mapping of the bounded area onto the group of specific polygons P 5 {p1,. . ., pn}, generated around a set of n distinct initial points S 5 {s1,. . .,sn}. Each polygon px is characterized by the two following features: • It is connected exactly to one point sx, where x 5 1. . .n. • Each point inside the polygon is closer to sx than to any other point from S. For the purposes of the present work points from the S set represent grains nuclei, while the polygon areas around these points, called Voronoi cellulars, represent final grains in the microstructure. Several algorithms used to determine Voronoi diagram can be distinguished. One of them, based on the Brutal Forces approach, points out bisectrixes between each pair of points and then eliminates bisectrixes, which do not satisfy Voronoi assumptions [67]. Nevertheless, this algorithm is very time consuming. Thus, the Fortune algorithm [67], which is much more efficient, however, is limited only to 2D applications is often used in practical applications. The solution to create Voronoi diagrams in 2D as well as in 3D space is a Delaunay triangulation. The general algorithm for triangulation starts by forming the super triangle enclosing all the points from set V that have to be triangulated. Then incrementally, a process of inserting the points p into the set V is performed. After every insertion step a search is made, to find the triangles whose circumcircles enclose p. Identified triangles are then deleted from the set. As a result an insertion polygon containing p is created. Edges between the vertices of the insertion polygon and p are inserted and form the new triangulation. After all the points are inserted, the Delaunay triangulation is created (Figure 6.1a).When all the centers of the circumcircles are connected, the Voronoi diagram is obtained (Figure 6.2b). Similar procedure is applied in 3D space. Examples of obtained results for different number of grains in 2D space are presented in Figure 6.2 while for 3D in Figure 6.3.

Trade off Between Accuracy and Efficiency

257

Figure 6.1 (a) Schematic idea of the Delaunay triangulation, and (b) Voronoi diagram.

Figure 6.2 Results of simulation for different number of grains (a) 10, and (b) 100 obtained from the Delaunay triangulation in 2D space.

As seen, this approach replicates polycrystalline character of metallic microstructure, however, obtained shapes of grains are still far from shapes of real grains observed under the optical microscopes. This is the main disadvantage of this approach. That is why authors decided to implement a grain growth algorithm based on the earlier described cellular automata method.

6.1.2 Cellular automata grain growth model The basic definitions of the cellular automata method that are necessary to create a robust numerical model are presented in Chapter 5. This method was adopted to develop a model of the grain growth in order to

258

Computational Materials Engineering

Figure 6.3 The digital microstructure containing 100 grains obtained from the Delaunay triangulation in 3D space.

create digital representation of microstructure morphology. The implemented algorithm for the purposes of this research is a simplified version of the CA grain growth algorithm. The first step in the CA method is to establish the discrete space composed of cellular automatons. As mentioned in the 2D space it will be a grid built of squares, whereas in 3D spaces this grid will be composed of cubes. In the next step of the algorithm a set of CA cells is selected randomly and then internal variable describing cells state is set to already grown. These cells represent grains nuclei. The second step of the algorithm is focused on grain growth. The transition rule for this stage is defined as follows: when a neighbor of a particular cell in the previous time step is in the already grown state, then this particular cell can also change its state. Particular grains grow with no restrictions until the impingement with other grains. After that they grow only in the empty space area where no grains are observed as seen in Figure 6.4. This process is performed until the entire space is filled with grains. Due to the fact that transition rules depend on the type of the neighborhood, several different neighborhoods presented in Chapter 5 are implemented to investigate differences between them. Examples of the obtained 2D and 3D microstructures using various CA neighborhoods are presented in Figure 6.5. The CA space size used in the 2D simulation was

259

Trade off Between Accuracy and Efficiency

(a)

(b)

(c)

Growing grains

Grains nuclei

Empty space

Figure 6.4 Stages of CA grain growth algorithm.

Figure 6.5 Microstructures obtained by (a) von Neumann, (b) Moore, (c) hexagonal random, and (d) random for 40 grains.

set to 500 3 500 px which is equal to 500 3 500 μm. In the 3D, the space size was set to 300 3 300 3 300 μm. Based on the space and required grain size the number of initial grain nuclei can be easily calculated. Due to their stochastic character, results obtained by the hexagonal random and random neighborhoods can represent the geometry of a real microstructure after annealing or static/metadynamic recrystallization.

260

Computational Materials Engineering

Figure 6.6 Microstructure with 10 grains generated (a) without and (b) with periodic boundary conditions.

That is why these two neighborhoods are recommended for generation of initial microstructures for further simulation of material behavior under loading conditions. Additionally, to obtain CA space/microstructure continuity, periodic boundary conditions are introduced. Using periodic boundary conditions no free edges in the CA space are present what assures space continuity. Microstructures generated with and without periodic boundary conditions are shown in Figure 6.6. Initial position of grain nuclei and type of neighborhood was exactly the same in both investigated cases. Similar approaches to generate microstructure and boundary conditions were implemented in the DMR model in 3D space. The equivalent of the hexagonal or pentagonal random neighborhoods is especially interesting in the 3D space because there are six possibilities in neighbor selection. Examples of obtained microstructures on the basis of two different types of neighborhoods are shown in Figure 6.7. The grain size distribution of the obtained 2D and 3D microstructures can be easily compared with the experimental results. If the generation of the microstructure does not provide the required grain size distribution, the process is performed one more time. An example of the grain size distribution obtained from the CA model is shown in Figure 6.8. It can be concluded that the CA approach with the appropriate neighborhood can very closely replicate, in a statistical manner, the geometry of a real microstructure after annealing or static/metadynamic recrystallization. The digital microstructure is then considered during numerical simulations of loading conditions as a statistical representation of a real microstructure.

261

Trade off Between Accuracy and Efficiency

Figure 6.7 Results of microstructure obtained by (a) Von Neumann, (b) Moore in 3D space. 160 140

Grain count

120 100 80 60 40 20 0 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 Grain diameter, µm

Figure 6.8 Distribution of the grain sizes in the digital microstructure representation.

Grain boundary CA cell outside the radius

CA cell within the radious Center of the spherical inclusion

Grains

Figure 6.9 Algorithm used to incorporate inclusions into the DMR model.

The algorithm can be easily extended to introduce other microstructure features, for example, inclusions. The inclusions are added using a developed algorithm illustrated in Figure 6.9. All CA cells belonging to the grains or grain boundaries lying within the specified radius R changes their state to represent inclusion. Spheres are added manually into the

262

Computational Materials Engineering

Figure 6.10 Initial microstructure with added inclusions.

microstructure, through user interface, thus by their proper combination any shape of the inclusion can be obtained as seen in Figure 6.10. The size and locations of these inclusions can be evaluated on the basis of experimental knowledge. Another simple modification to the presented algorithm can provide also digital representation of the two phase microstructure. In order to obtain representation of, for example, dual phase microstructures, selected grains are randomly assigned to second phase (Figure 6.11). However, the more complex the microstructure is investigated, the finite element mesh density increases to accurately represent the morphology of particular constituents. Thus, one of the approaches to reduce this computational effort is to apply significantly simplified geometry of microstructure morphology while maintaining its representativeness properties. The so-called statistically similar representation volume element (SSRVE) approach is presented next. The basic differences between both approaches are presented in Figure 6.12.

6.2 REDUCTION OF THE COMPUTATIONAL DOMAIN As it has been shown in the previous chapters, more advanced analysis of materials uses multiscale modeling techniques, where numerical

Figure 6.11 Examples of 2D digital microstructures representing dual phase materials obtained by modified cellular automata algorithm (a) separated grains, (b) binary form of microstructure, (c) generated specific finite element mesh.

Macro

Or

Micro

SS-RVE

DMR-RVE

Figure 6.12 Basic differences between DMR based on RVE and DMR based on SSRVE approaches.

264

Computational Materials Engineering

simulations are performed in many scales. FE model or alternative methods are used in the macroscale while microscale is modeled mainly by using FEM (FE2) or Cellular Automata (CAFE) [215] or other discrete methods. However, the extremely high computing costs of multiscale calculations are the main limiting factors for a wide application of that approach. As mentioned (Figure 6.12) application of the statistical representation of the microstructure is one of the possibilities of reduction of complexity of material representation in multiscale calculations. Introduction of the Statistically Similar or Statistically Equivalent Representative Volume Element (SSRVE/SERVE) [319,296,17] is one of the possibilities of simplifying the digital representation of the multiphase materials. It is based on an assumption that typical RVE can be further reduced into simplified form, which joined together periodically behaves in the same way as its larger equivalent. The general idea of the creation of the SSRVE is discussed in this section of the book. In general, the SSRVE is based on Non-Uniform Relational B-Splines representation and is determined by using optimization procedure, where objective function includes comparison of mechanical properties, shape coefficients, and statistical characteristics of real microstructure. The first of these elements requires application of finite element method allowing to simulate deformation of initial RVE and current SSRVE within elastic or elastic-plastic range. That requires a series of calculations and is time consuming. Therefore, an alternative approach allowing to replace time consuming FE solution with more efficient isogeometric analysis (IGA) was considered within the book as an example of possibilities provided by alternative computational approaches. This is of importance to have a broad knowledge on different available solutions, for example, boundary element method and meshfree methods, that may be more suited for particular problem. Appropriate selection of computational technique is also a way to deal with long computing times. The performance of IGA approach was analyzed and compared to the conventional FE-based methodology within the chapter. Additionally, possibilities of IGA implementation on heterogeneous hardware devices allowing to improve computational efficiency and decrease overall power consumption are also mentioned.

6.2.1 Statistical representation of microstructure 6.2.1.1 Idea of SSRVE As mentioned in Chapter 5, in the micro
macro modeling approach an RVE representing the underlying microstructure is attached at each Gauss

Trade off Between Accuracy and Efficiency

265

point of the macroscopic solution. The constitutive law describing material behavior in the macroscale is obtained by averaging the first Piola
Kirchoff stresses with respect to the RVE. The theoretical basis of the micro
macro modeling is well described in the scientific literature [319] and it is briefly repeated here. The focus is on the development of the simple SSRVE, which allow decreasing the computing costs and make micro
macro modeling approach more efficient. Considering micro-heterogeneous materials, the continuum mechanical properties at the macroscale are characterized by the morphology and by the properties of the particular constituents at the microscale. In this book the idea of the SSRVE is presented using two-phase materials. Since construction of the SSRVE requires that the second phase in the material occurs in the form of islands, which are separated from one another, DP steel microstructures composed of soft ferrite (70
80%), and hard martensite (20
30%) are considered. The theoretical basis of the method and an algorithm of the solution are presented for the DP steel, following the research described in [296]. In the analysis for DP steels only measures characterizing the hard martensite islands are considered. The material models of the individual constituents are assumed to be known. The description of the microstructure is based on statistical consideration as presented in the following section [33]. Description of a usual RVE is determined by the smallest possible sub domain, which is still able to represent the macroscopic behavior of the material. Although these RVEs are the smallest possible by definition, they still can be too complex for efficient calculations. Therefore, the construction of statistically similar RVEs, which are characterized by a lower complexity than the smallest possible substructure, was proposed in [319]. The basic idea is to replace RVE with arbitrary complex inclusion morphology by a periodic one composed of optimal unit cells, see Figure 6.13. This idea is applied in this book to the analysis of the DP steel microstructures. 6.2.1.2 Shape coefficients The parameters, describing fraction of different phases and their geometrical characteristics, form a group of the most important parameters, which, beyond statistical analysis of microstructure, is often applied. Ohser and Muecklich [251] proposed four basic parameters. Few more parameters, which are used in the image analysis, were added and described in the present book, as well. In consequence, the following thirteen parameters are considered.

266

Computational Materials Engineering

Figure 6.13 Schematic illustration of the basic concept of the SSRVE [43]—(a) RVE, and (b) periodically arranged SSRVE.

Volume fraction of hard constituent: ζ1 5

Vm V

(6.1)

where Vm is the total volume of hard constituent and V is the volume of the sample. The interfacial surface to volume ratio for hard constituent islands: ζ2 5

Sm Vm

(6.2)

where Sm is the interfacial surface. Another coefficient represents relation between the interfacial surface and the volume of hard constituent islands, but it is given in a dimensionless form. This coefficient is represented by dimensionless parameter defined as relative difference between pffiffiffiffiffiffiffiffiffiffiffiffiffi radii (RS
RV)/(RV) calculated from the surface Sm (RS 5 Sm =4π) and from the volume Vm pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi (RV 5 3 3Vm =4π): pffiffiffiffiffiffi Sm ffiffiffip ffiffiffiffiffiffip ζ3 5 p (6.3) ffiffiffiffiffiffiffi 3 6 3 4π 3 Vm 2 1 Roundness of the hard constituents is the next basic parameter, which characterizes the DP steel microstructure. Roundness is defined as the difference between volumes of the smallest enclosing ellipsoid and the largest enclosed ellipsoid: 4 ζ 4 5 πðABC 2 abcÞ 3

(6.4)

267

Trade off Between Accuracy and Efficiency

where A, B, C are the equatorial radii of the smallest enclosing ellipsoid and a, b, c are the equatorial radii of the largest enclosed ellipsoid. Ellipsoid fit is the degree of fitting to ellipsoid with the same area and length-width-height proportions: ζ5 5 2

ðx; y; zÞAvm :εm ðx; y; zÞ # 1 21 Vm

(6.5)

where vm is a set of points inside the hard constituent inclusion and εm is a set of points inside the ellipsoid. Ellipsoid diagonals are equal to Feret diameters calculated for particular grains transformed on the basis of Principal Component Analysis (PCA). Graphical interpretation of both coefficients describing roundness, for example, ζ 4 and ζ 5, is presented in Figure 6.14. Another parameter is the ratio between the minimum (r min) and maximum (r max) distance between contour and the center of gravity of an island: rmin ζ6 5 (6.6) rmax The following parameter is the border index, which is the ratio between the surface of the martensite constituent and the maximum surface of this constituent: ζ7 5

SLm SLmmax

(6.7)

where SLm and SLmmax are projections of the inclusion. Thus, SLm is the perimeter of the projection of the hard constituent and SLmmax is calculated as the perimeter of the smallest rectangle enclosing this projection.

Smax Sm Smin

νm εm

Roundness

Ellipsoid fit

Figure 6.14 Schematic illustration of the parameters characterizing degree of roundness.

268

Computational Materials Engineering

The curvature of the surface of the martensite constituent is considered in the following two parameters. The specific integral of the mean curvature is: ! ð 1 ζ8 5 min½κ 1 min½κ dSm (6.8) 2Vm Sm β β where β is the direction in the tangential plane and κ 5 κ (β) is the curvature. The specific integral of the total curvature is: ! ð 1 ζ9 5 min½κ min½κ dSm (6.9) Vm Sm β β The last two parameters provide statistical information concerning the degree of the curvature of the hard constituent inclusion. Additional four parameters are dedicated only for analysis of objects inside digital images and they have no equivalents in 3D case. Malinowska coefficient is defined as: Lm ζ 10 5 pffiffiffiffiffiffiffiffi 2 1 2 πSm

(6.10)

where Sm, Lm are area and perimeter of an analyzed island, respectively. Blair
Bliss coefficient is defined as: Sm ffi ζ 11 5 rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi (6.11) n P 2 2π ri i50

where n is the number of all pixels and ri is the distance between selected pixel and object’s center of gravity. Danielsson coefficient is calculated from the following equation: ζ 12 5 

Sm3 2 n P li

(6.12)

i50

where li is a minimal distance between a pixel and object’s contour. Finally, Haralick coefficient is given by a formula: vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u 2 n u P u di u i50 u (6.13) ζ 13 5 u P t n 2 n di 2 1 i50

269

Trade off Between Accuracy and Efficiency

Lm r1 sm

Malinowska

sm r2

Blair–Bliss

l1

sm

l2

Danielsson

d1

d2

Haralick

Figure 6.15 Schematic illustration of all the parameters dedicated to image analysis.

where n is the number of pixels belonging to the object’s contour and di is a distance between an analyzed pixel and object’s center of gravity. The graphical interpretation of some of the mentioned shape coefficients is presented in Figure 6.15. All these basic geometrical parameters were analyzed and evaluated in the work [296] and those having the strongest influence on the features of the microstructure were selected as presented in the following part of the chapter. Beyond this, the statistical measures of higher order are discussed briefly below, although they are of minor importance for typical martensite-ferrite microstructures of DP steels. 6.2.1.3 Statistical measures of higher order The statistical measures, which are used for microstructure characterization in the SSRVE for micro
macro modeling of two-phase materials are discussed in this section. Following [319], these measures are divided into n-point probability functions, spectral density and linear-path function. A number of measures proposed in the publication [319], as well as those used in the image analysis technique, are listed below and those used in the present book are explained wider. Since the parameters discussed in previous section are scalar values, they are not able to cover direction-dependent information, which is required for the representation of oriented inclusion leading to an overall anisotropy. Statistical measures of higher order are needed for the construction of the SSRVE for anisotropic microstructures. Schroeder et al. [319] introduced the following higher order statistical measures for microstructures: n-point probability functions, spectral density, and linear-path function. The latter parameter is described briefly below. This parameter describes the probability that a complete line segment a 5
a

1

a! 2 is located in the same phase, where a1 5 {x1,y1} and a2 5 {x2,y2} are coordinates of the ends of the line segment. Lu and Torquato, cited in [319], gave a general mathematical description of this

270

Computational Materials Engineering

measure for multiphase anisotropic materials. Simplified approach, which is applicable to DP steel microstructures, is presented below. Let D(i) denote the domain occupied by the selected phase i in the particular sample. The following modified indicator function is considered:  ðiÞ 1 if
a

1

a! ðiÞ




! 2 AD χ ð a1 a2 Þ 5 (6.14) 0 otherwise Function (6.14) checks if a whole line segment is located in the selected phase i. The general definition of a linear-path function requires computation of the ensemble average:




! ðiÞ




! ζ ðiÞ LP ð a1 a2 Þ 5 χ ð a1 a2 Þ

(6.15)

Procedure of numerical evaluation of ζ ðiÞ LP for a general microstructure is given in [319]. Since statistically homogeneous isotropic DP microstructures are analyzed as an example in the present book, several simplifications are introduced to this procedure. The origin of the line segment vanishes and the only dependency is on the length (b) and orientation of the line segment. (i) When b tends to zero, ζ ðiÞ LP tends to the volume fraction of the phase F . Further simplification is due to the fact that only martensite inclusions are analyzed in DP steels. Thus, in 2D binary images of statistically homogeneous microstructure with inclusions it takes the following form: 















! χðyÞ 5 1 if a1 ða1 1 bÞAD (6.16) 0 otherwise where b 5 {x2
x1, y2
y1} is the line segment vector. When evaluating the linear-path function for the periodic binary image of the two-phase microstructure a suitable template consisting of N 5 Nx 3 Ny pixels is defined. The pixel position in the binary image is defined by the pair (p,q). Length of the line segment vector in the binary image measured in pixels is b 5 {m,k}. For such a discrete case the linear-path function is: ζ lLP ðm; kÞ 5

pX M 21 qX K 21 1 χl ð~yÞ ðpM 2 pm ÞðqK 2 qk Þ p5pm q5qk

(6.17)

where pm 5 max½0; 2m; pM 5 min½Nx ; Nx 2m qk 5 max½0; 2k; qK 5 min½Ny ; Ny 2k

(6.18)

271

Trade off Between Accuracy and Efficiency

For periodic unit cells the linear-path function is: Nx X y 1 X χðbÞ Nx Ny p51 q51 N

ζ LP ðm; kÞ 5

(6.19)

Calculation of the linear-path function from Eqs.(6.17)
(6.19) is computationally expensive. The efficiency of the solution can be improved without noticeable loss of the accuracy. The following simplifications can be introduced: • The size of the template is reduced according to the maximum inclusion size. • Number of analyzed line segments is reduced and only certain set of line orientations is considered. • Not all positions of the original image are considered. Line segments are placed at random positions and this process is repeated Nr times. After these simplifications are introduced the linear-path function can be described as: ζ LP ðm; kÞ 5

Nr 1 X χðbj Þ Nr j51

(6.20)

The vector bj 5 {pj 1 m, qj 1 k}. The random positions of the coordinates have to fulfill following conditions: max½0; 2m # pj # min½Nx ; Nx 2 m 2 1 max½0; 2k # qj # min½Ny ; Ny 2 k 2 1

(6.21)

for random microstructures and: 0 # p j # Nx 2 1 0 # qj # Ny 2 1

(6.22)

for periodic unit cells. Presented statistical measures of higher order can be used to discover periodicity as well as different kinds of patterns in microstructure inclusions. 6.2.1.4 Sensitivity analysis of shape coefficients As mentioned in Chapter 2, sensitivity analysis is an effective tool to reduce the complexity of numerical approaches. Selection of the parameters, which are the most representative for the two-phase

272

Computational Materials Engineering

Table 6.1 The ranges of possible values of the shape coefficients Coefficient Symbol Units

Volume fraction Area/volume Area/volume Roundness Ellipsoid fit Contour to center Border index Mean curvature Total curvature Malinowska Blair
Bliss Danielsson Haralick

Steel A

ζ1 ζ2 ζ3 ζ4 ζ5 ζ6 ζ7 ζ8 ζ9 ζ 10 ζ 11 ζ 12 ζ 13

Range


1/m
m3


1/m2 1/m





Steel B

[0,1] [3/R,N] [0,N] [0,N] [21,1] [0,1] [1,N] [0,N] [0,N] [0,N] [0,1] [0,N] [0,1]

Steel C

Figure 6.16 Samples of microstructures obtained from DP steels with different phase fractions.

microstructure, is the prime task, which has to be performed to allow effective construction of the SSRVE. In the paper [296] this selection was made on the basis of the sensitivity analysis, which included evaluation of the importance of the parameters. This procedure is described below using DP steel microstructure as an example. The ranges of possible values of the shape coefficients for this steel are given in Table 6.1. Sensitivity analysis was performed on the basis of values of coefficients obtained for general set of dual phase steels characterized by different fractions of martensite. Typical examples of these microstructures are presented in Figure 6.16. Steel A is typical DP steel composed of 125 small separated islands without included matrix phase. The fraction of martensite phase in steel B exceeds 30%, thus it cannot be treated as typical DP steel,

273

Trade off Between Accuracy and Efficiency

Border index

Number of grains –1

Steel A Steel B Steel C 0

1 Coefficient value

2

3

Number of grains

Malinowska 40 35 30 25 20 15 10 5 0

40 35 30 25 20 15 10 5 0

Steel A Steel B Steel C 0

0.5

1

1.5

2

2.5

Coefficient value

Figure 6.17 Distributions of Malinowska (ζ 10) and border index (ζ 7) coefficients for selected microstructures characterized by singular value of the coefficient.

however, the shape coefficients can be still reliably calculated. In steel C the fraction of martensite was slightly higher than 50%. This caused very sophisticated shapes of analyzed hard constituent islands with many islands connected together, which highly influence reliability of coefficients based on perimeter-area ratio or on an object’s center of gravity. Analysis of calculated coefficients was performed twofold, that is, for one selected microstructure, and in comparison to other microstructures. In the first case, distributions of the calculated coefficients were estimated on the basis of obtained histograms. The analysis of distributions showed that some of the coefficients, for example, Malinowska or border index, are localized around one peak value. These peaks can be treated as expected values and they can be applied in optimization function as optimal values of particular coefficients (Figure 6.17). One-peak distributions are complemented by other types of plots, for example, distributions with more peaks, uniform distributions or monotonic distributions. The plot with more peaks would suggest that the coefficients are localized nearby two or more values, therefore SSRVE should be composed of two or more inclusions, characterized by particular expected values. The Haralick coefficient (ζ 13) is such an example in the case of steels A and B (Figure 6.18). It can be applied in optimization procedures for each inclusion of SSRVE separately. Besides strongly localized values, the coefficients, which are characterized by almost uniform distribution exist, for example, Blair
Bliss coefficient for steel C (Figure 6.19). Distributions of this coefficient are characterized by a specific range of values, where the coefficient has noticeable predominance. This information about uniform distribution can be passed to the cost function as a constraint for possible solutions.

274

Computational Materials Engineering

Haralick Number of grains

25 20 15 Steel A Steel B Steel C

10 5 0 0.85

0.9

0.95

1

Coefficient value

Figure 6.18 Distribution of the Haralick coefficient for selected microstructures characterized by double or triple peak values.

Blair-Bliss Number of grains

20 15 Steel A Steel B Steel C

10 5 0 0.5

0.6

0.7 0.8 0.9 Coefficient value

1

1.1

Figure 6.19 Distribution of the Blair
Bliss coefficient for selected microstructures.

Coefficients characterized by monotonic (ascending or descending) distributions are the least useful in the design of the SSRVE. Contour to center or ellipsoid fit coefficients for steel A or C, presented in Figure 6.20, are such examples. Independently of the number of inclusions inside SSRVE, these coefficients cannot be reliably incorporated into a cost function. Following the results of sensitivity analysis, the set of coefficients was chosen and applied in the cost function in the optimization procedure described in next section. 6.2.1.5 Construction of the SSRVE The process of SSRVE creation consists of the image analysis, calculation of shape coefficients, construction of a cost function for optimization

275

Trade off Between Accuracy and Efficiency

Contour to center

Ellipsoid fit 25

30 25 20 Steel A Steel B

15 10

Steel C

5 0 0

0.2

0.4

0.6

Coefficient value

0.8

1

Number of grains

Number of grains

35

20 15 10

Steel A Steel B

5

Steel C

0 0

0.2

0.4

0.6

0.8

1

1.2

Coefficient value

Figure 6.20 Distribution of the contour to center and ellipsoid fit coefficients for selected microstructures characterized by a uniquely described trend.

procedures, implementation of a proper optimization procedure, and selection of the most suitable results. Image analysis for the selected microstructure was performed first. The process of segmentation is very sophisticated and obtained results can be highly diversified even for the same input data and algorithm. This phenomenon depends mainly on the parameters established for the selected algorithm. Moreover, the automated assessment of results is very difficult, thus, it is hard to design a universal segmentation method able to work on various types of images. The methodology for analysis of DP microstructure images is proposed below. The procedure begins with the image pre-processing stage, which consists of the two supplementary approaches, that is, filtering and reconstruction. For the purposes of filtering the Dynamic Particles Algorithm (DPA) was designed and implemented [294]. The main advantage of the solution is its universality, which allows the application of the DPA method not only for images but also for signals and sophisticated multidimensional data. The basis of this method is in the definition of the particle, which can be treated as an N-dimensional object, which in the case of images is a set of 3D vectors (X-axis, Y-axis, grayscale value) related to particular pixel inside the analyzed image. The main idea of the DPA consists in the appropriate movement of each particle according to the following set of differential equations: 8 vi < m d~ vi 5 2 rVij 2 fc~ i dt (6.23) : d~ r 5~ v i dt where mi is mass; vi is the velocity of the i-th particle; Vij is the potential between particles i and j, dependent on the distance between them; and fc is a friction coefficient responsible mainly for the convergence of calculations.

276

(a)

Computational Materials Engineering

(b)

(c)

Figure 6.21 Results of image segmentation: (a) original DP steel image, (b) divided phases, (c) separated grains.

The magnitude of the force that causes the movement of the considered particle is reduced by the friction coefficient fc , 1. It was found that, according to the Newtonian laws of motion, if all of the pre-conditions would be fixed properly, the entire system would remain stable and convergent to the expected result. Additionally, the friction coefficient fc can be modified during the calculations of the algorithm, which influences the final smoothness of the results, as well as the stop criterion of the DPA. Image filtered by the described method is passed further to the process of image segmentation, which allows distinguishing each particular material grain as a separated object. The process of segmentation consists of two steps, that is, separation of phases and detailed analysis of grains inside each phase. Several images of the DP steel were analyzed and it was observed that there is a possibility to automatically find the grayscale threshold, which separates the bright and dark phases on the image. The value of such a threshold usually covers the first minimum or the first inflexion point located on the histogram on the left side of the peak [295]. Then, separated areas of different grains can be processed by detailed grain boundary detection algorithms to separate particular grains within each phase. At the moment, this process can be omitted, especially in the case of background phase, during the analysis for the purposes of the SSRVE creation. All of the earlier described shape coefficients do not take into account information about the matrix phase. Instead, the whole analysis is focused on the hard constituent phase. Calculation of shape coefficients is performed on the basis of the results of the segmentation, see typical results in Figure 6.21. Optimization procedure was performed next to find the best shape of hard constituent islands in the SSRVE. The selection of the cost function was crucial for the optimization procedure. A method for the construction of simple periodic structures for the special case of randomly

Trade off Between Accuracy and Efficiency

277

distributed circular inclusions with constant equal diameters was proposed in [280]. The positions of circular inclusions with given diameter were found by minimizing the objective function, which was defined as a square root error between spectral density of the periodic RVE and non-periodic real microstructure: vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u n "  2 # uX ζ 2ζ i iSSRVE Φ5t (6.24) wi ζi i51 where wi is weights;

n P i51

wi . 5 1; n is the number of coefficients; ζ i is the

average i-th shape coefficient of original microstructure; and ζ iSSRVE is the shape coefficient of the i-th SSRVE. At the first approach in the paper [296] the calculations were performed with the assumption that inside SSRVE there is only one inclusion. The implementation of the optimization procedure was based on genetic algorithm (GA), where chromosome is composed of m elements representing coordinates of crucial points determining SSRVE shape. These points were connected with spline functions forming smooth shape of the inclusion in the SSRVE. Calculation of the cost function was performed iteratively for each proposition of a new SSRVE shape. Few examples of calculations of the optimal SSRVE are presented in [296]. The example performed for the steel A (Figure 6.16) is presented below. The following parameters were assumed for optimization procedure: • Expected values obtained from 10-class histograms—phase fraction 5 26.3%, Malinowska coefficient 5 0.617, border index 5 0.855, maximum curvature 5 2.627. • Image parameters—original microstructure: width 5 640 px, height 5 427 px; SSRVE: width 5 height 5 300 px. • Optimization parameters—maximum number of iterations 5 20, minimum fitness 5 0.001, the cost functions composed of the following coefficients: phase fraction (w1 5 0.4), Malinowska coefficient (w2 5 0.2), border index (w3 5 0.2), max curvature (w4 5 0.2). The obtained results for different number of reference points are presented in Figure 6.22. The final value of the cost function for 5 reference points is 0.006085, for 10 points is 0.00098, for 15 points is 0.000961 and for 20 points is 0.000466. In the cases of SSRVE composed of 15 and 20 reference points, the optimization procedure stopped after 10 and 7 iterations, respectively.

278

Computational Materials Engineering

(a)

(b)

Figure 6.22 Example of original image (a) with related RVE (b) composed of periodically connected SSRVEs.

As mentioned, selection of suitable components of the objective function is crucial to obtain reliable SSRVE. That is why, presented approach can be additionally improved by application of multicriteria optimization objective function composed of, for example, shape coefficients, statistical measures and rheological behavior: vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi !ffi u k p k1l 3 X X X uX 2 Φ5t wi ζ 1 wi ϕ2 1 ws σ2 ðεj Þ (6.25) i

i51

sj

i2k

i5k11

s51

j51

ζi 5

ζ iRef 2 ζ iSSRVE ζ iRef

(6.26)

ϕi 5

ϕiRef 2 ϕiSSRVE ϕiRef

(6.27)

σsj 5

σsjRef 2 σsjSSRVE σsjRef

(6.28)

where Eq.(6.26) is a comparison of i-th shape coefficient, Eq.(6.27) is a comparison of statistical measures and Eq.(6.28) is a comparison of stresses obtained for different deformations setups of referential RVE microstructure and SSRVE (which depends on the number of dimensions of model), that is, uniform compression test without friction along each axis and pure shear test, wi are parameters weights, k is number of shape coefficients, l is number of statistical measures, s is number of rheological curves, p is number of iterations in numerical simulations. Such approach allows comparison of not only individual numbers but also the homogenized flow curves identified for RVE microstructure and SSRVE: I 1 ~εij 5 εij dS (6.29) S S

279

Trade off Between Accuracy and Efficiency

S,Mises (Avg: 75%) +1.420e+03 +1.308e+03 +1.195e+03 +1.083e+03 +9.707e+02 +8.583e+02 +7.460e+02 +6.337e+02 +5.213e+02 +4.090e+02 +2.967e+02 +1.843e+02 +7.200e+01

S,Mises (Avg: 75%) +1.420e+03 +1.304e+03 +1.188e+03 +1.072e+03 +9.562e+02 +8.403e+02 +7.243e+02 +6.084e+02 +4.924e+02 +3.765e+02 +2.605e+02 +1.446e+02 +2.861e+01

Figure 6.23 Example of strain distribution across the RVE and SSRVE obtained based on multicomponent objective function.

1 ~σij 5 S

I σij dS

(6.30)

S

Equations (6.29) and (6.30) are calculated in each iteration of numerical simulation (S is the area of the microstructure or SSRVE sample). This concept allows to fully capture the behavior of both the material flow under loading conditions as well as the appearance of individual phases. Example of obtained result is shown in Figure 6.23. The accuracy and reliability of the solution have to be always kept in mind when simplification of the model is applied. The final value of the cost function is a parameter, which represents accuracy of the optimization. This value may depend on applied optimization method, its parameters and randomness. The main sources of numerical errors are: • The procedures calculating perimeter and area of an analyzed object does not take into account real shape of an object but its digitalized approximation, for example, in the case of Malinowska coefficient, its value is sometimes lower than zero going out of range, what is caused by underestimated value of object’s perimeter in the numerator of the coefficient’s fraction. • Low reliability of calculation procedures estimating curvature of the object’s contour may cause problems in interpretation of results by using only singular value. Therefore, mean and total curvature parameters should be replaced with information about distribution of these parameters around the contour for more accurate description of all grains in original microstructure. • Dimensionality of coefficients’ units requires normalization of an image before calculation process starts.

280

Computational Materials Engineering

As presented identification of the SSRVE may also involve large number of finite element calculations, which is time consuming. Thus, to reduce the problem complexity, other alternative computational methods can be used. Examples of application of the Isogeometric Analysis to the SSRVE concept is presented as an example within the book.

6.2.2 Isogeometric analysis as alternative for the FE method in the SSRVE The main technologies that are used in the computer aided design (CAD) to represent exactly complex geometries are the B-splines and the Non-Uniform Rational B-Spline (NURBS), which were developed to improve designing process. The former approach only was considered in this book. There are many papers available in the scientific literature, which are focused on NURBS creation, applications and efficient algorithms for their evaluation [86,263,303]. Finite element method is now commonly used in various engineering applications to improve analyses of strains and stresses. This numerical technique enables finding approximate solutions to boundary value problems for partial differential equations. The developments of CAD and FEM software usually are independent of each other and, therefore, CAD geometry must be converted to the Analysis Suitable Geometry (ASG) to become the acceptable input for the FE model. To avoid numerical problems the small complicated details are omitted during creating ASG representation. The whole process can take up to 80% of the total time of analysis [147] and is commonly known as meshing or geometry discretization. Another approach was proposed by Hughes et al. [148], who stated that NURBS besides being used to describe the geometry, can also be used to construct finite approximations for analysis. This approach is called Isogeometric Analysis (IGA). It has been shown by various researchers that NURBS are better suited then FE discretization in application for many engineering problems. It is possible to obtain better accuracy from NURBS comparing to standard finite elements [28,29,32,371]. Many applications of IGA were done also for analysis of composite materials, for example, [205,47,71]. The differences between the standard FEM and IGA will be briefly discussed in application to the multiscale modeling of multiphase materials deformation in combination with SSRVE. 6.2.2.1 The general idea of the IGA Isogeometric analysis is a certain generalization of the FE analysis. For example, application of NURBS enables more precise geometric

Trade off Between Accuracy and Efficiency

281

representation of complex objects. Isogeometric analysis also simplifies mesh refinement procedures because the geometry is fixed at the coarsest level of refinement and remain unchanged throughout the refinement process. This eliminates geometrical errors and the necessity of linking the refinement procedure to a CAD representation of the geometry, as in classical FEA. Examples of evaluation of approximation and stability properties in the context of h-refinement will be presented within the chapter. Advantages of the IGA can be summarized as follows: • It enables precise geometric definition of complex engineering designs thus reducing errors caused by low-order, faceted geometric approximation by finite element. • It simplifies mesh refinement because even the coarsest model precisely represents the geometry. Thus, no link is necessary to the CAD geometry in order to refine the mesh, in contrast with the finite element method, in which each mesh represents a different approximation of the geometry. • It holds promise to simplify the mesh generation process, currently one of the most time consuming aspects of simulation. Possibility of application of the k-refinement process is another advantage of the IGA. This is unique feature of the IGA among geometrically flexible methodologies, which has been shown to possess significant accuracy and robustness properties, compared with the usual p-refinement utilized in FE methods [263,29]. k-refinement is a procedure in which the order of approximation is increased, as in the p-refinement, but continuity (i.e., smoothness) is likewise increased, in contrast to the p-method. In consequence, the IGA presents a unique combination of attributes that can be exploited on problems involving higher-order differential operators, namely, higherorder accuracy, robustness, two-and three-dimensional geometric flexibility, compact support, and, most importantly C1 and higher-order continuities. Mesh and geometry discretization. Differences in mesh generation for the FE and IGA approaches come off the NURBS basis. The mesh in IGA is directly defined by NURBS parameterization. The geometry domain discretization requires a node vector for each parametric direction: Ξξ 5 fξ 1 ; ξ 2 ; . . .; ξn1p11 g

(6.31)

Ξη 5 fη1 ; η2 ; . . .; ηn1q11 g

(6.32)

where p,q are the degree of basic functions in the two independent directions.

282

1.0 0.5 Y-Axis 0.0 –0.5

Y-Axis 0.0 –0.5 0.5

0.0 –0.5 X-Axis

–1.0

(b)

1.0

0.5

0.0 X-Axis

–0.5

–1.0

1.0

(a)

–1.0

–1.0

–0.5

Y-Axis 0.0

0.5

0.5

1.0

1.0

Computational Materials Engineering

–1.0

(c)

1.0

0.5

0.0

–0.5

–1.0

X-Axis

Figure 6.24 Control net (CAD) (a), IGA discretization (b), traditional FEM discretization (c).

Then new node vectors Ξ0ξ and Ξ0η are constructed by removing repeating values and as a result division of elements is obtained. These elements are not considered as new additional representations of the geometry (as it is the case of the FE method), but they are a part of the original mesh with its basis functions and control points. In consequence, the resulting element covers geometry exactly (Figure 6.24). This method can effectively reduce the error of geometric discretization and significantly improve the computational accuracy. Shape functions. Another major difference between IGA and traditional FE method is that the NURBS shape functions are not related to single element and its points but can also span over several elements. In FE method, shape functions are defined locally within one element. The degrees of freedom are called control variables instead of nodal variables. NURBS are constructed on the basis of B-splines, in which the control points and their basic functions influence the shape. The number of basis functions is equal to the number of control points (one basis function per control point) and the result of multiplication of control point position and its basis function determines contribution of this control point. B-spline basis functions are defined recursively based on the following equation and the knots vector, which represents coordinates in the parametric space:  1 ξi # ξ , ξ i11; p 5 0:Ni;0 ðξÞ 5 0 otherwise (6.33) ξ i1p11 2 ξ ξ 2 ξi Ni;p21 ðξÞ 1 Ni11;p21 ðξÞ p . 0:Ni;p ðξÞ 5 ξ i1p11 2 ξi11 ξi1p 2 ξi where p is the order of i-th basis function.

Trade off Between Accuracy and Efficiency

283

NURBS additionally introduces the weight of the of control points, which allows to specify the influence for each control point more precisely. It causes that the functions are defined more accurately for both curves and surfaces: Ni;p ðξÞwi p Ri ðξÞ 5 Pn ; i^51 Ni^;p ðξÞwi^

Ni;p ðξÞMj;q ðηÞwi;j p;q Ri;j ðξ; ηÞ 5 Pn Pm i^51 j^51 Ni^;p ðξÞMj^;q ðηÞwi^;j^ (6.34)

where p is the order of i-th basis function. 6.2.2.2 IGA mesh refinement There are several methods to refine NURBS in 2D and in 3D: insertion of additional nodes, increase in the order of the basis functions, or both. Inserting additional nodes is called h-refinement, increasing the order of approximation is referred to p-refinement and k-refinement considers both of mentioned methods simultaneously. In general refinement algorithms for IGA are less complicated then for standard FE solutions. Practical applications of the IGA are extensive and they are continuously increasing. Shape optimization used to the design of real three dimensional structural elements within the engineering community is definitely one of the most frequent applications. It is due to the earlier mentioned fact that the CAD models and numerical analysis models are separated from each other, and communication between them within the optimization process is costly. Another very common example of application of the IGA is within phase-field models. Traditionally the evolution of interface boundaries has been modeled using sharp interface models and the partial differential equations of the individual phases are coupled through interface boundary conditions, which poses major difficulties from the numerical point of view. As it was shown briefly in section 3.6.4 phase-field models provide an alternative description for phasetransition phenomena. The key idea is to replace sharp interfaces by thin transition regions where the interfacial forces are smoothly distributed. In consequence, the explicit front tracking is avoided. The properties of the IGA described at the beginning of this section open the way to application of this method to phase-field models. This application allows simulation of the higher-order operators using isogeometric analysis. Higher-order operators are encountered in biomedical applications and in many areas of engineering, such as, for example, liquid-liquid flows,

284

Computational Materials Engineering

liquid-vapor flows, emulsication, cancer growth, rotation-free thin shell theory, strain-gradient elastic and inelastic material models, dynamic crack propagation, etc. The simplicity of isogeometric analysis compared with many procedures that have been published in the literature is noteworthy. It is believed that it may be an effective procedure for solving higherorder differential equations on complex geometries. A noteworthy example of phase-field methodology is solution of the Cahn
Hilliard equation [115]. However, in the following section of this book we will focus on the application of the IGA to improve efficiency of the creation of the SSRVE for the multiphase materials. Example of simulations using this method is presented, as well. 6.2.2.3 Mesh refinement The h-refinement is a procedure, which is realized by node insertion. This procedure is well known in classical FE analysis. It enriches the basis by increasing the resolution of the parameter space. Initial node vector is given by Eq.(6.31). Based on this vector a new extended node vector can be introduced as: Ξ 5 fξ1 5 ξ1 ; ξ2 ; :::; ξn1m1p11 5 ξ n1p11 g

(6.35)

The condition that ΞAΞ has to be met. The application of the IGA uses the formulation of the FE method for the boundary value problem represented by equilibrium Eqs.(3.21). According to the FE analysis technique the new n 1 m basis functions are formed for new node vector (Ξ), and new n 1 m control points P 5 fP 1 ; P 2 ; :::; P n1m g are formed as linear combination of the original control points P 5 (P1, P1,. . .,Pn1m) by: 8 9 1; 1 # i # k 2 p; > > > > < ε 2 εi = ; k 2 p 1 1 # i # k; P 5 αi Pi 1 ð1 2 αi ÞPi21 ; αi 5 ε 2 ε i1p i > > > > : ; 0; k 1 1 # i # n 1 p 1 2: (6.36) As mentioned p-refinement is another mesh refinement method, which is well known from the FE analysis. In this method the degree of basis functions is increased. In the IGA the first step is dividing original NURBS (B-spline) to obtain separate Be´zier curves by node insertion based on Eq.(6.36). It raises the number of nodes and the polynomial degree. The procedure of increasing basis functions is applied to each

285

Trade off Between Accuracy and Efficiency

(a)

(b)

(c)

Figure 6.25 Basic control mesh and surface division (a), basic control mesh after knot insertion (b) and basic control mesh after raising of the interpolation degree (c).

Be´zier element and finally excess nodes are removed to obtain new NURBS, while the original continuity of basis functions is maintained. Figure 6.25 shows initial basic control mesh and surface division and basic control mesh after node insertion and after raising of the interpolation degree. Numerical simulations of 10% elastic deformation of the SSRVE for DP600 steel were performed and convergence of the solution was tested. This convergence was compared for FE and IGA approaches. After the first stimulation, the displacement values were checked in 400 points. Then refinement algorithm was performed and the results of simulation were obtained again for a new mesh. Displacements were evaluated again for the same points and the difference of results from these two simulations was compared. The following error coefficient was introduced: n pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi P ðΔxi Þ2 1 ðΔyi Þ2 e 5 i51 (6.37) n where n is the number of test points; Δxi, Δyi is the displacement difference in the x direction and in the y direction for point i before and after refinement. This procedure was repeated several times until a satisfactory small error values for both FEM and IGA solutions were obtained. Selected examples of the results are presented below. The four node elements were selected for FE simulations and Abaqus software was used as a refinement tool.

286

Computational Materials Engineering

(a) (c)

(b)

disp_y 0

disp_x 0.00123 0.001

–0.002 –0.004 0 –0.006 –0.008 –0.001 –0.00123

–0.01

Figure 6.26 Mesh with 161 nodes and 140 elements (a) and calculated displacements in the x (b) and y (c) directions.

(a)

(b)

(c) disp_x 0.00123 0.001

disp_y 0 –0.002 –0.004

0 –0.006 –0.008 –0.001 –0.00123

–0.01

Figure 6.27 Mesh with 951 nodes and 930 elements (a) and calculated displacements in the x (b) and y (c) directions.

During matrix assembling step, the Gaussian quadrature approximation was used for integration. Numerical simulations were performed for different numbers of elements to investigate and compare the convergence. The finite element mesh and the results of simulations of displacements are presented in Figure 6.26 for the mesh with 140 elements and in Figure 6.27 for the mesh with 930 elements. The same methodology was used for evaluation of the convergence for IGA simulations. Series of simulations for different number of elements was performed. It started from 16 elements for presented case study, which is the smallest number of elements needed to cover exactly the whole geometry. Before the next simulations the h-refinement was used in order to increase the number of elements and the selected results of simulations are presented in Figure 6.28 and Figure 6.29. The following parameters were used in the case 1: knots vectors 5 {0, 0, 0, 0.125, 0.25, 0.25, 0.375, 0.5, 0.5, 0.625, 0.75, 0.75, 0.875, 1, 1, 1}, 5 {0, 0, 0.5, 1, 1}, p 5 2, q 5 1 and 16 elements (Figure 5). These parameters for

287

Trade off Between Accuracy and Efficiency

(a) (b)

disp_x

(c)

0.000994

disp_y 0

0.0008 –0.002 0.0004

–0.004

0 –0.006

–0.0004

–0.008

–0.000725

–0.01

Figure 6.28 Isogeometrical representation for the case 1 (a) and calculated displacements in the x (b) and y (c) directions.

(a) (b)

(c) disp_x 0.00117 0.001

disp_y 0 –0.002 –0.004

0 –0.006 –0.008 –0.001 –0.00133

–0.01

Figure 6.29 Isogeometrical representation for the case 2 (a) and calculated displacements in the x (b) and y (c) directions.

the case 2 were as follows: knots vectors Ξξ 5 {0, 0, 0, 0.0625, 0.125, 0.1875, 0.25, 0.25, 0.3125, 0.375, 0.4375, 0.5, 0.5, 0.5625, 0.625, 0.6875, 0.75, 0.75, 0.8125, 0.875, 0.9375, 1, 1, 1}, Ξη 5 {0, 0, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1, 1}. Displacements calculated for p 5 2, q 5 1 and 128 elements are shown in Figure 6.30. Convergence of the two methods was compared next and results are shown in Figure 6.31. It is seen that FE requires much more elements to obtain similar error of the solution. As expected, numerical tests led to the conclusion that for considered problem the convergence of the IGA approach was better than for the standard FE. Obtained results were not surprising, and the reason of this is because IGA approach is devoid of discretization error and covers geometry exactly. Finally, the smaller number of elements causes that by using IGA it is possible to obtain a significant simplification of the computational domain. This contributes to a smaller final system of equations and enables decrease of the computing costs.

288

Computational Materials Engineering

disp_x 0.00117 0.001

disp_y

0

–0.002

–0.004 0 –0.006

–0.008 –0.001 –0.01

–0.00133

(a)

(b)

Figure 6.30 Displacements in the x (a) and y (b) directions calculated for p 5 2, q 5 1 and 128 elements.

1000

FEM IGA

Number of elements

800 600

400 200 0 0,0002

0,00016

0,00012

8E-005

4E-005

0

Figure 6.31 Convergence comparison for IGA and FE method.

To obtain even higher computational efficiency, an idea of parallelization for heterogeneous hardware architectures can also be applied (Figure 6.32). The idea of the solution is division of the tasks for main device and for additional computing devices. The main device would be responsible for storing global matrix of system of equations and scheduling work for other devices. The computing devices would be responsible for performing deformation simulations for each SSRVE.

289

Trade off Between Accuracy and Efficiency

Host machine, global system of equations

Compute devices (cards), performing SSRVE deformations for every Gauss point

Figure 6.32 Idea of parallelization of the IGA for heterogeneous hardware architectures.

Thus, as presented within the chapter, different methods of reduction of computational time can be applied simultaneously: model complexity reduction, selection of proper numerical approach as well as parallelization of the model.