Next Article in Journal
The Status of On-Board Hydrogen Storage in Fuel Cell Electric Vehicles
Next Article in Special Issue
Comparative Analysis of Various Hyperelastic Models and Element Types for Finite Element Analysis
Previous Article in Journal
A Laminated Spherical Tsunami Shelter with an Elastic Buffer Layer and Its Integrated Bulge Processing Method
Previous Article in Special Issue
Binary Integer Formulations for Task Allocation and Optimal Labor Cost in Small and Medium Apparel Manufacturing
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Unit-Load Approach for Reliability-Based Design Optimization of Linear Structures under Random Loads and Boundary Conditions

by
Robert James Haupin
and
Gene Jean-Win Hou
*
Department of Mechanical and Aerospace Engineering, Old Dominion University, Norfolk, VA 23454, USA
*
Author to whom correspondence should be addressed.
Designs 2023, 7(4), 96; https://doi.org/10.3390/designs7040096
Submission received: 19 May 2023 / Revised: 23 June 2023 / Accepted: 11 July 2023 / Published: 2 August 2023
(This article belongs to the Special Issue Design Sensitivity Analysis and Engineering Optimization)

Abstract

:
The low order Taylor’s series expansion was employed in this study to estimate the reliability indices of the failure criteria for reliability-based design optimization of a linear static structure subjected to random loads and boundary conditions. By taking the advantage of the linear superposition principle, only a few analyses of the structure subjected to unit-loads are needed through the entire optimization process to produce acceptable results. Two structural examples are presented in this study to illustrate the effectiveness of the proposed approach for reliability-based design optimization: one deals with a truss structure subjected to random multiple point constraints, and the other conducts shape design optimization of a plane stress problem subjected to random point loads. Both examples were formulated and solved by the finite element method. The first example used the penalty method to reformulate the multiple point constraints as external loads, while the second example introduced an approach to propagate the uncertainty linearly from the nodal displacement vector to the nodal von Mises stress vector. The final designs obtained from the reliability-based design optimization were validated through Monte Carlo simulation. This validation process was completed with only four unit-load analyses for the first example and two for the second example.

1. Introduction

Various uncertainties are inherent in engineering design processes with quantities such as tolerancing, material and mechanical properties, applied loads, and environmental conditions. To count these uncertainties in modeling and analysis, it is necessary to produce results with high confidence. For this reason, stochastic finite element methods have recently attracted attention in research and applications. Stochastic finite element methods can be broadly divided into three categories: sampling methods, the spectral stochastic finite element methods, and the perturbation finite element methods [1,2,3,4]. These methods have been extended to robust design optimization, fatigue failure, and structural health monitoring [5,6,7,8,9,10,11].
Sampling methods are non-inclusive, which include Monte Carlo simulation and smart sampling techniques [12]. In particular, Georgioudakis et al. [9] applied smart sampling techniques to compute the probability of fatigue failure caused by crack initiation and propagation. Zheng et al. [10,11] proposed an adjustable hybrid resampling technique to effectively identify the critical area to place sensors to monitor the health of a bridge structure.
On the other hand, the spectral and perturbation stochastic finite element methods are inclusive. The first step in these methods is to expand the input random fields and processes as a combination of random functions and deterministic base functions. The Karhunen–Loeve expansion is an example of such a process, which decouples a random field or a random process with a known distribution into stochastic functions and spatial or time-dependent functions [13,14,15,16].
In the spectral finite element method, the stochastic response function, which is a function of random variables, is expanded as a linear combination of orthogonal base functions of the random variables and the deterministic coefficients. The coefficients in the function expansion are integrals of the stochastic response function and each of the base functions. This integration can be time consuming, if the number of random variables involved in the problem increases. To alleviate this difficulty, polynomial dimension decomposition is introduced to sequentially decompose the stochastic function into component functions with an increasing number of random variables [6,17]. These component functions can be further expended by the orthonormal functions with respect to the given probability density function. The deterministic coefficients in the expanded component functions are again evaluated through the weighted integration. However, the number of required integrations in this case is much less than the original expansion process. Recently, Song et al. [18] employed polynomial dimension decomposition directly to integrate the coefficients in the function expansion. In this case, the polynomial dimension decomposition was introduced to express the integrand in the integral in terms of the component functions with an increasing number of random variables. In short, Rahman [16] applied polynomial dimension decomposition directly to the original stochastic response function while Song et al. [18] applied it to the product of the stochastic response function and the orthonormal basis functions introduced in the original expansion. Both approaches achieved the computational efficiency for problems with multiple variables compared with the original expansion process. Kaintura et al. [17] offered the tensor reduction method for further dimension reduction. The dimension decomposition method has been used to support reliability-based design optimization [6].
Taylor’s series expansion in terms of random variables has been used for a while by practicing engineers to set up the standard for design safety [19]. In the perturbation finite element methods, Taylor’s series expansion is used to expand the finite element matrix equation into a polynomial of random variables. The coefficients of this expansion are the derivatives with respect to the random variables of any component in the matrix equation including the stiffness matrix, the mass matrix, the displacement vector, and the force vector. These derivatives are evaluated at the mean values of the random variables. The perturbation finite element methods are relatively easier to implement than the spectral finite element methods as the former can be conveniently extended to broad applications such as vibration [20,21], heat transfer [22], and dynamic analysis [23].
Most applications of the perturbation stochastic finite element method (PSFEM) are limited themselves up to second-order derivatives, which produce acceptable results for cases with covariances of the random variables of less than 0.15. Kaminski [24] extended PSFEM to include terms with an order of derivatives higher than two. As expected, more accurate results were achieved. Wu et al. [25] replaced analytical derivatives in PSFEM by finite differencing to ease the complexity in computation. Kaminski [26] presented high-order Taylor’s series expansion of strain and stress in terms of random material properties. However, the reported results of most of the stochastic finite element methods are nodal displacement, natural frequency, and the critical number of cycles in fatigue analysis. Only a very few studies have reported the nodal von Mises stresses as output. This could be because the von Mises is a nonlinear function of random variables. The randomness of the Young’s modulus and the Poisson’s ratio will further increase such nonlinearity [27]. Furthermore, the state of stresses in the hp finite element method is discontinuous across the boundary between a pair of adjacent elements. It is necessary to use a stress recovery scheme to smooth the stress field globally in the whole domain to find the nodal von Mises stress [14,28] or to smooth the stress field locally in the sub-element surrounding the node of concern [29]. Rahman and Rao [30] did report the statistics of the state stresses at nodes as the output, which was conducted by using the meshless finite element method.
As previously mentioned, the input uncertainties of most of the works carried out by stochastic finite element methods including both spectral and perturbation methods are material properties and domain geometries. Only a few have dealt with loading and the boundary conditions as input uncertainty [7,22]. Zhao et al. [7] conducted robust topology optimization of a linear static structure under loading uncertainty. The random loads in [7] were represented by a linear combination of the random amplitudes and their associated deterministic load vectors. The random amplitudes alone could model the uncertainties in both the magnitude and direction of the applied random load. The analysis and sensitivity analysis of the matrix state equation in the paper were then computed as a combination of the linear structural analyses, each of which was subjected to one of the deterministic load vectors. The objective and the constraint functions in the optimization formulation [7] were expressed in terms of the second-order and the fourth-order products of random variables. Their means and variances of products of random variables were evaluated by Monte Carlo simulation, which does not involve any structural analysis. Xiu and Karniadakis [22] applied both spectral and perturbation finite element methods to study the solution statistics of the viscous Burgers’ equation subjected to uncertainty in the boundary conditions. They concluded that the perturbation finite element method was unable to produce an accurate solution with large variations.
The goal of this study was to conduct reliability index-based design optimization of a linear elastic structure subjected to random loads and boundary conditions. The constraints of concern were formulated in terms of the associated reliability indices. In most reliability design optimization, the constraint reliability indices are calculated based upon the reliability index approach or performance measurement approach [31,32,33,34]. However, in this study, the reliability indices were directly formulated in terms of the associated means and standard deviations, which were computed based upon the low-order Taylor’s series expansion. A unit-load approach was employed in this study throughout the entire design optimization process to achieve computational efficiency by taking advantage of the linear nature of the structure. The unit-load approach referred to in this study, with a given stiffness matrix K, conducted all related analyses with the unit-load equation, K U i = I i , where the load vector I i is all zero except for its ith component, which was set as 1.
The Materials and Methods section presents a detailed process to formulate a reliability-based design optimization problem, which includes a general finite element equation for a linear static structure subjected to random load and multipoint boundary conditions, and one-side and interval probabilities to model the uncertainties associated with the equality and inequality constraints. The constraint equations of concern comprised the nodal displacement vector, the nodal reactions, and the nodal von Mises stresses. Among these three responses, the nodal displacement vector is the direct solution of the finite element equation, which leads to the computation of the nodal reaction and nodal von Mises stresses. The first- and second-order derivatives of nodal displacements, reactions, and von Mises stresses with respect to random loads and boundary conditions are also presented in this section, which are required to compute the reliability indices of the concerned constraints through Taylor’s series expansion.
The Results section presents two examples to demonstrate the use of the means and standard deviations, which were approximated by Taylor’s series expansion, to compute the reliability indices of the constraints for reliability-based design optimization. The means and standard deviations of the constraints calculated in these examples were validated with those obtained by Monte Carlo simulation. The final optimization solutions of the reliability-based design examples were also be validated by Monte Carlo simulation. Our concluding remarks are summarized in the final section of the discussion.
The first example was to minimize the weight of a two-bar truss structure while allowing the randomness of the nodal displacement at the free end to fall into an allowable bound. Besides two standard single point constraints, the state equation of the truss example was subjected to two multiple-point constraints. The coefficients of the multiple point constraints were set to be random. The design problem had six design variables in total, which were all deterministic. Two of the six design variables were related to the cross-sectional areas of the truss members and the rest were the standard deviations of four random variables related to the coefficients of two multiple point constraints. These standard deviations were treated as the design variables to control the level of the randomness of the design. The penalty method was applied here to explicitly incorporate the multiple-point constraints into the state equation of the free truss structure. The augmented state equation is now equipped with a random stiffness matrix and a random external force, which can be solved by four unit-load solutions for a given set of design variables.
The second example aimed to find the best shape and location of a rectangular cutoff at the base of a mooring bracket to minimize the weight. The mooring bracket was modeled by a plane stress problem subjected to a point load described by two random variables. The domain of the problem was discretized into four-node quadrilateral elements. The constraint set was made of nodal displacements, reactions, and von Mises stresses. It should be noted that the nodal displacements and reactions were the direct output of the state equation, but not the nodal von Mises stresses. Moreover, the nodal von Mises stresses were nonlinear in terms of the nodal displacements. Traditionally, the latter are computed in three steps in the finite element methods. The first step is to compute the state of the stresses at the interior integration points within each element in terms of the nodal displacements of the selected element. The von Mises stresses at the same integration points can then be computed in the second step within each element. In the third step, a global stress recovery scheme is introduced to extrapolate the von Mises stresses from the integration points within each element to the nodal points of the entire domain. In this study, the sequence of the last two steps was switched to extend the linearity of uncertainty propagation from the state of stresses at the integration points within each element to nodal points of the entire domain. The latter could then be used to compute the von Mises stresses at the nodal points directly.
In summary, the unit-load approach is the key element in this study, which enables the reliability indices of any performance function to be computed efficiently. The first example used the penalty method to convert the random multiple boundary conditions imposed upon the truss structure to the random load. This conversion allows the reliability index of any response function of the given truss structure to be computed by the unit-load approach. The second example applied a new approach to propagate the uncertainty smoothly from the nodal displacements that resulted from the unit-load analysis to the nodal von Mises, which is a nonlinear function of nodal displacements. Both examples demonstrated the computational efficiency of using the unit-load approach for reliability-based design optimization.

2. Materials and Methods

A 3-step formulation process is presented in this section to propagate the uncertainties in the loads and boundary conditions from the finite element state equation to the reliability-based design optimization problem. Each step is described sequentially as a subsection as follows. The first step is to use the penalty method to construct a state equation that includes the randomness explicitly as part of the equation. The second subsection focuses on finding the first- and second-order derivatives of the structural responses with respect to the input random variables. These derivatives enable the use of the Taylor’s series expansion to estimate the means and standard deviations of the response functions of concern. These estimated means and standard deviations can then be directly applied to calculate the reliability indices of the response constraints. With the availability of these reliability indices, the third subsection builds a design optimization problem with uncertainty.

2.1. Linear Static Equation with Random Loads and Boundary Conditions

The static finite element equation of a free structure is represented as
K q = f
where the n × 1   vectors, f and q , are the global force and the displacement, respectively, while K is the n × n   global stiffness matrix, which is non-invertible at this moment. To ensure that it produces a unique solution, Equation (1) must include proper boundary conditions. The finite element equation is then modified as
K q = f + R
where the boundary reaction force, R , is installed to enforce the boundary conditions upon the involved boundary degrees of freedom, q B . Two types of linear boundary conditions, called the single point constraint and the multiple point constraint, are commonly considered in finite element analysis. Assume that the kth boundary conditions in analysis, C k , is a typical linear, two-point constraint that involves two degrees of freedom, q k p and q k p , which can be expressed as
C k = β k p q k p + β k q q k q β k 0 = 0
where the user-provided constants β k p , β k q , and β k 0 are treated as random variables in this study. A single point constraint is a special case of a multi-point constraint, Equation (3), in which the coefficient, q k p , can be set to zero. Note that this multi-point constraint in Equation (3) can be put in a vector form as
C k = v k T q β k 0 = 0
where v k is a n × 1   vector with two non-zero values, β k q and β k q , located respectively at the degrees of freedom, p and q. Equation (4) can be extended to count for all boundary constraints in the analysis as
C = V T q β 0 = 0
where C is an m × 1 vector and V is a n × m matrix with m being the total number of constraints. Note that the kth column of V is the vector, v k , in Equation (4), with two non-zero coefficients, β k p and β k q .
The penalty method can be used to solve Equation (1) constrained by Equation (5). This is carried out by introducing a large penalty coefficient, μ , to transfer Equation (1) into a new form
K + μ V V T q = f + μ V β 0
Comparing Equation (6) with Equation (2), the boundary condition of Equation (5) has now been included as part of the forcing term. Specifically, the reaction force, R, in Equation (2) can now be identified as
R = μ V C = μ V V T q β 0 = μ V V T q + μ V β 0
The reaction force in this case enforces the constrained degrees of freedom to satisfy the required multiple point constraint, Equation (6). Note that the penalty method does not produce a perfect solution that fully accomplishes C = 0 . However, the error in C in the penalty method facilitates the calculation of the reaction force.
As for computing the nodal von Mises stresses, it is not as straight forward as computing the nodal displacements and reactions. This is because the stresses are not continuous across the elements in the finite element method. In a deterministic finite element analysis process, the nodal von Mises stresses are computed in three steps. In the first step, the state of the stresses is computed at the interior integration points within an element. Once the state of the stresses becomes available, the von Mises stress at the same location can be calculated. This is followed by the use of a smoothing method to obtain the nodal von Mises stresses by curve-fitting the von Mises stresses calculated at the interior points throughout the entire domain.
Specifically, the state of the stresses at every location in an element can be calculated based upon the known nodal displacement field q . For example, one can relate the state of the stresses of a plane stress problem at an interior integration point i to the displacement field in a specified finite element, e, as
σ e = σ x σ y τ x y = D B ξ i , η i q e
where matrix D is a function of the material properties and matrix B is composed of the derivatives of shape function evaluated at ξ i , η i , the natural coordinates of the location of the interior integration point i [35]. Note that the natural coordinates, ξ ,   η , are non-dimensional and are independent of the shape of a finite element. Once the state of the stresses at the integration points becomes available, the von Mises stress can be calculated at the same interior integration point as
σ V M e 2 = σ x 2 σ x σ y + σ y 2 + 3 τ x y 2 = σ e T A σ e
where the symmetric matrix A is given by
A = 1 1 / 2 0 1 / 2 1 0 0 0 3
This is followed by employing a smoothing process to map the von Mises from the interior integration points to the elemental nodal points. The least-squares fit method suggested by Chandrupatla and Belegundu [35] was employed in this study for smoothing. Like the displacement field in an element, the von Mises stress field in an individual element is also approximated by a linear combination of shape functions and nodal von Mises stresses, Σ e , as
σ V M e ξ , η = N e T ξ , η Σ e
The best-fitting nodal von Mises stresses are those that minimize the difference, collected through the whole domain, between the von Mises stresses evaluated at the given interior integration points and those evaluated by the approximation, Equation (10). The process results in a matrix equation,
M Σ = σ V M
where the vector, Σ , is unknown, which is the collection of the von Mises stresses at the nodal points, the vector σ V M is known, which are those computed at the interior integration points, and the matrix M , is the assembly of element shape functions evaluated at the interior integration points over the whole domain as
M = e N e ξ i , η i N e T ξ i , η i
Note that the matrix, M , is independent of the element shape or element material properties for an iso-parametric element.
Once the input load vector, f ,   and the boundary conditions, C = 0 , are considered random, the output structural responses in the finite element analysis, the nodal displacements, reactions, and von Mises stresses will also become random. The challenge now is to quantify the uncertainties of these output responses with measurable terms. To this end, the low-order Taylor’s series expansion is employed to approximate the mean and the variance of any response function, g x i , for a given set of random variables, x i ,   as
E g = μ g g μ x 1 , μ x 2 , , μ x n + 1 2 i = 1 n g , x i s x i 2
V a r g = s g 2 = i = 1 n g , x i s x i 2
where μ x i and s x i are the mean and the standard deviation of a random variable x i and g , x i   a n d   g , x i   are the first- and the second-order derivatives of the function, g x i .   The latter are all evaluated at means, μ x i . Assume that the random variables are statistically independent to each other, and that the response functions have sufficient regularity. The above approximations provide acceptable accuracy for a covariance of less than 0.15.
To take advantage of Equations (12) and (13), it is necessary to develop a method that can efficiently evaluate the first- and second-order derivatives of the response functions resulting from the finite element analysis. Specifically, in this study, we needed to compute the first- and second-order derivatives of the structural responses, nodal displacements, reactions, and von Mises stresses with respect to the random variables in f , V , and β 0 . The following subsection demonstrates that these derivatives can be found with a few unit-load solutions. Note that it is straight forward to derive different orders of derivatives of the nodal displacements and reactions, but not the nodal von Mises stresses, as the latter involves additional uncertainty propagation.

2.2. Derivatives with Respect to Random Loads and Boundary Conditions

This section provides a detailed derivation to find the high-order derivatives of the state equation and its output of a linear structure with respect to the random external forces and the random multiple point constraints. Note that the solution q of the state equation, Equation (6), is linear in terms of the load vector, f and β 0 . Therefore, the order of the derivatives of Equation (6) higher than one with respect to any component in f and β 0   will be zero. However, deriving the design derivatives of Equation (6) with respect to the rest of coefficients, β k p and β k q , in the k t h multipoint constraint will be more involved. This is because matrix V on the right of Equation (6) and V V T on the left are the linear and nonlinear functions of β k p and β k q , respectively. Consequently, the displacement q is now nonlinear in β k p and β k q .
To start the derivation, set the vector U i , being the solution of the linear structure under the unit-load, I i , as
K + μ V V T U i = I i
where the subscript i specifically denotes the non-zero degree of freedom, the vector, I i , has only one non-zero value of 1 at the ith component, and V   defined in Equation (5) is a function of β k p and β k q . Note that U i is in fact the derivative of q with respect to the ith component of the external force vector f i . That is, U i = q , f i . For simplicity, the solutions of the finite element matrix equation of Equation (14) due to any unit-load imposed at different degrees of freedom such as U i , U k p , or U k q will be referred to as a unit-load displacement in the rest of this paper.
One can now differentiate Equation (6) with respect to β k 0 , a variable on the right-hand side of the kth constraint, to obtain
K + μ V V T q , β k 0 = μ v k = μ β k p I k p + μ β k q I k q
or
q , β k 0 = μ β k p U k p + μ β k q U k q
where U k p is the solution to Equation (14) with a unit-load placed at the degree of freedom, kp, while U k q is associated with a unit-load at the degree of freedom, kq.
The next step is to differentiate Equation (6) with respect to random coefficients in the constraint equation, β k p and β k q . For simplicity, only β k p is considered as an independent variable in the following derivation. The first-order derivative of Equation (6) is now obtained as
K + μ V V T q , β k p = μ V V T , β k q q + μ V β 0 , β k q
where the matrix product of the vector, V , which is a collection of v k , defined in Equation (4)
V V T = k = 1 m v k v k T
and the first- and second-order terms on the right-side of Equation (16) can be expressed in detail as follows:
V V T , β k q q = v k , β k p v k T + v k v , β k p T q = 2 β k p q k p + β k q q k q I k p β k q q k p I k q
and
V β 0 , β k q = β k 0 I k p
With the help of the above equations, Equation (16) can be reorganized as
K + μ V V T q , β k p = μ 2 β k p q k p + β k q q k q + β k 0 I k p μ β k q q k p I k q
The solution of the above equation, which is the first-order displacement derivative with respect to β k p , is now obtained as a linear combination of two unit-load solutions as
q , β k p = μ 2 β k p q k p + β k q q k q + β k 0 U k p μ β k q q k p U k q
The second-order displacement derivative with respect to β k p can be carried out in a similar matter, which again yields a linear combination of the same two unit-load solutions as
q , β k p = 2 μ 2 β k p q k p , β k p + β k q q k q , β k p + q k p U k p 2 μ β k q q k p , β k p U k q
The same process can be further extended to obtain the higher-order derivatives of q as
q , β k p ( n ) = n μ 2 β k p q k p , β k p n 1 + β k q q k q , β k p n 1 + n 1 q k p , β k p n 2 U k p n μ β k q q k p , β k p n 1 U k q
where the order of differentiation, n, is higher than 2. In short, the derivatives of the displacement vector with respect to β k q in different orders can be obtained in terms of two unit-load displacements, U k p and U k q .
Next, the derivatives of the reaction force vector, R, can be found by directly differentiating Equation (7). For example, the derivative of R with respect to a force component, f i , will be linear in q , f i  as
R , f i = μ V V T q , f i = μ V V T U i
while the derivative of R with respect to β k 0   is found to be
R , β k 0 = μ V V T q , β k 0 + μ v k
The derivatives of R with respect to the coefficients of the constraint, β k p and β k q , can also be found by the direct differentiation of Equation (7). Again, for illustration, one may select β k p as an independent variable to obtain the following relations,
R , β k p = μ 2 β k p q k p + β k q q k q + β k 0 I k p μ β k q q k p I k q μ V V T q , β k p
and
R , β k p = 2 μ 2 β k p q k p , β k p + β k q q k q , β k p + q k p I k p 2 μ β k q q k p , β k p I k q μ V V T q , β k p
Furthermore, the high-order derivatives of the reaction force can also be found as
R , β k p n = n μ 2 β k p q k p , β k p n 1 + β k q q k q , β k p n 1 + n 1 q k p , β k p n 2 I k p n μ β k q q k p , β k p n 1 I k q μ V V T q , β k p n
Note that the first two coefficients of I k p and I k q on the right-hand side of R , β k p are the same as those of U k p and U k q of q , β k p   in Equation (18). The same similarity was also observed between R , β k p and q , β k p in Equation (19) and their higher-order derivatives.
In summary, the nth derivatives of the displacement vector with respect to either β k p or β k q as the independent variable can be briefly expressed as, based upon Equations (18)–(20),
q , β k p n = a n U k p + b n U k q
where U k p and U k q are the unit-load solutions of Equation (14), and a n and b n are scalar functions of random variables, β k p or β k q , and the lower-order derivatives, q , β k p n 1 , and   q , β k p n 2 , which should be available at the time when q , β k p n is computed. The derivatives of the constraint reaction force vector can also be expanded in terms of I k p and I k q in a similar matter. Specifically, one has the following relation,
R , β k p n = a n I k p + b n I k q μ V V T q , β k p n
Once the displacement field is treated as random, its randomness will be propagated to the von Mises stresses at the interior integration points first by Equation (9), and then to the nodal points through the smoothing process of Equation (11). In the deterministic environment, the von Mises stresses are usually calculated first at the interior points in each element, based upon Equation (9). These are then collected through the entire domain and mapped to the nodal points based upon Equation (11). To take advantage of the linear nature of the unit-load approach, the original sequence to calculate the nodal von Mises stresses was reversed in this study. The state of the stresses and its derivatives are first calculated as before at the interior integration points,
σ e , β k p n = σ x , β k p n σ y , β k p n σ x y , β k p n = D B q e , β k p n = a n D B U e , k p + b n D B U e , k q
where σ e , β k p n consists of three components of the stresses at an interior point in an element. Each of the stress components is then collected and mapped separately to the nodes in the structural by using Equation (9). As an example, one can assemble a vector, σ x e , β k p n , which collects all axial stresses σ x , β k p n reported at the interior points. The value of the axial stress σ x , β k p n can be recast as
σ x , β k p n = a n D 1 B U e , k p + b n D 1 B U e , k q
where D 1 is the first row of D. It is noted that the stress, σ x , β k p n , is a linear combination of unit-loads. Thus, the smoothing process can be conducted separately to map σ x e , β k p n from interior points to nodes as
S x , β k p n = a n S x , k p n + b n S x , k q n
where M S x , k p n = D 1 B U k p and M S x , k q n = D 1 B U k q are solved for the states of the stresses at the nodal points produced by the unit-loads, U k p and U k q , respectively. The M matrix described here is the same one described in Equation (11). The same process can be extended to find the normal stress S y , β k p n and shear stress T x y , β k p n . As a result, the nth order derivatives of the state of the stresses at a nodal point, S , β k p n , can now be expressed as a linear combination of the nth-order derivatives of the nodal state of the stresses produced by the unit-0loads, respectively, as
S , β k p n = S x , β k p n S y , β k p n T x y , β k p n = a n S k p n + b n S k q n
It should be noted that the above equations are derived based upon the assumption that the random variable of concern, β k p , is associated with only one constraint, C k = 0 .   The same derivation process can be further extended to the cases where the random variables of concern involve more than one constraint. Example 1 in the Results section demonstrates a case when a random variable is involved in two multipoint constraints.
Once the derivatives of all stress components at the nodes become available, one can use Equations (12) and (13) to find the mean and standard deviation of the nodal von Mises stresses. As an example, set β k p as a random variable with the mean and standard deviation, μ β k p and σ β k p . The first- and second-order derivatives of the von Mises stress at the nodal point of concern can now be calculated by differentiating Equation (9) at each nodal point as
σ V M , β k p = 1 σ V M S T A S , β k p
and
σ V M = 1 σ V M S , β k p T A S , β k p + S T A S , β k p 1 σ V M 3 S T A S , β k p 2
where σ V M   in Equations (28) and (29) is the nodal von Mises stresses node of concern, and S , β k p and S , β k p are the first- and second-order derivatives of nodal state of stresses as defined in Equations (22) and (23) with n = 1 and n = 2 , respectively. All terms in the right-hand sides of Equations (28) and (29) are evaluated at the mean of the random variable, μ β k p . Once the first- and the second-order derivatives of the nodal von Mises stress become available, one can use Equations (12) and (13) to compute its mean and standard deviation.
In summary, it has been demonstrated in this subsection that only a few unit-load solutions are required to compute any order of derivatives of the nodal displacements, reactions, and von Mises stresses with respect to the random point loads and boundary conditions. As a result, the unit-load approach benefits the use of the Taylor series expansion to efficiently calculate the means and standard deviations of any nodal response function. Such benefits can be further extended to conduct reliability-based design optimization in which the design gradients of the reliability indices can be obtained analytically by differentiating the few unit-load solutions.
The next subsection will elaborate on the design optimization formulation for the reliability-based design of a linear structure subjected to random loads and boundary conditions, and on the design gradient computation by directly differentiating the unit-load solutions. This is followed by Section 3, which presents two examples to demonstrate the processes in the reliability analysis, design sensitivity analysis, and design optimization.

2.3. Reliability-Based Design Optimization with Random Loads and Boundary Conditions

Reliability-based design optimization has been widely accepted as a method that balances the optimal designs with the confidence of probabilistic design constraints. The method incorporates statistical measures into design constraints by developing reliability indices to quantify uncertainties in terms of the levels of reliability.
An optimum design problem aims to identify the best set of design variables, b , that can minimize the objective, f, and satisfy the given and inequality constraints, g i 0 , and equality ones, h j = 0 . Specifically, the mathematical expression of a generalized optimum design problem can be stated in a deterministic form as follows,
m i n b R k   f b s u b j e c t e d   t o : g i b , q 0 , i = 1,2 , m h j b , q = 0 , j = 1,2 , n b L b b U ,
In most engineering design problems, the constraints usually involve limits on nodal displacements, nodal stresses, and nodal reactions. These limiting factors are assigned to guard against the failure of a structure to sustain applied loads. In this study, the applied loads and the boundary conditions were treated as random. Their uncertainty will propagate through the structure to its stress, deflection, and reaction, which make up the constraint equations defined in Problem (P1). Therefore, both the equality and inequality constraints stated in Problem (P1) must incorporate these uncertainties with measurable terms. For instance, the success of a design against a stress inequality constraint at a given node i, g i = σ V M i σ y d   0 , must now be measured by a targeted level of failure probability,   P f 0
        P g i P g i 0       P g 0                   i = 1,2 , n d  
where n d   is the total number of element nodes, and σ y d   a n d   σ V M i   represent the yield stress and the von Mises stress, respectively. The probability failure of the stress constraint,   P g i , can also be stated by a complementary term, called reliability,
      R g i = 1 P g i ,                       i = 1,2 , n d
If the randomness follows the normal distribution, the above reliability can be calculated by the CDF of a single variable, called the reliability index, β g i , as
R g i = Φ β g i
which is a ratio between the mean and the standard deviation of all possible values of g i . That is,
β g i = μ g i s g i = μ V M i μ y d s V M i 2 + s y d 2
The inequality constraint, Equation (26), can then be recast as a single value constraint of the reliability indices,
β g 0 β g i 0             i = 1,2 , n d
where the value of the reliability index, β g 0 , is associated with the targeted reliability,
β g 0 = Φ 1 1   P g 0
As for the equality constraints in Problem (P1), its uncertainty can be controlled by imposing tolerances. For example, one can require the constraint, h j = 0 , fall into a range of a specified bilateral tolerance, τ h j , with a pre-determined probability, P h j as
P τ h j h j τ h j P h j
ssuming that the distribution is normal, one may set the mean value of the constraint as zero, μ h j = 0 , and the standard deviation of s t d h j as follows, based upon the required tolerance, τ h j and the probability, P h j , as
  s t d h j = τ h j β P h j / 2
where the reliability index, β P h j / 2 , is determined based upon the specified probability, P h j , as
β P h j / 2 = Φ 1 1 1   P h j 2
For instance, the constraint that falls into the range 0.02   h j 0.02 with 99.73% probability can be replaced by a pair of constraints. One is an equality constraint on its mean and the other is an inequality constraint on its standard deviation as
μ h j = 0 s t d h j 0.0067
Note that in this case, β P h j / 2 = 3 or h j = 0 ± 0.02 , follows the format of natural tolerance.
In summary, with specified probabilities of failure, the deterministic Problem (P1) can now be reformulated in terms of reliability indices as
m i n b R k   f b subjected   to :       β g 0 β g i μ g i , s g i , 0 , i = 1,2 , m μ h j = 0 , s t d h j τ h j β P h j / 2   , j = 1,2 , n b L b b U ,
All variables in Problem (P2) are deterministic including the means and the standard deviations of the input random variables and the constraints. The latter can be calculated based upon Equations (12) and (13). It is worthwhile mentioning that the magnitude of the mean value calculated by Equation (12) is most likely dominated by the first term, and not the second term as the former is a function of the mean values while the second term is the sum of the square terms of the standard deviations. On the other hand, the magnitude of the standard deviations, as indicated in Equation (13), is a sum of the square terms of the standard deviations. Thus, from a designer point of view, it will be more effective to accomplish the constraints on the mean values such as, μ h j = 0 by changing the mean values of the design variables, while it will be more effective to change the standard deviations to accomplish the constraints on reliability such as s t d h j 0.0067 in Equation (32). In Section 3, Example 1 will take into consideration the standard deviations of the random variables as design variables to satisfy the standard deviation requirements.
To effectively support gradient-based design optimization, it is preferable to have the first-order design sensitivity, with respect to the design variables, computed efficiently for response functions involved in the design consideration. In a structural reliability-based design optimization problem, the probability of failure involves the means and the standard deviations of the nodal displacement, reaction forces, and von Mises stresses. The design derivatives of these means and standard deviations with respect to the structural design variables can be calculated based upon Equations (12) and (13). These design variables can be either random or deterministic.
Set b l as a deterministic variable that is related to the material properties and dimensions of the structure of concern. The design derivatives of the mean and the standard deviation of a response function, g x i ,   are then obtained as
d μ g d b l d g μ x 1 , μ x 2 , , μ x n d b l + 1 2 i = 1 n d g , x i d b l s x i 2
d S g d b l = i = 1 n d g , x i d b l s x i 2 2 S g
As presented in Section 2.2, the response functions, g x i , which include the nodal displacements, reactions, and von Mises stresses of a linear structure, and their derivatives, g , x i and g , x i , evaluated at the means, can be recursively expressed as a linear combination of the unit-load solutions, denoted by U k p and U k q . Therefore, their first-order design derivatives can be obtained as a linear combination of the design derivatives of these unit-load solutions. This is also be conducted by directly differentiating Equation (14) as
K + μ V V T d U i d b l = d K d b l U i μ d V d b l V T U i + V d V T d b l U i
Note that the mean value of any random variable, say x k , can be selected as a design variable (i.e., μ x k b l ) . Since function g and its derivatives,   g , x i and g , x i , in Equations (35) and (36), are all evaluated at the means of the random variables, one can follow the standard procedure to differentiate g , g , x i , and g , x i with respect to μ x k , as per Equations (35) and (36) with respect to b l . In the case where the standard deviation of any random variable, S x k , is treated as a design variable, Equations (35) and (36) can now be simplified as
d μ g d S x k = g , x k s x i
and
d S g d S x k = g , x k 2 s x k S g
This is again because g ,   g , x i , and g , x i are evaluated at the means of all random variables, and consequently, they are independent of the standard deviations.

3. Results

Two examples of RBDO are presented in this section to demonstrate the use of the unit-load solutions in Taylor’s series expansion to estimate the means and the standard deviations of the response functions computed by the finite element method. The first example was a truss problem subjected to random multipoint constraints. It incorporates an equality constraint and includes standard deviations of random variables as design variables. The second example conducted shape optimization of a mooring bracket subjected to two random point loads. The bracket was modeled as a plane stress problem, which was discretized with quadrilateral elements. The unit-load approach presented in Section 2 was employed here to compute the constraints, which were formulated in terms of the reliability indices of the nodal displacements, reactions, and von Mises stresses. Monte Carlo simulation was used later to validate the results of the Taylor series expansion and the reliability-based design optimization.

3.1. Example 1–Random Multipoint Constraints

A horizontal rigid bar with length l 3 was hinged at the left end and subjected to a vertical point load P at the free end, as shown in Figure 1. The problem was taken from [35]. The bar was supported by two vertical truss bars, Bar 1 and Bar 2. They had the same lengths, while their Young’s modulus and cross-sectional area were not the same. The length denoted by L was set as 36 inches. The Young’s modulus and the cross-sectional area of Bar 1, E 1 ,   a n d   A 1 were set as 30 × 10 6 psi and 1 in2, and those of Bar 2, E 2   a n d   A 2 , were set as 10 × 10 6 psi and 1.25 in2, respectively. The distances from the hinge point to Bar 1 and to Bar 2 are denoted by l 1 and l 2 , respectively. The lengths, l 1 , l 2 , l 3 , and the point load P were considered as normal random variables in this study. With the initial covariances being 0.15, they can be given by l 1 ~ N 15,2.25 ,     l 2 ~ N 27,4.05 ,     l 3 ~ N 36,5.40 ,   and P ~ N 15,000 ,   2250 .   The goal of this design optimization exercise was to find the lightest design that can control the deflection at the free end of the rigid bar be in the range of 0.448 δ m a x 0.452 , with 99.73% of probability. To this end, the cross-sectional areas of the vertical bars, A 1 , A 2 ,   as well as the standard deviations of l 1 , l 2 , and l 3 were considered as the design variables. The solution of this design problem is presented in three sub-sections. The problem statement for deterministic design optimization is stated in Section 3.1.1, which involved only one equality constraint. The second section, Section 3.1.2, uses the second-order Taylor’s series expansion to compute the mean and standard deviations of the nodal displacements and reactions. The last section, Section 3.1.3, presents the results of several reliability-based design optimization runs with different levels of reliability of constraint satisfaction.

3.1.1. Deterministic Finite Element State Equation

The global stiffness and the external load vector of the two-member truss structure are stated as
K 0 = s 1 0 0 s 2 s 1 0 s 1 0 0 0 s 2 0 s 1 0 0 0 s 2 0 0 0 s 2 0 0 0 0   and   f 0 = 0 0 0 0 P
where the initial values of s 1 and s 2  are given by  s 1 = E 1 A 1 L = 0.833 × 10 5 , s 2 = E 2 A 2 L = 0.347 × 10 5 ,  and the point load P is taken at its mean, P = 15 × 10 3 . The problem is subjected to two single point constraints at Nodes 3 and 4,
q 3 = q 4 = 0
and two multiple point constraints
l 3 q 1 l 1 q 5 = 0 l 3 q 2 l 2 q 5 = 0
The above multipoint constraints are imposed to ensure that Nodes 1, 2, and 5 maintain a straight line, even after truss deformation to comply with the rigidness assumption of the horizontal bar.
One may implement the single point constraints first by eliminating rows and columns 3 and 4 to obtain the revised K and f as
K = s 1 0 0 0 s 2 0 0 0 0   and   f = 0 0 P
Next, one can implement the multipoint constraints by the penalty method. This leads to a set of three equations as
K + μ v 1 v 1 T + μ v 2 v 2 T q = f
where v 1 and v 2 are the coefficients of the given multiple point constraints; v 1 T = l 3 , 0 , l 1 and v 2 T = 0 , l 3 , l 2 . More specifically, Equation (39) is detailed as follows for a unique solution of q T = q 1 , q 2 , q 5
s 1 + μ l 3 2 0 μ l 1 l 3 0 s 2 + μ l 3 2 μ l 2 l 3 μ l 1 l 3 μ l 2 l 3 μ l 1 2 + l 2 2 q 1 q 2 q 5 = 0 0 P
along with the reaction force vector at Nodes 1, 2, and 5,
R = μ l 3 2 q 1 l 1 l 3 q 5 l 3 2 q 2 l 2 l 3 q 5 l 1 l 3 q 1 l 2 l 3 q 2 + l 1 2 + l 2 2 q 5

3.1.2. Means and Standard Deviations of Displacement and Reactions

The means and standard deviations of the nodal displacement vector, q T = q 1 , q 2 , q 5 , can be calculated based upon Equations (12) and (13), resulting from the Taylor’s series expansion, as
μ q = q ¯ + 1 2 q , l 1 σ l 1 2 + q , l 2 σ l 2 2 + q , l 3 σ l 3 2 + q , P σ P 2
and
V a r q = σ q 2 = q , l 1 σ l 1 2 + q , l 2 σ l 2 2 + q , l 3 σ l 3 2 + q , P σ P 2
where μ q and σ q are the mean and standard deviation vectors of q , and q ¯ denotes the displacement solution of Equation (39) with all random variables evaluated at their means.
The next task is to derive the first- and second-order derivatives presented in Equations (41) and (42) and evaluate them at the means of random variables. The first- and the second-order derivatives of the displacement vector with respect to the load force, P , is straight forward as
q , P = U 5
and
q , P = 0
where U 5 is the unit-load displacement vector at Node 5. The first- and second-order derivatives of the displacement vector with respect to the random variables, l 1 and l 2 , can be directly derived by using Equations (18) and (19), as they appear only in one of the constraint equations. The resultant equations are listed below
q , l 1 = μ 2 l 1 q 5 l 3 q 1 U 5 + μ l 3 q 5 U 1 q , l 2 = μ 2 l 2 q 5 l 3 q 2 U 5 + μ l 3 q 5 U 2
and
q , l 1 l 1 = 2 μ 2 l 1 q 5 l 3 q 1 + q 5 U 5 + 2 μ l 3 q 5 U 1 q , l 2 l 2 = 2 μ 2 l 2 q 5 l 3 q 2 + q 5 U 5 + μ l 3 q 5 U 2
where U 1 and U 2 are the unit-load displacements with non-zero force imposed at Nodes 1 and 2, respectively.
As for the random variable, l 3 , it is presented in both multipoint constraints, which is different from other random variables, l 1 and l 2 , as each of the latter is presented in only one multipoint constraint. Consequently, a revision of Equations (18) and (19) is required to complete the differentiation as expressed below,
q , l 3 = μ 2 l 3 q 1 l 3 q 5 U 1 μ 2 l 3 q 2 l 2 q 5 U 2 μ l 1 q 1 + l 2 q 2 U 5 q , l 3 l 3 = 2 μ 2 l 3 q 1 l 3 q 5 + q 1 U 1 2 μ 2 l 3 q 2 l 2 q 5 + q 2 U 2 2 μ l 1 q 1 + l 2 q 2 U 5
The first- and second-order derivatives of the reaction vectors, R , with respect to the random variables, l 1 , l 2 and l 3 , can also be calculated based upon Equations (21) and (22), where R is defined by Equation (40).Again, some revision is required to calculate the derivatives with respect to l 3 , which is involved in both multipoint constraints simultaneously. The equations of the first- and second-order derivatives of R   with respect to l 1 , l 2 , and l 3 are listed in Appendix A for reference.
Once the derivatives of the displacement and the reaction force vectors with respect to the random variables at their means become available, one can follow Equations (12) and (13) to find the mean and standard deviation of these vectors. For validation, let the variables, l 1 , l 2 , l 3 ,   P be normally distributed and statistically independent with covariances of 0.15. Specifically, at the initial design, these are given by l 1 ~ N 15 ,   2.25 ,   l 2 ~ N 27 ,   4.05 , l 3 ~ N 36 ,   5.40 , and P ~ N 15,000 ,   2250 . Table 1 and Table 2 compare the means and standard deviations of the displacement and the reaction vector calculated by the Taylor’s series expansion presented in this section with those by the Monte Carlo simulation with samples of 106. The maximal error of the displacement vector was observed at the standard deviation of the displacement at the free end of the horizonal bar, Node 5, while the maximal error of the reaction vector was the standard deviation at Node 1. With reduced variations, the errors in standard deviations were reduced, as demonstrated in Table 3, for the displacement vector. It should be noted that only three unit-load analyses for U 1 , U 2 , and U 5 were required to complete the entire validation process presented in Table 1, Table 2 and Table 3.

3.1.3. Reliability-Based Design Optimization

The goal of the design optimization in Example 1 was to find the lightest truss bars that could control the deflection at Node 5, the free end of the rigid bar falling into the range 0.43,0.47 , with different levels of reliability. To achieve this goal, the concerned design variables included the cross-sectional areas of the truss members, which are deterministic, and the standard deviations of the input random variables. In particular, the mathematical expression of the design optimization problem is expressed in terms of statistical terms as
m i n A 1 ,   A 2 ,     s t d l 1 , s t d l 2 , s t d l 3 , s t d p f = A 1 + A 2     subject   to :   μ q 5 = 0.45 s t d q 5 0.02 β P h j / 2
where the cross-sectional areas of the vertical bars A 1 , A 2   as well as the standard deviations of l 1 , l 2 , and l 3 , denoted as s t d l 1 , s t d l 2 , and s t d l 3 , are considered as the design variables. The standard deviation of the point load, P, s t d p ,   can also be included as a design variable. The reliability index, β P h j / 2 , in the above constraint, is associated with a 0.02 tolerance of the mean value μ q 5 , as defined by Equations (32)–(34).
Five different cases of design optimization formulation (P3) were solved by the MATLAB built-in function, fmincon. The required gradients were supplied by the finite difference method. Their optimal solutions are summarized in Table 4. Case 1 was a deterministic design optimization in which only two variables, A 1 and A 2 , were treated as design variables, and the force and the boundary conditions were deterministic. Next, problem (P3) was solved to achieve 95% of reliability, with the standard deviation of P set as a constant, 750. The problem failed to converge. However, the problem could converge once the targeted reliability was reduced to 60%. The associated results are listed in Case 2. It is shown that the cross-sectional area, A 1 ,   approached its lower limit. The problem as then resolved by including the standard deviation of the pressure, s t d p , as an additional design variable. The results of optimization runs for different ranges of P h j described in Equation (31) are listed in Cases 3 to 5 in Table 4. The results show that the free end deflection, q 5 , could reach 99.5% of confidence or fall into the range of 0.43,0.47 .   Finally, the last column in Table 4 indicates that the reliability levels of the optimal solutions computed by the Taylor’s series expansion approach were matched very well with those by the Monte Carlo simulation with 10 6   samples. Note that only three unit-load displacement vectors, U 1 , U 2 , and U 5 were required in this example to complete this Monte Carlo simulation task at any optimal solution of concern.

3.2. Example 2–Shape Optimization of a Mooring Bracket Subject to Random Loads

The example problem was a bracket that is used for mooring small boats. The dimension of the mooring bracket is shown in Figure 2. The random point load was imposed at the contact, Point P, on the side of the hole, caused by the pressure of the mooring clip. The shape of the rectangular cutoff at the base of the mooring bracket was the concern of the design. The goal of this design example was to find the best cutoff shape that minimized the weight of the mooring bracket, while maintaining a satisfactory level of reliability for constraints on the displacement, von Mises stress, and reaction forces. The shape of the cutoff was measured by three variables: the left edge, x L , the height, x H , and the right edge, x R . Their best values were scaled between their initial values, (1.14, 0.13, 2.8) and their upper bounds (1.5, 0.31, 3.07) by introducing three design variables ranging between 0 and 1 as
x L = 1.50 + b 1 1.14 1.50 x H = 0.13 + b 2 0.31 0.13 x R = 2.80 + b 3 3.07 2.80
The initial finite element model of the mooring bracket is displayed in Figure 3, which was made of 238 quadrilateral elements and 284 total nodes, of which 13 were along the fixed base boundary. The two components of the applied load at points P, P x , and P y were assumed to be independent of each other and normally distributed. The maximal allowable nodal displacement, von Mises stress, and reaction force were also assumed to be normal, which were set as N 0.0013 ,   0.00013 i n , N 3.0 E 5 ,   3.0 E 4 p s i , and N 290,29 l b s , respectively.
To validate the use of the Taylor’s series expansion for the reliability analysis of this example problem, the results of the unit-load-based approach described in Section 2 were compared with those generated by the Monte Carlo simulation with 15,000 samples. The output responses of concern in this example were the reliability indices of the magnitudes of the nodal displacement, nodal von Mises stresses, and nodal boundary reactions. Two independent random variables were considered in this study, P x and P y . The mean values of these two load components were set as 750 lbs and 600 lbs, respectively. The results of the maximal displacement and the maximal von Mises stresses at Point P and the maximal reaction at Point R are tabulated in Table 5, which were computed for the values of the variances set at 10%, 15%, and 20% for the loads P x and P y . It was noted that increasing the value of COV in random loads would increase more in the standard deviations of the structural responses than in the mean values, as the former are linear in the standard deviations of input variables while the latter are quadratic. It is worthwhile mentioning here that all the results listed in Table 5 were completed by using only two unit-load analyses.
Next, following problem statement (P2), the shape optimization of the given mooring bracket is now modeled as a standard design problem with inequality constraints expressed in terms of the reliability indices of nodal displacements, β q i ,   stresses,   β s t j , and boundary reactions,   β r k :
m i n b 1 , b 2 , b 3 V o l u m e
subject to
β q 0 β q 0 ,         i = 1,2 , n β s t j β s t 0 ,         j = 1,2 , n β r k β r o ,         k = 1,2 , m
where the lowest bounds of the required reliabilities, β q 0 ,   β s t 0 ,     a n d   β r o , were set at 2.0, which corresponded to 97.7% of reliability, Furthermore, the total nodal number, n, in the above problem statement was found to be 284, while the total boundary node, m, was 13, based upon the initial finite element mesh in Figure 3. In the initial design, all of the design variables were set at zero. The corresponding reliability index distributions of the displacement field and the stress field are displayed in Figure 4 and Figure 5. The volume of the initial design was 0.6213 in3. However, it comes along with a violation on the boundary reaction at point Q, as indicated in Figure 2. The values of the mean and standard deviation of this reaction force were 203.7 lbs and 19.12 lbs, which resulted in a reliability index of 1.908, less than the required value of 2.
To conduct the shape optimization, the MATLAB built-in function, fmincon, was used to solve the optimization problem. The required gradients were computed by the finite difference method. The final design had one design variable hitting the upper bound and two constraints hitting the lower bound. Specifically, the design variables of the final design were 1.0 ,   0.8247 ,   0.3604 with a volume of 0.5412 in3. One of the tight constraints was associated with the reliability index of the nodal von Mises at Point P and the other was the nodal boundary reaction force at Point R. Both had a reliability index of 2. The Monte Carlo simulation with 15,000 samples was used to compute the reliability indices of all constraints imposed upon the optimal design. The maximal reliability index was 2.0008, very close to the upper bound set by the optimization formulation. This maximal reliability index was associated with the maximal reaction force at Point R. The corresponding mesh and the reliability index distributions of the displacement field and the stress field of the optimal design are displayed in Figure 6, Figure 7 and Figure 8.

4. Conclusions

Deterministic design problems are, in general, not practical, as many do not count variables with inherent uncertainty into the formulation such as mechanical properties, support conditions, and external loadings. Implementing a reliability-based design can help encapsulate many of the uncertainties and propagate the stochastic variables into the structural design process to produce more reliable and confident results. Unfortunately, incorporating uncertainties into a design is often computationally expensive, especially when they involve an iterative optimization solution process. This study considered the design of a linear structure subjected to the random point loads and boundary conditions. For its simplicity in implementation, the Taylor series expansion was applied to find accurate statistical measures of the structural responses, which included the nodal displacements, the reaction forces, and the nodal von Mises stresses. To take advantage of the linear nature of structural problems, a unit-load approach was applied in this study. For a given design, this approach can generate reliability indices of all the nodal values of the structural responses by using only a few unit-load displacement vectors. These unit-loaded displacement vectors can be repeatedly used to find any order of derivatives required in the Taylor series expansion with respect to random variables. Two design optimization examples were presented in this study. The first example was a truss design problem subjected to randomness in the multipoint constraints. The penalty method was used to convert the random multipoint constraints to random loads. Consequently, only four unit-load analyses were required in each design iteration of the reliability-based optimization process. The second example dealt with the nodal von Mises stresses in the shape optimization of a plane stress problem under random point loads. A new approach was implemented here to propagate the uncertainty linearly from the nodal displacements to the nodal von Mises stresses, so only two unit-load equations were solved in each design optimization iteration to compute the required reliability indices. The results of the optimization were validated by Monte Carlo simulation with 106 samples for Example 1 and 15,000 for Example 2. The entire validation process, however, was completed with only four unit-load equations for Example 1 and two for Example 2.
These examples demonstrate that the method of unit-loads is very effective and accurate to handle uncertainty and support reliability index analysis as well as design optimization for a wide variety of applications. The applications of the current study focused upon linear static analysis. It is expected that the proposed unit-load approach can be extended to linear dynamic problems, which will be a topic for future study.

Author Contributions

Conceptualization, R.J.H. and G.J.-W.H.; methodology, R.J.H. and G.J.-W.H.; software, R.J.H. and G.J.-W.H.; validation, G.J.-W.H.; formal analysis, R.J.H.; investigation, R.J.H. and G.J.-W.H.; resources, R.J.H. and G.J.-W.H.; data curation, R.J.H.; writing—original draft preparation, R.J.H.; writing—review and editing, G.J.-W.H.; visualization, R.J.H.; supervision, G.J.-W.H.; project admin-istration, G.J.-W.H. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

The first- and second-order derivatives of the reactions at Nodes 1, 2, and 5 of Example 1 are listed below for reference.
R , l 1 = μ l 3 2 q 1 , l 1 l 1 l 3 q 5 , l 1 l 3 q 5 l 3 2 q 2 , l 1 l 2 l 3 q 5 , l 1 l 3 q 5 l 1 l 3 q 1 , l 1 l 2 l 3 q 2 , l 1 + l 1 2 + l 2 2 q 5 , l 1 l 3 q 1 + 2 l 1 q 5 R , l 2 = μ l 3 2 q 1 , l 2 l 1 l 3 q 5 , l 2 l 3 q 5 l 3 2 q 2 , l 2 l 2 l 3 q 5 , l 2 l 3 q 5 l 1 l 3 q 1 , l 2 l 2 l 3 q 2 , l 2 + l 1 2 + l 2 2 q 5 , l 2 l 3 q 2 + 2 l 2 q 5 R , l 3 = μ l 3 2 q 1 , l 3 l 1 l 3 q 5 , l 3 + 2 l 3 q 1 l 1 q 5 l 3 2 q 2 , l 3 l 2 l 3 q 5 , l 3 + 2 l 3 q 2 l 2 q 5 l 1 l 3 q 1 , l 3 l 2 l 3 q 2 , l 3 + l 1 2 + l 2 2 q 5 , l 3 l 1 q 1 l 2 q 2
while the second-order derivatives are
R , l 1 l 1 = μ l 3 2 q 1 , l 1 l 1 l 1 l 3 q 5 , l 1 l 1 2 l 3 q 5 , l 1 l 3 2 q 2 , l 1 l 1 l 2 l 3 q 5 , l 1 l 1 l 3 q 5 , l 1 l 1 l 3 q 1 , l 1 l 1 l 2 l 3 q 2 , l 1 l 1 + l 1 2 + l 2 2 q 5 , l 1 l 1 2 l 3 q 1 , l 1 + 4 l 1 q 5 , l 1 + 2 q 5 R , l 2 l 2 = μ l 3 2 q 1 , l 2 l 2 l 1 l 3 q 5 , l 2 l 2 l 3 q 5 , l 2 l 3 2 q 2 , l 2 l 2 l 2 l 3 q 5 , l 2 l 2 2 l 3 q 5 , l 2 l 1 l 3 q 1 , l 2 l 2 l 2 l 3 q 2 , l 2 l 2 + l 1 2 + l 2 2 q 5 , l 2 l 2 2 l 3 q 2 , l 2 + 4 l 2 q 5 , l 2 R , l 3 l 3 = μ l 3 2 q 1 , l 3 l 3 l 1 l 3 q 5 , l 3 l 3 + 4 l 3 q 1 , l 3 2 l 1 q 5 , l 3 l 3 2 q 2 , l 3 l 3 l 2 l 3 q 5 , l 3 l 3 + 4 l 3 q 2 , l 3 2 l 2 q 5 , l 3 l 1 l 3 q 1 , l 3 l 3 l 2 l 3 q 2 , l 3 l 3 + l 1 2 + l 2 2 q 5 , l 3 l 3 2 l 1 q 1 , l 3 2 l 2 q 2 , l 3
The means and the standard deviations of the reaction force vector, R T = R 1 , R 2 , R 5 , can be calculated based upon Equations (12) and (13), resulting from the Taylor’s series expansion as
μ R = R ¯ + 1 2 R , l 1 σ l 1 2 + R , l 2 σ l 2 2 + R , l 3 σ l 3 2 + R , P σ P 2
and
V a r R = σ R 2 = R , l 1 σ l 1 2 + R , l 2 σ l 2 2 + R , l 3 σ l 3 2 + R , P σ P 2
where μ R and σ R are the mean and the standard deviation vectors of R , and R ¯ is computed based upon Equation (7), in which all random variables are evaluated at the means.

References

  1. Stefanou, G. The Stochastic Finite Element Method: Past, Present and Future. Comput. Methods Appl. Mech. Eng. 2009, 10, 1031–1051. [Google Scholar] [CrossRef] [Green Version]
  2. Gunzburger, M.D.; Webster, C.G.; Zhang, G. Stochastic Finite Element Methods for Partial Differential Equations with Random Input Data. Acta Numer. 2014, 23, 521–650. [Google Scholar] [CrossRef]
  3. Arregui-Mena, J.D.; Margetts, L.; Mummery, P.M. Practical Applications of the Stochastic Finite Element Method. Arch. Comput. Methods Eng. 2016, 23, 171–190. [Google Scholar] [CrossRef]
  4. Aldosary, M.; Wang, J.; Li, C. Structural Reliability and Stochastic Finite Element Methods: State-of-the-Art Review and Evidence-based Comparison. Eng. Comput. 2018, 35, 2166–2214. [Google Scholar] [CrossRef] [Green Version]
  5. Xu, H.; Rahman, S. Decomposition Methods for Structural Reliability Analysis. Probabilistic Eng. Mech. 2005, 20, 239–250. [Google Scholar] [CrossRef]
  6. Rahman, S. A Polynomial Dimensional Decomposition (PDD) for Stochastic Computing. Int. J. Numer. Methods Eng. 2008, 76, 2091–2116. [Google Scholar] [CrossRef]
  7. Zhao, J.; Wang, C. Robust Topology Optimization under Loading Uncertainty based upon Linear Elastic Theory and Orthogonal Diagonalization of Symmetric Matrices. Comput. Methods Appl. Mech. Eng. 2014, 273, 204–218. [Google Scholar] [CrossRef]
  8. Guo, T.; Frangopol, D.M.; Chen, Y. Fatigue Reliability Assessment of Steel Bridge Details Integrating Weight-in-Motion Data and Probabilistic Finite Element Analysis. Comput. Struct. 2012, 112–113, 245–257. [Google Scholar] [CrossRef]
  9. Georgioudakis, M.; Lagaros, N.D.; Papdrakakis, M. Probabilistic Shape Design Optimization of Structural Components under Fatigue. Comput. Struct. 2017, 182, 252–266. [Google Scholar] [CrossRef]
  10. Zheng, W.; Shen, J. Adjustable Hybrid Resampling Approach to Computationally Efficient Probabilistic Inference of Structural Damage Based on Vibration Measurements. J. Civ. Struct. Health Monit. 2016, 6, 153–173. [Google Scholar] [CrossRef]
  11. Zheng, W.; Qian, F.; Shen, J.; Xiao, F. Mitigating Effects of Temperature Variations through Probabilistic-based Machine Learning for Vibration-based Bridge Scour Detection. J. Civ. Struct. Health Monit. 2020, 10, 957–972. [Google Scholar] [CrossRef]
  12. Haldar, A.; Mahadevan, S. Probability, Reliability and Statistical Methods in Engineering Design; John Wiley & Sons: Hoboken, NJ, USA, 2000. [Google Scholar]
  13. Ghanem, R.G.; Spanos, P.D. Stochastic Finite Elements: A Spectral Approach, Revised ed.; Dover Publication: Mineola, NY, USA, 2003. [Google Scholar]
  14. Lucor, D.; Su, C.-H.; Karniadakis, G.E. Generalized Polynomial Chaos and Random Oscillators. Int. J. Numer. Methods Eng. 2004, 60, 571–596. [Google Scholar] [CrossRef]
  15. Sudret, B. Polynomial Chaos Expansions and Stochastic Finite Element Methods. In Risk and Reliability in Geotechnical Engineering; Phoon, K.-K., Ching, J., Eds.; CRC Press: Boca Raton, FL, USA, 2014; Chapter 6; pp. 265–300. [Google Scholar]
  16. Xiu, D.; Karniadakis, G.E. The Weiner-Askey Polynomial Chaos for Stochastic Differential Equations. SIAM J. Sci. Comput. 2002, 24, 619–644. [Google Scholar] [CrossRef]
  17. Kaintura, A.; Dhaene, T.; Spina, D. Review of Polynomial Chaos-Based Methods for Uncertainty Quantification in Modern Integrated Circuits. Electronics 2018, 7, 30. [Google Scholar] [CrossRef] [Green Version]
  18. Son, J.; Du, D.; Du, Y. Modified Polynomial Chaos Expansion for Efficient Uncertainty Quantification in Biological Systems. Appl. Mech. 2020, 1, 153–173. [Google Scholar] [CrossRef]
  19. Budynas, R.G.; Nisbett, J.K. Shigley’s Mechanical Engineering Design, 10th ed.; McGraw Hill Education, 2012 Penn Plaza: New York, NY, USA, 2015. [Google Scholar]
  20. Kaminski, M. Stochastic Perturbation Approach to Engineering Structure Vibration by the Finite Difference Method. J. Sound Vib. 2002, 251, 651–670. [Google Scholar] [CrossRef]
  21. Ghosh, D.; Ghanem, R.G.; Red-Horse, J. Analysis of Eigenvalues and Modal Interaction of Stochastic Systems. AIAA J. 2005, 43, 2196–2201. [Google Scholar] [CrossRef]
  22. Xiu, D.; Karniadakis, G.E. Supersensitivity due to Uncertain Boundary Conditions. Int. J. Numer. Methods Eng. 2004, 61, 2114–2138. [Google Scholar] [CrossRef]
  23. Cavdar, O.; Bayraktar, A.; Cavdar, A.; Adanur, S. Perturbation Based Stochastic Finite Element Analysis of the Structural Systems with Composite Sections under Earthquake Forces. Steel Compos. Struct. 2008, 8, 129–144. [Google Scholar] [CrossRef]
  24. Kaminski, M. On Generalized Stochastic Perturbation-Based Finite-Element Method. Commun. Numer. Methods Eng. 2006, 22, 23–31. [Google Scholar] [CrossRef]
  25. Wu, F.; Gao, Q.; Xu, X.-M.; Zhong, W.-X. A Modified Computational Scheme for the Stochastic Perturbation Finite Element Method. Lat. Am. J. Solids Struct. 2015, 12, 2480–2505. [Google Scholar] [CrossRef] [Green Version]
  26. Kaminski, M. On the Dual Iterative Stochastic Perturbation-Based Finite Element Method in Solid Mechanics with Gaussian Uncertainties. Int. J. Numer. Methods Eng. 2015, 104, 1038–1060. [Google Scholar] [CrossRef]
  27. Ghosh, D.; Farhad, C. Strain and Stress Computations in Stochastic Finite Element Methods. Int. J. Numer. Methods Eng. 2007, 74, 1219–1239. [Google Scholar] [CrossRef]
  28. Varz, M., Jr.; Munoz-Rojas, P.A.; Fillippini, G. On the Accuracy of Nodal Stress Computation in Plane Elasticity Using Finite element Volumes and Finite Elements. Comput. Struct. 2009, 87, 1044–1057. [Google Scholar]
  29. Jia, Y.; Bergheau, J.-M.; Leblond, J.-B.; Roux, J.-C.; Bouchaoui, R.; Gallee, S.; Brosse, A. A New Nodal-Integration-Based Finite Element Method for the Numerical Simulation of Welding Process. Metals 2020, 10, 1386. [Google Scholar] [CrossRef]
  30. Rahman, S.; Rao, B.N. A Perturbation Method for Stochastic Meshless Analysis in Elastostatics. Int. J. Numer. Methods Eng. 2001, 50, 1969–1991. [Google Scholar] [CrossRef] [Green Version]
  31. Chirataksanakul, A.; Mahadevan, S. First-Order Approximation Methods in Reliability-Based Design Optimization. J. Mech. Des. 2005, 127, 851–857. [Google Scholar] [CrossRef]
  32. Neves Careiro, G.D.; Antoniao, C.C. Reliability-based Robust Design Optimization with the Reliability Index Approach Applied to Composite Laminate Structures. Compos. Struct. 2019, 209, 844–855. [Google Scholar] [CrossRef]
  33. Wang, X.; Meng, Z.; Cheng, C.; Long, K.; Li, J. Reliability-based Design Optimization of Material Orientation and Structural Topology of Fiber-reinforced Composite Structure under Load Uncertainty. Compos. Struct. 2022, 291, 115537. [Google Scholar] [CrossRef]
  34. Kharmanda, G.; Gowid, S.; Mahdi, E.; Shokry, A. Efficient System Reliability-Based Design Optimization Study for Replaced Hip Prosthesis Using New Optimized Anisotropic Bone Formulations. Materials 2020, 13, 362. [Google Scholar] [CrossRef] [Green Version]
  35. Chandrupatla, T.R.; Belegundu, A.D. Introduction to Finite Elements in Engineering, 5th ed.; Cambridge University Press: Cambridge, MA, USA, 2022. [Google Scholar]
Figure 1. Truss structure subjected to multipoint constraints.
Figure 1. Truss structure subjected to multipoint constraints.
Designs 07 00096 g001
Figure 2. The initial design of the mooring bracket.
Figure 2. The initial design of the mooring bracket.
Designs 07 00096 g002
Figure 3. Initial mesh of the mooring bracket.
Figure 3. Initial mesh of the mooring bracket.
Designs 07 00096 g003
Figure 4. Reliability index distribution of the displacements in the initial mooring bracket.
Figure 4. Reliability index distribution of the displacements in the initial mooring bracket.
Designs 07 00096 g004
Figure 5. Reliability index distribution of the stresses in the initial mooring bracket.
Figure 5. Reliability index distribution of the stresses in the initial mooring bracket.
Designs 07 00096 g005
Figure 6. Optimal mesh of the mooring bracket.
Figure 6. Optimal mesh of the mooring bracket.
Designs 07 00096 g006
Figure 7. Reliability index distribution of displacements in the optimal mooring bracket.
Figure 7. Reliability index distribution of displacements in the optimal mooring bracket.
Designs 07 00096 g007
Figure 8. Reliability index distribution of stresses in the optimal mooring bracket.
Figure 8. Reliability index distribution of stresses in the optimal mooring bracket.
Designs 07 00096 g008
Table 1. Mean and standard deviation of nodal displacement with all covariance equal to 0.15.
Table 1. Mean and standard deviation of nodal displacement with all covariance equal to 0.15.
Nodal No.Mean—MCMean—Equation (40)Std—MCStd—Equation (41)
10.18470.18470.05160.0504
20.33030.33030.08370.0823
50.46250.46180.19000.1757
Table 2. Mean and standard deviation of nodal MPC reaction forces × 10 3 with all covariance equal to 0.15.
Table 2. Mean and standard deviation of nodal MPC reaction forces × 10 3 with all covariance equal to 0.15.
Nodal No.Mean—MCMean—Equation (40)Std—MCStd—Equation (41)
115.38615.3884.3024.201
211.46311.4612.9042.856
5−15.001−15.0002.2502.250
Table 3. Comparison of the differences in standard deviations of displacements for different covariances in all random variables.
Table 3. Comparison of the differences in standard deviations of displacements for different covariances in all random variables.
Nodal No.Covariance = 0.15Covariance = 0.10Covariance = 0.05
MCEquation (41)MCEquation (41)MCEquation (41)
10.05160.05040.03400.03360.01690.0168
20.08360.08230.05530.05490.02750.0274
50.18980.17570.12110.11720.05910.0586
Table 4. Optimal designs of Example 1 with different P h j .
Table 4. Optimal designs of Example 1 with different P h j .
CaseObjective
Weight
Design VariablesConstraint
A 1 A 2 μ l 1 15 μ l 2 = 27 μ l 3 = 6 μ P = 1.5   ×   10 4 μ q 5 = 0.45 P h j % MC%
s t d l 1 s t d l 2 s t d l 3 s t d p s t d q 5
12.13370.00372.130000000
22.13360.00052.1330.48250.02340.3091750.000.023860.059.99
32.14060.20831.9410.52730.30150.4124640.160.023860.059.96
42.14960.21681.9320.47680.17820.2201204.560.010295.095.00
5 2.13750.05392.0830.52810.11300.1604150.910.007199.599.48
Table 5. Mean and standard deviations of the structural responses for different COV values.
Table 5. Mean and standard deviations of the structural responses for different COV values.
ResponsesInput COVTaylor’s Series ExpansionMonte Carlo Simulation
MeanStandard Dev.MeanStandard Dev.
Displacement × 10 4 at Point P10%8.16970.74418.18650.7415
15%8.16971.11648.20711.1074
20%8.16961.48888.23911.4691
V M stress × 10 5
at Point P
10%2.28520.18942.28530.1893
15%2.28530.28412.28520.2839
20%2.28540.37882.28580.3789
BC reactions × 10 2 at Point R10%−1.60800.1138−1.60810.1138
15%−1.60800.1708−1.60780.1706
20%−1.60800.2277−1.60810.2278
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Haupin, R.J.; Hou, G.J.-W. A Unit-Load Approach for Reliability-Based Design Optimization of Linear Structures under Random Loads and Boundary Conditions. Designs 2023, 7, 96. https://doi.org/10.3390/designs7040096

AMA Style

Haupin RJ, Hou GJ-W. A Unit-Load Approach for Reliability-Based Design Optimization of Linear Structures under Random Loads and Boundary Conditions. Designs. 2023; 7(4):96. https://doi.org/10.3390/designs7040096

Chicago/Turabian Style

Haupin, Robert James, and Gene Jean-Win Hou. 2023. "A Unit-Load Approach for Reliability-Based Design Optimization of Linear Structures under Random Loads and Boundary Conditions" Designs 7, no. 4: 96. https://doi.org/10.3390/designs7040096

APA Style

Haupin, R. J., & Hou, G. J. -W. (2023). A Unit-Load Approach for Reliability-Based Design Optimization of Linear Structures under Random Loads and Boundary Conditions. Designs, 7(4), 96. https://doi.org/10.3390/designs7040096

Article Metrics

Back to TopTop