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
where the
vectors,
and
, are the global force and the displacement, respectively, while
is the
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
where the boundary reaction force,
, is installed to enforce the boundary conditions upon the involved boundary degrees of freedom,
. 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,
, is a typical linear, two-point constraint that involves two degrees of freedom,
and
, which can be expressed as
where the user-provided constants
,
, and
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,
, can be set to zero. Note that this multi-point constraint in Equation (3) can be put in a vector form as
where
is a
vector with two non-zero values,
and
, 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
where
is an
vector and
is a
matrix with
being the total number of constraints. Note that the
kth column of
V is the vector,
, in Equation (4), with two non-zero coefficients,
and
.
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
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
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 . 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
. For example, one can relate the state of the stresses of a plane stress problem at an interior integration point
to the displacement field in a specified finite element,
e, as
where matrix
D is a function of the material properties and matrix
B is composed of the derivatives of shape function evaluated at
, the natural coordinates of the location of the interior integration point
[
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
where the symmetric matrix
A is given by
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,
, as
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,
where the vector,
, is unknown, which is the collection of the von Mises stresses at the nodal points, the vector
is known, which are those computed at the interior integration points, and the matrix
, is the assembly of element shape functions evaluated at the interior integration points over the whole domain as
Note that the matrix, , is independent of the element shape or element material properties for an iso-parametric element.
Once the input load vector,
and the boundary conditions,
, 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,
, for a given set of random variables,
as
where
and
are the mean and the standard deviation of a random variable
and
are the first- and the second-order derivatives of the function,
The latter are all evaluated at means,
. 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 , and 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 of the state equation, Equation (6), is linear in terms of the load vector, and . Therefore, the order of the derivatives of Equation (6) higher than one with respect to any component in and will be zero. However, deriving the design derivatives of Equation (6) with respect to the rest of coefficients, and , in the multipoint constraint will be more involved. This is because matrix on the right of Equation (6) and on the left are the linear and nonlinear functions of and , respectively. Consequently, the displacement is now nonlinear in and .
To start the derivation, set the vector
, being the solution of the linear structure under the unit-load,
, as
where the subscript
i specifically denotes the non-zero degree of freedom, the vector,
, has only one non-zero value of 1 at the
ith component, and
defined in Equation (5) is a function of
and
. Note that
is in fact the derivative of
with respect to the
ith component of the external force vector
. That is,
. 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
,
, or
will be referred to as a unit-load displacement in the rest of this paper.
One can now differentiate Equation (6) with respect to
, a variable on the right-hand side of the
kth constraint, to obtain
or
where
is the solution to Equation (14) with a unit-load placed at the degree of freedom,
kp, while
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,
and
. For simplicity, only
is considered as an independent variable in the following derivation. The first-order derivative of Equation (6) is now obtained as
where the matrix product of the vector,
, which is a collection of
, defined in Equation (4)
and the first- and second-order terms on the right-side of Equation (16) can be expressed in detail as follows:
and
With the help of the above equations, Equation (16) can be reorganized as
The solution of the above equation, which is the first-order displacement derivative with respect to
, is now obtained as a linear combination of two unit-load solutions as
The second-order displacement derivative with respect to
can be carried out in a similar matter, which again yields a linear combination of the same two unit-load solutions as
The same process can be further extended to obtain the higher-order derivatives of
as
where the order of differentiation,
n, is higher than 2. In short, the derivatives of the displacement vector with respect to
in different orders can be obtained in terms of two unit-load displacements,
and
.
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,
, will be linear in
as
while the derivative of
R with respect to
is found to be
The derivatives of
R with respect to the coefficients of the constraint,
and
, can also be found by the direct differentiation of Equation (7). Again, for illustration, one may select
as an independent variable to obtain the following relations,
and
Furthermore, the high-order derivatives of the reaction force can also be found as
Note that the first two coefficients of and on the right-hand side of are the same as those of and of in Equation (18). The same similarity was also observed between and in Equation (19) and their higher-order derivatives.
In summary, the
nth derivatives of the displacement vector with respect to either
or
as the independent variable can be briefly expressed as, based upon Equations (18)–(20),
where
and
are the unit-load solutions of Equation (14), and
and
are scalar functions of random variables,
or
, and the lower-order derivatives,
, and
, which should be available at the time when
is computed. The derivatives of the constraint reaction force vector can also be expanded in terms of
and
in a similar matter. Specifically, one has the following relation,
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,
where
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,
, which collects all axial stresses
reported at the interior points. The value of the axial stress
can be recast as
where
is the first row of
D. It is noted that the stress,
is a linear combination of unit-loads. Thus, the smoothing process can be conducted separately to map
from interior points to nodes as
where
and
are solved for the states of the stresses at the nodal points produced by the unit-loads,
and
, 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
and shear stress
. As a result, the
nth order derivatives of the state of the stresses at a nodal point,
, 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
It should be noted that the above equations are derived based upon the assumption that the random variable of concern, , is associated with only one constraint, 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
as a random variable with the mean and standard deviation,
and
. 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
and
where
in Equations (28) and (29) is the nodal von Mises stresses node of concern, and
and
are the first- and second-order derivatives of nodal state of stresses as defined in Equations (22) and (23) with
and
, respectively. All terms in the right-hand sides of Equations (28) and (29) are evaluated at the mean of the random variable,
. 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,
that can minimize the objective,
f, and satisfy the given and inequality constraints,
, and equality ones,
. Specifically, the mathematical expression of a generalized optimum design problem can be stated in a deterministic form as follows,
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,
, must now be measured by a targeted level of failure probability,
where
is the total number of element nodes, and
represent the yield stress and the von Mises stress, respectively. The probability failure of the stress constraint,
, can also be stated by a complementary term, called reliability,
If the randomness follows the normal distribution, the above reliability can be calculated by the CDF of a single variable, called the reliability index,
, as
which is a ratio between the mean and the standard deviation of all possible values of
. That is,
The inequality constraint, Equation (26), can then be recast as a single value constraint of the reliability indices,
where the value of the reliability index,
, is associated with the targeted reliability,
As for the equality constraints in Problem (P1), its uncertainty can be controlled by imposing tolerances. For example, one can require the constraint,
fall into a range of a specified bilateral tolerance,
, with a pre-determined probability,
as
ssuming that the distribution is normal, one may set the mean value of the constraint as zero,
, and the standard deviation of
as follows, based upon the required tolerance,
and the probability,
, as
where the reliability index,
, is determined based upon the specified probability,
, as
For instance, the constraint that falls into the range
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
Note that in this case, or 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
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,
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
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
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,
are then obtained as
As presented in
Section 2.2, the response functions,
, which include the nodal displacements, reactions, and von Mises stresses of a linear structure, and their derivatives,
and
, evaluated at the means, can be recursively expressed as a linear combination of the unit-load solutions, denoted by
and
. 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
Note that the mean value of any random variable, say
, can be selected as a design variable (i.e.,
. Since function
and its derivatives,
and
, in Equations (35) and (36), are all evaluated at the means of the random variables, one can follow the standard procedure to differentiate
,
, and
with respect to
, as per Equations (35) and (36) with respect to
. In the case where the standard deviation of any random variable,
, is treated as a design variable, Equations (35) and (36) can now be simplified as
and
This is again because , and are evaluated at the means of all random variables, and consequently, they are independent of the standard deviations.