Bending of Orthotropic Scalene Triangle Plates: Finite Difference Modeling

1. Introduction

In various branches of technology (construction, mechanical engineering, aircraft building, and shipbuilding), orthotropic plates are widely used, in which the physical and mechanical characteristics of materials are symmetrical about three mutually perpendicular planes (wood, plywood, reinforced concrete). Many designs also have orthotropic properties, for example, corrugated and ribbed plates. A lot of work is devoted to the study of orthotropic plates. Recent works include [1,2,3].

A study [4] showed that the general anisotropic property and the variable in-plane stiffness property inherent in composites create significant challenges in the static and dynamic analysis of these new types of variable stiffness composite structures. Therefore, it is important to develop mechanical models and calculation methods for composite lamellar-shell structures of variable stiffness to understand their complex mechanical mechanisms and promote their further application in aerospace engineering.

Publications devoted to triangular orthotropic plates appear much less frequently. In one of the first such publications [5], large amplitude oscillations of right triangular anisotropic plates were analyzed based on von Kármán governing equations. Utilizing the Bubnov-Galerkin procedure, a nonlinear second-order differential equation for the unknown time function was drawn. The equation was solved in terms of Jacobian elliptic functions.

Vibration analysis and multi-objective optimization of stiffened triangular plate was done in [6]. The triangular plate included stiffeners (ribs), which were parallel to each other and parallel to one edge of the triangle. The governing equation of transversal deflection of the plate was obtained by considering the effects of orthotropic characteristics and external excitation. The ordinary differential equation for the system's time response was obtained using the Bubnov-Galerkin method. In the next step, a multi-objective optimization was carried out considering two conflicting objective functions, i.e., maximizing the system's nonlinearity and minimizing the amplitude of vibration. Four decision variables were considered, including the plate's thickness, geometry, and the distance of the stiffeners. Finally, the effects of different parameters on the optimal solutions and the distribution of decision variables were investigated.

Buckling analysis of general triangular anisotropic plates in [7] and [8] with different boundary conditions subjected to combined in-plane loads was considered. Solutions for plate buckling were obtained using the Rayleigh-Ritz method combined with a variational formulation. The numerical results were obtained for various triangular geometries with isotropic and anisotropic material properties. The effect of transverse-shear deformation was studied for different triangular geometries. The results confirm the importance of including the effect of transverse-shear deformation in the buckling analysis of composite plates. These results were later used in [9] and [10] for optimal design of composite grid-stiffened aircraft panels subjected to global and local buckling constraints. The local buckling of aircraft skin segments is assessed with material anisotropy and transverse shear flexibility. The local buckling of stiffener segments was also assessed.

An accurate and simplified solution was provided in [11] for the free vibration problem of simply supported thin general triangular plates. The proposed method applies to thin plates with linear boundaries regardless of their geometrical shapes. Results were compared with previously published data for the isosceles and general triangles. Good agreements were reported. Although the paper deals only with simple support conditions, the article claims that any combination of classical boundary conditions can be handled with or without complicating factors.

The p-Ritz method operates with mathematically complete two-dimensional polynomial functions, and boundary polynomial equations raised to appropriate powers are used to approximate the displacements. In the [12], the p-Ritz method was used to derive the governing eigenvalue equation for the buckling behavior of triangular plates with both translational and rotational elastic edge constraints. The influence of elastic edge supports on the buckling factors for triangular plates of various vertex angles (aspect ratios) and boundary conditions was examined. The buckling solutions for isosceles and right-angled triangular plates with elastic edge constraints were presented.

The unilateral buckling behavior of point-restrained triangular plates was studied [13]. The Rayleigh-Ritz method with polynomial-based shape functions was used to investigate the unilateral buckling behavior of triangular plates with various loading combinations, aspect ratios, boundary conditions, and point restraint configurations. Polynomials and tensionless foundations modeled the displacement functions and restraining medium, respectively. The results are obtained for different boundary conditions, aspect ratios, and various in-plane compressive and shear loadings. Convergence and comparison studies were undertaken to confirm the validity and precision of the solution method.

Free and forced multi-frequency vibrations of stiffened triangular plates with the stiffeners were investigated [14]. The governing motion equation for a triangular plate was developed based on the von Kármán theory. The nonlinear ordinary differential equation of the system using the Bubnov-Galerkin approach was obtained. Closed-form expressions for the free undamped and large-amplitude vibration of an orthotropic triangular elastic plate were presented using two analytical methods, namely, the energy balance method [15] and the variational approach. It is demonstrated that those two methods are straightforward and reliable techniques for solving those nonlinear differential equations.

A nonlinear vibration of a triangular shape plate with several stiffeners was studied [16]. The Bubnov-Galerkin method obtained the ordinary differential equation for the system's time response. A genetic-based multi-objective optimization was performed to find the geometry and locations of the plate's stiffeners.

Three-dimensional structures in the form of triangular plates of thin and medium thickness with homogeneous and inhomogeneous characteristics of materials (isotropic and orthotropic) were considered [17]. The finite difference method was applied using a grid of scalene triangles as a resolving method. Based on a numerical algorithm and the author's Fortran programs, the leading parameters were obtained from transverse and in-plane loads under various boundary conditions (in bending, buckling, and free oscillations).

In [18], free vibration characteristics of moderately thick composite materials triangular plates under multi-point support boundary conditions were analyzed. An improved Fourier series method was used. Energy equations are established based on the first-order shear deformation theory. The Rayleigh-Ritz technique is adopted to solve unknown coefficients of energy equations with the multi-point support boundary conditions. This study aimed to simulate real engineering structures with multi-point support. It was shown that the method has higher accuracy than the finite element method.

In-plane vibration of arbitrary laminated triangular plates with elastic boundary conditions was studied [19] by the Chebyshev-Ritz method. The coordinate transformation mapped the arbitrarily shaped triangular laminated plate into a square plate to facilitate energy calculation. The displacement functions of the square plate after transformation were expressed as two-dimensional Chebyshev polynomials multiplied by coefficients. The arbitrary elastic boundary conditions of the plate were obtained by changing the stiffness values of each spring with artificial virtual spring technology. The in-plane free vibration characteristics of the triangular laminated plate under different boundary conditions were calculated. The accuracy of this method was verified by comparing finite element results with experimental results. The comparison shows that the present method has good convergence and satisfactory actuarial accuracy.

The bending of cantilever triangular plates at the same inclination angles of the side edges to the base was investigated [20]. A finite difference method was applied. Combining the results was applied to solve the problem of an acute angle at the top of the plate. The results of calculating a cantilever bar of variable bending stiffness were combined with similar results of calculating a triangular plate supported along the contour using a reduction factor. This study's theoretical provisions and applied results could be partially used to investigate the bending of orthotropic scalene triangle plates supported along the edges.

The frequency and mode of free vibrations of thin isotropic triangular plates with a central hole for different boundary conditions were investigated [21]. The finite-element method was used. Some plates' topology of vibration modes was compared to square plates with hinged and clamped edges. The numerical values of the natural frequencies and modes of triangular plates are in good agreement with the experimental results.

The simplistic superposition method was used for analytic free vibration solutions of right triangular plates [22]. The problem was solved by dividing it into three subproblems that were solved by the symplectic techniques via imposing the variable separation on the Hamiltonian-system-based governing equation and symplectic eigen expansion. Then, the analytic frequency and mode shape solutions were obtained by requiring the equivalence between the original problem and the superposition. The comparison with the numerical results for the right triangular plates confirmed the approach's convergence and accuracy. For modeling channel flows, approaches based on the numerical solution of two-dimensional Saint-Venant equations have shown high efficiency and sufficient accuracy. The derivation of these equations, numerical integration algorithms, and examples of calculations are given, for example, in [5,6,7,8]. However, the Saint-Venant model is open-ended for deformable channels and requires addition, for instance, with an equation or system of equations to find variables in the time and space of bottom marks.

This paper aimed to develop and adapt the complex triangular mesh for the finite difference method. The main application is to receive the stress-strain state of orthotropic scalene triangle plates, as it wasn't researched enough. The objective is to analyze the stress-strain state of orthotropic triangle plates under different loads and boundary conditions. The other objective is the parametric analysis of orthotropic coefficients in boundary conditions and the influence of base angles on strength parameters and load-bearing capacity.

2. Methods

The object of study is a transversely bent triangular plate made of orthotropic material fixed along the edges of the plate under the action of a uniformly distributed load. Figure 1 shows the calculation scheme for the triangular plate.

The configuration of the plate is determined by three parameters: a this is the base of the triangle (plate); α, β are the angles of inclination of the side edges to the base (in the general case, α is not equal to β).

The study used the well-known initial differential equation to bend orthotropic plates [24,25,26]. The equation has the form:

$$D_1\frac{{\partial }^{4}w}{\partial x^4}+2D_3\frac{{\partial }^{4}w}{\partial x^4\partial y^2}+2D_2\frac{{\partial }^{4}w}{\partial y^4}=q,$$

(1)

where w=w(x,y) is the deflection function,

q is the intensity of the stationary transverse distributed load in the plane of the plate,

$D_1,D_2,D_3,D_k$ are bending and torsional stiffness of the orthotropic material,

x, y are Cartesian coordinates.

The flexural and torsional stiffness of an orthotropic material is determined as follows:

$$D_2=\frac{{Е}_2{t}^3}{12\left(1-{\nu }_1{\nu }_2\right)};D_3=D_1\cdot {\nu }_2+2D_k;D_k=G{t}^3/12.$$

(2)

Here:

$E_1$,${Е}_2$ are the elastic modulus of the material in two mutually perpendicular planes;

${\nu }_1$, ${\nu }_2$ are Poisson's ratios in two mutually perpendicular planes;

G is the shear modulus;

$D_k$ is the torsional bending stiffness;

t is the thickness of the plate.

Next, the orthotropy coefficients are introduced:

$${\alpha }_0=\frac{2\cdot D_3}{D_1},{\beta }_0=\frac{D_2}{D_1}.$$

(3)

Equation (1), taking into account the accepted notations (2) and (3), has the form:

$$\frac{{\partial }^{4}w}{\partial x^4}+2{\alpha }_0\frac{{\partial }^{4}w}{\partial x^4\partial y^2}+{\beta }_0\frac{{\partial }^{4}w}{\partial y^4}=\frac{q}{D_1}.$$

(4)

Bending moments and torque, taking into account the orthotropy of the structural material, have the form:

$$M_x=-D_1\left({\nu }_2\frac{{\partial }^{2}w}{\partial y^2}+{\alpha }_0\frac{{\partial }^{2}w}{\partial x^2}\right);M_y=-D_1\left({\nu }_1\frac{{\partial }^{2}w}{\partial x^2}+{\beta }_0\frac{{\partial }^{2}w}{\partial y^2}\right);M_{хy}=-D_k\frac{{\partial }^{2}w}{\partial x\partial y}.$$

(5)

The boundary conditions of plates supported along their edges are written in the same way as for bending isotropic plates (with hinged or rigidly fixed edges of the plate):

$$w_i=0;M_{n_j}=0;$$

(6)

where

j is edge numbers,

$n_j$ is the directions of the normal to the corresponding edges of the plate,

$M_{n_j}$ is bending moment in the normal direction

Table 1. Fractional composition of sand in the model.

Fraction diameter (mm)	1-0.5	0.5-0.25	0.25-0.1	<0.1
Content (%)	0.2	31.9	67.7	0.2

Table 1. Fractional composition of sand in the model.

Fraction diameter (mm)	1-0.5	0.5-0.25	0.25-0.1	<0.1
Content (%)	0.2	31.9	67.7	0.2

The cross-sectional shape of the model closely resembled half the cross-section of a trape.

3. Results

It is known that the finite difference method can be a preferable choice for problems on a uniform grid. The finite difference method with a mesh of scalene triangles was used to apply equation (4) to the bending calculation of a triangular orthotropic plate (Figure 1). The grid of scalene triangles fits well into the oblique contour of the triangular plates. Figure 2 shows a fragment of a grid of scalene triangles:

For the i^th node of a grid of scalene triangles, equation (4) will take the form (without taking into account boundary conditions), which is given in [17]:

$$\begin{array}{l}{\psi }_1{w}_i+{\psi }_2\left(w_o+w_r\right)+{\psi }_3\left(w_P+w_{1S}\right)+{\psi }_4\left(w_q+w_t\right)\\ +{\psi }_5\left(w_n+w_f\right)+{\psi }_6\left(w_b+w_h\right)+{\psi }_7\left(w_d+w_k\right)\\ +{\psi }_8\left(w_m+w_e\right)+{\psi }_9\left(w_a+w_g\right)+{\psi }_{10}\left(w_c+w_j\right)=\frac{qh{y}^4}{D_1}.\end{array}$$

(7)

In this case, mesh parameters ${\psi }_1,{\psi }_2,{\psi }_3,{\psi }_4,{\psi }_5,{\psi }_6,{\psi }_7,{\psi }_8,{\psi }_9,{\psi }_{10}$are defined as follows:

$$\begin{array}{l}{\psi }_1=6{С}^4+{\alpha }_0{С}^2\left(-6AB+4\right)\\ +{\beta }_0\left[4{(AB-1)}^2+2{(AB)}^2+2A^2+2B^2\right];\\ {\psi }_2=-4{С}^4+{\alpha }_0{С}^2\left(4AB-2\right)+{\beta }_0\left[-4AB\left(AB-1\right)+2AB\right];\\ {\psi }_3={\alpha }_0{С}^2\left(B-2A\right)+{\beta }_0\left[4A\left(AB-1\right)-2A^{2}B\right];\\ {\psi }_5={\alpha }_0{С}^2+2{\beta }_0{A}^{2}B;{\psi }_6=2{\beta }_{0}AB;{\psi }_7={\alpha }_0{С}^{2}B-2{\beta }_{0}A{B}^2;\\ {\psi }_8={С}^4-{\alpha }_0{C}^{2}AB+{\beta }_0{(AB)}^2; {\psi }_9={\beta }_0{A}^2; {\psi }_{10}={\beta }_0{B}^2.\end{array}$$

(8)

Angles $\alpha$ and $\beta$ define parameters A, B, C and U:

$$\begin{gathered} A=\frac{\sin \beta \cos \alpha }{\sin (\alpha +\beta )};В=\frac{\sin \alpha \cos \beta }{\sin (\alpha +\beta )}; \\ С=\frac{\sin \alpha \sin \beta }{\sin (\alpha +\beta )};U=C^2-AB; \end{gathered}$$

(9)

The equation in the particular case for an isotropic plate for the i-th grid node is also written according to [17] in the following form:

$$\begin{array}{l}{\phi }_1{w}_1+{\phi }_2\left(w_o+w_r\right)+{\phi }_3\left(w_p+w_s\right)+{\phi }_4\left(w_q+w_t\right)\\ +{\phi }_5\left(w_n+w_f\right)+{\phi }_6\left(w_b+w_h\right)+{\phi }_7\left(w_d+w_k\right)\\ +{\phi }_8\left(w_m+w_e\right)+{\phi }_9\left(w_a+w_g\right)+{\phi }_{10}\left(w_c+w_j\right)=\frac{qh{y}^4}{D_1},\end{array}$$

(10)

where

$$\begin{array}{l}{\phi }_1=2\left[2{(U+1)}^2+A^2+B^2+U^2\right];\\ {\phi }_2=-2\left[2(U\left(U+1\right)-AB\right];{\phi }_3=-2\left[2A\left(U+1\right)-BU\right];\\ {\phi }_4=-2\left[2B\left(U+1\right)-AU\right];{\phi }_5=2AU;\\ \hspace{1em}{\phi }_6=2AB;{\phi }_7=2BU;{\phi }_8=U^2;{\phi }_9=A^2;{\phi }_{10}={В}^2.\end{array}$$

(11)

Thus, equations (7) and (10) differ only in expressions (8) and (11), respectively.

Appropriate boundary conditions at the edges of the plate must accompany the original equation (1). For continuously supported sides of the plate along the perimeter (Figure 1), the boundary conditions have the form:

$$w_i=0;\partial w_i/\partial n_j=0;{\partial }^2{w}_i/\partial n_j^2=0,$$

(12)

where $w_i$ is the deflection at the i-th node of the plate contour;

$n_j$ is normals to the edges of the plate (Figure 1);

j=1,2 are the plate edge numbers.

Boundary conditions (11) for the i-th node of the triangular mesh are written in finite differences in group form as follows:

For the left edge of the plate (Figure 3a):

$$w_q=w_i=w_t=0;\left({\psi }_8{w}_0+0.5{\psi }_7{w}_p\right)={\gamma }_1\left({\psi }_8{w}_r+0.5{\psi }_7{w}_S\right).$$

(13)

For the right edge of the plate (Figure 3b):

$$w_q=w_i=w_S=0;\left({\psi }_8{w}_r+0.5{\psi }_5{w}_q\right)={\delta }_1\left({\psi }_8{w}_o+0.5{\psi }_5{w}_t\right).$$

(14)

For the base of the plate (Figure 3c):

$$w_o=w_i=w_r=0;\left(0.5{\psi }_7{w}_S+0.5{\psi }_5{w}_t\right)={\theta }_1\left(0.5{\psi }_7{w}_p+0.5{\psi }_5{w}_t\right).$$

(15)

The coefficients in expressions (13) – (15) took the following values:

$\left({\gamma }_1={\delta }_1={\theta }_1\right)=-1.0$ for hinge-supporting edges of the plate;

$\left({\gamma }_1={\delta }_1={\theta }_1\right)=+1.0$ for hard plate edges.

Boundary conditions (13) – (15) modify the main finite-difference equation for orthotropic triangular plates (7) and take the form of seven standard finite-difference equations for seven characteristic mesh nodes. Figure 4 illustrates the numbering of typical grid equations.

With the number of mesh partitions (N=8) with continuous support of the edges along the perimeter of the triangular plate, we obtain 21 design (intra-circuit) nodes. Figure 5 shows the numbering of some calculation nodes.

Finite-difference equations for all 21 computational nodes of the triangular mesh (Figure 5) constitute a system of linear algebraic equations of the 21st order, which in matrix form has the form:

$${\mathit{D}}_u\cdot \overrightarrow{\mathit{w}}=\overrightarrow{\mathit{P}.},$$

(16)

where ${\mathit{D}}_u$is a square matrix of the 21^st order,

$\overrightarrow{\mathit{w}}=\left|\left\{w_1,w_2,\mathrm{...},w_{21}\right\}\right|$ is the vector of the desired deflection,

$\overrightarrow{\mathit{P}}=\left|\left\{P_1,P_2,\mathrm{...},P_{21}\right\}\right|$ is the vector of free terms that considers the transverse load acting on the surface.

With a uniform distribution of the intensity of the transverse distributed load q over the plate is:

$$P_i=q{h}_y^4/D_1,$$

(17)

where $D_1$ is taken according to expression (3); $h_y=AC/8$ is the grid step along the axis (Figure 5).

Thus, the algorithm proposed by the authors for calculating orthotropic plates consists in using expressions (2), (3), (7)-(9), (13)-(15) to implement equation (16).

Using the SMath Studio program [27,28], several research problems were solved based on expression (16). The influence of:

- type of boundary conditions (using coefficients${\gamma }_1,{\delta }_1,{\theta }_1$),

- geometry of the plate (while setting the values of the angles α and β at the base of the plate), ${\alpha }_0 and {\beta }_0$- orthotropy coefficients ${\alpha }_0 and {\beta }_0$, determined by the formula (3), by the amount of nodal deflections ($w_i=1,2,\mathrm{..},21$).

Figure 6 and Figure 7 show the dependences of deflections at grid nodes δ=1, 6, 8, 10, 17, 21 (Figure 5) for equilateral (α=β=60°) triangular orthotropic plates (at α₀-=0.66; 1; 1.5; 2.0 and β₀=0.5; 1; 1.5; 2.0) for rigid and simply supported edges with the number of mesh partitions N=8.

The calculations were conducted using (4) and (5). A rectangular grid of 1800 cells measuring 0.1 m × 0.2 m was constructed for a tray with 18.0 m × 2.0 m dimensions. At the entrance boundary corresponding to the first row of cells, it was assumed that the bottom was not eroded, which matched the entrance section of the tray reinforced with a cement crust in the laboratory experiment.

Table 2 and Table 3 present the results of numerical calculations of equilateral triangular plates for various methods of pinning the edges.

$\mathrm{Deflections}\left(w_i\cdot {10}^4\frac{D_1}{q{a}^4}\right)$

According to Table 2 and Table 3, the following can be concluded.

1)

With an increase in the orthotropy coefficient α₀ (Figure 6) for different β₀ and different boundary conditions, the deflections at the nodes monotonically decrease.

2)

With an increase in the orthotropy coefficient β₀ (Figure 7) for constant values (α₀=1 ≈ const), the deflections at the plate nodes change according to a complex dependence.

3)

For an equilateral isotropic triangular plate (α₀=2, β₀=1), there are three axes of symmetry along its three medians. In this case, deflections are equal at nodes 1, 6, and 21.

4)

There is only one axis of elastic symmetry for equilateral orthotropic triangular plates. This axis of symmetry is the y-axis. In this case, the deflections are equal only at nodes 6 and 21 (Figure 6 and 7, Table 2 and Table 3).

This problem was solved for two partition numbers of the sides of the plate as 4 and 8, for assessing the mesh convergence. The results were as follows:

1)

When the number of partitions is N = 4, the deflection in the middle of the triangle median $w_4^a$ equals to 5.96.

2)

When the number of partitions is N = 8, the deflection in the middle of the triangle median $w_8^a$ equals 4.97.

The resulting convergence graph (Figure 8) is obtained using Richardson's extrapolation formula [29].

To assess the reliability of the results presented, a finite element calculation of an equilateral triangular plate with hinged support along the perimeter was performed to determine the deflection in node 8 (Figure 5) for $E_1=1.5\cdot {10}^{11}$ Pa, $E_2=1.0\cdot {10}^{11}$ Pa, ${\nu }_1=0.1$, ${\nu }_2=0.2$, $t=2.5\cdot {10}^{-3}$m, $q=1.0$kPa. The following value was obtained: $u=2.584\cdot {10}^{-3}$m. When switching to dimensionless coefficients, the deflection in the node with ${\alpha }_0=0.66$, ${\beta }_0=2$ is:

$$w_8=\frac{D_1\cdot u}{q\cdot a^4}\cdot {10}^4=5.15.$$

(18)

The uncertainty with the table value is $\delta =\left(5.15-4.97\right)/4.97\cdot 100\%=3.6\%$, which indicates the acceptable accuracy of the described calculation algorithm based on the finite difference method.

4. Conclusions

The object of study is a transversely bent triangular plate made of orthotropic material, fixed along the edges of the plate under the action of a uniformly distributed load. An algorithm for simple engineering calculation of triangular orthotropic plates based on the finite difference method is presented, which makes it possible to perform manual or machine calculations in variant (or preliminary) design with sufficient accuracy. At the same time, the technical problem of boundary conditions for the edges of triangular plates was solved by the idea of paired exclusion of deflections of contour nodes, which arises when using a mesh of scalene triangles. When using the well-known rectangular grid, this problem does not arise.

In the proposed research method, the original 4th-order differential equations with variable orthotropy coefficients based on finite-difference approximation are replaced by grid equations using a grid of scalene triangles. As a result, a mathematical transition occurs to a system of linear algebraic equations, the solution of which results in deflections at the nodes of the applied mesh, and then the corresponding internal forces are calculated.

The results obtained allow us to draw the following conclusions.

• The main and standard finite-difference equations for a mesh of scalene triangles (Figure 2) are obtained, taking into account changes in the parameters of the geometry of the plates (angles α and β), boundary conditions (coefficients γ₁, δ₁, θ₁), orthotropy of the material (coefficients α₀, β₀), number of grid partitions N.
• Boundary conditions for a mesh of scalene triangles are written in the original group (paired) form, which makes it possible to take into account the presence of oblique edges of irregular (non-equilateral) triangular plates;
• Calculation of the deflection of equilateral (α=β=60°) orthotropic triangular plates (Equation 16) The deflection of equilateral (α=β=60°) orthotropic triangular plates was studied with the number of mesh partitions N=8, providing sufficient engineering accuracy.

The results of this study expand the range of problems to be solved for triangular orthotropic plates with their geometric asymmetry (oblique, scalene plates), which is the basis for further theoretical and practical developments in the mechanics of thin-walled structures.