Fatigue Crack Growth Life Maximization

This fatigue crack growth life maximization designs a structure such that the most cycles are required for the crack to grow from \(a_0\) (strating crack length) to \(a_{\text{end}}\) (final crack length). The theory behind the algorithm is explained in at Fatigue Crack Growth Life Maximization The crack path must be kown before running the optimization algorithms An example as how to use the optimization is shown in an example optimization example.py

Density Constraints
Load Cases
Finite Element Solvers
Optimization Module
Plotting Module

Density Constraints 

Constraints class used to specify the density constraints of the topology optimisation problem. It contains functions for minimum and maximum element density in the upcomming iteration and the magnitude of the volume constraint function itself of the current design. This version of the code is meant for the fatigue live maximization.

Bram Lagerweij Aerospace Structures and Materials Department TU Delft 2018

class src_FatigueLive.constraints.DensityConstraint(nelx, nely, move, volume_frac, density_min=0.0, density_max=1.0)

This object relates to the constraints used in this optimization. It can be used for the MMA updatescheme to derive what the limit is for all element densities at every itteration. The class itself is not changed by the itterations.

Parameters

nelx (int) – Number of elements in x direction.
nely (int) – Number of elements in y direction.
move (float) – Maximum change in density of an element over 1 itteration.
volume_frac (float) – Maximum volume that can be filled with material.
volume_derivative (2D array size(1, nelx*nely)) – Sensityvity of the density constraint to the density in each element.
density_min (float (optional)) – Minumum density, set at 0.0 if not specified.
density_max (float (optional)) – Maximum density, set at 0.0 if not specified.

nelx

Number of elements in x direction.

Type: int

nely

Number of elements in y direction.

Type: int

move

Maximum change in density of an element over 1 itteration.

Type: float

volume_frac

Maximum volume that can be filled with material.

Type: float

volume_derivative

Sensityvity of the density constraint to the density in each element.

Type: 2D array size(1, nelx*nely)

density_min

Minumum density, set at 0.0 if not specified.

Type: float, optional

density_max

Maximum density, set at 0.0 if not specified.

Type: float, optional

current_volconstrain(x)

Calculates the current magnitude of the volume constraint funcion:

\[V_{\text{constraint}} = \frac{\sum v_e X_e}{ V_{\max}}-1\]

Parameters: x (2D array size(nely, nelx)) – Density distribution of this itteration.
Returns: curvol – Curent value of the density constraint function.
Return type: float

xmax(x)

This function calculates the maximum density value of all ellements of this itteration.

Parameters: x (2D array size(nely, nelx)) – Density distribution of this itteration.
Returns: xmax – Maximum density values of this itteration after updating.
Return type: 2D array size(nely, nelx)

xmin(x)

This function calculates the minimum density value of all ellements of this itteration.

Parameters: x (2D array size(nely, nelx)) – Density distribution of this itteration.
Returns: xmin – Minimum density values of this itteration for the update scheme.
Return type: 2D array size(nely, nelx)

Load Cases 

This file contains the Load class that allows the generation of an object that contains geometric, mesh, loads and boundary conditions that belong to the load case. This version of the code is meant for the fatigue live maximization.

Bram Lagerweij Aerospace Structures and Materials Department TU Delft 2018

Parent Load Case

class src_FatigueLive.loads.Load(nelx, nely, young, Emin, poisson, ext_stiff, hoe)

Load parent class that contains the basic functions used in all load cases. This class and its children do contain information about the load case considered in the optimisation. The load case consists of the mesh, the loads, and the boundaries conditions. The class is constructed such that new load cases can be generated simply by adding a child and changing the function related to the geometry, loads and boundaries.

Parameters

nelx (int) – Number of elements in x direction.
nely (int) – Number of elements in y direction.
young (float) – Young’s modulus of the materials.
Emin (float) – Artificial Young’s modulus of the material to ensure a stable FEA. It is used in the SIMP based material model.
poisson (float) – Poisson ration of the material.
ext_stiff (float) – Extra stiffness to be added to global stiffness matrix. Due to interactions with mechanisms outside design domain.
hoe (dict) – Dictionary with for every cracklength the x end y element locations that need to be enriched.

nelx

Number of elements in x direction.

Type: int

nely

Number of elements in y direction.

Type: int

dim

Amount of dimensions considered in the problem, set at 2.

Type: int

edof

Dictionary containing list with all elements and their degree of freedom numbers for all crack_lengths, str(length) is the key.

Type: dict

x_list

Dictionary with a 1D list that contains the x indices of all degrees of freedom for all crack lengths, str(length) is the key.

Type: dict

y_list

Dictionary with a 1D list that contains the y indices of all degrees of freedom for all crack lengths, str(length) is the key.

Type: dict

num_dofs

Amount of degrees of freedom.

Type: int

young

Young’s modulus of the materials.

Type: float

Emin

Artificial Young’s modulus of the material to ensure a stable FEA. It is used in the SIMP based material model.

Type: float

poisson

Poisson ration of the material.

Type: float

k_list

Dictionary containing a list for every crack length, these lists contain the element stiffness matrices of full density for all elements, str(length) is the key.

Type: dict

kmin_list

Dictionary containing a list for every crack length, these lists contain the empty element stiffness matrices for all elements, str(length) is the key.

Type: list len(nelx*nely)

ext_stiff

Extra stiffness to be added to global stiffness matrix. Due to interactions with mechanisms outside design domain.

Type: float

alldofs()

Returns a list with all degrees of freedom.

Returns: all – List with numbers from 0 to the maximum degree of freedom number.
Return type: 1-D list

edofcalc(hoe)

Generates an array with the position of the nodes of each element in the global stiffness matrix. This takes the Higher Order Elements in account.

Parameters

hoe (list) – A list containing the x and y location of the higher order elemens for this crack length.

Returns

edof (2-D list size(nelx*nely, # degrees of freedom per element)) – The list with all elements and their degree of freedom numbers.
x_list (1-D array) – The list with the x indices of all elements to be inserted into the global stiffness matrix.
y_list (1-D array) – The list with the y indices of all elements to be inserted into the global stiffness matrix.
num_dofs (int) – The amount of degrees of freedom.

fixdofs(length_i)

Returns a list with indices that are fixed by the boundary conditions.

Parameters: length_i (int) – Length of the crack for the current mesh
Returns: fix – List with all the numbers of fixed degrees of freedom. This list is empty in this parent class.
Return type: 1-D list

force()

Returns an 1D array, the force vector of the loading condition.

Returns: f – Empty force vector.
Return type: 1-D column array length covering all degrees of freedom

freedofs(length_i)

Returns a list of arr indices that are not fixed

Parameters: length_i (int) – Length of the crack for the current mesh
Returns: free – List containing all elements of all dofs except those that appear in the freedos list.
Return type: 1-D list

import_stiffness(elementtype, E, nu)

This function imports a matrix from a csv file that has variables to the material properties. The correct material properties are added.

Parameters

elementtype (str) – Describes what .csv file should be used for the import.
E (float) – Young’s modulus of the material.
nu (float) – Poisson’s ratio of the material.

Returns

lk – Element stiffness matrix

Return type

array size(dofs, dofs)

kiloc()

The location of the stress intensity factor KI can be found at the second last index.

Returns: l – Zeros except for the second last index.
Return type: 1-D column array length covering all degrees of freedom

lk(E, nu, hoe)

Generates a list with all element stiffness matrices. It differentiates between the element types used.

Parameters

E (float) – Young’s modulus of the material.
nu (float) – Poisson’s ratio of the material.

Returns

k – Returns a list with all local stiffness matrices.

Return type

list len(nelx*nely)

node(elx, ely)

Calculates the topleft node number of the requested element. Does not toke Higher Order Elements in account.

Parameters

elx (int) – X position of the considered element.
ely (int) – Y position of the considered element.

Returns

topleft – The node number of the top left node.

Return type

int

nodes(elx, ely)

Calculates all node numbers of the requested element. Does not take Higher Order Elements in account.

Parameters

elx (int) – X position of the considered element.
ely (int) – Y position of the considered element.

Returns

n0 (int) – The node number of the bottom left node.
n1 (int) – The node number of the bottom right node.
n2 (int) – The node number of the top left node.
n3 (int) – The node number of the top right node.

passive()

Returns three lists containing the location and magnitude of fixed density values

Returns

elx (1-D list) – X coordinates of all passive elements, empty for the parent class.
ely (1-D list) – Y coordinates of all passive elements, empty for the parent class.
values (1-D list) – Density values of all passive elements, empty for the parent class.
fix_ele (1-D list) – List with all element numbers that are allowed to change.

reset_Kij(): Resets the global variable Kij. This is necessary as function import_stiffness will not clean up its local variables itself.

Child Load Cases

class src_FatigueLive.loads.EdgeCrack(nelx, nely, crack_length, young, Emin, poisson, ext_stiff)

Bases: src_FatigueLive.loads.Load

This child class of Load class represents the symmetric top half of an edge crack. The crack is positioned to the bottom left and propagates towards the right. Special elements are placed around the crack tip. The plate is subjected to a distributed tensile load (\(\sigma=1\)) on the top.

For a perfectly flat plate analytical expressions for K_I are known. 2

The stress intensity factors calculated can be be interperted in two ways:

Without scaling. This means that all elements have a size of 2 length units.
With scaling, comparison to reality should be based upon.

\[K^{\text{Real}} = K^{\text{FEA}}(\sigma=1) \sigma^{\text{Real}} \sqrt{\frac{a^{\text{Real}}}{2a^{\text{FEA}}}}\]

where \(a^{\text{FEA}}\) is the cracklength in number of elements.

Parameters

nelx (int) – Number of elements in x direction.
crack_length (array) – An array containing all crack lengths considered.
young (float) – Young’s modulus of the materials.
Emin (float) – Artificial Young’s modulus of the material to ensure a stable FEA. It is used in the SIMP based material model.
poisson (float) – Poisson ration of the material.
ext_stiff (float) – Extra stiffness to be added to global stiffness matrix. Due to interactions with mechanisms outside design domain.

crack_length

Is the amount of elements that the crack is long.

Type: int

hoe

List containing the x end y element locations that need to be enriched.

Type: list len(2)

References

2: Tada, H., Paris, P., & Irwin, G. (2000). “Part II 2.10-2.12 The Single Edge Notch Test Specimen”, The stress analysis of cracks handbook (3rd ed.). New York: ASME Press, pp:52-54.

fixdofs(length_i)

The boundary conditions limit y-translation at the bottom of the design space (due to symetry) and x-translations at the top (due to the clamps)

Parameters: length_i (int) – Length of the crack for the current mesh
Returns: fix – List with all the numbers of fixed degrees of freedom.
Return type: 1-D list

force(length_i)

The top of the design space is pulled upwards by 1MPa. This means that the nodal forces are 2 upwards, except for the top corners they have a load of 1 only.

Parameters: length_i (int) – Length of the crack for the current mesh
Returns: f – Force vector.
Return type: 1-D column array length covering all degrees of freedom

passive()

Returns three lists containing the location and magnitude of fixed density values. The elements around the crack tip are fixed at a density of one.

Returns

elx (1-D list) – X coordinates of all passive elements, empty for the parent class.
ely (1-D list) – Y coordinates of all passive elements, empty for the parent class.
values (1-D list) – Density values of all passive elements, empty for the parent class.
fix_ele (1-D list) – List with all element numbers that are allowed to change.

class src_FatigueLive.loads.CompactTension(nelx, crack_length, young, Emin, poisson, ext_stiff)

Bases: src_FatigueLive.loads.Load

This child class of Load class represents the symmetric top half of an compact tension specimen. The crack is positioned to the bottom left and propagates towards the right. Special elements are placed around the crack tip. The plate is subjected to upwards load of one. The design follows the ASTM standard. 3

For a perfectly flat plate analytical expressions for K_I do exist. 4

The stress intensity factors calculated can be be interpreted in two ways: 1. Without scaling. This means that all elements have a size of 2 length units. 2. With scaling, comparison to reality should be based upon.

\[K^{\text{Real}} = K^{\text{FEA}}(F=1) F^{\text{Real}} \sqrt{\frac{2W^{\text{FEA}}}{W^{\text{Real}}}}\]

where \(W^{\text{FEA}}\) is the width in number of elements.

Parameters

nelx (int) – Number of elements in x direction.
crack_length (array) – An array containing all crack lengths considered.
young (float) – Young’s modulus of the materials.
Emin (float) – Artificial Young’s modulus of the material to ensure a stable FEA. It is used in the SIMP based material model.
poisson (float) – Poisson ration of the material.
ext_stiff (float) – Extra stiffness to be added to global stiffness matrix. Due to interactions with mechanisms outside design domain.

nely

Number of y elements, this is now a function of nelx.

Type: int

crack_length

Is for all cracks considered the crack_length.

Type: array

References

3: ASTM Standard E647-15e1, “Standard Test Method for Measurement of Fatigue Crack Growth Rates,” ASTM Book of Standards, vol. 0.30.1, 2015.
4: Tada, H., Paris, P., & Irwin, G. (2000). “Part II 2.19-2.21 The Compact Tension Test Specimen”, The stress analysis of cracks handbook (3rd ed.). New York: ASME Press, pp:61-63.

fixdofs(length_i)

The bottom of the design space is fixed in y direction (due to symmetry around the x axis). While at the location that the load is introduced x translations are constraint.

Parameters: length_i (int) – Length of the crack for the current mesh
Returns: fix – List with all the numbers of fixed degrees of freedom.
Return type: 1-D list

force(length_i)

The ASTM standard requires the force to be located approx. 1/5 of nelx and at 0.195 * nely from the top.

Parameters: length_i (int) – Length of the crack for the current mesh
Returns: f – Force vector
Return type: 1-D column array length covering all degrees of freedom

passive()

Returns three lists containing the location and magnitude of fixed density values. The elements around the crack tip are fixed at a density of one.

Returns

elx (1-D list) – X coordinates of all passive elements, empty for the parent class.
ely (1-D list) – Y coordinates of all passive elements, empty for the parent class.
values (1-D list) – Density values of all passive elements, empty for the parent class.
fix_ele (1-D list) – List with all element numbers that are allowed to change.

Finite Element Solvers 

Finite element solvers for the displacement from stiffness matrix, force and adjoint vector. This version of the code is meant for the fatigue crack growth maximization.

Bram Lagerweij Aerospace Structures and Materials Department TU Delft 2018

Parent Solver

class src_FatigueLive.fesolvers.FESolver(verbose=False)

This parent FEA class can only assemble the global stiffness matrix and exclude all fixed degrees of freedom from it. This function, gk_freedofs is used in all FEA solvers classes. The displace function is not implemented in this parrent class as it does not contain a solver for the linear problem.

Parameters: verbose (bool, optional) – False if the FEA should not print updates

verbose

False if the FEA should not print updates.

Type: bool

displace(load, x, penal, length)

FE solver based upon the sparse SciPy solver that uses umfpack.

Parameters

load (object, child of the Loads class) – The loadcase(s) considerd for this optimisation problem.
x (2-D array size(nely, nelx)) – Current density distribution.
penal (float) – Material model penalisation (SIMP).
length (int) – Length of the current crack conciderd.

Returns

u (1-D column array shape(max(edof), 1)) – The displacement vector.
lamba (1-D column array shape(max(edof), 1)) – Adjoint equation solution.

gk_freedofs(load, x, penal, length)

Generates the global stiffness matrix with deleted fixed degrees of freedom. It generates a list with stiffness values and their x and y indices in the global stiffness matrix. Some combination of x and y appear multiple times as the degree of freedom might apear in multiple elements of the FEA. The SciPy coo_matrix function adds them up at the background. At the location of the force introduction and displacement output an external stiffness is added due to stability reasons.

Parameters

load (object, child of the Loads class) – The loadcase(s) considerd for this optimisation problem.
x (2-D array size(nely, nelx)) – Current density distribution.
penal (float) – Material model penalisation (SIMP).
length (int) – Length of the current crack conciderd.

Returns

k – Global stiffness matrix without fixed degrees of freedom.

Return type

2-D sparse csc matrix

Child Solvers

class src_FatigueLive.fesolvers.CvxFEA(verbose=False)

Bases: src_FatigueLive.fesolvers.FESolver

This parent FEA class can assemble the global stiffness matrix and solve the FE problem with a Supernodal Sparse Cholesky Factorization. It solves for both the equalibrium and adjoint problem.

verbose

False if the FEA should not print updates.

Type: bool

displace(load, x, penal, length)

FE solver based upon a Supernodal Sparse Cholesky Factorization. It requires the instalation of the cvx module. It solves both the FEA equalibrium and adjoint problems. 1

Parameters

load (object, child of the Loads class) – The loadcase(s) considerd for this optimisation problem.
x (2-D array size(nely, nelx)) – Current density distribution.
penal (float) – Material model penalisation (SIMP).
length (int) – Length of the current crack conciderd.

Returns

u (1-D column array shape(max(edof), 1)) – The displacement vector.
lamba (1-D column array shape(max(edof), 1)) – Adjoint equation solution.

References

1: Y. Chen, T. A. Davis, W. W. Hager, S. Rajamanickam, “Algorithm 887: CHOLMOD, Supernodal Sparse Cholesky Factorization and Update/Downdate”, ACM Transactions on Mathematical Software, 35(3), 22:1-22:14, 2008.

class src_FatigueLive.fesolvers.CGFEA(verbose=False)

Bases: src_FatigueLive.fesolvers.FESolver

This parent FEA class can assemble the global stiffness matrix and solve the FE problem with a sparse solver based upon a preconditioned conjugate gradient solver. The preconditioning is based upon the inverse of the diagonal of the stiffness matrix.

Recomendations

Make the tolerance change over the iterations, low accuracy is required for first itteration, more accuracy for the later ones.
Add more advanced preconditioner.
Add gpu accerelation.

verbose

False if the FEA should not print updates.

Type: bool

ufree_old

Displacement field of previous iteration for every crack length, the keys are the related cracklengths.

Type: dict

lambafree_old

Ajoint equation result of previos iteration for every crack length, the keys are the related cracklengths.

Type: array len(freedofs)

displace(load, x, penal, length)

FE solver based upon the sparse SciPy solver that uses a preconditioned conjugate gradient solver, preconditioning is based upon the inverse of the diagonal of the stiffness matrix. Currently the relative tolerance is hardcoded as 1e-3.

Parameters

load (object, child of the Loads class) – The loadcase(s) considerd for this optimisation problem.
x (2-D array size(nely, nelx)) – Current density distribution.
penal (float) – Material model penalisation (SIMP).
length (int) – Length of the current crack conciderd.

Returns

u (1-D array len(max(edof)+1)) – Displacement of all degrees of freedom
lamba (1-D column array shape(max(edof), 1)) – Adjoint equation solution.

Optimization Module 

Topology Optimization class that handles the itterations, objective functions, filters and update scheme. It requires to call upon a constraint, load case and FE solver classes. This version of the code is meant for the fatigue live maximization.

Bram Lagerweij Aerospace Structures and Materials Department TU Delft 2018

class src_FatigueLive.topopt.Topopt(constraint, load, fesolver, weights, C, m, verbose=False, x0_loc=None)

This is the optimisation object itself. It contains the initialisation of the density distribution.

Parameters

constraint (object of DensityConstraint class) – The constraints for this optimization problem.
load (object, child of the Loads class) – The loadcase(s) considerd for this optimisation problem.
fesolver (object, child of the CSCStiffnessMatrix class) – The finite element solver.
weights (array length load cases) – The weight given to each of the load cases.
C (float) – Multiplication part of Paris-Erdogan law.
m (float) – Power part of Paris-Erdogan law.
verbose (bool) – Printing itteration results.
x0_loc (str) – Set initial design with numpy ‘.npy’ file location.

constraint

The constraints for this optimization problem.

Type: object of DensityConstraint class

load

The loadcase(s) considerd for this optimisation problem.

Type: object, child of the Loads class

fesolver

The finite element solver.

Type: object, child of the CSCStiffnessMatrix class

verbose

Printing itteration results.

Type: bool

itr

Number of iterations performed

Type: int

weights

The weight given to each of the load cases.

Type: array length load cases

C

Multiplication part of Paris-Erdogan law.

Type: float

m

Power part of Paris-Erdogan law.

Type: float

free_ele

All element nubers that ar allowed to change.

Type: 1-D list

x

Array containing the current densities of every element.

Type: 2-D array size(nely, nelx)

xold1

Flattend density distribution one iteration ago.

Type: 1D array len(nelx*nely)

xold2

Flattend density distribution two iteration ago.

Type: 1D array len(nelx*nely)

low

Column vector with the lower asymptotes, calculated and used in the MMA subproblem of the previous itteration.

Type: 1D array len(nelx*nely)

upp

Column vector with the lower asymptotes, calculated and used in the MMA subproblem of the previous itteration.

Type: 1D array len(nelx*nely)

densityfilt(rmin, filt)

Filters with a normalized convolution on the densities with a radius of rmin if:

>>> filt=='density'

The relusting geometry retains passive elements.

Parameters

rmin (float) – Filter size.
filt (str) – The filter type that is selected, either ‘sensitivity’ or ‘density’.

Returns

xf – Filterd density distribution.

Return type

2-D array size(nely, nelx)

iter(penal, rmin, filt)

This funcion performs one itteration of the topology optimisation problem. It

loads the constraints,
calculates the stiffness matrices,
executes the density filter,
executes the FEA solver,
calls upon the displacment objective and its sensitivity calculation,
executes the sensitivity filter,
executes the MMA update scheme,
and finaly updates density distribution (design).

Parameters

penal (float) – Material model penalisation (SIMP).
rmin (float) – Filter size.
filt (str) – The filter type that is selected, either ‘sensitivity’ or ‘density’.

Returns

change (float) – Largest difference between the new and old density distribution.
volcon (float) – Amount of volume of this itteration.
N (float) – Fatigue live in cycles of the crack.
Obj (float) – Objective in weigted cycles.

kicalc(x, u, lamba, penal, length)

This fuction calculates displacement of the objective node and its sensitivity to the densities.

Parameters

x (2-D array size(nely, nelx)) – Possibly filterd density distribution.
u (1-D array size(max(edof), 1)) – Displacement of all degrees of freedom.
lamba (2-D array size(max(edof), 1)) –
ke (2-D array size(8, 8)) – Element stiffness matrix with full density.
penal (float) – Material model penalisation (SIMP).
length (int) – Length of the crack conciderd.

Returns

ki (float) – Displacement objective.
dki (2-D array size(nely, nelx)) – Displacement objective sensitivity to density changes.

layout(penal, rmin, delta, loopy, filt, history=False)

Solves the topology optimisation problem by looping over the iter function.

Parameters

penal (float) – Material model penalisation (SIMP).
rmin (float) – Filter size.
delta (float) – Convergence is roached when delta > change.
loopy (int) – Amount of iteration allowed.
filt (str) – The filter type that is selected, either ‘sensitivity’ or ‘density’.
history (bool) – Do the intermediate results need to be stored.

Returns

xf (array size(nely, nelx)) – Density distribution resulting from the optimisation.
xf_history (list of arrays len(itterations size(nely, nelx))) – List with the density distributions of all itterations, None when history != True.
ki (array) – Stress intensity factor at each crack length increment.

mma(m, n, itr, xval, xmin, xmax, xold1, xold2, f0val, df0dx, fval, dfdx, low, upp, a0, a, c, d)

This function mmasub performs one MMA-iteration, aimed at solving the nonlinear programming problem:

\[\begin{split}\min & f_0(x) & + a_0z + \sum_{i=1}^m \left(c_iy_i + \frac{1}{2}d_iy_i^2\right) \\ &\text{s.t. }& f_i(x) - a_iz - y_i \leq 0 \hspace{1cm} & i \in \{1,2,\dots,m \} \\ & & x_{\min} \geq x_j \geq x_{\max} & j \in \{1,2,\dots,n \} \\ & & y_i \leq 0 & i \in \{1,2,\dots,m \} \\ & & z \geq 0\end{split}\]

Parameters

m (int) – The number of general constraints.
n (int) – The number of variables \(x_j\).
itr (int) – Current iteration number (=1 the first time mmasub is called).
xval (1-D array len(n)) – Vector with the current values of the variables \(x_j\).
xmin (1-D array len(n)) – Vector with the lower bounds for the variables \(x_j\).
xmax (1-D array len(n)) – Vector with the upper bounds for the variables \(x_j\).
xold1 (1-D array len (n)) – xval, one iteration ago when iter>1, zero othewise.
xold2 (1-D array len (n)) – xval, two iteration ago when iter>2, zero othewise.
f0val (float) – The value of the objective function \(f_0\) at xval.
df0dx (1-D array len(n)) – Vector with the derivatives of the objective function \(f_0\) with respect to the variables \(x_j\), calculated at xval.
fval (1-D array len(m)) – Vector with the values of the constraint functions \(f_i\), calculated at xval.
dfdx (2-D array size(m x n)) – (m x n)-matrix with the derivatives of the constraint functions \(f_i\). with respect to the variables \(x_j\), calculated at xval.
low (1-D array len(n)) – Vector with the lower asymptotes from the previous iteration (provided that iter>1).
upp (1-D array len(n)) – Vector with the upper asymptotes from the previous iteration (provided that iter>1).
a0 (float) – The constants \(a_0\) in the term \(a_0 z\).
a (1-D array len(m)) – Vector with the constants \(a_i1 in the terms :math:\).
c (1-D array len(m)) – Vector with the constants \(c_i\) in the terms \(c_i*y_i\).
d (1-D array len(m)) – Vector with the constants \(d_i\) in the terms \(0.5d_i (y_i)^2\).

Returns

xmma (1-D array len(n)) – Column vector with the optimal values of the variables \(x_j\) in the current MMA subproblem.
low (1-D array len(n)) – Column vector with the lower asymptotes, calculated and used in the current MMA subproblem.
upp (1-D array len(n)) – Column vector with the upper asymptotes, calculated and used in the current MMA subproblem.
Version September 2007 (and a small change August 2008)
Krister Svanberg <krille@math.kth.se>
Department of Mathematics KTH, SE-10044 Stockholm, Sweden.
Translated to python 3 by A.J.J. Lagerweij TU Delft June 2018

sensitivityfilt(x, dOi, rmin, filt)

Filters with a normalized convolution on the sensitivity with a radius of rmin if:

>>> filt=='sensitivity'

Parameters

x (2-D array size(nely, nelx)) – Current density ditribution.
dOi (2-D array size(nely, nelx) – Objective sensitivity to density changes.
rmin (float) – Filter size.
filt (str) – The filter type that is selected, either ‘sensitivity’ or ‘density’.

Returns

dOif – Filterd sensitivity distribution.

Return type

2-D array size(nely, nelx)

solvemma(m, n, epsimin, low, upp, alfa, beta, p0, q0, P, Q, a0, a, b, c, d)

This function solves the MMA subproblem with a primal-dual Newton method:

\[\begin{split}\min &\sum_{j-1}^n& \left( \frac{p_{0j}^{(k)}}{U_j^{(k)} - x_j} + \frac{q_{0j}^{(k)}}{x_j - L_j^{(k)}} \right) + a_0z + \sum_{i=1}^m \left(c_iy_i + \frac{1}{2}d_iy_i^2\right) \\ &\text{s.t. }& \sum_{j-1}^n \left(\frac{p_{ij}^{(k)}}{U_j^{(k)} - x_j} + \frac{q_{ij}^{(k)}}{x_j - L_j^{(k)}} \right) - a_iz - y_i \leq b_i \\ & & \alpha_j \geq x_j \geq \beta_j\\ & & z \geq 0\end{split}\]

Returns: x – Column vector with the optimal values of the variables x_j in the current MMA subproblem.
Return type: 1-D array len(n)

Plotting Module 

Plotting the simulated TopOpt geometry with boundery conditions and loads. This version of the code is meant for mixed element types problems for the fatigue live maximization.

Bram Lagerweij Aerospace Structures and Materials Department TU Delft 2018

class src_FatigueLive.plotting.Plot(load, directory, title=None)

This class contains functions that allows the visualisation of the TopOpt algorithem. It can print the density distribution, the boundary conditions and the forces.

Parameters

load (object, child of the Loads class) – The loadcase(s) considerd for this optimisation problem.
directory (str) – Relative directory that the results should be saved in.
title (str, optional) – Title of the plot, optionaly.

nelx

Number of elements in x direction.

Type: int

nely

Number of elements in y direction.

Type: int

fig

An empty figure of size nelx/10 and nely/10 inch.

Type: matplotlib.pyplot figure

ax

The axis system that belongs to fig.

Type: matplotlib.pyplot axis

images

This list contains all density distributions that need to be plotted.

Type: 1-D list with imshow objects

directory

Location where the results need to be saved.

Type: str

add(x, animated=False)

Adding a plot of the density distribution to the figure.

Parameters

x (2-D array size(nely, nelx)) – The density distribution.
animated (bool) – An animated figure is genereted when history = True.

boundary(load)

Plotting the boundary conditions.

Parameters: load (object, child of the Loads class) – The loadcase(s) considerd for this optimisation problem.

find(dof)

This function returns the location, x,y of any degree of freedom by corresponding it with the edof array.

Parameters

dof (int) – Degree of freedom number of unknown location.

Returns

x (float) – x location of the dof.
y (float) – y location of the dof.

loading(load)

Plotting the loading conditions.

Parameters: load (object, child of the Loads class) – The loadcase(s) considerd for this optimisation problem.

save(filename, fps=10)

Saving an plot in svg or mp4 format, depending on the length of the images list. The FasterFFMpegWriter is used when videos are generated. These videos are encoded with a hardware accelerated h264 codec with the .mp4 file format. Other codecs and encoders can be set within the function itself.

Parameters

filename (str) – Name of the file, excluding the file exstension.
fps (int) – Amount of frames per second if the plots are animations.

saveXYZ(x, x_size, thickness=1)

This function allows the export of the density distribution as a point cloud. This can be used to create .stl files in the following steps:

Open meshlab and ‘import mesh’ on all .xyz files.
Use ‘Per Vertex Normal Fnction’ on all point clouds.
- bot with [nx, ny, nz] = [ 0, 0,-1]
- top with [nx, ny, nz] = [ 0, 0, 1]
- x- with [nx, ny, nz] = [-1, 0, 0]
- x+ with [nx, ny, nz] = [ 1, 0, 0]
- y- with [nx, ny, nz] = [ 0,-1, 0]
- y+ with [nx, ny, nz] = [ 0, 1, 0]

3. Apply the ‘Screened Poisson Surface Reconstruction’ filter with the option of ‘Merge all visible layers’ as True

Parameters

x (2-D array) – Density array.
x_size (float) – X dimension of the mesh.
thickness (foat) – Thickness of the mesh.

show(): Showing the plot in a window.

class src_FatigueLive.plotting.FasterFFMpegWriter(**kwargs)

Bases: matplotlib.animation.FFMpegWriter

FFMpeg-pipe writer bypassing figure.savefig. To improof saving speed

classmethod bin_path(): Return the binary path to the commandline tool used by a specific subclass. This is a class method so that the tool can be looked for before making a particular MovieWriter subclass available.

cleanup(): Clean-up and collect the process used to write the movie file.

finish(): Finish any processing for writing the movie.

frame_size: A tuple (width, height) in pixels of a movie frame.

grab_frame(**savefig_kwargs)

Grab the image information from the figure and save as a movie frame.

Doesn’t use savefig to be faster: savefig_kwargs will be ignored.

classmethod isAvailable(): Check to see if a MovieWriter subclass is actually available.

saving(fig, outfile, dpi, *args, **kwargs)

Context manager to facilitate writing the movie file.

*args, **kw are any parameters that should be passed to setup.

setup(fig, outfile, dpi=None)

Perform setup for writing the movie file.

Parameters

fig (~matplotlib.figure.Figure) – The figure object that contains the information for frames
outfile (str) – The filename of the resulting movie file
dpi (int, optional) – The DPI (or resolution) for the file. This controls the size in pixels of the resulting movie file. Default is fig.dpi.

Fatigue Crack Growth Life Maximization

Density Constraints

Load Cases

Parent Load Case

Child Load Cases

Finite Element Solvers

Parent Solver

Child Solvers

Optimization Module

Plotting Module

Density Constraints 

Load Cases 

Finite Element Solvers 

Optimization Module 

Plotting Module 