Efficient strong Unified Formulation for stress analysis of non-prismatic beam structures

The Unified Formulation (UF) has gained attention as a powerful tool for efficient design of structural components. Due to the inherent flexibility of its kinematics representation, arbitrary shape functions can be selected in different dimensions to achieve a high‐fidelity characterisation of structural response under load. Despite this merit, the classical isoparametric description of UF limits the application to prismatic structures. The weak‐form anisoparametric approach adopted to overcome this limitation in a recent work by Patni et al. proves to be versatile yet computationally challenging owing to the expensive computation of its UF stiffness matrix by means of full volume integrals. We propose a strong‐form anisoparametric UF (SUF) based on the Serendipity Lagrange Expansion (SLE) cross‐sectional finite element and differential quadrature beam element. The main objective of the SUF is to achieve an efficient computation of the UF stiffness matrix by restricting Gauss operations to the variable cross‐sections of non‐prismatic structures in a discrete sense, thus eliminating the need for full volume integrals. When assessed against weak‐form based UF, ABAQUS FE and analytical solutions, the static analysis of non‐prismatic beam‐like structures under different loads by the SUF is shown to be accurate, numerically stable, and computationally more efficient than state‐of‐the‐art methods.


Introduction
The applications of variable cross-section, non-prismatic beam-like structures is common in engineering fields due to weight and efficiency requirements that necessitate the use of tailored geometry for functional purposes, for example, in mechanical and aeronautical industries. For such variable geometry structures, accurate stress analyses for optimisation involve expensive computational procedures that can characterise global and localised deformations and their far field effects, especially regarding structural systems with non-homogenous material or load distributions. According to studies by Balduzzi et al. [2,3], cross-sectional variation generates a non-trivial stress distribution in non-prismatic structures partly due to the coupled axial and shear-bending effects.
A crucial aspect of structural optimisation involves development of mathematical models that are simultaneously accurate and computationally efficient. Obtaining 3D finite element (FE) based solution that can accurately capture mechanical responses like in-plane bending, stretching or through-thickness transverse stresses poses computational challenges, thus necessitating the adoption of simplified models. Over the last few decades, researchers have focused on developing tools for modelling non-prismatic beam-like structures. In this regard, several finite element models have been proposed in the literature for analysis of tapered structures. For example, Yang and Yau [4] and Kim and Kim [5] investigated, respectively, the instability and lateraltorsional buckling and vibration of doubly symmetric tapered I-beam while Zang and Tong [6] performed linear stability analysis of webtapered beams with doubly symmetric I-section using a potential energy based method. Other examples of FE based works on nonprismatic structures can be found in [7][8][9][10][11][12].
In the context of analytical approach, Yuan et al. [13] obtained the lateral-buckling load of steel web-tapered T-section cantilevers under a uniformly distributed load or concentrated load at the free edge. Benyamina [14] proposed a nonlinear formula for lateral buckling analysis of simply supported tapered beams with doubly symmetric cross-sections. Hodges et al. developed a variational asymptotic method (VAM) [15] for stress analysis of linearly tapered beams. However, the VAM-based beam theory is limited to problems for which elasticity solutions exist. Furthermore, Balduzzi et al. [2] derived a set of ordinary differential equations (ODEs) which permit the generation of analytical and numerical solutions to study the effect of nonprismatic geometry on beam's behaviour. It was established that the shear stress distribution of the cross-section is non-trivial as it is significantly influenced by generalised stresses as well as the beam geometry. Balduzzi et al. [3] further employed a simplified constitutive relation and stress representation to present a Timoshenko-like beam model for multi-layer non-prismatic structures in which non-trivial shear stress distribution is resolved by means of Jourawsky theory [16]. Due to the extremely simple kinematics, Balduzzi et al.'s model [3] is, unfortunately, not suited for capturing boundary effects or phenomena around discontinuities. In a recent development, Vilar et al. [17] improved the recovery of transverse normal stresses of linearly tapered and general asymmetric beams by accounting for the contributions of internal force derivatives.
For practical purposes, the choice of numerical tool for analysing non-prismatic structures is considerably dictated by the requirements for generality and versatility of the method on the one hand and the constraints to satisfy efficient design on the other hand. In this regard, the theory of Unified Formulation (UF) proposed in [18,19] for 3D analysis of engineering structures shows promise. In the context of beam structures, the UF displacement fields are approximated as a product of predetermined cross-sectional deformation often defined by hierarchical expansion functions over the cross-section, and unknown 1D displacement fields. So far, many researchers have exploited the merits of the UF theory for different structural analyses ranging from 3D stress predictions [20,21], buckling and postbuckling [22][23][24][25], free vibration [26][27][28] and failure [29,30] in a computationally efficient manner without compromising numerical accuracy. Remarkably, these UF-based models have been implemented either by considering a weak-form solution in a finite element framework or adopting a strong-form solution based on a high-order collocation strategy such as the differential quadrature method (DQM). It is worth noting that the DQM proposed by Bellman and Casti [31] has been widely applied for many mechanics problem [20,21,24,25,[32][33][34][35][36] while the recently proposed inverse counterpart known as inverse differential quadrature method (iDQM) proposed by Ojo et al. [37] and Trinh et al. [38] shows promise, due to its characteristics of spectral accuracy and fast convergence.
Notwithstanding the merits of UF for structural analysis, the classical isoparametric representation is generally limited to applications involving prismatic structures, thereby largely excluding complex variable geometric configurations characteristic of real-life structures. According to [39,40], non-prismatic effect in tapered structures can be captured by employing a strategy that takes the tapered plane as the cross-section and treat the beam width as the longitudinal axis. However, this approach is limited by the high aspect ratio of the cross-sectional elements which tend to compromise the fundamental kinematics of the beam leading to numerical inaccuracies. To overcome this limitation, Patni et al. [1] employed a 3D shape function to describe the structural geometry separately from the shape functions that capture the beam kinematics in a FE framework. In essence, two distinct mesh configurations are required to capture the geometry and structural kinematics (see Fig. 1). This approach is general application-wise, as it allows derivation of elements of arbitrary geometry. However, the consequence of such anisoparametric representation is the substitution of the split volume integrals required to compute the stiffness matrix with full volume integrals. From the computational point of view, this approach can increase the number of loop operations for computing the stiffness matrix by many orders depending on the total number of Gauss points considered in the volume to achieve satisfactory accuracy. Consequently, this approach poses some challenges for efficient design of structural components.
In line with the investigation in [20], where comparisons between strong, high-order, and low-order weak UFs for prismatic structures were examined, it was noted that the high-order Weak Unified Formulation (WUF) (henceforth referred to as WUF unless stated otherwise) which employs a single high-order beam element in a FE framework offers an alternative strategy to improve the efficiency of low-order WUF-FE scheme that uses several beam elements without compromising the numerical accuracy or stability of the system. On the other hand, it was established in [20] that the strong UF (SUF) based on DQM discretisation guarantees spectral convergence and numerical stability of the system solutions, similarly to WUF, in addition to employing less Gauss operations to compute the global stiffness matrix. In this work, we develop an efficient SUF based on the anisoparametric configuration proposed by Patni et al. [1] to investigate static deformation of homogeneous and non-homogeneous nonprismatic structures. The main idea of the SUF lies in the preservation of split integrals for prismatic structures in a manner that allows 1D strong-form equilibrium equations to be derived based on variational principles. As such, this strategy restricts the 2D Gauss operations (as per surface integrals) to the variable cross-section in a discrete sense. Considering the computational merits of the WUF over loworder WUF as noted in [20], an anisoparametric WUF-FE is further proposed as a control solution with which the numerical performance of the SUF is assessed, since both formulations have comparable convergence. Except for using a single high-order beam element against several low-order beam elements, the WUF proposed here is similar in all regards to the low-order WUF proposed by Patni et al. [1]. The rest of the paper is presented as follows: Fig. 1. Geometric and kinematic representations of non-prismatic beam-like structure [1]. Section 2 explains the preliminaries of anisoparametric Unified Formulation for beams and also describes the theoretical framework of high-order WUF. In Section 3, strong-form equilibrium equations for non-prismatic structures is derived together with the discretisation in a DQM framework. Furthermore, the details surrounding computation of the stiffness matrices of SUF and WUF are explained. Section 4 is dedicated to numerical examples to demonstrate the performance of the proposed SUF while Section 5 highlights the computational aspects SUF vis-à-vis WUF subject to varying non-prismatic effects. Finally, Section 6 contains the concluding remarks of the paper.

Anisoparametric Unified Formulation for beams
Consider an arbitrary non-prismatic beam-like structure of length L with a varying cross-sectional area A as shown in Fig. 2. The y-axis is defined along the beam axis while the z-and x-axes define the crosssection.
According to the theory of UF, the three-dimensional displacement fields of the beam can be expressed as where F s are the Serendipity Lagrange Expansion (SLE) shape functions which capture the cross-sectional deformation of the beam structure while u s is the 1D vector of displacement fields dependent on y. The termsu x , u y and u z are the 1D displacement components in the x-, yand z-directions respectively. The M variable represents the number of expansion terms and the repeated subscript s indicates summation.
In line with the strategy adopted in [41], the cross-section of the beam is discretised with four-node rectangular SLE elements and enrichment of the displacement field within the subdomain is accomplished by increasing the order of the SLE elements arbitrarily. Readers are referred to [41] for details of the SLE formulation. Generalised stress and strain components follow from the UF, assuming small deformations where E T ¼ E xx E yy E zz E yz E xz E xy ½ and D is the linear differential operator expressed as D T ¼ @ @x 0 0 0 @ @z @ @y 0 @ @y 0 @ @z 0 @ @x 0 0 @ @z @ @y @ @x The stress is recovered by using the constitutive relation where S T ¼ S xx S yy S zz S yz S xz S xy ½ and the transformed constitutive matrix C for a general material composition is expressed as C αβ are the components of transformed constitutive matrix whose derivation is detailed in [42].

Strain energy
Based on the principle of virtual work (PVD), neglecting inertial effects, the balance of external and internal virtual work for a static analysis can be expressed as in which δL int is the virtual strain energy contribution while δL ext is the virtual work done by external loads. The virtual variation of strain energy in terms of displacement is expressed in a compact form as

3D geometric mapping
To describe the geometry of the structure, the global coordinate vector, x ∈ R 3 , is interpolated by a 3D shape function N 3D through the relation While noting Eq. (10), the first-order derivative of the displacement fields of general non-prismatic beam-like structures in the context of UF can be realised as Coordinate representation of non-prismatic beam-like structure [1]. To assess the performance of the proposed SUF, WUF is develop in this study. This approach allows for a reasonable comparison of the computational accuracies and efficiencies of the strong and weak formulations. Adopting the high-order Lagrange polynomial basis, L, as shape functions in the axial dimension for the WUF, the displacement fields for the WUF can be reproduced in light of Eq. (11), as Furthermore, the virtual strain energy equation can now be expressed as Substituting for the stresses in terms of displacement fields leads to the compact form Based on Eqs. (6) and (15), the generalised governing equations of the WUF reads

1D strong Unified Formulation
To obtain a strong-form of Eq. (7), we consider the explicit form of the virtual strain energy as which is rewritten in the local coordinate system as where V x∪y∪z ð Þ∈ R 3 and Ω α∪ζ∪β ð Þ∈ R 3 are undeformed volumes of the generic element in the global and local coordinate systems, respectively. Applying partial integration by parts to Eq. (19) after substituting for displacement variables while noting that the 3D Jacobian, J, is a function of ζ leads to where b J ¼ @ @ζ J j j and Γ b is the cross-section at the beam boundary. Eq. (20) in compact form reads where S ¼ CD τ and the operator D τ is the derived version of D τ after partial integration by parts and C ¼ Eq. (21) can further be expressed in terms of 1D variation along the beam axis ζ as where and Γ is the variable cross-section along the beam axis.
It should be noted that the stiffness operator K τs is a function of ζ since the Jacobian J is a function of ζ. In view of the principle of variational calculus, Eq. (24) can be rewritten as where q τ is the external load due to line and volume forces while p τ is the external load due to surface traction. Equation (25) represents the one-dimensional strong-form of the governing differential equations and boundary conditions of generic beam model which can be expanded to any order s for static analysis of non-prismatic beam structures. Moreover, K τs and Π τs are 3 Â 3 matrix operators whose explicit expression are given in the Appendix and are precomputed by means of Gaussian integration as in finite element method while the solution of the 1D displacement field u s is achieved by means of DQM.

Discretisation of the strong-form equation
A differential quadrature (DQ) approach is employed for discretising the strong-form Eq. (25) in which for a computational domain, a high-order, non-zero polynomial is used to approximate the displacement field and its derivatives as a linear combination of all functional values in the domain. Furthermore, to foster optimum accuracy and avoid numerical instabilities arising from the beam boundaries, a Lagrange basis polynomial combined with a Chebyshev nodal distribution along the beam is adopted, similar to the WUF. Accordingly, in line with the procedures outlined in [43], a continuous variable u ζ ð Þ and its nth order derivative can be approximated as where L i ζ ð Þ and a n ð Þ ij are weighting coefficients of u ζ ð Þ and its nth derivative defined by some sets of base polynomials. Readers are referred to [43] for explicit computation of L i y ð Þ and a n ð Þ ij that guarantees optimum accuracy and numerical stability. In light of Eq. (26), Eq. (25) can be generalised as After assembly of the fundamental nucleus of the SUF, the global system of equations is now given as

Computation of K ijτs andK ijτs
For non-prismatic UF, the computation of K ijτs andK ijτs contributes significantly to the overall efficiency of SUF and WUF, since, unlike the prismatic kinematic representation, K ijτs andK ijτs are continuous functions of the beam axial dimension. To assess the computational costs of computing K ijτs andK ijτs which are expressed as where axially varying cross-sectional moment parameters are expressed as To compute the cross-sectional moment parameters, 2D Gauss quadrature operations are required at each discrete point i along the beam. It should be noted that components of the cross-sectional parameters E ij and their derivatives b E ij at each discrete point i of the beam can be pre-computed within the 2D Gauss quadrature loops before expansion over the free indices i and j and summation over the dummy index k. Therefore, the computational cost of computing the SUF fundamental nucleus K ijτs at a discrete beam node i is x gp Â z gp loops, x gp and z gp being the number of Gauss points in the x-and z-directions, respectively. After computing K ijτs at N discrete points, the total cost of computing the global stiffness matrix K is In the case of WUF, the components ofK ijτs xx explicitly reads where axially varying cross-sectional moment parameters are expressed as The computational cost of computing the fundamental nucleusK ijτs at a discrete nodal point i is x gp Â y gp Â z gp loops, x gp , y gp and z gp being the number of Gauss points in the x-, y-and z-directions, respectively. Unlike the SUF, the cross-sectional moment parameters at each nodal point i must be computed within x gp Â y gp Â z gp loops plus additional N loops required to expand the moment parameters along the axial dimension. Therefore, the total computational cost of computing the WUF global stiffness matrix K ∼ is x gp Â y gp Â z gp Â N. Comparatively, the computational cost of SUF and WUF for prismatic and nonprismatic UF representations is summarised in Table 1.
Clearly, according to Table 1, the extension of WUF to nonprismatic structures increases the computational cost much more significantly than for SUF.

Numerical examples
The numerical accuracy and computational efficiency of the proposed SUF model is assessed in this section, where static analysis of homogeneous and nonhomogeneous non-prismatic beam-like structures under different type of loads is investigated. The validation and the computational merits of SUF are evaluated by comparisons with WUF, ABAQUS FE (A-FE) and an analytical solution proposed in the literature. Since SUF and WUF models for prismatic structures have comparable spectral convergence, according to [20], an Nth order Chebyshev-biased 1D DQ and FE element (see Fig. 3) combined with a fifth-order SLE element are used for structural discretisation of the SUF and WUF models respectively, while a 27-node 3D brick element is used for the geometric discretisation (see Fig. 1). For this reason, high-order accuracy is guaranteed for both SUF and WUF models. The computational time of the UF-based models (i.e., SUF and WUF) implemented in MATLAB on a 64-bit, core i7 16 GB RAM, 3.20 GHz machine is significantly influenced by the computation and inversion of the stiffness matrices K and K ∼ . Therefore, to evaluate the efficiency of the SUF model, computational times to implement these models are reported for different examples considered in this work.
To allow for comparisons with the analytical solution on the one hand, and to enable the simulation of plane strain/plane stress elements based examples in ABAQUS due to limited computational resources at the authors' disposal, the condition of plane strain is invoked for the UF-based models. To realise this condition, material coefficients C 12 , C 13 , C 14 , C 15 , C 16 C 24 , C 25 , C 26 C 34 , C 35 , C 36 C 45 , C 46 , and C 56 in Eq. (5) are set to zero. The robustness of the SUF model is then assessed by investigating the numerical stability of the strong-form solutions subject to variation of geometric parameters that influence the asymptotic stabilities of the structures. In this regard, the WUF solution is taken as the control solution, which is considered numerically stable, since the WUF stiffness matrix is symmetric and positive definite.

Homogeneous structures
The performance of SUF for modelling homogeneous linearly tapered and curved structures is evaluated in the section. The SUF estimates are benchmarked against WUF and A-FE solutions to assess the computational gains.

Asymmetric tapered beam under transverse load (example 1)
An asymmetric beam-like structure with geometric and material properties shown in Fig. 4 is subjected to distributed transverse force F z applied at the top and bottom surfaces of the beam while the end of the beam at y ¼ 0 is clamped. To obtain a converged solution for this problem, the SUF and WUF models, respectively, employ 41-node 1D DQ and FE elements with a fifth-order SLE element to characterise the beam's kinematics while only one 27-node 3D brick element is used for the geometric mesh. The converged A-FE solution is accomplished with 139,800 CPS8R 8-node plane stress reduced integration elements which amounts to a total of 559,200 DOFs. According to Fig. 5, SUF stress predictions agree excellently with the WUF as well as A-FE solutions. In all cases, the beam's bending response is well captured with accurate axial and shear stress predictions. Moreover, the normal traction at the beam's top and bottom surfaces are characterised accurately. Due to the tapered feature, the structure experiences shear traction at the top surface while the straight bottom is free from shear traction. This effect is well accounted for by the SUF model as its prediction agrees well with WUF and the A-FE solution. Thus, the proposed SUF model demonstrates excellent accuracy in this regard.
According to Table 2, although the SUF and WUF models employ the same number of DOFs to attain convergence, the computational times to obtain the SUF and WUF solutions differ significantly, with the SUF model attaining up to 84% in time efficiency. Instructively, as noted in section 3, the WUF model requires at least one order of Gauss operations more than the SUF model to compute the stiffness matrix, thus significantly increasing the time requirement of the WUF model.

Arched beam under linearly distributed load (example 2)
In this example, we examine the capability of SUF to capture the deformation of a non-prismatic beam with a curved geometric feature. An arched beam clamped at one end and free at the other is subjected a to linearly distributed shear traction as shown in Fig. 6. The converged SUF and WUF models consist of 51-node 1D DQ and FE elements, respectively, and a fifth-order SLE element to characterise the crosssectional deformation leading to 3519 DOFs. In addition, a 27-node 3D brick element is employed for geometric discretisation of the structure. Concerning the A-FE solution, 145,140 CPS8R elements comprising 440,379 nodes and 880,758 DOFs are used to achieve good convergence. According to Fig. 7, the accuracies of shear and transverse normal stress estimates by SUF and WUF are comparable, and both models agree satisfactorily with A-FE based solutions. As the structure is under shear load at the upper straight edge, the bottom arched edge bears compressive shear stress to counteract the shear load, thus satisfying static equilibrium. Furthermore, the straight top surface of the beam is free from normal traction while the bottom is induced with normal traction due to the non-prismatic effect imposed by the arched bottom surface. The SUF model satisfactorily predicts these effects with good agreement with WUF and A-FE solutions. Table 2 shows that, although SUF and WUF models are the same in terms of number of DOFs required for satisfactory convergence, the proposed SUF model demonstrates superior efficiency over the WUF model with up to 92% gain in computational time. Obviously, the 3D Gauss operations employed for WUF contribute significantly to the computational cost. Therefore, SUF proves computationally more efficient and is therefore promising for modelling arched beam structures.

Non-homogeneous structures
The SUF model is now used to investigate cases of throughthickness structural discontinuities induced by material or load distributions. Structural discontinuity may promote shear failure at the interface of the discontinuity. The induced interfacial shear is further complicated by the non-prismatic effect which can stimulate 3D-like stress-fields at the interface. Hence the need for accurate stress predictions of non-homogeneous structures.

Bi-symmetric tapered beam under discontinuous body force (example 3)
This example is indicative of general discontinuous load conditions in beam-like structures which may contain structural discontinuities Fig. 4. Cantilevered asymmetric tapered beam, like welds or cracks that stimulate different bending response at different regions of the structure. The system which consists of a beam with material and geometric properties described in Fig. 8 is subjected, separately, to normal surface and shear body forces. Although the beam is homogeneous material-wise, the discontinuously distributed loads require that the structure is partitioned into three different regions (layers). For the UF-based models (i.e. SUF and WUF), the distributed loads are applied as a body force at the middle layer. In total, 10,431 DOFs is required for the UF-based models consisting of a structural mesh of 61-node 1D DQ and FE elements and one fifth-order SLE element per layer, and a geometric mesh of a single 27-node 3D brick element. The fine structural mesh requirement for the UF-based models stems from the need to adequately capture the continuity of shear and transverse stresses at the interfaces of the transversely and shear loaded structures, respectively. Furthermore, ABAQUS FE simulation consisting of 181,760 4-noded CPS4R plane stress reduced integration elements with 1,094,426 DOFs is performed to validate the solution of the system. According to Figs. 9 and 10, it is evident that the SUF stress predictions at different regions along the beam agree excellently with WUF, and both UF-based models are in satisfactory agreement with A-FE solutions. The structure under normal body force experiences global bending without interfacial slip at the point of discontinuity as evidenced by the continuous distribution of through-thickness transverse stress. On the other hand, the shear body force induces local bending in the structure which is manifested by interfacial slip at the point of discontinuity. These through-thickness piecewise, continuously varying shear stresses, and transverse normal stresses, respectively, for the transversely and shear loaded structures are accurately captured.
In terms of DOFs, the SUF and WUF models are comparable, and demonstrate superior efficiency over A-FE solutions. However, according to Table 2, the SUF model shows improved computational time efficiency over the WUF model with up to 91% time savings due to less Gauss operations required to obtain the numerical solution.

Multi-layer bi-symmetric tapered beam under shear load (example 4)
This example follows [3], where an analytical solution was developed for multi-layer non-prismatic structures that are being employed in different engineering fields due to optimised structural response. According to Fig. 11, the clamped-free multi-layer tapered beam-like structure is subjected to a shear traction at the beam's free edge. Due to the tapered feature, beam deformation is influenced by geometric coupling effects characteristic of 3D structures. The structure is further complicated by the sandwich-like material configuration which typically affects interfacial behaviour due to different orders of deformation experienced by different layers.
Converged solutions for both SUF and WUF models are accomplished with 8721 DOFs consisting of a 51-node 1D DQ and FE element coupled with a fifth-order SLE element per layer. According to Fig. 12, stress predictions by SUF and WUF models computed at the middle and close to the free edge of the beam agree satisfactorily with the analytical solution proposed in [3]. Beam bending response, described by the linear piecewise continuous and the quadratic piecewise continuous variation, respectively, of through-thickness axial and shear stresses, is well captured by the UF-based models, except in the core where minor discrepancies are observed for the shear stress distribution close to the free edge. As usual, the SUF model shows promise in computational time efficiency over the WUF model with almost 90% gain in computational time savings, according to Table 2.   Fig. 9. Through-thickness distributions of (a) shear stress, (b) transverse normal stress of the beam with transversely loaded body force.

Computational aspects of SUF
According to previous studies by the authors [20], SUF and WUF performance for prismatic structures are comparable in terms of spectral accuracy, computational efficiency and numerical stability (measured by the condition number of the stiffness matrix), with WUF slightly edging SUF in numerical stability due to symmetry and positive definiteness of the stiffness matrix. With the introduction of non-prismatic effect, it is clear from studies described in the previous section that the computational time efficiency of SUF is superior to WUF due to less Gauss operations required. Notwithstanding this merit, it is important to assess the numerical stability of the SUF model subject to varying non-prismatic parameters that influence the asymptotic stability of the structure. This is especially necessary considering that the SUF stiffness matrix is nonsymmetric and non-positive definite with implications for numerical stability of the system solution. In comparison, the WUF solution is taken as the control solution that is numerically stable due to symmetric and positive definite properties of the stiffness matrix.
Considering examples 3 and 4, the non-prismatic parameters defined by the midspan-thickness ratio (h 0 /h m ) and the beam tapering ratio (h 0 /h L ), respectively, are varied and the conditioning of the stiffness matrix is computed with the rcond MATLAB command. According to Table 3, compared with the prismatic structure (h 0 /h L = 1), it is clear that the WUF solution for the tapered beam (in example 4) is shown to be more numerically stable than the SUF solution of the same structure, since the WUF model preserves the stability of the solution for increasing tapering ratio while the stability of the SUF solution deteriorates slightly on increasing the taper angle. In spite of this observation, Fig. 13 shows that SUF is able to excellently preserve the accuracy of the stress predictions for increasing tapering ratio showing the robustness of the SUF model. With respect to the arched beam (example 3), it is clear that the non-prismatic effect can influ-ence the stability of SUF and WUF solutions, with WUF showing a slightly improved stability over SUF. Nonetheless, the accuracy of the predictions by SUF and WUF compare satisfactorily for increasing midspan thickness ratio (see Fig. 14). Therefore, SUF is numerically robust to handle different degrees of non-prismatic effects.

Conclusion
Strong Unified Formulation (SUF) has been developed for static analysis of non-prismatic structures in the present study. Employing an anisoparametric representation of the 3D structure, the SUF structural mesh consists of a 2D SLE-based finite element in a layerwise framework and a high-order 1D differential quadrature element while the geometric mesh of the structure is characterised by a 27-node 3D brick element. A high-order weak Unified Formulation (WUF) with similar anisoparametric configuration based on a high-order 1D FE element has been further developed as a control solution to assess the computational merits of the proposed SUF. From the examples performed in this study, the UF-based models show significant computational savings in DOFs over ABAQUS FE solutions. Furthermore, the SUF demonstrates excellent numerical performance in terms of accurate stress predictions of non-prismatic beam-like structures, as SUF estimates agree satisfactorily with WUF, ABAQUS FE and analytical solutions. Instructively, the computational requirement of the SUF is characterised by a stiffness matrix computed via a 2D Gauss operations at N discrete points while WUF uses a 3D Gauss operations at N discrete points. This computational disparity significantly affects the efficiency of the UF-based models with SUF showing a high promise in computational time savings up to 92% compared to WUF-based solutions even though both SUF and WUF models use the same number of DOFs to attain convergence. Further analysis reveals that the SUF is numerically robust to model different degrees of non-prismatic effects as SUF solutions of systems with varying non-prismatic parameters are numerically stable to preserve the accuracies of the stress predictions.

Data availability statement
The raw/processed data required to reproduce these findings cannot be shared at this time as the data also forms part of an ongoing study.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.     h 0 is the beam thickness at y = 0, h L is the beam thickness at y = L, and h m is the beam thickness at y = L/2.