Dynamic Analysis of Beams with Interval Parameters

Abstract

The present study deals with the transient interval analysis of a shallow beam having uncertainty in structural parameters viz. mass density and applied load. It is quite difficult to obtain information regarding the exact values of these parameters in several practical situations. Use of precise (deterministic) values of structural parameters in such a situation leads to erroneous results as the mathematical model built using deterministic structural parameters does not account for the uncertainty present in the system. In the present work, uncertainty present in the system is represented by interval parameters. In the research work carried out in the past quarter century, several methods were developed to model structural response of uncertain structural systems subjected to static loads under conditions of linear elasticity. The partial differential equations of motion of a Euler-Bernoulli beam are solved using Finite difference and finite element methods under conditions of linear elasticity. The resulting interval equations are solved using search and direct methods. Further, direct optimization approach is used to compute the bounds of displacement. The applicability and effectiveness of presented methods is demonstrated by solving example problems.

1. Introduction

Analysis of vibration of structures such as wings of aero planes, tall buildings and bridges is very important to ensure safe design. The issue of structural safety is of utmost importance in the field of structural engineering. This is to ensure the safety and security of structures so that unexpected and disastrous failures causing significant loss of life and property can be avoided. It is very important to conduct safety studies to avoid failure during the lifetime of a structure [1]. Given the tremendous strides made in enhancement of computational capabilities, it has now become possible to prepare detailed numerical models of structures for analysis and design. These mathematical models are built based on several structural parameters that may be subject to uncertainty. Thus, it is essential to build mathematical models that represent uncertainty in structural parameters so that the mathematical models can predict uncertainty in the behavior of structural systems. This brings out the importance of risk/reliability analysis of structural systems.

All structures possess physical and geometrical uncertainties due to physical imperfections, model inaccuracies, and system complexities [2]. Lack/scarcity of data is the main reason that causes uncertainty in the parameters of a structural system. For example, concrete undergoes progressive deterioration due to age and steel is subject to corrosion with the passage of time. These phenomena result in changes in the values of structural parameters and cause uncertainty. Also, uncertainty is associated with the determination of external loads acting on structures.

To begin with, the concept of interval uncertainty is discussed [3]. To describe the uncertainty of a structural parameter $\tilde{r}$, it is necessary to know the error in measurement $\Delta r=\left|\tilde{r}-r\right|$. Usually, the upper bound of this error $\Delta$ is known such that $\Delta r\le \Delta$. It is possible to represent all the possible values that represent the uncertainty in the parameter $\tilde{r}$ as an interval value $r=\left[\tilde{r}-\Delta ,\tilde{r}+\Delta \right]$. It is guaranteed that this interval value contains the exact value of $\tilde{r}$. However, it is to be noted that the interval values of parameters must be carefully handled using the rules of interval algebra so that a result that can guarantee to bound an exact solution can be obtained.

Initial work in interval algebra was published by Moore [3] in the year 1966. However, interval methods were applied to the problems of structural analysis only during the past two decades. This period saw major advances in the areas of design based on the concepts of structural reliability. Uncertainty was modelled under conditions of limited availability of data. Several researchers developed methods to use interval variables to incorporate uncertainty into response of structure. Right from the development of Interval Finite Element Methods in the mid-1990s, numerous researchers have explored various facets of this approach. Key contributors include Alefeld and Herzberger [4], Rao and Sawyer [5], Nakagiri and Yoshikawa [6], Rao and Berke [7], Koyluoglu and Elishakoff [8], Rao and Chen [9], and Muhanna and Mullen [10].

Researchers have particularly focused on two major challenges: obtaining solutions for the resulting linear interval equations with practical bounds on the structural response to minimize overestimation and deriving realistic bounds on quantities that are derived from structural response. For instance, when structural response involves displacement, the derived variables might include forces or stresses, which depend on displacements. Achieving tight bounds on these derived variables is quite challenging due to their dependency on the primary variables, which are already overestimated.

In the past decade, significant advancements have been made in applying interval FEM to problems of structural mechanics. Researchers such as Zhang [11], Muhanna, Mullen and Rama Rao [12] have made notable contributions. These studies have primarily focused on computing structural responses under uncertainty in loads, stiffness, and element cross-sectional areas. Neumaier and Pownuk [13] introduced a notable improvement in IFEM formulation for truss problems. They presented an iterative method for computing tight bounds on the solutions of linear interval systems, achieving a frequently small amount of overestimation. Despite these improvements, the two-step approach in their formulation still results in additional overestimation when evaluating derived quantities. Skalna et al. [14] presented various methods for solving systems of fuzzy equations in structural mechanics. Skalna [15] presented a direct method to solve interval linear systems with non-affine dependencies. Zou et al. [16] worked on the prediction of bounds of Vehicle-Induced Bridge Responses using the Interval analysis method. Xiao, Fedele and Muhanna [17] presented an interval-based approach for parameter identification in structural static problems. Muhanna and Shahi [18] worked on the problem of inclusion of uncertainty in boundary conditions. Misraji et al. [19] worked on the linear structural systems subject to stochastic Gaussian loading. Xiao et al. [20] presented a new decomposition of the stiffness matrix that yields a significant reduction of overestimation. The IFEM equations are solved utilizing a new variant of the iterative enclosure method to obtain the tightest guaranteed enclosure of the exact solution. Mullen and Muhanna [21] investigated the effect of multiplicative uncertainty (present at both element and group levels in the elements of structure) on structural behavior. Pownuk [22] used adaptive approximation for the interval solution of the equations of dynamics. Callens et al. [23] presented a computationally efficient technique for interval analysis and linear models using the Multilevel Quasi-Monte Carlo framework. Betancourt and Muhanna [24] developed a novel interval deep neural network (DINN) for providing reliable input data to mechanical models from uncertain data.

Thus, over the past three decades, various methods have been developed to model the structural response of uncertain systems under static loads within the framework of linear elasticity. However, research on the response of such systems to dynamic loads remains limited. This paper aims to address this gap by predicting the transient structural response of uncertain structural systems. The concept of uncertainty is discussed first, followed by a brief description of research work carried out by several researchers in applications of interval methods to structural analysis, followed by a brief description of the work carried out by the authors.

The present research investigates the transient dynamic behavior of a Euler-Bernoulli beam with uncertain structural parameters, specifically mass density and applied load. Typically, mathematical models rely on fixed values for parameters, but in practice, obtaining precise information about these values can be challenging or impossible. Using deterministic values in such cases can lead to inaccurate results, as the model fails to account for inherent uncertainties. To address this, various methods can be employed to model uncertainty, including interval parameters, which are particularly useful when only upper and lower bounds are available. This study represents uncertainty using interval parameters, ensuring that the dynamic response of the system reflects this uncertainty. While previous research has focused on static loads and linear elasticity, this paper contributes to the understanding of transient responses in uncertain structural systems under dynamic loads and non-linear elasticity. The analysis uses finite difference and finite element methods to solve the partial differential equation for the motion of a straight elastic beam. To compute the time-history response under uncertainty, interval dynamic beam equations are solved using search, direct, and optimization methods. The effectiveness of these methods is demonstrated through numerical examples with interval-valued mass density and dynamic loading, and the results are presented.

In the present work, vibration analysis of a Euler-Bernoulli beam under conditions of linear elasticity with interval mass density and interval load is considered. The response of the beam with interval structural parameters is obtained by the solution of interval dynamic beam equations are solved using three methods viz. the search method, direct method, and direct optimization. The effectiveness of these methods is demonstrated through examples involving beams with Young’s modulus, mass density and applied load being interval quantities. The response of the beam is computed in the example problems, with results presented to illustrate the methods’ applicability and effectiveness.

When employing the first two approaches, specifically search and direct methods [14,15], the beam’s governing equation is discretized using finite differences in both time and space. This discretization yields a set of equations with interval coefficients, which are then solved to obtain the beam’s displacement at each time step. These methods leverage the monotonic relationship between structural response and parameters, enabling the determination of response bounds at interval endpoints. Monotonicity can be verified using Taylor series or interval methods [14,15]. While methods exploiting monotonicity are efficient for large problems, they may yield conservative estimates. Alternatively, parametric interval linear systems can be solved to compute the solution without relying on monotonicity [15]. These methods provide a guaranteed enclosure of the solution, but their applicability is limited by factors such as the degree of uncertainty. In the third method, the structure is modelled using the Finite Element analysis. Wilson-method and optimization approach are used to compute the structural response [25,26].

The paper has the structure outlined as follows. Concepts of interval uncertainty are introduced in Section 2. Mathematical theory related to the forced vibration of a beam is presented in Section 3. The process of discretization of the beam using finite differences is described in Section 4. The process of discretization of the beam and implementation of finite element method to the vibration problem is outlined in Section 5. Section 6 provides the description of methods for the solution of interval linear structural systems. Section 7 presents the Wilson $\theta$ method for calculating the transient vibration response of a structure to sudden loading. Numerical examples are solved in Section 8. Conclusions are given in Section 9.

2. Interval Uncertainty

In case only a limited knowledge of value of a uncertain parameter is available, then this value can be defined by an interval variable follows:

$$r_i\in \left[{\tilde{r}}_i-\Delta r_i,{\tilde{r}}_i+\Delta r_i\right]=\left[{\underset{\_}{r}}_i,{\overline{r}}_i\right]={\mathit{r}}_i$$

(1)

where ${\stackrel{⌢}{r}}_i$ is taken as an approximation of the exact value of uncertain parameter $r_i$. Also, $\Delta r_i$ is an estimate of error. In this paper, interval parameters are represented by bold face font.

Let us consider an quantity y to be related to parameters $r$ by a function $z=g\left(r\right)$, then the process of computation of the outcome, assuming that $r$ varies within $\mathit{r}$, is comparable to computing solution given by:

$$z_S=\left\{z:z=g\left(r\right),r\in \mathit{r}\right\}$$

(2)

Here $z_S$ represents a set of possible solutions of the interval analysis. This set cannot be described by a hypercube or an interval. We can interpret that $z_S$ comprises all vectors $z$ that can be computed from the application of function $g$ on all vectors $r$ bounded by $\mathit{r}$.

It is often difficult to provide an accurate description of the solution set $z_S$. Typically, an interval vector $x^{*}\in z_S$, termed outer enclosure to the solution, is computed instead. The objective is to make $x^{*}$ as much narrow as feasible. This will result in a hull solution that denotes the narrowest interval vector containing $z_S$. It is also possible to compute the inner enclosure to the solution. The inner enclosure to the solution is defined as an interval enclosed by the hull of the solution. However, the use of straightforward interval arithmetic results in overestimation due to the problem of dependency [1]. The excess width of the solution can be decreased by keeping track of the dependency of intermediate results on input data. There are several research papers that successfully handle the dependency problem to obtain narrow bounds to the solution [14,15].

3. Forced Vibration

The mathematical model of the forced vibration of a Euler-Bernoulli beam (Figure 1) is given by

$${\displaystyle \frac{{𝜕}^2}{𝜕x^2}}\left(EI{\displaystyle \frac{{𝜕}^{2}w}{𝜕x^2}}\right)=q-\rho A{\displaystyle \frac{{𝜕}^{2}w}{𝜕t^2}}$$

(3)

where $q$ is an external load, material properties of the beam are given by mass density $\rho$ and Young’s modulus $E$. $A$ is the cross-sectional area, $I$ is the moment of inertia. In Equation (3) displacement $w$ depends solely on variables $x$ and $t$. This model represents the behavior of shallow beams under small transverse deformations. For a uniform beam with constant value of flexural rigidity $EI$, Equation (3) can be simplified as:

$$EI{\displaystyle \frac{{𝜕}^{4}w}{𝜕x^4}}=q-\rho A{\displaystyle \frac{{𝜕}^{2}w}{𝜕t^2}}$$

(4)

Vibration of a beam is an initial-boundary value problem that involves specification of both boundary and initial conditions to compute displacement $w\left(x,t\right)$. The problem requires four boundary conditions and two initial conditions since the second order derivative with respect to time and fourth order derivative w.r.t space are present in Equation (4).

$$\begin{array}{ll}w\left(0,t\right)=w\left(L,t\right)=0& {\displaystyle \frac{{𝜕}^{2}w}{𝜕x^2}}\left(0,t\right)={\displaystyle \frac{{𝜕}^{2}w}{𝜕x^2}}\left(L,t\right)=0\\ w_0\left(x\right)=w\left(x,0\right)& v\left(x,0\right)={\displaystyle \frac{𝜕w}{𝜕t}}\left(x,0\right)=v_0\left(x\right)\end{array}$$

(5)

Considering ends of the beam as simply supported, the transverse displacements at the end points are zero. This is expressed as $w\left(0,t\right)=w\left(L,t\right)=0$ for $0\le t\le T$. Owing to the given end conditions, bending moments at endpoints are zero. Thus,

$$EI{\displaystyle \frac{{𝜕}^{2}w}{𝜕x^2}}\left(0,t\right)=M_x\left(0,t\right) \mathrm{and} EI{\displaystyle \frac{{𝜕}^{2}w}{𝜕x^2}}\left(L,t\right)=M_x\left(L,t\right)$$

and thus

$$\mathrm{For} 0\le t\le T,{\displaystyle \frac{{𝜕}^{2}w}{𝜕x^2}}\left(0,t\right)={\displaystyle \frac{{𝜕}^{2}w}{𝜕x^2}}\left(L,t\right)=0$$

For $t$ is equal to zero, both initial displacement and initial velocity are zero. This means:

$$\mathrm{for} 0\le x\le L, w_0\left(x\right)=0 \mathrm{and} v_0\left(x\right)=0.$$

4. Modelling of the Beam Using Finite Difference Method

The procedure for solving the problem of dynamic beam vibrations using the finite difference method is outlined in the following steps [14,15]:

(a)

The beam is discretized in space into elements as shown in Figure 2 into $N$ elements with $N+1$ nodes. Discretization is space domain is denoted by index $i$. The number of points in the space domain are $N+1$ such that $1\le i\le N+1$

(b)

Time domain $0\le t\le T$ is discretized into $M$ number of time steps $\delta t$. Discretization in the time domain is indicated by index $j$. The number of points in the time domain are $M+1$ such that $1\le j\le M+1$.

(c)

So, the total number of points in the grid formed in the domain of space-time is $\left(N+1\right)\left(M+1\right)$.

(d)

Equation (4) is discretized at the point $\left(i,j+1\right)$ in the above grid as follows:

$${\left(EI{\displaystyle \frac{{𝜕}^{4}w}{𝜕x^4}}\right)}_{i,j+1}=q_{i,j+1}-{\left(\rho A{\displaystyle \frac{{𝜕}^{2}w}{𝜕t^2}}\right)}_{i,j+1}$$

(6)

(e)

Finite differences are substituted in Equation (6) to obtain:

$$\begin{array}{l}{E}_{i,j+1}{I}_{i,j+1}\left({\displaystyle \frac{w_{i+2,j+1}-4w_{i+1,j+1}+6w_{i,j+1}-4w_{i-1,j+1}+w_{i-2,j+1}}{\Delta x^4}}\right)\\ +\left({\displaystyle \frac{{\rho }_{i+1.j+1}{A}_{i+1.j+1}}{\Delta t^2}}\right)w_{i,j+1}=q_{i,j+1}-{\rho }_{i.j+1}{A}_{i.j+1}\left({\displaystyle \frac{2w_{i,j}-w_{i,j-1}}{\Delta t^2}}\right)\end{array}$$

(7)

(f)

Equation (7) represents the set of finite difference equations that are applicable over the points of the grid specified in step 3.

(g)

On similar lines, initial and boundary conditions can also be discretized as follows:

$$\{\begin{array}{c}{w}_{0,j+1}=0\\ w_{0,j+1}-2w_{1,j+1}+w_{2,j+1}=0\\ \begin{array}{l}{E}_{i,j+1}{I}_{i,j+1}{\displaystyle \frac{w_{i+2,j+1}-4w_{i+1,j+1}+6w_{i,j+1}-4w_{i-1,j+1}+w_{i-2,j+1}}{\Delta x^4}}\\ +{\displaystyle \frac{{\rho }_{i+1.j+1}{A}_{i+1.j+1}}{\Delta t^2}}{w}_{i,j+1}=q_{i,j+1}-{\rho }_{i.j+1}{A}_{i.j+1}{\displaystyle \frac{2w_{i,j}-w_{i,j-1}}{\Delta t^2}}\end{array}\\ w_{n-2,j+1}-2w_{n-1,j+1}+w_{n,j+1}=0\\ w_{n,j+1}=0;\\ w_{i,0}=w_i^{*}\\ w_{i,1}=w_{i,0}+v_i^{*}\Delta t\end{array}$$

(8)

It shall be noted that $w_{i,j+1}=w\left(r_{i,j+1}\right)$ where $r_{i,j+1}=\left(E_{i,j+1},{\rho }_{i,j+1},q_{i,j+1}\right)$.

(h)

Finite difference equations given by Equation (7) are applied at each point on the grid formed in the space-time domain, subject to the boundary conditions specified in Equation (8) in finite difference form.

(i)

Thus, the process of discretization does reduce the problem of calculating the transient displacement of the structure to the problem of obtaining solution to a set of parametric linear equations.

(j)

These linear interval equations need to be solved, considering parametric uncertainty. The process is described in Section 6.

(k)

To enhance the precision of the finite difference method, a finite difference of order > 3 is applied to the time step.

$${\left({\displaystyle \frac{{𝜕}^{2}w}{𝜕t^2}}\right)}_{i,j}={\displaystyle \frac{2w_{i,j}-5w_{i,j-1}+4w_{i,j-2}-w_{i,j-3}}{\Delta t^2}}$$

(9)

5. Modelling of the Beam Using Finite Element Method

The procedure for solving the problem of dynamic beam vibrations using the finite element method is outlined in the following steps [27]:

(a)

The beam is discretized in space into elements as shown in Figure 2 into $N$ elements with $N+1$ nodes. Length of each element is $l_e$

(b)

The stiffness matrix of each element is given by

$$\left[k^{(e)}\right]=\left[\begin{array}{cccc}{\displaystyle \frac{12EI}{l_e^3}}& {\displaystyle \frac{6EI}{l_e^2}}& -{\displaystyle \frac{12EI}{l_e^3}}& {\displaystyle \frac{6EI}{l_e^2}}\\ {\displaystyle \frac{6EI}{l_e^2}}& {\displaystyle \frac{4EI}{l_e}}& -{\displaystyle \frac{6EI}{l_e^2}}& {\displaystyle \frac{2EI}{l_e}}\\ -{\displaystyle \frac{12EI}{l_e^3}}& -{\displaystyle \frac{6EI}{l_e^2}}& {\displaystyle \frac{12EI}{l_e^3}}& -{\displaystyle \frac{6EI}{l_e^2}}\\ {\displaystyle \frac{6EI}{l_e^2}}& {\displaystyle \frac{2EI}{l_e}}& -{\displaystyle \frac{6EI}{l_e^2}}& {\displaystyle \frac{4EI}{l_e}}\end{array}\right]$$

(10)

(c)

Lumped mass matrix for each element is given by

$$\left[m^{(e)}\right]=diag\lfloor \begin{array}{cccc}{\displaystyle \frac{\overline{m}{l}_e}{2}}& 0& {\displaystyle \frac{\overline{m}{l}_e}{2}}& 0\end{array}\rfloor$$

(11)

(d)

The stiffness matrix $\left[K\right]$ and mass matrix $\left[M\right]$ of the structure are assembled from the element stiffness and mass matrices $\left[k^{(e)}\right]$ and $\left[m^{(e)}\right]$ respectively.

(e)

Damping matrix $\left[C\right]$ of the structure is expressed as a linear combination of $\left[K\right]$ and $\left[M\right]$, considering Rayleigh damping. The procedure is explained in Section 7.

(f)

Element force vector, due to loads carried by the element, is also computed as follows:

$$\left\{f^{(e)}\right\}={\lfloor \begin{array}{cccc}-{\displaystyle \frac{q{l}_e}{2}}& -{\displaystyle \frac{q{l}_e^2}{12}}& -{\displaystyle \frac{q{l}_e}{2}}& {\displaystyle \frac{q{l}_e^2}{12}}\end{array}\rfloor }^T$$

(12)

(g)

The force vector acting on the structure $\left\{F\left(t\right)\right\}$ is computed by the assembly of element force vectors. In addition, any loads acting at the nodes of the structure are also added to the structure force vector.

(h)

For the given beam, the stiffness and mass matrices and force vector of the structure are assembled.

(i)

Boundary conditions are applied on the resulting system of matrix equations.

The mass, stiffness and damping matrices and the force vector of the structure obtained following the steps above define the matrix equation for the vibration of the multi-degree of freedom (MDOF) system for which solution is obtained using Wilson-$\theta$ method. The detailed procedure is outlined in Section 7.

6. Solution of Parametric Linear Equations

There are two main approaches for addressing interval uncertainty. The interval arithmetic strategy seeks to enclose the interval result exactly within a hypercubic boundary from the outside, relying on the computation of outer bounds. In contrast, the global optimization method focuses on computing an inner approximation. While interval arithmetic methods are generally less computationally demanding, the resulting interval solution gets overestimated due to the problem of dependency. Conversely, although optimization-based approaches are more computationally intensive and time-consuming, they offer more acceptable solutions for real-world engineering problems.

6.1. Solution of Interval Equations

The interval solution is a function of well-defined parameter combinations, which are explicitly stated.

$${\underset{\_}{w}}_{i,j}=w_{i,j}\left(p_{i,j}^{\min }\right), {\overline{w}}_{i,j}=w_{i,j}\left(p_{i,j}^{\max }\right)$$

(13)

In the continuous case, one can write,

$$\begin{array}{l}\underset{\_}{w}\left(x,t\right)=w_{i,j}\left(x,t,r^{\min }\left(x,t\right)\right), \\ \overline{w}\left(x,t\right)=w_{i,j}\left(x,t,r^{\max }\left(x,t\right)\right)\end{array}$$

(14)

The interval solution sometimes depends only on a given combination of parameters [9] for a domain given by $D_{\alpha }\in \left[0,L\right]\times \left[0,T\right]$

$$\begin{array}{l}\underset{\_}{w}\left(x,t\right)=w_{i,j}\left(x,t,r_{\alpha }^{\min }\right), \\ \overline{w}\left(x,t\right)=w_{i,j}\left(x,t,r_{\alpha }^{\max }\right)\end{array}$$

(15)

Then one can compute the interval solution exactly using finite combinations of parameters $r_1^{\min },r_1^{\max },\dots r_{\alpha }^{\min },r_{\alpha }^{\max },\dots r_p^{\min },r_p^{\max }$ where $D_1\cup D_2\dots \cup D_p=\left[0,L\right]\times \left[0,T\right]$.

6.1.1. Search Method

Search method is used to compute the solution. The method attempts to solve interval linear system of equations that correspond to the specified combination of parameters. Thus, every interval parameter ${\mathit{r}}_i$ is substituted by:

$${\mathit{r}}_i\approx \left\{r_{i,1},\dots r_{i,k}\right\}$$

(16)

An interval $\mathit{r}=\left[{\mathit{r}}_1,{\mathit{r}}_2,\dots ,{\mathit{r}}_m\right]$ is represented by the discrete set of points given by:

$$\mathit{r}=\left[{\mathit{r}}_1,{\mathit{r}}_2,\dots ,{\mathit{r}}_m\right]\approx \left\{\left(r_{1,i1},\dots r_{1,ip}\right):0\le i_1,\dots i_p\le p\right\}$$

(17)

The equation contains $p^m$ elements, where $m$ represents the number of interval parameters and $p$ signifies the number of points within each interval. Notably, when $p$ is equal to 2, the method becomes equivalent to the end-points combination method. In this scenario, the solution is calculated as follows:

$${\underset{\_}{w}}_{i,j}\approx \min \left\{w_{i,j}\left(p_{1,i1},\dots ,p_{m,ip}\right):\left(p_{1,i1},\dots ,p_{m,ip}\right)\in {\mathrm{P}}_{m,p}\right\}$$

(18)

$${\overline{w}}_{i,j}\approx =\max \left\{w_{i,j}\left(p_{1,i1},\dots ,p_{m,ip}\right):\left(p_{1,i1},\dots ,p_{m,ip}\right)\in {\mathrm{P}}_{m,p}\right\}$$

(19)

The search method lets us compute approximate values of $p_{i,j}^{\min .search},p_{i,j}^{\max .search}$ such that

$${\underset{\_}{w}}_{i,j}\approx w_{i,j}\left(r_{i,j}^{\min ,search}\right),{\overline{w}}_{i,j}\approx w_{i,j}\left(r_{i,j}^{\max ,search}\right)$$

(20)

According to numerical experiments,

$$r_{i,j}^{\min }\approx r_{i,j}^{\min ,search},r_{i,j}^{\max }\approx r_{i,j}^{\max ,search}$$

(21)

which implies that Search method can compute precise values of $p_{i,j}^{\min },p_{i,j}^{\max }$.

6.1.2. Direct Method

To validate the results computed through the methods discussed, the direct method (DM) developed by Skalna [15] is employed to solve parametric interval linear systems. Affine arithmetic is utilized to manage nonlinear dependencies.

7. Wilson-θ Method

The matrix system of equations developed in Section 5 are utilized in the implementation of Wilson-$\theta$ method. Vibration of the beam is modelled using a multi-degree of freedom (MDOF) system defined as

$$\left[M\right]\left\{\ddot{w}\right\}+\left[C\right]\left\{\dot{w}\right\}+\left[K\right]\left\{w\right\}=\left\{F\left(t\right)\right\}$$

(22)

The stiffness matrix $\left[K\right]$, mass matrix $\left[M\right]$ and force vector $\left\{F\left(t\right)\right\}$ of the structure given in Equation (22) are computed using the procedure outlined in Section 5. The damping matrix C is computed as

$$\left[C\right]={\alpha }_0\left[M\right]+{\alpha }_1\left[K\right]$$

(23)

Here, Rayleigh damping with damping ratios ${\xi }_1$ and ${\xi }_2$ in the first and second modes of vibration. These modes correspond to frequencies ${\omega }_1$ and ${\omega }_2$ respectively. The coefficients ${\alpha }_0$ and ${\alpha }_1$ are computed as shown below:

$$\left[\begin{array}{c}{\alpha }_0\\ {\alpha }_1\end{array}\right]={\displaystyle \frac{1}{2}}\left[\begin{array}{cc}{\displaystyle \frac{1}{{\omega }_1}}& {\omega }_1\\ {\displaystyle \frac{1}{{\omega }_2}}& {\omega }_2\end{array}\right]\left[\begin{array}{c}{\xi }_1\\ {\xi }_2\end{array}\right]$$

(24)

The problem of transient vibration of beam is solved using Wilson-θ method which follows implicit integration [23]. This method computes transient response of a multi degree of freedom systems by adopting a stepwise integration procedure. The Wilson-θ method considers acceleration to vary linearly across the time interval $\left[t,t+\theta \Delta t\right]$, where $\theta \ge 1$ and $\Delta t$ is the timestep. This method turns out to be stable for θ ≥ 1.38, as per the work done by Wilson.

Optimization Approach

In this approach, Interval uncertainty is present in mass density and Young’s modulus as well as applied load. An optimization scheme is used to compute solution to the resulting interval multi-degree of freedom system. This is accomplished by using the function fmincon that belongs optimization toolbox of MATLAB [24]. The function computes a constrained minimum of an objective function by the solution of:

$$\left\{x\right\}=fmincon\left(objfunc,\left\{x_0\right\},\left[A\right],\left[B\right],\left[A_{eq}\right],\left\{b_{eq}\right\},\left\{lb\right\},\left\{ub\right\}\right)$$

(25)

subject to the inequality constraints

$$\left[A\right]\left\{x\right\}\le \left\{b\right\}$$

(26)

and equality constraints

$$\left[A_{eq}\right]\left\{x\right\}=\left\{b_{eq}\right\}$$

(27)

and bounds

$$\left\{lb\right\}\le \left\{x\right\}\le \left\{ub\right\}$$

(28)

where $x_0$ is the starting point for search.

Equation (28) describes bounds on the variables $\left\{x\right\}$ to find a solution subject to $\left\{lb\right\}\le \left\{x\right\}\le \left\{ub\right\}$.

The interval solution is computed as a function of time by extending the Wilson-θ method. This is done by formulating objfunc to yield a unique output that describes the deterministic value of transient nodal displacement. Optimization of displacement response is carried out to compute bounds for the displacement response at every time step $\Delta t$ for 0 ≤ t ≤ tmax. The normalized uncertainties associated with mass and stiffness and load terms are represented by normalized intervals quantities ${\mathit{r}}_1$, ${\mathit{r}}_2$ and ${\mathit{r}}_3$ respectively. The upper and lower bounds of these interval parameters form the vertices of an n-dimensional hypercube ${\mathit{r}}^I$. A point $\left(r_1,r_2,r_3\right)$ lying inside the bounds $\left\{lb\right\}$ and $\left\{ub\right\}$ (that are described in Equation (25)) are defined as

$$\left\{lb\right\}={\lfloor \begin{array}{ccc}{\underset{\_}{r}}_1& {\underset{\_}{r}}_2& {\underset{\_}{r}}_3\end{array}\rfloor }^T,\left\{ub\right\}={\lfloor \begin{array}{ccc}{\overline{r}}_1& {\overline{r}}_2& {\overline{r}}_3\end{array}\rfloor }^T$$

(29)

Equation (22) is recast in interval parametric form as

$${\mathit{r}}_1\left[M\right]\left\{\ddot{w}\right\}+\left[C\right]\left\{\dot{w}\right\}+{\mathit{r}}_2\left[K\right]\left\{w\right\}={\mathit{r}}_3\left\{F\left(t\right)\right\}$$

(30)

where parameters ${\mathit{r}}_i$ are defined as

$${\mathit{r}}_i=\left[{\underset{\_}{r}}_i,{\overline{r}}_i\right],\left(i=1,2,3\right)$$

(31)

The objective function is evaluated at a point $r$ given by coordinates $\left(r_1,r_2,r_3\right)$ within the hypercube ${\mathit{r}}^I$ which forms the search domain. Thus, the deterministic procedure is rendered as an interval procedure using an approach built on global optimization. In this approach, the lower bound and upper bound of interval displacement ${\mathit{w}}_{\mathit{m}}$ at a particular node $m$ are computed, taking into consideration that the uncertain parameters $r$ can vary within the interval ${\mathit{r}}^I$. This interval ${\mathit{w}}_{\mathit{m}}$ displacement is computed by a minimization and a maximization over the uncertainty interval ${\mathit{p}}^I$.

$${\mathit{w}}_m=\left[\underset{p\in p^I}{\min }\left(w_m\right),\underset{p\in p^I}{\max }\left(w_m\right)\right]$$

(32)

This is carried out by calculating the displacement vector $\left\{w\right\}$ at a distance x along the beam at time $t$, using the equation given below, by applying Wilson-θ method. This equation is deterministic in nature.

$$p_1\left[M\right]\left\{\ddot{w}\right\}+\left\{C\right\}\left\{\dot{w}\right\}+p_2\left[K\right]\left\{w\right\}=p_3\left\{F\left(t\right)\right\}$$

(33)

The displacement response’s time history is calculated by determining the displacement bounds at each time step. To compute the displacement at a specific time, all preceding time steps must be evaluated. However, to mitigate the computational burden of local optimization with fmincon, a database stores function evaluations of the objective function, allowing the optimizer to reuse them in subsequent optimizations. This approach enables efficient optimization initialization from the optimal function value found thus far, either the minimum or maximum. The optimization procedures are performed sequentially, from the lowest to the highest level, ensuring a comprehensive search for the optimal solution.

8. Numerical Examples

Numerical examples are solved to illustrate the applicability of the methods described in the paper. The first example involves a beam subjected to a uniform load distribution along its full span. The study presents findings for two different cases of uncertainty, highlighting the impact of uncertain parameters on the beam’s response.

Example 1.

The beam presented in Figure 1 is considered with a distributed load of 2.5 kN per meter for a time of $9\times {10}^{-3}$ s. The span of the beam is 4 m, cross-sectional area $0.01 m^2$, area moment of inertia $8:333\times {10}^{-6} m^4$. Young’s modulus and mass density of the material of the beam are 200 GPa and 7850 kg/m³ respectively. The interval uncertainties in mass density and load are taken respectively as ±0.5% and ±20% from their mean values. One should note that the uncertainty of mass density of the beam is usually small. However, the uncertainty of dynamic load acting of the beam can be quite large. This gives two interval parameters $\mathit{r}=\left(r_1,r_2\right)=\left(\rho ,q\right)$. The beam is discretized into 20 divisions ($n_x$) so that each division $\Delta x$ is equal to $L/n_x$, number of time steps $nt$ is 100 and time step $\Delta t$ is 0.0015 s. The load is applied for a duration of 0.009 s. Figure 3 gives bounds of the interval displacement obtained using the search method when m equals 2 and k equals 7.

The results of Direct and Search Methods are compared in Figure 4. The solid black line in the figure represents the deterministic solution and coincides with the solution obtained using FEM. The inner and outer bounds of the solution shown in gray color represent the interval solution obtained using direct and search methods. Results obtained using the Direct method coincide with those obtained using Search method. This is because the structural behavior of the beam under vibration is monotonic in nature.

Example 2.

A beam cited as the first example is taken up for analysis another time. The beam is subjected to a sudden distributed load of 5 kilo Newtons per meter run over the entire span for a duration of 0.4 s. 5% damping is taken to be present. In the first case, mass density and Young’s modulus are taken to have a small uncertainty of ±0.5% about their mean values whereas load is taken to be deterministic i.e., ${\mathit{r}}_1={\mathit{r}}_2=\left[0.995,1.005\right]$ and ${\mathit{r}}_3=\left[1.0,1.0\right]$. This is because variation of mass density and Young’s modulus of a structure can be quite small as per the observations of other researchers [14].

The transient displacement of the midspan of the beam is calculated following the method given in Section 6.1. Figure 5 presents plot of displacement at mid-span. The figure represents an envelope of the vibration containing the lower and upper bounds of vibration response of the beam. It is observed that the plot is dark near the center and becomes gray towards the upper and lower bounds. The uncertainty of the response increases greatly with the increase in time.

Figure 6 illustrates the second case where Young’s modulus and mass density are precisely defined, whereas the load exhibits a ±20% interval uncertainty around its mean value, introducing a range of possible values. Figure 6 shows the plot of mid-span displacement that corresponds to ${\mathit{r}}_1={\mathit{r}}_2=\left[1.0,1.0\right]$, ${\mathit{r}}_3=\left[0.8,1.2\right]$ with a damping of 5 percent. The figure depicts an envelope encapsulating the upper and lower bounds of the beam’s vibration response. A notable observation from Figure 5 is the temporal shift in the response peaks, attributable to the escalating uncertainty in displacement over time. This phenomenon arises from the propagation of uncertainties in mass and stiffness, culminating in a broader range of uncertainty in the beam’s natural frequencies. In contrast, Figure 6 exhibits no such peak shifting, as the deterministic nature of stiffness and mass ensures that the natural frequencies remain invariant with time. This observation aligns with established findings in the literature [25].

Example 3.

Consider a cantilever beam of span $4 m$ as shown in Figure 7. The Young’s modulus of the material of the beam is 2 × 10⁶ N/m² and the cross-sectional area is 0.01 m². The beam is discretized 5 elements each of length 1 m. Mass matrix is computed using the lumped mass approach. The uncertainties of Young’s modulus and mass density are ±5 percent about their respective mean values. Rayleigh damping of 5 percent is considered. The beam is subjected to a sudden load of 10 kN in the vertically downward direction at its tip for a duration of 0.04 s. Now, the bounds of structural response for the vertical displacement of the tip of the cantilever beam are computed. Figure 8 shows the time history response of the bounds of vertical displacement of the tip of the cantilever beam. Here also shifting of peaks at various levels of uncertainty is noticed, similar to Figure 5.

9. Conclusions

This paper presents three approaches for the solution of vibration problems of beams with interval uncertainty. In the first example, the structure is discretized using the Finite difference approach and solution is obtained using search method and direct method. First, the resulting parametric interval equations are solved using the direct method that gives guaranteed bounds to the solution. Direct method can eliminate the overestimation from the interval computations. Search method obtains bounds to the interval solution by considering the solutions at the vertices of the hypercube formed by the endpoints of the interval parameters. Owing to the monotonic nature of the beam vibration problem, it is observed that in the case of vibration of beam under conditions of linear elasticity, the interval displacement only depends on the combinations of end-points of interval parameters. Thus, the results obtained using two methods coincide.

In the second example, a beam with interval values of load, stiffness and mass density is considered. The structure is discretized the Finite Element approach. An objective function that corresponds to the maximum displacement of the beam is taken up for optimization subjected to constraints. The MATLAB optimization toolbox is used for constrained optimization in the search space corresponding to the hypercube formed by the endpoints of the interval parameter. The lower and upper bounds of resulting interval values of displacement are obtained at each time step and are used to make a time history plot of structural response at several levels of uncertainty. It is seen that the process of obtaining results at each time step using Optimization approach takes time but yields satisfactory results even with wide intervals of input parameters.