A Novel Algorithm for the Determination of Walker Damage in Loaded Disc Springs

Abstract

In this paper, a novel algorithm for the determination of Walker damage in loaded disc springs is presented. The algorithm takes a 3D-scan of a disc spring, measured residual stresses, material parameters, and spring loads as inputs. It outputs a distribution of Walker damage over the surface area of the input disc spring. As the algorithm allows a fully automated determination of the Walker damage, it can be used by disc spring manufacturers to reduce the working time spent on this task by specialized engineers significantly. Compared to spreadsheet applications using analytical formulas and finite element models using idealized geometry, this approach offers a superior description of the stress states in disc springs.

1. Introduction

The fatigue behavior of disc springs can be described by models ranging widely in complexity. Simple models like the one described in [1,2] can be implemented using spreadsheet applications, e.g., Microsoft Excel. As models become more sophisticated, the expense of building and evaluating the model rises. With the latest step in creating more complex models, the introduction of scanned geometries, superposed residual stresses, and the Walker damage parameter, the need for a novel algorithm for the determination of Walker damage in loaded disc springs has arisen.

The algorithm described in this paper is implemented in the Spring_stack Python module, which can also be used for other applications. For more information on the simulation of single disc springs without a 3D-scanned geometry or with multiple springs in one assembly, see [3,4].

The algorithm is the first published algorithm to build and evaluate a finite element model from a 3D-scanned geometry without user interaction. It allows the user to obtain an understanding of the influence of geometric deviations that are present in a batch of disc springs on the lifetime of individual disc springs under cyclic loading. This would be impossible without the algorithm because manually conducting a finite element simulation for each disc spring is prohibitively expensive.

Introductions to the mathematical description of disc springs and to the Walker damage parameter are given in Section 2. In Section 3, an algorithm used to describe Walker damage at the surface of a 3D-scanned disc spring is presented. Its implementation is described in Section 4. An example application is presented in Section 5. We recommend that readers without programming experience read Section 5 before Section 4.

2. Related Research

2.1. Mechanical Behavior of Disc Springs

The first mention of conical disc springs with a rectangular cross-section in the literature was Belleville’s patent in 1861 [5,6,7]. The first formulas to compute the characteristic of disc springs, as well as the stresses present in disc springs under load were published by Almen and Laszlo in 1936. These formulas are still in use in today’s standards [1,2], supplemented by a friction formulation published by Curti and Montanini [8]. There is a variety of other approaches for the analytical assessment of disc springs, each bringing its own advantages, such as improved accuracy [9,10,11,12,13,14,15,16,17,18,19], improved simplicity [20], applicability for different geometries [21,22,23,24,25,26,27,28], consideration of new material laws [29,30,31], offering a new concept for friction [32], or allowing the computation of resonance frequencies [33].

Analytical assessment of disc springs allows a direct understanding of the mechanical nature of disc springs and even the direct identification of links between geometric features and stresses. It also offers a sufficiently precise description of disc springs for applications like fatigue design according to EN 16983 and EN 16984 [1,2] at low computational costs. They can also be automated easily, which makes these approaches attractive for optimization algorithms like [34,35] and for analytical models of systems containing multiple components like [36,37].

Since the 1980s [38], disc springs have also been assessed using finite element analysis (FEA). While FEA is expensive in licensing, computational cost, and training, it allows models to be adapted to new, similar problems quickly. The first FEA models were used to exceed the accuracy of the analytical approaches available at that time [39,40]. FEA is still in use for the verification of new analytical models. FEA has extensively been used to describe residual stresses and changes in the characteristics created by plastic deformation and creep effects [41,42,43,44,45,46]. It has also been used to describe disc springs with complex geometries [47,48,49,50,51] or complex load cases [52], as well as for disc springs made from new materials [48,53,54,55,56,57] and to describe the behavior of assemblies containing disc springs [58,59,60].

Today, 2D models of disc springs solve quickly. An automated tool for the processing of 2D simulations of disc springs was implemented as early as 2000 [61]. FEA has been included in the curricula of most engineering programs. Therefore, graduates with basic FEA skills are available to companies. Disc springs have been numerically simulated by free finite element software [62]. These effects have led to the literature being split roughly in half between analytical formulas and FEA. The authors believe this distribution to roughly stay the same in the future because both approaches offer unique advantages and because of a legacy effect. The legacy effect is created by fatigue tests being evaluated using a certain method and the raw data like material laws and spring geometries not being documented. Furthermore, engineers are more experienced in the design of springs using analytical formulas, and this experience can only partially be used for the design of springs using FEA.

2.2. The Walker Damage Parameter

The Walker damage parameter [63] helps to compare the fatigue behavior of materials under different mean stresses. Compared with other approaches to the characterization of mean stress effects in steels [64,65,66], the Walker damage parameter shows a superior lifetime prediction [67]. In this paper, the stress-based approach is described. For the strain-based approach, see [68].

The Walker damage parameter $P_{\mathrm{Walker}}$ is computed from the maximum stress ${\sigma }_{\max }$, the stress amplitude ${\sigma }_{\mathrm{amp}}$, and the Walker exponent $\gamma$:

$$P_{\mathrm{Walker}}={\sigma }_{\max }^{1-\gamma }{\sigma }_{\mathrm{amp}}^{\gamma }$$

(1)

The Walker exponent $\gamma$ is a material parameter that can be identified by fitting fatigue-curves with different mean stresses or R-ratios. A high Walker exponent implies a low sensitivity to mean stress effects, while a low Walker exponent implies a high sensitivity. If the Walker exponent $\gamma$ is fixed to 0.5, the stress Walker approach is equivalent to the stress Smith–Watson–Topper [66] approach.

The Walker equation requires ${\sigma }_{\max }$ and ${\sigma }_{\mathrm{amp}}$ to be scalar values. At the surface of the disc spring, the stress state is two-dimensional; therefore, an equivalent stress must be computed. Here, the von Mises equivalent stress is used.

For proportional loads in Quadrant 1, the computation of the maximum stress ${\sigma }_{\max }$ and the stress amplitude ${\sigma }_{\mathrm{amp}}$ from the minimum and maximum principal stresses ${\sigma }_{\mathrm{I},\max }$, ${\sigma }_{\mathrm{II},\max }$, ${\sigma }_{\mathrm{I},\min }$, and ${\sigma }_{\mathrm{II},\min }$ is obvious.

$${\sigma }_{\min }=\sqrt{\frac{{\sigma }_{\mathrm{I},\min }^2+{\left({\sigma }_{\mathrm{I},\min }-{\sigma }_{\mathrm{II},\min }\right)}^2+{\sigma }_{\mathrm{II},\min }^2}{2}}$$

(2)

$${\sigma }_{\max }=\sqrt{\frac{{\sigma }_{\mathrm{I},\max }^2+{\left({\sigma }_{\mathrm{I},\max }-{\sigma }_{\mathrm{II},\max }\right)}^2+{\sigma }_{\mathrm{II},\max }^2}{2}}$$

(3)

$${\sigma }_{\mathrm{amp}}=0.5·\left({\sigma }_{\max }-{\sigma }_{\min }\right)$$

(4)

The stress state in disc springs however is non-proportional, especially in shot peened specimens. An example why the formulas given for the proportional case may give inconsistent results for the non-proportional case is depicted in Figure 1. Given the stress states ${\sigma }_{\mathrm{a}}$, ${\sigma }_{\mathrm{b}}$, and ${\sigma }_{\mathrm{c}}$ in Quadrant 1, the stress amplitude between ${\sigma }_{\mathrm{a}}$ and ${\sigma }_{\mathrm{c}}$ is equal to the stress amplitude between ${\sigma }_{\mathrm{b}}$ and ${\sigma }_{\mathrm{c}}$ if the von Mises equivalent stress in ${\sigma }_{\mathrm{a}}$ is equal to the von Mises equivalent stress in ${\sigma }_{\mathrm{b}}$. The same applies for mean stresses. In reality, the damage inflicted by cycling between ${\sigma }_{\mathrm{b}}$ and ${\sigma }_{\mathrm{c}}$ is smaller than the damage inflicted by cycling between ${\sigma }_{\mathrm{a}}$ and ${\sigma }_{\mathrm{c}}$.

To avoid the described misrepresentation of damage inflicted, the modified Manson–McKnight method [69] is used. Instead of computing two scalar stress states ${\sigma }_{\max }$ and ${\sigma }_{\min }$ and deducing a scalar amplitude and a scalar mean stress, a two-dimensional stress amplitude and a two-dimensional mean stress are deduced from a two-dimensional maximum and minimum stress state:

$${\sigma }_{\mathrm{I},\mathrm{amp}}=0.5·\left({\sigma }_{\mathrm{I},\max }-{\sigma }_{\mathrm{I},\min }\right)$$

(5)

$${\sigma }_{\mathrm{II},\mathrm{amp}}=0.5·\left({\sigma }_{\mathrm{II},\max }-{\sigma }_{\mathrm{II},\min }\right)$$

(6)

$${\sigma }_{\mathrm{I},\mathrm{mean}}=0.5·\left({\sigma }_{\mathrm{I},\max }+{\sigma }_{\mathrm{I},\min }\right)$$

(7)

$${\sigma }_{\mathrm{II},\mathrm{mean}}=0.5·\left({\sigma }_{\mathrm{II},\max }+{\sigma }_{\mathrm{II},\min }\right)$$

(8)

These pairs of stress components are converted into von Mises equivalent stresses:

$${\sigma }_{\mathrm{amp}}=\sqrt{\frac{{\sigma }_{\mathrm{I},\mathrm{amp}}^2+{\left({\sigma }_{\mathrm{I},\mathrm{amp}}-{\sigma }_{\mathrm{II},\mathrm{amp}}\right)}^2+{\sigma }_{\mathrm{II},\mathrm{amp}}^2}{2}}$$

(9)

$${\sigma }_{\mathrm{mean}}=\frac{{\sigma }_{\mathrm{I},\mathrm{mean}}-{\sigma }_{\mathrm{III},\mathrm{mean}}}{{\sigma }_{\mathrm{I},\mathrm{mean}}+{\sigma }_{\mathrm{III},\mathrm{mean}}}·\sqrt{\frac{{\sigma }_{\mathrm{I},\mathrm{mean}}^2+{\left({\sigma }_{\mathrm{I},\mathrm{mean}}-{\sigma }_{\mathrm{II},\mathrm{mean}}\right)}^2+{\sigma }_{\mathrm{II},\mathrm{mean}}^2}{2}}$$

(10)

The maximum stress ${\sigma }_{\max }$ is defined with a lower bound to avoid complex numbers as Walker damage parameters.

$${\sigma }_{\max }=max\left({\sigma }_{\mathrm{mean}}+{\sigma }_{\mathrm{amp}},0\right)$$

(11)

3. The Algorithm

3.1. General Concept

The algorithm starts with a surface mesh (given as an .STL file) of a disc spring placed randomly in space. After aligning the disc spring approximately symmetrically around the y-axis, it imports the geometry information into a finite element application. In our implementation of the algorithm, the finite element application was used. Through the Abaqus Scripting interface [70,71], it builds a model around the geometry data, using additional inputs like loads and friction coefficients provided by the user. It generates an output database file using the Abaqus/Standard Solver. Utilizing Abaqus/Viewer functionalities via the Scripting interface, a field of residual stresses is generated from input data. Aggregated minimum and maximum stress fields are computed by superposing the computed load stress field and the residual stress field. A Walker damage parameter field is computed from the minimum and maximum stress fields. The Walker damage parameter at the surface of the disc spring is computed and exported as tabular data.

3.2. Architecture

The software architecture used is called a pipeline in programming terms and resembles a production line: an instance of the class Spring_stack is passed through a series of different processing stations. The processing stations gradually transform the instance from input data to a solved finite element model with post-processed output fields and from there to easily readable tabular data and a graphical representation thereof. In programming terms, these processing stations are called methods. The flow of data is depicted in Figure 2.

3.3. User-provided Inputs

To conduct the inquiry described above, the algorithm requires several inputs. These inputs are provided by the user. They are comprised of an .STL file containing the triangle surface mesh of a scanned disc spring and one .JSON file for the description of the simulated physical situation and the numerical configuration of the simulation, respectively.

The .JSON file describing the simulated physical situation contains an array of parts defining the geometry, the Walker exponent $\gamma$, Young’s modulus, the Poisson number, and a volumetric mass density (for numerical stability) for each. The geometry of scanned disc springs is given by means of a link to an .STL file. Plasticity may be defined for disc springs in the present implementation. However, this is not in the scope of this paper. Furthermore, the friction coefficients between springs stacked in parallel and springs stacked in series, as well as between springs and the plates and the pillar are committed. The axial load is committed as an array of paths, giving the displacements of one of the plates at the end of each time step. Measured residual stresses in the radial and tangential direction are committed for two points with different coordinates. These two points should be selected carefully by the user and ideally describe a linear model representing more measurements. Furthermore, a target value for the number of triangles of the coarsened triangle surface mesh is committed.

The .JSON file describing the numeric inputs contains several naming definitions, element types, solver options, contact and friction formulations, computing resource allowances, and flags determining whether stresses etc. are to be written to input files. Different modeling conventions for plates and pillars may be used in the current implementation. However, this is outside the scope of this paper.

4. Implementation

4.1. Aligning the Disc Spring

In order to build a model around a part, the alignment of the part in space must be known. For example, to apply axial loads, the axial direction of the loaded part must be known. Because it makes the rest of the algorithm much easier, the input disc spring is preprocessed so it always has the same orientation. By convention, this position is in the point of origin, closely axially symmetric to the y-axis.

In the alignment process, the mesh ${\mathcal{M}}_{\mathrm{start}}$ is processed as a set of vertices ${\mathit{V}}_{\mathrm{start}}=\{{\mathit{v}}_1,\cdots ,{\mathit{v}}_v\}$ and a set of triangles $\mathit{e}=\{{\mathit{e}}_1,\cdots ,{\mathit{e}}_e\}$ connecting these vertices. The vertices ${\mathit{v}}_{i,\mathrm{start}}=(x_{i,\mathrm{start}},y_{i,\mathrm{start}},z_{i,\mathrm{start}})$ are defined in a Cartesian coordinate system $(x,y,z)$. The point set is aligned by a translation $\mathit{t}\left(t_x,t_y,t_z,\mathit{V}\right)$ along and a rotation $\mathit{A}\left({\theta }_x,{\theta }_y,{\theta }_z\right)$ around the principal axes.

$${\mathit{V}}_{\mathrm{aligned}}\left(\mathit{A},\mathit{t},{\mathit{V}}_{\mathrm{start}}\right)=\mathit{A}\left({\theta }_x,{\theta }_y,{\theta }_z\right)·\mathit{t}\left(t_x,t_y,t_z,{\mathit{V}}_{\mathrm{start}}\right)$$

(12)

The edges are defined as connectors between vertices; therefore, a geometric transformation of any mesh $\mathcal{M}$ is fully described by a geometric transformation of its vertices $\mathit{V}$. The resulting mesh ${\mathcal{M}}_{\mathrm{aligned}}$ is defined by the vertices ${\mathit{V}}_{\mathrm{aligned}}$ and the edges $\mathit{e}$.

A feasible combination of translation $\mathit{t}$ and rotation $\mathit{A}$ is found by solving an optimization problem.

$$\underset{{\mathit{\theta }}_{x,y,z},{\mathit{t}}_{x,y,z}}{min}f\left({\mathit{V}}_{\mathrm{aligned}}\left(\mathit{A},\mathit{t},{\mathit{V}}_{\mathrm{start}}\right)\right)$$

(13)

To compute the objective function f, the transformed point cloud ${\mathit{V}}_{\mathrm{aligned}}$ is projected into the $xy$-plane, producing ${\mathit{V}}_{xy}=\left\{{\mathit{v}}_{xy,1},\cdots ,{\mathit{v}}_{xy,v}\right\}$ by rotation around the y-axis. ${\mathit{V}}_{xy}$ is computed using the coordinates of ${\mathit{V}}_{\mathrm{aligned}}$:

$${\mathit{v}}_{xy,i}=\left(\begin{array}{c}\sqrt{x_{i,\mathrm{aligned}}^2+y_{i,\mathrm{aligned}}^2}\\ z_{i,\mathrm{aligned}}\end{array}\right)$$

(14)

The objective function f is defined as the area A of the convex hull of ${\mathit{V}}_{xy}$.

$$f=A\left(Conv\left({\mathit{V}}_{xy}\right)\right)$$

(15)

The convex hull is computed using the Quickhull algorithm [72,73]. The objective function f is invariant to rotations around the y-axis. Therefore, ${\theta }_y$ is fixed to zero. To further simplify the optimization problem, it is assumed that the centroid of the aligned mesh $c\left({\mathcal{M}}_{\mathrm{aligned}}\right)$ lies in the point of origin. Although this assumption is applicable only for perfectly symmetric spring geometries, we assumed that for springs possessing minor deviations regarding their symmetry, the assumption holds true, too.

Based on this assumption, $\mathit{t}$ is defined as translating the volumetric centroid c of the mesh ${\mathcal{M}}_{\mathrm{start}}$ into the point of origin.

$$c\left({\mathcal{M}}_{\mathrm{translated}}\right)=\left(0,0,0\right)$$

(16)

The computational cost for finding $\mathit{t}$ is low because the volumetric centroid of any closed mesh $\mathcal{M}$ can be computed inexpensively.

In the following rotation $\mathit{A}$, the mesh ${\mathcal{M}}_{\mathrm{translated}}$ is only rotated around the principal axes. As the volumetric centroid lies on all three principal axes, it is invariant under said rotation $\mathit{A}$.

$$c\left({\mathcal{M}}_{\mathrm{aligned}}\right)=c\left({\mathcal{M}}_{\mathrm{translated}}\right)=\left(0,0,0\right)$$

(17)

By fixing the translation additionally to the rotation around the y-axis, the dimensionality of the optimization problem is reduced from six to two.

$$\underset{{\mathit{\theta }}_{x,z}}{min}f\left({\mathit{V}}_{\mathrm{aligned}}\left(\mathit{A},\mathit{t},{\mathit{V}}_{\mathrm{start}}\right)\right)$$

(18)

The parameters ${\theta }_x$ and ${\theta }_z$ are calculated using the L-BFGS-B method [73,74,75].

The resulting mesh ${\mathcal{M}}_{\mathrm{aligned}}$ may be aligned upside down. This being the case, it can simply be rotated by 180 degrees around the x-axis. To further ease the following steps, it is translated along the y-axis, so its lowest point is in the $xz$-plane.

$$min\left(y_{i,\mathrm{postprocessed}}\right)=0$$

(19)

The postprocessed, aligned mesh ${\mathcal{M}}_{\mathrm{postprocessed}}$ is exported as an Abaqus input file (.inp).

4.2. Building the Finite Element Model

The creation of the model is structured as a pipeline inside the pipeline depicted in Figure 2. The sub-pipeline is depicted in Figure 3. The allotment of tasks to processing stations (methods) along the pipeline is based on the allotment of functionalities to modules in the Abaqus/CAE environment. Each of the processing stations is customized according to input data. The following paragraphs describe the individual processing stations.

For the simulation of a single disc spring, four components are required: the disc spring itself, a guide pin, and two flat plates. The geometry of the disc spring is read from the input file generated from the 3D-scan data. The guide pin is defined as an idealized cylinder with the same height as the disc spring, and the plates are defined as idealized planes.

The plate and the guide pin part instances are defined as analytic rigid surfaces and therefore do not need to be meshed. The part instance representing the disc spring is imported as a triangle surface mesh and is converted into a volume tetrahedron mesh by the free mesher implemented in Abaqus/CAE. Mesh size cannot be controlled actively, but only by changing the mesh size of the imported surface mesh. The surface mesh size is adjusted by exporting very fine surface meshes from the 3D-scanning software GOM Scan and increasing the mesh size using the quadric edge collapse decimation algorithm [76]. In this application, quadric edge collapse is especially suitable because it reduces the mesh size at the flat surfaces, where the mesh may be more coarse, and keeps it nearly constant at the edges, where a fine mesh is needed.

In the approach described in this paper, the finite element analysis utilizes a purely elastic material law. Residual stresses are incorporated by superposition of the measured stresses committed by the user. This is implemented by building two axisymmetric field outputs, one for ${\sigma }_{\mathrm{residual},\tan }$ and one for ${\sigma }_{\mathrm{residual},\mathrm{rad}}$, in Abaqus/CAE. These output fields will later be used for the computation of ${\sigma }_{\mathrm{I},\min }$, ${\sigma }_{\mathrm{I},\max }$, ${\sigma }_{\mathrm{II},\min }$, and ${\sigma }_{\mathrm{II},\max }$ according to Equations (24) to (27).

To introduce inertia and improve convergence, a volumetric mass density is defined for the disc spring. Because the spring was aligned beforehand, the assembly is straightforward. All four parts are joined in an assembly with a single coordinate system. The lower plate, the disc spring, and the guide pin are already in place. The upper plate is positioned in its place by a translation along the y-axis.

The lower plate and the guide pin are fixed in space, and loads are applied to the upper plate. To improve convergence, a small gravitational load is defined. No other loads are applied directly to the disc. Forces are transmitted via contacts. Surface to surface contacts are defined between the disc spring and each of the other part instances.

The Abaqus Scripting interface includes an option to customize output requests. This option is used in the algorithm to request stress and coordinate outputs. However, coordinates cannot be requested at integration points by means of the Scripting interface. Coordinate data in the integration points are necessary to compute residual stresses in the integration points.

A job instance is created to generate an input file. The input file is manipulated to create an output request for coordinate information in the integration points, and a second job instance pointing to the new input file is created. This way, a job is created that is identical to the first job except that it includes an output request for coordinate information in the integration points. The second job is converted into a system of partial differential equations and solved using Abaqus/Standard.

4.3. Computing Stresses and Walker Damage Parameters

After solving the system of partial differential equations for the displacements of the nodes in the model, Abaqus/Standard computes stresses and coordinates in the integration points at different points in step-time, which refer to the minimum and the maximum load. In this section, the computation of Walker damage parameters at the surface of the disc spring between Edge II and Edge III (see Figure 4) is described. Here, we may assume a plane stress state with a dominating tension component:

$${\sigma }_{\mathrm{load},\mathrm{II},\min }=0\vee {\sigma }_{\mathrm{load},\mathrm{III},\min }=0$$

(20)

$${\sigma }_{\mathrm{load},\mathrm{II},\max }=0\vee {\sigma }_{\mathrm{load},\mathrm{III},\max }=0$$

(21)

For integration points close to the surface, this is approximately true. The assumption becomes more realistic with a finer mesh and is true with an infinitesimal mesh size. The minimal and maximal stress components ${\sigma }_{\mathrm{load},\mathrm{II},\min }^{\prime }$ and ${\sigma }_{\mathrm{load},\mathrm{II},\max }^{\prime }$ are redefined:

$${\sigma }_{\mathrm{load},\mathrm{II},\min }^{\prime }={\sigma }_{\mathrm{load},\mathrm{II},\min }+{\sigma }_{\mathrm{load},\mathrm{III},\min }$$

(22)

$${\sigma }_{\mathrm{load},\mathrm{II},\max }^{\prime }={\sigma }_{\mathrm{load},\mathrm{II},\max }+{\sigma }_{\mathrm{load},\mathrm{III},\max }$$

(23)

In the surface of disc springs, ${\sigma }_{\mathrm{load},\mathrm{I}}$ points in the tangential direction and ${\sigma }_{\mathrm{load},\mathrm{II}}^{\prime }$ points in the radial direction. Residual stresses are measured in the tangential direction (${\sigma }_{\mathrm{residual},\tan }$) and in the radial direction (${\sigma }_{\mathrm{residual},\mathrm{rad}}$). The input stresses for the computation of the Walker damage parameters are defined as:

$${\sigma }_{\mathrm{I},\max }={\sigma }_{\mathrm{load},\mathrm{I},\max }+{\sigma }_{\mathrm{residual},\tan }$$

(24)

$${\sigma }_{\mathrm{II},\max }={\sigma }_{\mathrm{load},\mathrm{II},\max }^{\prime }+{\sigma }_{\mathrm{residual},\mathrm{rad}}$$

(25)

$${\sigma }_{\mathrm{I},\min }={\sigma }_{\mathrm{load},\mathrm{I},\min }+{\sigma }_{\mathrm{residual},\tan }$$

(26)

$${\sigma }_{\mathrm{II},\min }={\sigma }_{\mathrm{load},\mathrm{II},\min }^{\prime }+{\sigma }_{\mathrm{residual},\mathrm{rad}}$$

(27)

Walker damage parameters are computed according to Equations (1) to (11). For each surface triangle, a surface area $A_i$ and a Walker damage parameter $P_{\mathrm{Walker},i}$ are computed. The surface Walker damage parameter is computed by averaging the Walker damage parameters of the neighboring nodes. These data pairs are saved as the set S.

4.4. Extracting Tabular Data

The set S is too rich in information to be captured holistically without extraction and/or condensation. This is done by transforming the information into tabular data. Therefore, it is condensed into human readable tabular data. Each element of the set is defined as a Dirac delta function:

$$f_i\left(P_{\mathrm{Walker}}\right)=A_i·\delta \left(P_{\mathrm{Walker}}-P_{\mathrm{Walker},i}\right)$$

(28)

The accumulated surface area $A_{\mathrm{acc}}$ assigned to a given Walker damage parameter is computed as:

$$A_{\mathrm{acc}}\left(P_{\mathrm{Walker}}\right)={\int }_{P_{\mathrm{Walker}}}^{\infty }\sum _i{f}_i\left(P_{\mathrm{Walker}}\right)\mathrm{d}{P}_{\mathrm{Walker}}$$

(29)

It describes the size of the surface area of the disc spring with a Walker damage equal to or greater than the given $P_{\mathrm{Walker}}$. A graphical representation [77] of the accumulated surface area is created.

5. Example Application

In this section, an example application of the algorithm presented in Section 3 is given. A single disc spring was modelled. The surface mesh of the disc spring under investigation is presented in Figure 5. The mesh was obtained using the commercially available 3D-scanning device GOM ATOS and the software GOM Scan. The number of surface triangles was already reduced for the displayed mesh, from 253,904 to 50,152. Convergence studies with different loads showed this mesh to be a good compromise between computational cost (about seven hours of CPU time on a i7-9800X and 16 GB of RAM vs. about 110 h of CPU time and 115 GB of RAM for the 253,904 surface triangle model) and accuracy (no significant bias of the Walker damage-surface area plot) [4].

Additionally, an idealized geometry derived from the scanned data is presented in Figure 4. Edges I to IV are labelled for the reader’s orientation. A major difference between both geometries is that the surfaces between Edges I and IV, as well as II and III of the idealized geometry are straight, while those of the 3D-scanned disc spring are curved. This is also visible in Figure 6. Of course, the 3D-scanned geometry also was not perfectly symmetric.

In Figure 7, the geometry according to the standard [1] is presented. It differs from the idealized geometry in having sharp edges and all angles between faces being 90°. This simplified geometry is usually utilized for the analysis of disc springs, regardless of whether analytical formulas or finite element models are used.

The triangle surface mesh was aligned and imported into Abaqus/CAE. Afterwards, it was converted into a tetrahedron volume mesh using the generateMesh method included in Abaqus Scripting [71]. The resulting linear tetrahedron C3D4 mesh was converted into a quadratic C3D10 mesh. Modified quadratic tetrahedrons C3D10M offer improved contact behaviour [78]; however, due to poor mesh quality, C3D10 tetrahedrons performed better in our experience. The elements were assigned a Young’s modulus of 206,000 MPa, a Poisson ratio of 0.3, and a volumetric mass density of 8.05 g/mm${}^3$.

The volumetrically meshed disc spring was incorporated into an assembly. The load cycle was implemented in four steps. To obtain good convergence behavior, the steps were defined as dynamic steps. The implicit solver Abaqus/Standard was used. The purpose of Step 1 was to apply the lower load. Step 2 was to make sure there was no dynamic influence on the computed stresses. Step 3 was to apply the higher load. Step 4 was, again, implemented to eliminate dynamic influences. Step times for Steps 1 to 4 were 10 s, 1 s, 10 s. and 1 s. The boundary conditions were applied to the upper plate as displacements in the axial direction at a reference point; see Figure 6. The prescribed displacements were 0.425 mm for Steps 1 and 2 and 1.19 mm for Steps 3 and 4. All other degrees of freedom of the reference points were fixed to zero. For the lower plate, all degrees of freedom were fixed to zero.

All contact formulations used in this model were defined as surface-to-surface contacts with a finite sliding penalty formulation and a friction coefficient of 0.01.

A job was created and committed to the solver Abaqus/Standard. The resulting output database was loaded. The resulting highest tensile load stresses in a cross-section are displayed in Figure 8. The upper half of the disc spring was loaded compressively. This is why disc springs in general break between Edges II and III.

Two field outputs representing residual stresses in the tangential and radial direction were generated based on user input and the initial coordinates of the integration points (which were requested as outputs earlier). Since the measured residual stresses for the spring under investigation are confidential, the used inputs values were not the result of a measurement. They were however realistic for disc springs. The residual stresses between Edges II and III were approximated by a linear function, ignoring non-symmetric effects. The computed residual stresses were obviously wrong anywhere else. This does not matter here because the Walker damage parameter is a measure used to predict fracture. Fracture is initiated by cracks, which normally initiate from the surface between Edges II and III [79]. Cracks originating from between Edges I and IV are usually caused by too low preloading forces. The bias in the Walker damage parameter introduced outside the surface between Edges II and III is non-conservative. Therefore a false positive for fracture in this area can be ruled out. The computed output fields between Edges II and III are displayed in Figure 9. The significantly higher residual stresses in the tangential direction compared to the radial direction are normal in disc springs because disc springs are overloaded in production to prevent plastic deformation in use, to increase lifetime, and to decrease creep effects [80,81,82,83,84,85].

The output fields resulting from the finite element simulation describing stresses after Step 2 and the output fields describing the residual stresses were added up to compute output fields describing ${\sigma }_{\mathrm{I},\max }$ and ${\sigma }_{\mathrm{II},\max }$. The calculation followed Equations (24) and (25). Field outputs describing ${\sigma }_{\mathrm{I},\min }$ and ${\sigma }_{\mathrm{II},\min }$ were computed following Equations (26) and (27).

Based on these, field outputs describing the minimum and maximum equivalent stress, ${\sigma }_{\min }$ and ${\sigma }_{\max }$, were computed according to Equations (2) and (3); see Figure 10. Especially in the minimum equivalent von Mises stress visualization, the contact line between the disc spring and the lower plate can be identified by a circle of locations with high compressive stresses.

These output fields were used to compute the output fields representing the stress amplitude ${\sigma }_{\mathrm{amp}}$ according to Equation (4) and finally the Walker damage parameter $P_{\mathrm{Walker}}$ according to Equation (1); see Figure 11. The Walker exponent $\gamma =0.5$ was used, which makes the Walker damage parameter equivalent to the Smith–Watson–Topper damage parameter. Compared to the other fields, the stress amplitude field was very smooth. The reason for this is that the elastic deformation of the spring partially compensated for the small asymmetries that were present. For higher load increments, the additional elastic stresses were therefore distributed more homogeneously.

The algorithm created a list of all surface triangles and computed an average Walker damage parameter $P_{\mathrm{Walker},i}$, as well as a surface area $A_i$ for each triangle. The graph of the accumulated surface area over the Walker damage parameter was created; see Figure 12. From this particular graph, the user can for example extract the surface area where the Walker damage parameter is over 750 MPa, which is 4.9 mm${}^2$. About half of the surface had a Walker damage parameter of zero because stresses there were purely compressive; see Figure 12. This surface corresponds to the dark blue parts of the surface in Figure 8. As can be seen on the graph on the right, the function starts to jump at high stresses. This is because Walker damage parameters were averaged over surface triangles. The part of the graph at very high stresses exists purely because of numerical singularities and therefore does not correspond to the fatigue behaviour of the disc spring.

The Walker damage parameter is not directly accessible through experiments. To evaluate the quality of the finite element model, a characteristic derived from a similar model (only boundary conditions were changed) was compared to a characteristic obtained in an experiment; see Figure 13. As is customary for disc springs, the deflection was normalized over the deflection at which the disc spring lies flat on the ground. They agree well; especially, the correct representation of the progressive behaviour at the very start of the experiment has only been achieved by models directly implementing 3D-scanned geometry. To our knowledge, all published models directly implementing 3D-scanned geometry have been created using the algorithm presented in Section 3. The numerical characteristic is somewhat stiffer than the experimental one. This may be due to a misrepresentation of Young’s modulus, which was set to the normative default of 206,000 MPa; however, tensile tests on specimens from the same batch of material and a similar heat treatment did not show a sufficient deviation in Young’s modulus to use a lower value.

6. Summary

A novel algorithm as implemented in the Spring_stack module was presented in this paper. The algorithm receives geometry data, residual stresses, material parameters, load cases, and further inputs and builds an FE model based on these. It evaluates the FE model after solving with respect to the Walker damage inflicted locally using a Manson–McKnight approach. In a post-processing step, the accumulated surface area as a function of the Walker damage parameter is computed. An example application of the algorithm was presented.