A domain decomposition solver for spectral stochastic finite element: an approximate sparse expansion approach

Stochastic finite element: spectral approach

In the response analysis of a structural system with stochasticity in the model parameters, the equilibrium equation is usually expressed as

where $$ \mathbf{K}(\theta) $$ is the stochastic stiffness matrix related to the stochastic process $$ \omega(x, \theta ) $$ of the system, $$ \mathbf{u}(\theta) $$ is the random response vector, and $$ \mathbf{f} $$ is the external forcing vector.

In one particular formulation of the SFEM, according to the KL of the stochastic process $$ E(x, \theta) $$, the stochastic stiffness matrix, with the KL expansion^[1], can be given as

where $$ \mathbf{K}_{0} $$ is average stiffness, $$ \xi_{i} $$ is a set of standard Gaussian random variables, and $$ \{\mathbf{K}_{i} \} $$ are deterministic, which are functions of the eigenvalues and eigenfunctions; we called them stochastic part stiffness matrices in this work.

The resulting vector $$ \mathbf{u}(\theta) $$ of Equation (1), which is expanded along a PCE, is calculated by

where $$ \{\psi_{j} \} $$ are multidimensional orthogonal polynomials in $$ \{\xi_{j} \} $$, and each $$ \mathbf{u}_{j} $$ denotes a vector of deterministic coefficients. The number of terms $$ M $$ in Equation (3) depends on: (i) the number N of the uncorrelated random variables $$ \{\xi_{j} \} $$ used to describe the uncertainty; and (ii) the highest degree $$ p $$ of the $$ \{\psi_{j} \} $$ orthogonal polynomials. This term is calculated using the following formula^[25]:

Substituting Equations (2) and (3) into (1) and forcing the residual to be orthogonal to the approximating space spanned by the PC $$ \{\psi_{j} \} $$ yields the following system of linear equations^[36]:

These equations can be assembled into a matrix of size $$ M\cdot n \times M \cdot n $$:

where $$ \mathbf{K}^{(j, k)}= \sum\limits_{i=0}^{N} c_{i j k} \mathbf{K}_{i} \text {, and } c_{i j k}=\langle\xi_{i} \psi_{j} \psi_{k}\rangle $$

It is evident that the augmented system is an equation that couples the determinism of finite elements with the randomness of structures. The number of rows and columns of its coefficient matrix are the product of: (1) the spatial degrees of freedom (DOF) after meshing the structure; and (2) the number of expansion terms for random responses using the PC method. When the structural grid is subdivided densely and multiple terms are expanded for random responses, the size of the augmented system can become quite large, which exceeds the computational capacity of current computers. Ultimately, this fact makes it impossible for SSFE to analyze the stochastic response of large-scale structures. Therefore, it is quite necessary to develop efficient algorithms for the augmented system.

Domain decomposition approach for SSFE

When dealing with large-scale problems, DDM is one of the best candidates due to its proven performance and scalability. In this method, the computational domain $$ \Omega $$ is partitioned into $$ N_S $$ non-overlapping subdomains $$ \Omega=\bigcup_{s=1}^{N_{s}}\left\{\Omega_{s}\right\} $$; the random response vector of every typical subdomain s, $$ \mathbf{u}^{s}(\theta) $$, is divided into interior response vector $$ \mathbf{u}_{I}^{s}(\theta) $$ corresponding to the interior nodes and interface response vector $$ \mathbf{u}_{\Gamma}^{s}(\theta) $$ corresponding to interface nodes (shared by two or more adjacent subdomains)^[26,27].

According to this partition, the DDM numbers, first the subdomains and last the interface, that is, let $$ \mathbf{u}(\theta)=\{\mathbf{u}_{I}^{1}(\theta) \quad \mathbf{u}_{I}^{2}(\theta) \cdots \mathbf{u}_{I}^{N}(\theta) \quad \mathbf{u}_{\Gamma}(\theta)\} $$, at this time, the vector before any base after PCE of the random response vector $$ \mathbf{u}(\theta) $$ can be expressed as $$ \mathbf{u}_{i}=\{\mathbf{u}_{i, I}^{1} \quad \mathbf{u}_{i, I}^{2} \cdots \mathbf{u}_{i, I}^{N} \quad \mathbf{u}_{i, \Gamma} $$}. The new order of the vector of random response vector leads to an "arrow" pattern of every stochastic part stiffness matrix $$ \mathbf{K}_{i} $$ in Equation (2) as:

where the restriction operators $$ \mathbf{B}^{s} $$ acts as a scatter or gather operator between global and local components of the deterministic interface solution vectors^[31], expressed as $$ \mathbf{u}_{\Gamma}(\theta)=\left(\mathbf{B}^{s}\right)^{T} \mathbf{u}_{\Gamma}^{s}(\theta)\; {{\text {and }}} \mathbf{u}_{\Gamma}^{s}(\theta)=\mathbf{B}^{s} \mathbf{u}_{\Gamma}(\theta) $$.

Rearranging the augmented system Equation (5) based on the order of $$ \mathbf{u}(\theta) $$ will lead to an "arrow" pattern of the coefficient matrix, as given in

where the expansion restriction operators $$ \mathbf{B}_{S}^{s} $$ is obtained by the restriction operators $$ \mathbf{B}^{s} $$ as $$ \mathbf{B}_{S}^{s}=\operatorname{blockdiag}\left(\mathbf{B}^{s}, \ldots, \mathbf{B}^{s}\right), \; \mathbf{u}_{I}^{s}=\left\{\begin{array}{llll} \mathbf{u}_{0, I}^{s} & \mathbf{u}_{1, I}^{s} & \ldots & \mathbf{u}_{M-1, I}^{s} \end{array}\right\}^{T} $$ for $$ \alpha, \beta \in\{I, \Gamma\} $$, and the solution vector is split into $$ \mathbf{u}_{\Gamma}=\left\{\begin{array}{cccc} \mathbf{u}_{0, \Gamma} & \mathbf{u}_{1, \Gamma} & \ldots & \mathbf{u}_{M-1, \Gamma} \end{array}\right\}^{T} $$; $$ \mathbf{K}_{\alpha, \beta}^{s} $$ is a small augmented matrix and can be expressed as a Kronecker product^[37], determined by

where $$ \mathbf{A}_{i} $$ is a stochastic Galerkin matrix of size $$ M\times M $$, and $$ \mathbf{A}_{i}^{j, k}= < \xi_i\psi_j\psi_k> $$; the operator $$ \Theta(\mathbf{A}, \mathbf{B}) $$ denotes the Kronecker product of matrices $$ \mathbf{A} $$ and $$ \mathbf{B} $$.

Then, using Gaussian elimination technique, Equation (7) becomes a global extended e-SC system, as given in

where the e-SC matrix $$ \overline{\mathbf{K}} $$ and the e-SC loads $$ \overline{\mathbf{f}} $$ can be computed by

By solving the e-SC system Equation (9), the global interface solution coefficients $$ \mathbf{u}_{\Gamma} $$ are obtained; the local interior solution coefficient $$ \mathbf{u}_{I}^{s} $$ on each subdomain can be obtained independently as:

As given in Equation (9), the key idea of DDM is to transform the original super-large-scale augmented system into a smaller e-SC system so as to reduce the complexity for calculating the stochastic response of structures. However, even if DDM has improved the efficiency of analyzing the stochastic response of structures through the reduction of complexity, since the coefficient matrix is dense and has lost its structure similar to that of the augmented matrix, a suitable preconditioner has been lacking to date. The ultimate result is that the analysis of stochastic response in large-scale structures is not efficient enough.

Approximate sparse expansion of subdomain-level augmented matrix

In this section, we develop an approximate sparse expansion of the subdomain-level augmented matrix $$ \mathbf{K}_{II}^{s} $$ in Equation (10). Through such an expansion, the subdomain-level augmented matrix can be approximately expressed as a block diagonal matrix and a sparse matrix. This is because the random stiffness matrix at the subdomain level can be decomposed into a form of multiplication of a deterministic stiffness matrix and an approximately diagonal matrix. Since the expansion is established, the current dilemma of being unable to establish an effective preconditioner for the e-SC system will be overcome. Additionally, the multiplication steps between the e-SC matrix and an arbitrary vector will also be optimized.

In order to establish the approximate sparse expansion, we first develop an improved mean-based approximation. Considering the original mean-based approximation $$ \tilde{\mathbf{K}}^{(0)}=\Theta\left(\mathbf{I}, \mathbf{K}_{0, I I}^{s}\right) $$, where $$ \Theta $$($$ \mathbf{A} $$, $$ \mathbf{B} $$) denotes the Kronecker product of matrices $$ \mathbf{A} $$ and $$ \mathbf{B} $$, it is easy to find the inverse^[37], because it is a diagonal matrix, given as $$ \left(\tilde{\mathbf{K}}^{(0)}\right)^{-1}=\Theta\left(\mathbf{I}, \left(\mathbf{K}_{0, I I}^{s}\right)^{-1}\right) $$. However, with the increasing variance of the random field, the effect of this diagonal precondition on improving computational efficiency is getting worse. For this reason, if the subdomain-level augmented matrix is expanded based on the mean-based approximation, as $$ \mathbf{K}_{II}^{s}=\Theta(\mathbf{I}, \mathbf{K}_{0, I I}^{s})\mathbf{R}_{II}^{s} $$, it is only necessary to find an approximate value $$ \tilde{\mathbf{R}}_{II}^{s} $$ for the relation matrix $$ \mathbf{R}_{II}^s $$ of $$ \mathbf{K}_{II}^s $$ in order to obtain a suitable approximation for $$ \mathbf{K}_{II}^s $$.

In order to establish the approximate value $$ \tilde{\mathbf{R}}_{II}^{s} $$, we first consider the approximate sparsity of the random stiffness matrix at the subdomain level. Based on our findings about the properties of random stiffness matrices recently, the expansion term of the random stiffness matrix of the s-number subdomain $$ \mathbf{K}_{i}^s $$ ($$ i=0, 1, 2, \cdots, m $$) can be approximated as a product form, obtained as

where $$ \hat{\mathbf{R}}_{i}^{s} $$ is a diagonal matrix, an approximate value of the relation matrix $$ \mathbf{R}_{i}^{s}=\left(\mathbf{K}_{0}^{s}\right)^{-1} \mathbf{K}_{i}^{s} $$ for $$ \mathbf{K}_{i}^{s} $$. The number of terms $$ \mathbf{K}_{i}^{s}, m $$, depends on the order $$ d $$ of PCE of the stochastic process $$ E(x, \theta) $$ and number $$ N $$ of the KL expansion of the underlying Gaussian field $$ \omega(x, \theta) $$ expressed as $$ m=(N+d)!/(N !d !) $$.

The value of $$ \hat{\mathbf{R}}_{i}^{s} $$ is related to the results of KL expansion of random field at the subdomain; that is, the value of this matrix can be calculated through theoretical computation as:

This result further indicates that if the matrix $$ \mathbf{K}_{i}^{s} $$ is rewritten as a block matrix, as formulated in

Then, its sub-matrices can also be expressed in this approximate form. The submatri $$ \mathbf{K}_{i, I I}^{s} $$ of $$ \mathbf{K}_{i}^{s} $$ can be approximately expressed as

It means that any matrix $$ \mathbf{K}_{i, I I}^{j} $$ can be expressed as the product of the mean term and the approximate diagonal matrix, and the diagonal elements of this approximate diagonal matrix can be approximately obtained through theoretical calculation.

Then, we establish the approximate relation matrix $$ \tilde{\mathbf{R}}_{II}^{s} $$, considering the small augmented matrix $$ \mathbf{K}_{II}^{s} $$ according to the Kronecker product form given in Equation (8); if it is left multiplied by a block diagonal matrix $$ \Theta(\mathbf{I}, (\mathbf{K}_{0, I I}^{j})^{-1}) $$, its relation matrix will be established as $$ \mathbf{R}_{I I}^{s}= \sum_{i=0}^{N} \Theta\left(\mathbf{A}_{i}, \mathbf{R}_{i, I I}^{s}\right) $$. In the relational matrix, each submatrix $$ \mathbf{R}_{i, II}^{s} $$ will be an approximate diagonal matrix.

The submatrix of the relation matrix $$ \mathbf{R}_{i, II}^{s} $$ obtained by direct calculation is dense, and a large number of elements are quite close to 0 in it; these elements should be deleted and the relation matrix should be transformed into an approximate sparse structure with only large elements. Therefore, we establish a function $$ \ell(\mathbf{R}) $$ that converts a dense relationship matrix into an approximately sparse relationship matrix. In the approximate sparse function $$ \ell(\mathbf{R}) $$, firstly, the matrix is rewritten as a set, as:

where in this equation, $$ \vartheta_{j k} $$ is the element in the matrix $$ \mathbf{R} $$.

Then, according to the magnitude of these elements, they are sorted, and a certain amount of the big elements are retained. The preserved elements constitute another set $$ \tilde{\Upsilon} $$, which can be represented as a sparse matrix to obtain the approximate sparse form of the relation matrix $$ \tilde{\mathbf{R}} $$. Using the Frobenius norm to ensure the accuracy of results through the verification errors, that is, to let the Frobenius norm of the errors low enough as $$ ||\mathbf{R}-\tilde{\mathbf{R}}|| \leq \varepsilon $$, where the norm can be calculated by $$ ||\mathbf{A}||_{F}=\sqrt{\operatorname{tr}\left(\mathbf{A}^{T} \mathbf{A}\right)} $$ and $$ \operatorname{tr}(\mathbf{A})=\sum_{i=1}^{n}[\mathbf{A}]_{i, i} $$ is the trace of a matrix.

The above steps represent the implementation process of approximate sparsity, and we define the steps as a function, denoted as

Then, the approximate sparse relationship matrix of the sub-domain augmented matrix will be established through the following steps as given in Algorism 1 (in the APPENDIX). Now, all the elements where each node is located are defined as a set $$ \Xi_{v} $$; $$ v $$ denotes the number of one element, and let $$ \Omega(\Xi_{v}) $$ be the total domain of the elements set $$ \Xi_{v} $$.

Through Algorism 1, a series of relation matrices of random stiffness matrices at the subdomain level are obtained as $$ \tilde{\mathbf{R}}_{i, II}^{s} $$, and the approximate sparse relation matrix of the augmented matrix at this level is assembled as $$ \tilde{\mathbf{R}}_{I I}^{s}=\sum_{i=0}^{m-1} \Theta\left(\mathbf{A}_{i}, \tilde{\mathbf{R}}_{i, I I}^{s}\right) $$. According to the approximate sparse relation matrix, the approximate sparse matrix of the augmented matrix at this level can be obtained as

As given in Equation (18), approximate sparse expansion $$ \mathbf{M}^{s} $$ of $$ \mathbf{K}_{II}^{s} $$ is a product of a block diagonal matrix and a sparse matrix with a very small condition number; as an improvement over the mean-based approximation, its inverse matrix is easily established, and it is sufficiently approximate to $$ \mathbf{K}_{II}^{s} $$ under any conditions. Through this approximation, an approximate value of the e-SC matrix can be also derived, and its preconditioner can be also established. At the same time, a new method can be established to multiply the e-SC matrix with any vector, and using the method, it is possible to avoid performing direct operations with relatively dense e-SC matrices.

Approximate sparse preconditioner of e-SC system

Once the approximate sparse expansion of the subdomain-level augmented matrix is established, as determined by Equation (18), the approximate sparse preconditioner of e-SC can be established in this section. The preconditioner is expressed in the form of multiplication of a block diagonal matrix with a sparse matrix called the approximate relation matrix. The preconditioner is quite close to the e-SC matrix, which makes the iteration of PCG converge within a relatively small number of iteration steps. Furthermore, because the inverse matrix of the approximation relation matrix can be easily established, the computational complexity of each iteration step is considerably low. These facts make the solving efficiency of the e-SC system extremely high and ultimately substantially improve the efficiency of structural stochastic response analysis.

Approximate relation matrix for e-SC matrix

In order to efficiently calculate the approximate values of the relation matrix $$ \mathbf{R} $$, we further develop an algorism based on the two-level domain decomposition. By virtue of the two-level domain decomposition, in this algorism, the inverse of the m-SC matrix is translated into a form of Cholesky decomposition, meaning that only several small-scale matrices need to be inverse to obtain the inverse of the m-SC matrix. Meanwhile, the product of the inverse matrix of m-SC and the submatrix of the e-SC matrix is transformed into a series of small-scale matrix multiplications that can be performed efficiently in parallel. Thus, the problem of obtaining the approximate relation matrix in Section 3.2.1 has been solved.

(1) Computation of the inverse of m-SC matrix

The division of two-level subdomains is based on the division of one-level subdomains, and several adjacent first-level subdomains form a second-level subdomain. The order of e-SC matrix solutions is rearranged by two-level division; then, each sub-matrix has an arrow-shaped block structure according to the element distribution law of its coefficient matrix.

On the basis of finite element mesh, the structure is divided into a certain number $$ N_{s} $$ of subdomains; each subdomain contains a certain number of elements shown in Figure 1A, which is called the first-level domain decomposition. Then, according to the first-level domain decomposition, the structure is divided into fewer subdomains; each subdomain contains a certain number of the first-level subdomains, as shown in Figure 1B.

Figure 1. The two-level domain decomposition of a square plate ($$ \bullet $$ （red） denotes the interface nodes, and the interior nodes are omitted in this figure). (A) the first-level domain decomposition; (B) the second-level domain decomposition.

After the structure is divided into two levels, the order of Schur complementary system solutions is adjusted according to the results of the division, so that its coefficient matrix becomes a block matrix with an arrow structure element distribution in each submatrix. In order to explain this structure more clearly, firstly, the element distribution of the mean interface stiffness matrix $$ \overline{\mathbf{K}}_{0} $$ after two-level division and adjustment of the order of solutions is considered. The m-SC matrix comes from the SC process of a temporary deterministic problem, denoted as

(32) $$ \begin{equation} \begin{aligned} \mathbf{K}_{0}\mathbf{v}=\mathbf{f} \end{aligned} \end{equation} $$

According to the deterministic domain decomposition, the SC process of Equation (32) will get the following m-SC system, as formulated in

(33) $$ \begin{equation} \begin{aligned} \overline{\mathbf{K}}_{0}\mathbf{v}_{\Gamma}=\overline{\mathbf{f}} \end{aligned} \end{equation} $$

By the two-level domain decomposition, the m-SC matrix $$ \overline{\mathbf{K}}_{0} $$ can be expressed by the subdomain-level m-SC matrices and Boolean matrices as $$ \overline{\mathbf{K}}_{0}=\sum_{s=1}^{N_{s}}\left(\mathbf{B}^{s}\right)^{T} \overline{\mathbf{K}}_{0}^{s} \mathbf{B}^{s} $$. It is obvious that, the matrix is a series of small-scale matrices mapped and added by Boolean matrices; therefore, there are a large number of zero elements in the matrix, and the distribution of zero elements satisfies

(34) $$ \begin{equation} \begin{aligned} \overline{\mathbf{K}}_{0}^{(i j)}\left\{\begin{array}{ll} \neq 0 & \text { if } \mathbf{v}_{\Gamma}^{d} \text { and } \mathbf{v}_{\Gamma}^{j} \text { in the same subdomain } \\ =0 & \text { if } \mathbf{v}_{\Gamma}^{i} \text { and } \mathbf{v}_{\Gamma}^{j} \text { in different subdomains } \end{array}\right. \end{aligned} \end{equation} $$

According to the distribution law of zero elements given in Equation (34) and the division of the second-level domain decomposition shown in Figure 1, we reorder the solutions of the interface $$ \mathbf{v}_{\Gamma} $$ and block them as

(35) $$ \begin{equation} \begin{aligned} \mathbf{v}_{\Gamma}=\left\{\begin{array}{lllll} \mathbf{v}_{\Gamma}^{1} & \mathbf{v}_{\Gamma}^{2} & \ldots & \mathbf{v}_{\Gamma}^{L_{t}} & \mathbf{v}_{\mathrm{X}} \end{array}\right\} \end{aligned} \end{equation} $$

where $$ \mathbf{v}_{\Gamma}^{t} $$ denotes the second-level interior solution, and $$ \mathbf{v}_{\mathrm{X}} $$ denotes the second-level interface solution.

After rearranging the order of SC matrices according to the new order of interface solution vectors derived in Equation (35), based on the law of element distribution in Equation (34), the matrix will become the following arrow shape expressed as

(36) $$ \begin{equation} \begin{aligned} \overline{\mathbf{K}}_{0}=\left[\begin{array}{ccc|c} \overline{\mathbf{K}}_{0, I I}^{(1)} & \cdots & 0 & \overline{\mathbf{K}}_{0, I \Gamma}^{(1)}\left(\mathbf{B}_{L}^{(1)}\right)^{T} \\ \vdots & \ddots & \vdots & \vdots \\ 0 & \cdots & \overline{\mathbf{K}}_{0, I I}^{L_{t}} & \overline{\mathbf{K}}_{0, I \Gamma}^{L_{t}}\left(\mathbf{B}_{L}^{\left(L_{t}\right)}\right)^{T} \\ \hline \mathbf{B}_{L}^{(1)} \overline{\mathbf{K}}_{0, \Gamma I}^{(1)} & \cdots & \mathbf{B}_{L}^{\left(L_{t}\right)} \overline{\mathbf{K}}_{0, \Gamma I}^{L_{t}} & \sum\limits_{t=1}^{\left(L_{t}\right)} \mathbf{B}_{L}^{(t)} \overline{\mathbf{K}}_{0, \Gamma \Gamma}^{(t)}\left(\mathbf{B}_{L}^{(t)}\right)^{T} \end{array}\right] \\ \end{aligned} \end{equation} $$

where $$ \mathbf{B}_{L}^{(t)} $$ acts as a scatter or gather operator between global and local components of the deterministic second-level interface solution vectors as

(37) $$ \begin{equation} \begin{aligned} \mathbf{v}_{\mathrm{X}}^{t}=\mathbf{B}_{L}^{(t)} \mathbf{v}_{\mathrm{X}}, \quad \mathbf{v}_{\mathrm{X}}=\left(\mathbf{B}_{L}^{(t)}\right)^{T} \mathbf{v}_{\mathrm{X}}^{t} \end{aligned} \end{equation} $$

The form of the block matrix given in Equation (36) is quite easy to find the inverse. Because it can be easily expanded into a simpler form by block Cholesky decomposition, denoted as

(38) $$ \begin{equation} \begin{aligned} \overline{\mathbf{K}}_{0}=\left[\begin{array}{cc} \mathbf{I}_{I} & \mathbf{0} \\ \overline{\mathbf{K}}_{0, \Gamma I} \overline{\mathbf{K}}_{0, I I}^{-1} & \mathbf{I}_{\Gamma} \end{array}\right]\left[\begin{array}{cc} \overline{\mathbf{K}}_{0, I I} & \mathbf{0} \\ \mathbf{0} & \mathbf{S} \end{array}\right]\left[\begin{array}{cc} \mathbf{I}_{I} & \overline{\mathbf{K}}_{0, I I}^{-1} \overline{\mathbf{K}}_{0, I \Gamma} \\ \mathbf{0} & \mathbf{I}_{\Gamma} \end{array}\right] \\ \end{aligned} \end{equation} $$

where $$ \overline{\mathbf{K}}_{0, \Gamma I}=\left[\begin{array}{llll} \mathbf{B}_{L}^{(1)} \overline{\mathbf{K}}_{0, \Gamma I}^{(1)} & \mathbf{B}_{L}^{(2)} \overline{\mathbf{K}}_{0, \Gamma I}^{(2)} & \ldots & \mathbf{B}_{L}^{\left(L_{t}\right)} \overline{\mathbf{K}}_{0, \Gamma I}^{\left(L_{t}\right)} \end{array}\right] $$, $$ \mathbf{I}_{I} $$ and $$ \mathbf{I}_{\Gamma} $$ denote the identity matrix, $$ \overline{\mathbf{K}}_{0, II}=\text{blockdiag}(\overline{\mathbf{K}}_{0, II}^{(1)}, \overline{\mathbf{K}}_{0, II}^{(2)}, \ldots, \overline{\mathbf{K}}_{0, II}^{(L_{t})}) $$ is a block diagonal matrix, $$ \overline{\mathbf{K}}_{0, I\Gamma }=\overline{\mathbf{K}}_{0, \Gamma I}^{T} $$, and $$ \mathbf{S} $$ is the second-level SC matrix which can be obtained by

(39) $$ \begin{equation} \begin{aligned} \mathbf{S}= \sum\nolimits_{t=1}^{L_{t}} \mathbf{B}_{L}^{(t)}\left(\overline{\mathbf{K}}_{0, \Gamma \Gamma}^{t}-\overline{\mathbf{K}}_{0, \Gamma I}^{t}\left(\overline{\mathbf{K}}_{0, II}^{t}\right)^{-1} \overline{\mathbf{K}}_{0, I \Gamma}^{t}\right)\left(\mathbf{B}_{L}^{(t)}\right)^{T} \end{aligned} \end{equation} $$

It can be easily observed that the inverse of this matrix can be expressed as

(40) $$ \begin{equation} \begin{aligned} \overline{\mathbf{K}}_{0}^{-1}=\left[\begin{array}{cc} \mathbf{I}_{I} & -\overline{\mathbf{K}}_{0, I I}^{-1} \overline{\mathbf{K}}_{0, I \Gamma} \\ \mathbf{0} & \mathbf{I}_{\Gamma} \end{array}\right]\left[\begin{array}{cc} \overline{\mathbf{K}}_{0, I I}^{-1} & \mathbf{0} \\ \mathbf{0} & \mathbf{S}^{-1} \end{array}\right]\left[\begin{array}{cc} \mathbf{I}_{I} & \mathbf{0} \\ -\overline{\mathbf{K}}_{0, \Gamma I} \overline{\mathbf{K}}_{0, I I}^{-1} & \mathbf{I}_{\Gamma} \end{array}\right] \end{aligned} \end{equation} $$

where only a series of internal stiffness matrices for second-level sub-domains $$ \overline{\mathbf{K}}_{0, II} $$ and the second-level SC matrix $$ \mathbf{S} $$ needs to be inversed. Usually, the second-level SC matrix will be fairly small-scale, which can be directly inversed, because after two-level domain decomposition, there will be quite few second-level interface nodes. Especially for structures with many small sections, such as frame structures and power transmission tower structures, it is quite easy to find a two-level subdomain division scheme with very few second-level interface nodes.

Once the inverse of the m-SC matrix is obtained as Equation (40), the multiplication of the inverse of the m-SC matrix and a vector (or matrix) can be transformed into a more efficient form. Taking multiplication of the matrix and a vector $$ \mathbf{a}=\overline{\mathbf{K}}_{0}^{-1} \mathbf{b} $$ as an example, the multiplication will be transformed into three steps: (1) the SC process of the right-hand side vector; (2) the solution of the SC system to obtain the interface solution; and (3) substituting the interface solution into the calculation of the internal solution, as determined by

(41) $$ \begin{equation} \begin{aligned} \left\{\begin{array}{l} \mathbf{h}_{\Gamma}=\mathbf{b}_{\Gamma}-\overline{\mathbf{K}}_{0, \Gamma I} \overline{\mathbf{K}}_{0, I I}^{-1} \mathbf{b}_{I} \quad\quad (1)\\ \mathbf{a}_{\Gamma}=\mathbf{S}^{-1} \mathbf{h}_{\Gamma} \quad\quad\quad \quad\quad\quad \quad\; (2)\\ \mathbf{a}_{I}=\overline{\mathbf{K}}_{0, I I}^{-1}\left(\mathbf{b}_{I}-\overline{\mathbf{K}}_{0, I \Gamma} \mathbf{a}_{\Gamma}\right) \quad\; (3) \end{array}\right. \end{aligned} \end{equation} $$

where only multiplication of the small-scale matrix and vector and inverse of the small-scale matrix exist, and clearly, it is more efficient to multiply than computing the inverse matrix and then multiplying it with the vector directly.

(2) Determination of relation matrix

Based on our previous research findings, it has been concluded that the e-SC matrix can be represented in a two-level block form, where each sub-block has the same distribution of elements as the determined SC matrix, obtained as

(42) $$ \begin{equation} \begin{aligned} \overline{\mathbf{K}}=\left[\begin{array}{cccc} \sum\nolimits_{s=1}^{N_{s}} \mathbf{B}^{s} \overline{\mathbf{K}}_{0, 0}^{s}\left(\mathbf{B}^{s}\right)^{T} & \sum\nolimits_{s=1}^{N_{s}} \mathbf{B}^{s} \overline{\mathbf{K}}_{0, 1}^{s}\left(\mathbf{B}^{s}\right)^{T} & \cdots & \sum\nolimits_{s=1}^{N_{s}} \mathbf{B}^{s} \overline{\mathbf{K}}_{0, M-1}^{s}\left(\mathbf{B}^{s}\right)^{T} \\ \sum\nolimits_{s=1}^{N_{s}} \mathbf{B}^{s} \overline{\mathbf{K}}_{1, 0}^{s}\left(\mathbf{B}^{s}\right)^{T} & \sum\nolimits_{s=1}^{N_{s}} \mathbf{B}^{s} \overline{\mathbf{K}}_{1, 1}^{s}\left(\mathbf{B}^{s}\right)^{T} & \cdots & \sum\nolimits_{s=1}^{N_{s}} \mathbf{B}^{s} \overline{\mathbf{K}}_{1, M-1}^{s}\left(\mathbf{B}^{s}\right)^{T} \\ \vdots & \vdots & \ddots & \vdots \\ \sum\nolimits_{s=1}^{N_{s}} \mathbf{B}^{s} \overline{\mathbf{K}}_{M-1, 0}^{s}\left(\mathbf{B}^{s}\right)^{T} & \sum\nolimits_{s=1}^{N_{s}} \mathbf{B}^{s} \overline{\mathbf{K}}_{M-1, 1}^{s}\left(\mathbf{B}^{s}\right)^{T} & \cdots & \sum\nolimits_{s=1}^{N_{s}} \mathbf{B}^{s} \overline{\mathbf{K}}_{M-1, M-1}^{s}\left(\mathbf{B}^{s}\right)^{T} \end{array}\right] \end{aligned} \end{equation} $$

where each submatrix can be transformed into this arrow-shaped representation akin to Equation (36) through the two-level region decomposition, as derived in

(43) $$ \begin{equation} \begin{aligned} \overline{\mathbf{K}}_{j k}=\left[\begin{array}{ccc|c} \overline{\mathbf{K}}_{j k, I I}^{(1)} & \cdots & 0 & \overline{\mathbf{K}}_{j k, I \Gamma}^{(1)}\left(\mathbf{B}_{L}^{(1)}\right)^{T} \\ \vdots & \ddots & \vdots & \vdots \\ 0 & \cdots & \overline{\mathbf{K}}_{j k, I I}^{\left(L_{i}\right)} & \overline{\mathbf{K}}_{j k, I \Gamma}^{\left(L_{i}\right)}\left(\mathbf{B}_{L}^{\left(L_{i}\right)}\right)^{T} \\ \hline \mathbf{B}_{L}^{(1)} \overline{\mathbf{K}}_{j k, I \Gamma}^{(1)} & \cdots & \mathbf{B}_{L}^{(L_{i})} \overline{\mathbf{K}}_{j k, I \Gamma}^{\left(L_{i}\right)} & \sum\limits_{i=1}^{\left(L_{i}\right)} \mathbf{B}_{L}^{(i)} \overline{\mathbf{K}}_{j k, \Gamma \Gamma}^{(i)}\left(\mathbf{B}_{L}^{(i)}\right)^{T} \end{array}\right] \end{aligned} \end{equation} $$

It is noteworthy that substituting matrices of the same number of rows as $$ \overline{\mathbf{K}}_{0} $$ for $$ \mathbf{a} $$ and $$ \mathbf{b} $$ also yields a valid equation. Therefore, Equation (41) can also be applied to the construction of the relationship matrix $$ \mathbf{R}_{j k}=\overline{\mathbf{K}}_{0}^{-1}\overline{\mathbf{K}}_{jk} $$. The obtained relation matrix through this approach will be block-structured as

(44) $$ \begin{equation} \begin{aligned} \mathbf{R}_{j k}=\left[\begin{array}{cccc} \mathbf{R}_{j k}^{(1, 1)} & \mathbf{R}_{j k}^{(1, 2)} & \cdots & \mathbf{R}_{j k}^{\left(1, L_{i}+1\right)} \\ \mathbf{R}_{j k}^{(2, 1)} & \mathbf{R}_{j k}^{(2, 2)} & \cdots & \mathbf{R}_{j k}^{\left(2, L_{i}+1\right)} \\ \vdots & \vdots & \ddots & \vdots \\ \mathbf{R}_{j k}^{\left(L_{i}+1, 1\right)} & \mathbf{R}_{j k}^{\left(L_{i}+1, 2\right)} & \cdots & \mathbf{R}_{j k}^{\left(L_{i}+1, L_{i}+1\right)} \end{array}\right] \end{aligned} \end{equation} $$

Directly approximating sparseness of the relationship matrix would consume a considerable amount of time due to the lowest computational complexity of sorting, which is at least $$ \text{O}(n\text{log}\; n) $$. Therefore, it is possible to take advantage of the fact that the sub-matrices of the relationship matrix can be calculated independently and perform two rounds of sparse approximation. For any column of the relationship matrix, according to the above steps, the interface solution is obtained through steps (1) and (2) in Equation (41) first, and then this interface solution is used to determine the interior solution of the subdomain by step (3) in Equation (41).

Since the inverse of second-level SC matrix $$ \mathbf{S} $$ will be used repeatedly in calculations, it is advisable to compute and save the inverse $$ \mathbf{S}^{-1} $$ at the beginning of the analysis. Firstly, considering any column $$ t $$ in the first $$ L_t $$ columns, through the calculation of steps (1) and (2) in Equation (41), the value of the interface solution can be obtained as

(45) $$ \begin{equation} \begin{aligned} \mathbf{R}_{j k}^{\left(L_{i}+1, t\right)}=\mathbf{S}^{-1}\left(\mathbf{B}_{L}^{(t)} \overline{\mathbf{K}}_{j k, I \Gamma}^{(t)}-\overline{\mathbf{K}}_{0, I \Gamma}^{(t)}\left(\mathbf{B}_{L}^{(t)}\right)^{T}\left(\overline{\mathbf{K}}_{0, I I}^{(t)}\right)^{-1} \overline{\mathbf{K}}_{j k, I I}^{(t)}\right) \end{aligned} \end{equation} $$

where a series of matrix multiplication problems is involved, including the multiplication of Boolean matrix and other matrices. In order to improve the calculation efficiency of this step, we change the original calculation order, Using the properties of Boolean matrix, $$ \left(\mathbf{B}_{L}^{(t)}\right)^{T}\left(\overline{\mathbf{K}}_{0, I I}^{(t)}\right)^{-1} $$ is calculated first to reduce the scale of the matrix, then other steps are calculated, so that each operation is a multiplication between small-scale matrices, then the calculation time is reduced. By substituting the calculation result of Equation (45) into step (3) of (41), the rest solutions can be obtained as

(46) $$ \begin{equation} \begin{aligned} \mathbf{R}_{j k}^{(o, t)}=\left(\overline{\mathbf{K}}_{0, I I}^{(t)}\right)^{-1}\left(\delta_{o}^{t} \overline{\mathbf{K}}_{j k, I I}^{(t)}-\overline{\mathbf{K}}_{0, I \Gamma}^{\left(L_{i}\right)}\left(\mathbf{B}_{L}^{\left(L_{i}\right)}\right)^{T} \mathbf{R}_{j k}^{\left(L_{i}+1, t\right)}\right) \end{aligned} \end{equation} $$

Similar to formula a, there is also a Boolean matrix in the calculation of Equation (46), so although what needs to be calculated is the multiplication of $$ \overline{\mathbf{K}}_{0, I \Gamma}^{(L_{i})}{(\mathbf{B}_{L}^{(L_{i})})}^{T} $$ and $$ \mathbf{R}_{j k}^{\left(L_{i}+1, t\right)} $$, in practical calculation, in order to reduce the calculation cost, we choose to calculate $$ {\left(\mathbf{B}_{L}^{\left(L_{i}\right)}\right)}^{T}\mathbf{R}_{j k}^{\left(L_{i}+1, t\right)} $$ first, and then multiply $$ \overline{\mathbf{K}}_{0, I \Gamma}^{\left(L_{i}\right)} $$ left to get the final result. For every $$ o=1, 2, \cdots, L_{t} $$, $$ \mathbf{R}_{j k}^{(o, t)} $$ can be calculated independently, so in order to reduce the storage cost as much as possible, when each $$ \mathbf{R}_{j k}^{(o, t)} $$ is calculated, it is approximately sparse before calculating the next item $$ \mathbf{R}_{jk}^{(o+1, t)} $$.

For the last column of the relationship matrix in Equation (44), the interface solution $$ \mathbf{R}_{j k}^{\left(L_{i}+1, L_{i}+1\right)} $$ can be obtained by adopting the same calculation method as

(47) $$ \begin{equation} \begin{aligned} \mathbf{R}_{j k}^{\left(L_{i}+1, L_{i}+1\right)}=\mathbf{S}^{-1}\left( \sum\limits_{s=1}^{(L_{i})}\left(\mathbf{B}_{L}^{(s)} \overline{\mathbf{K}}_{j k, \Gamma \Gamma}^{(s)}\left(\mathbf{B}_{L}^{(s)}\right)^{T}-\mathbf{K}_{0, \Gamma I}^{s}\left(\mathbf{K}_{0, I I}^{s}\right)^{-1} \overline{\mathbf{K}}_{j k, I \Gamma}^{(s)}\left(\mathbf{B}_{L}^{(s)}\right)^{T}\right)\right) \end{aligned} \end{equation} $$

Bring the calculation results into the sub-formula (3) of Equation (41), the rest of the solution can be obtained as

(48) $$ \begin{equation} \begin{aligned} \mathbf{R}_{j k}^{\left(s, L_{i}+1\right)}=\left(\overline{\mathbf{K}}_{0, I I}^{(s)}\right)^{-1}\left(\overline{\mathbf{K}}_{j k, I \Gamma}^{(s)}\left(\mathbf{B}_{L}^{(s)}\right)^{T}-\overline{\mathbf{K}}_{0, I \Gamma}^{\left(L_{1}\right)}\left(\mathbf{B}_{L}^{\left(L_{i}\right)}\right)^{T} \mathbf{R}_{j k}^{\left(L_{i}+1, L_{i}+1\right)}\right) \end{aligned} \end{equation} $$

It can be observed that in Equation (47), the same contents $$ \mathbf{K}_{0, \Gamma I}^{s}\left(\mathbf{K}_{0, I I}^{s}\right)^{-1} $$ will be calculated every time the submatrix of the relationship matrix is established, so these repeated calculations can be calculated in advance to avoid unnecessary repetition.

During the calculations, according to the specified order of computation, Equations (48) and (46) can be written in an approximate sparse form first as $$ \tilde{\mathbf{R}}_{j k}^{\left(s, L_{i}+1\right)}=\ell\left(\mathbf{R}_{j k}^{\left(s, L_{i}+1\right)}\right) \text { and } \tilde{\mathbf{R}}_{j k}^{(o, t)}=\ell\left(\mathbf{R}_{j k}^{(o, t)}\right) $$; then, Equations (47) and (45) are expressed as $$ \tilde{\mathbf{R}}_{j k}^{\left(L_{i}+1, L_{i}+1\right)}=\ell\left(\mathbf{R}_{j k}^{\left(L_{i}+1, L_{i}+1\right)}\right) \text { and } \tilde{\mathbf{R}}_{j k}^{\left(L_{i}+1, t\right)}=\ell\left(\mathbf{R}_{j k}^{\left(L_{i}+1, t\right)}\right) $$. According to the above calculation results, the relation matrix can be roughly approximated as $$ \tilde{\mathbf{R}}_{j k}^{\prime} $$, as formulated in

(49) $$ \begin{equation} \begin{aligned} \tilde{\mathbf{R}}_{j k}^{\prime}=\left[\begin{array}{cccc} \tilde{\mathbf{R}}_{j k}^{(1, 1)} & \tilde{\mathbf{R}}_{j k}^{(1, 2)} & \cdots & \tilde{\mathbf{R}}_{j k}^{\left(1, L_{i}+1\right)} \\ \tilde{\mathbf{R}}_{j k}^{(2, 1)} & \tilde{\mathbf{R}}_{j k}^{(2, 2)} & \cdots & \tilde{\mathbf{R}}_{j k}^{\left(2, L_{i}+1\right)} \\ \vdots & \vdots & \ddots & \vdots \\ \tilde{\mathbf{R}}_{j k}^{\left(L_{i}+1, 1\right)} & \tilde{\mathbf{R}}_{j k}^{\left(L_{i}+1, 2\right)} & \cdots & \tilde{\mathbf{R}}_{j k}^{\left(L_{i}+1, L_{i}+1\right)} \end{array}\right] \end{aligned} \end{equation} $$

Further approximate sparsity processing is performed on the matrix $$ \tilde{\mathbf{R}}_{jk}^{\prime} $$ to obtain the final approximate sparse relation matrix $$ \tilde{\mathbf{R}}_{j k}=\ell\left(\mathbf{R}_{j k}^{\prime}\right) $$. Assembling these submatrices results, the approximate sparse relation matrix $$ \tilde{\mathbf{R}} $$ for e-SC.

In the calculation process, all the computations involved are operations between small matrices and the approximate sparse processing of small matrices. Since the relation matrix of each submatrix is established by calculating the same amount of content separately, in the actual calculation, this part will be calculated in parallel, and the algorithm for finally calculating the relationship matrix of an e-SC matrix is expressed as Algorithm 4 (in the APPENDIX), where we first choose the row and column numbers of the eigen relation matrices by step 5; then, for every eigen relation matrix, which is split into a block matrix akin to Equation (44), the submatrices in the bottom row are first obtained by steps 6 and 8; using the result, the other row can be calculated by steps 9 and 11; steps 12 16 involve creating the approximate value of these submatrices, then use these submatrices to assembly the approximate eigen relation matrix by steps 17 and 18, finally create the approximate relation matrix by step 20.

Once the approximate relation matrix $$ \tilde{\mathbf{R}} $$ of the e-SC matrix is established by Algorism 4, it is possible to obtain the preconditioner for solving the e-SC system, expressed as

(50) $$ \begin{equation} \begin{aligned} \mathbf{M}=\Theta\left(\mathbf{I}, \overline{\mathbf{K}}_{0}\right) \tilde{\mathbf{R}} \end{aligned} \end{equation} $$

which denotes the proposed approximate sparse preconditioner that is quite close to the e-SC matrix; therefore, using this matrix as the preconditioner for PCG will lead to convergence within a relatively small number of iterations. On the other hand, the matrix is a product of a diagonal matrix and a sparse matrix with very little condition number, which makes the solution of the preconditioning equation relatively easy. Ultimately, adopting this preconditioner can make PCG more efficient in solving the e-SC system, which greatly improves the solution efficiency of the random response of the structure.

Approximate sparse approach for DDM-based SSFE

According to the derivation of Sections 3.1 and 3.2, we develop an approximate sparse approach for DDM-based SSFE; the proposed method targets the step of solving e-SC systems with PCG. An efficient algorism for the multiplication of an e-SC matrix and an arbitrary vector, and an approximate sparse preconditioner are established to improve the computing efficiency of PCG, thus making the structural random response more efficient. The following steps show the new computational framework based on the approximate sparse approach for DDM-based SSFE in this work:

$$ \mathbf{step\; 1} $$ Describe the stochastic process $$ E(x, \theta) $$ through KL expansion (or PCE, depending on the process is Gaussian or non-Gaussian). Calculate the matrix $$ \mathbf{A_i} $$ according to $$ c_{ijk} $$.

$$ \mathbf{step\; 2} $$ Mesh the structure and divide two-level subdomains, establish the first-level Boolean matrix $$ \mathbf{B}^{s} $$ and second-level Boolean matrix $$ \mathbf{B}_{L}^{s} $$.

$$ \mathbf{step\; 3} $$ For $$ s=1, 2, 3, \cdots, N_s $$ and $$ i=0, 1, 2, \cdots, N $$, establish the submatrix of subdomain-level stochastic stiffness matrix $$ \mathbf{K}_{i}^{s} $$, including $$ \mathbf{K}_{i, II}^{s}, \mathbf{K}_{i, I \Gamma}^{s}, \mathbf{K}_{i, \Gamma I}^{s} $$, $$ \mathbf{K}_{i, \Gamma \Gamma}^{s} $$.

$$ \mathbf{step\; 4} $$ Calculate approximate relation matrix $$ \tilde{\mathbf{R}}_{i, II}^{s} $$ of $$ \mathbf{K}_{i, II}^{s} $$ according to Algorism 1

$$ \mathbf{step\; 5} $$ Calculate Second-level m-SC matrix $$ \mathbf{S} $$ and inverse matrix $$ \mathbf{S}^{-1} $$ by Equation (39).

$$ \mathbf{step\; 6} $$ Establish the eigen set $$ \mathfrak{J} $$ by Equation (27).

$$ \mathbf{step\; 7} $$ Calculate the approximate relation matrix $$ \tilde{\mathbf{R}} $$ for the e-SC matrix $$ \overline{\mathbf{K}} $$ by Algorism 4.

$$ \mathbf{step\; 8} $$ Parallelly calculate the interface solution $$ \mathbf{u}_\Gamma $$ according to the proposed solution for the e-SC system Algorism 5 (in the APPENDIX).

$$ \mathbf{step\; 9} $$ Obtain the local interior solution $$ \mathbf{u}_I^s $$ by Equation (11) parallelly.

In the established algorithm, PCG in step 8 is used to solve the e-SC system; different from the traditional solution, we propose a new calculation scheme based on the approximate sparse stiffness matrix for the two most critical steps in PCG: the e-SC matrix and temporary vector multiplication and the solution of preconditioned equations, to solve the problems that the existing solution methods must save large dense matrices and their calculation efficiency are low.

In step 4, according to the derivation in Section 3.1, the approximate relation matrix of expansion term of the random stiffness matrix in subdomain $$ \tilde{\mathbf{R}}_{i, II}^{s} $$ is established. Then, the multiplication of the e-SC matrix and temporary vector (steps 2 and 7 in Algorism 5) can adopt the steps of Algorism 3 for parallel calculation. In Algorism 3, the direct multiplication of e-SC and vectors at the subdomain level will no longer be used, but an algorithm based on the calculation formulated in Equation (10) will be established instead. The subdomain augmented matrix and vector multiplication in Algorism 3 is further transformed into the multiplication given in Algorism 1; then, the preconditioner for solving any subdomain-level augmented system $$ \mathbf{K}_{II}^{s}\mathbf{a}=\mathbf{b} $$ by PCG is built as $$ \mathbf{M}^{s}=\Theta\left(\mathbf{I}, \mathbf{K}_{0, I I}^{s}\right) \tilde{\mathbf{R}}_{I I}^{s} $$ through the approximate relation matrix $$ \tilde{\mathbf{R}}_{i, II}^{s} $$ of $$ \mathbf{K}_{i, II}^{s} $$; this preconditioned matrix is quite close to the subdomain augmented matrix $$ \mathbf{K}_{II}^{s} $$, and because it is composed of a block diagonal matrix and a sparse matrix with low condition number, the preconditioned equation $$ \mathbf{M}^{s}\mathbf{a}=\mathbf{b} $$ is quite easy to solve. By this means, compared with the traditional methods, the proposed new method of multiplying e-SC matrices and temporary vectors, $$ \mathbf{b}=\overline{\mathbf{K}}\mathbf{a} $$, does not need to save large dense matrices and has quite high computational efficiency.

On the basis of Section 3.1, after further derivation, the proximate relationship matrix of e-SC matrix $$ \tilde{\mathbf{R}} $$ is established, and steps 5-7 show the specific calculation process. By creating the eigen set $$ \mathfrak{J} $$, the number of needed multiplication steps between the inverse of the m-SC matrix and submatrix of an e-SC matrix will be reduced, then by the two-level domain decomposition given in Algorism 4, the calculation process of every submatrix of relation matrix is transformed into a series of operations between small matrices, which greatly improves the calculation efficiency while avoiding the direct participation of large dense matrices in the operation. After the calculation, the preconditioner for solving the e-SC is obtained as $$ \mathbf{M}=\Theta(\mathbf{I}, \overline{\mathbf{K}}_{0})\tilde{\mathbf{R}} $$; this matrix ensures the similarity with the original e-SC matrix $$ \overline{\mathbf{K}} $$, and the high calculational efficiency of solving the preconditioned equation on steps 3 and 12 of Algorism 5, because it is composed of block diagonal matrix and a sparse matrix with low condition number.

To sum up, by rewriting the step of the e-SC matrix and the vector multiplication based on the approximate sparse random stiffness matrix, we avoid the large dense matrix directly participating in the operation and improve the calculation efficiency. At the same time, based on the approximate sparse random stiffness matrix, the relationship matrix of the e-SC matrix is further established to create a preconditioning matrix. Using this matrix as the preconditioner for solving the e-SC matrix by PCG can greatly improve the calculation efficiency by controlling the iteration steps. Finally, an efficient framework for structural random response analysis is established.

Stochastic response of a thin square plate

As a test case, we have considered a stochastic response analysis of a thin square plate. The dimensions of this plate are $$ L_x=L_y=1 $$; thickness is 0.1, and Poisson's ratio is 0.3. One of its edges is a fixed end, receiving uniformly distributed load with a size of $$ p=1 $$, as shown in Figure 2A. The elastic modulus of the thin square plate $$ E $$ is modeled as a Gaussian field $$ E(x, \theta) $$ with a mean value of $$ \rm{\mu} $$ = 1.0; the exponential covariance function of the field is expressed as

Figure 2. Problem setting(A) and typical finite element mesh with 49 subdomains (B).

where standard deviation is $$ \sigma=0.2 $$; the correlation lengths of $$ x $$ and $$ y $$ directions are $$ b_x=b_y=1 $$.

The physical domain with unit-square geometry is discretized using an unstructured finite element mesh. After partitioning the mesh, the number of nodes is 5, 041, and the count of Quadrilateral elements totals 4, 900; the number of subdomains involves 16, 25, 36, 49, and 64; the mesh and the division of subdomains are shown in Figure 2B illustrating the situation of 49 subdomains.

The random field $$ E(x, \theta) $$ is expressed by the KL expansion; in this example, the terms of expansion are $$ N $$ = 2, as $$ E(x, \theta)=E_{0}+\sum_{i=1}^{N} E_{i} \xi_{i} $$; the result is shown in the following Figure 3.

Figure 3. KL expansion result of the random field $$ E(x, \theta) $$: Coefficients before the standard Gaussian random variables $$ \xi_{i} $$. (A) $$ E_{1} $$; (B) $$ E_{2} $$.

The response is expressed using sixth order PCE ($$ M $$ = 28) leading to a linear system of order 282296. In this work, the traditional DDM, the DDM based on approximate sparse expansion proposed (the PCGM is considered to have converged when the error satisfies: $$ ||\overline{\mathbf{K}} \mathbf{u}-\overline{\mathbf{f}}\| /\|\overline{\mathbf{f}}\| \leq 10^{-6} $$) and 50, 000 times of Monte Carlo method are used to calculate the random response of structures. Firstly, the result of the stochastic response calculated through the method proposed is shown in Figure 4.

Figure 4. The solution process: Coefficients before the first 4 PC basis. (A) $$ u_{0} $$; (B) $$ u_{1} $$; (C) $$ u_{0} $$; (D) $$ u_{1} $$.

In order to verify the accuracy of the calculated results more clearly, we consider the node No. 2485 (point A in Figure 4A), which has a large average response, and draw the probability density function (PDF) of its displacement in the y direction according to the calculation results of different methods shown as Figure 5.

Figure 5. Stochastic response of y-direction displacement of node No. 2485: (A) the PDF and (B) tail of the PDF.

From the PDF shown in the figure, it can be seen that the random response results in the y-direction of node number 2485 obtained by the proposed method are almost consistent with the traditional SSFE method. Moreover, due to the maximum mean displacement in the y-direction of this node, it can be considered that the random response results of other nodes obtained by the proposed method also meet the accuracy requirements. Therefore, it indicates that the proposed approach has the same accuracy as the traditional SSFE method and can meet the accuracy requirements when the response PCE order is sufficient.

Next, the distribution law of the coefficient matrix elements of the e-SC system is verified. Firstly, the element distribution of the augmented matrix at the subdomain level is verified; the relation matrix $$ \mathbf{R}_{i, I I}^{s}=\left(\mathbf{K}_{0, I I}^{s}\right)^{-1} \mathbf{K}_{i, I I}^{s} $$ is obtained by left multiplying the inverse matrix of the average stiffness matrix of the subdomain level; then, Figure 6 shows the distribution of elements of $$ \mathbf{R}_{1, I I}^{1} $$ and $$ \mathbf{R}_{2, I I}^{1} $$; in the figures, the $$ x $$-axis and $$ y $$-axis represent the row number and column number of the matrix elements, while the z-axis represents the absolute value of the matrix element.

Figure 6. Element distribution of relation matrix of random partial stiffness matrix at the subdomain level: (A) $$ \mathbf{R}_{1, I I}^{1} $$ and (B) $$ \mathbf{R}_{2, I I}^{1} $$.

Figure 6 presents that for the relation matrices of the subdomain-level stochastic stiffness matrices, the non-diagonal elements are almost zero, while the values of the diagonal elements are relatively large compared to the non-diagonal elements, which confirms our findings.

Based on the above verification results, we further verify whether the properties of the relationship matrix of the e-SC matrix are consistent with the conclusion of proof in Section 3.2. Firstly, it is necessary to verify the distribution of the elements in the relation matrix $$ \mathbf{R}=\left(\Theta\left(\mathbf{I}, \overline{\mathbf{K}}_{0}\right)\right)^{-1} \overline{\mathbf{K}} $$. As shown in Figure 7, we choose some of the submatrices of the relation matrix, and the distribution of their elements was charted. Similar to Figure 6, in Figure 7, the $$ x $$-axis and $$ y $$-axis represent the row number and column number of the matrix elements, while the $$ z $$-axis represents the absolute value of the matrix element.

Figure 7. Element distribution of some relation matrix’s submatrices: (A) the first row, second column $$ \mathbf{R}_{1,2} $$ and (B) the first row, third column $$ \mathbf{R}_{1,3} $$.

It is evident that the sub-matrices of the relation matrix of the e-SC matrix also exhibit similar patterns of element distribution; that is, the non-diagonal elements are almost zero, while the values of the diagonal elements are relatively large compared to the non-diagonal elements.

Therefore, it can be said that, as demonstrated in this article, the distribution pattern of the elements in the relation matrix has been confirmed. Thus, the approximate sparse preconditioner established in Section 3.2 is approximately equal to the e-SC matrix, thus greatly improving the efficiency of solving the e-SC system.

On this basis, we study the magnitude of improvement in computational efficiency provided by the method proposed in our research. We compared the computational time of the traditional DDM-based SSFE and the proposed method, including the time of the SC process, the direct CG solver and the PCG solver with the proposed preconditioner of the e-SC system, as shown in Table 1.

As demonstrated in Table 1, it can be obviously seen that compared with the traditional DDM-based SSFE, which uses the direct CG solver to solve the e-SC system, PCG with the proposed approximate sparse preconditioner is highly efficient, and as a result, the efficiency of stochastic response analysis of structures is greatly improved. Especially when there are a large number of subdomains and the solving time of the e-SC system becomes a higher proportion, this improvement in computational efficiency will become more evident from the fact shown in Table 1.

Stochastic response of a frame structure

The second example features the stochastic response analysis of a frame structure, where we pay greater attention to its computational efficiency under different circumstances, including varying standard deviations and correlation lengths, than in Section 4.1.

As shown in the following Figure 8, the structure we consider in this example is a five-layer frame structure with a layer height of 3 m, the distance between each column is 2 m, the cross-sectional widths of the beams and columns are 0.3 m, thickness is 0.1, and Poisson's ratio is 0.25. The structure is subjected to the self-weight $$ G $$ with planar density $$ \rho=7.85 \times 10^{3}\; \mathrm{Kg} / \mathrm{m}^{2} $$ and a lateral uniformly distributed load P with P = 700 N/m. The elastic modulus of the material in this frame structure is considered as a gaussian random field with a mean of $$ \bar{E}=2.06 \times 10^{5}\; \mathrm{Mpa} $$ and a standard deviation of $$ \sigma=0.2 \bar{E} $$, the correlation length of the field is $$ b_{x}=b_{y}=15.3\; \mathrm{\; m} $$, and the exponential covariance function is the same as given in Equation (51). The finite element mesh of the structure consists of 9, 552 elements and 11, 855 nodes, and 65 subdomains are generated after subdivision.

Figure 8. The geometry and subjected load of the frame structure.

The KL expansion is used to express the random field; in this example, the number of terms $$ N $$ used in the truncated KL is 2, the maximum polynomial order of the PCE is 6, the total terms of the PC basis of the solution are $$ M $$ = 28, and the augmented system to be solved will have a coefficient matrix of $$ 663, 880 \times 663, 880 $$. Using our approximate sparse expansion-based domain decomposition solver, we get the solution of the augmented system, and Figure 9 shows the coefficients before the first PC basis of the solution process of the x-direction and y-direction.

Figure 9. The coefficient before the first PC basis of solution process: (A) Component in the x-direction (B) Component in the y-direction.

Then, we use PDF of the y-direction displacement at point A and the x-direction displacement at point B to describe the calculated results (points A and B are shown in Figure 9; A is the point with the largest average displacement in the y-direction shown in Figure 9A, while B is the point with the largest average displacement in the x-direction shown in Figure 9B) and compare them with the results of 50000 Monte-Carlo simulations, as shown in Figure 10.

Figure 10. The PDF of (A) x-direction displacement of point A and (B) y-direction displacement of point B.

It is obvious that the solution process calculated by the proposed approach is similar to the MC method in Figure 10, indicating that the proposed one can meet the accuracy requirements when the PCE order of the response is sufficient. Then, we are concerned about the numerical scalability of the improvement of computational efficiency, that is, the computational efficiency improvement of the proposed method under different conditions (the PCGM is considered to have converged when the error satisfies: $$ (||\overline{\mathbf{K}} \mathbf{u}-\overline{\mathbf{f}}\| /\|\overline{\mathbf{f}}\| \leq 10^{-8}) $$). So, we first consider the different standard deviation $$ \sigma $$ of the elastic modulus $$ E $$, taking values of $$ 0.1 \bar{E}, 0.15 \bar{E}, 0.2 \bar{E} $$ and $$ 0.25 \bar{E} $$, respectively. Then, we consider the different correlation lengths of $$ b $$ ($$ b=b_x=b_y $$), taking values of $$ 0.2b , 0.4b , 0.6b , 0.8b $$ and $$ b $$, respectively. The computation time statistics are documented in the following Table 2 and Table 3:

From the computational time presented in Table 2 and Table 3, it is evident that the proposed approximate sparse preconditioner exhibits high efficiency; moreover, the degree of improvement does not show significant variations as the standard deviation and the correlation length change, thus indicating its numerical scalability.

Authors’ contributions

Conceptualization, methodology, analysis and writing: Luo B

Software coding, numerical analysis and validation: Cao W

Methodology, supervision and review: Dai H, Zhou Z

Availability of data and materials

Not applicable.

Financial support and sponsorship

This research was supported by a Grant from the National Natural Science Foundation of China (Project 12272109). This support is gratefully acknowledged.

Conflicts of interest

All authors declared that there are no conflicts of interest

Ethical approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Copyright

© The Author(s) 2024.