Multiwall Rectangular Plates under Transverse Pressure—A Non-Linear Experimental and Numerical Study

Large deflection of rectangular plates under transverse pressure is described by Föppl–von Kármán equations, which have only approximated solutions. One of these methods is the separation into a small deflection plate and a thin membrane described by a simple third order polynomial expression. The present study presents an analysis to obtain analytical expressions for its coefficients by using the plate’s elastic properties and dimensions. To validate the non-linear relationship between the pressure and the lateral displacement of the multiwall plate, a vacuum chamber loading test is used to measure the plate’s response, with a large number of plates and length–width combinations. In addition, to further validate the analytical expressions, several finite element analyses (FEA) were performed. It has been found that the polynomial expression fairly describes the measured and calculated deflections. This method allows the prediction of plate deflections under pressure as soon as the elastic properties and the dimensions are known.


Introduction
The problem of large deflection of plates has attracted much attention since the end of the 19th century, due to its technical importance. Unfortunately, this problem has been found to be difficult to solve. The difficulties were raised at an early stage with Föppl equations for large deflections of membranes [1], where no closed-form solution was found for the rectangular membrane case. In 1910, this equation set was enhanced by Theodor von Kármán [2] to include the bending resistance of plates.
The Föppl-von Kármán equation set has challenged many researchers over the years. Nevertheless, only approximated solutions were developed, most of which are rather difficult to implement.
It was August Föppl himself who suggested an approximate approach. This approach is mentioned in Timoshenko [3] (p. 423) footnote 1, which mentions Föppl's "Drang und Zwang" [4] (p. 345). The approach is that the transverse distributed load q on the plate can be separated into two parts: q = q 1 + q 2 . The first part q 1 is balanced by the plate's bending and shearing resistance, which are calculated through the plate small deflection linear theory. The second part q 2 is balanced by the large deflections in-plane membrane forces only. Using the mid-point deflection w, this approximation is written as: The plate's small deflection coefficient A has been calculated in many previous studies, with most of them using a summation of the Fourier series. The large deflection coefficient B, however, has no exact solution, as it is ascribed to the difficult Föppl's membrane problem [1]. This expression (1) for the load-deflection behavior is quoted by many sources, such as Timoshenko [3] (p. 424) for square isotropic plates, Ugural [5] (p. 358), Wang and El-Sheikh [6] (p. 816), as well as many others.
One of them is Riber [7], who had suggested a "combined analytical solution" (see p. 71 in [7]) to find the constants in Equation (1). He used an energy method to obtain rather complex expressions for the coefficients A and B (see Equation (1)). He also presented simplified expressions for the constant B but with some internal inconsistencies. Riber [7] assumed in-plane immovable edges, which is not the case to be presented in the present study-see the BCs (Boundary Conditions) discussion later in the paper. Nevertheless, his B equation has inspired the presentation of a better expression for the coefficient B later in the present study.
Awrejcewicz et al. [8] present many efficient numerical methods that can be used to calculate specific rectangular plates without orthotropic and transverse shear behaviors.
Battaglia et al. [9] analyze orthotropic membranes and show that the load is proportional to w 3 , where the proportion coefficient can be compared to our cases.
Maier-Schneider et al. [10] found the membrane proportion coefficient with both the improved analytic energy method and the finite element method with a good agreement.
Niyogi [11] solves the simply supported orthotropic plate problem with an approximate Galerkin-Bubnov procedure. The resulting q = Aw + Bw 3 expression is verified for the isotropic plate only. The edge in-plane BCs are immovable.
Wang et al. [12] also solve the static behavior of flat glass plate and present loaddeflection data graphically. No explicit mathematical expression is given for that, although it can be extracted from the graph. Here, the in-plane BCs are totally immovable.
During the literature survey, several sources were found referring to tests and calculation methods of plates' large deflection. In order to compare the results of these papers, it was necessary to normalize the various data to a common comparable structure. The structure was a thin square isotropic plate with movable edges and an evenly distributed transverse load. The coefficients of Equation (1), A and B, were calculated considering the plate dimensions and the material properties. The result of this comparison has shown the considerable variability of the coefficient B. This variability was unexpected since most of the data was based on real laboratory tests that should respond in a similar way. This may demonstrate the fact that it is not easy to correctly measure this property. A full description of the comparison with a suggested explanation is presented in Appendix B.

Multiwall Plates
The structure of multiwall plates consists of two thin face sheets separated by an internal structure of ribs and walls. The plate is usually produced by extrusion, in which a melted material is pressed through a die with the required shape. The materials used are aluminum and various plastics. The result is a thick, endless plate with a fixed cross-section shape along the extrusion direction and width according to the equipment size. The plate is then cut to the desired length.
The plate is made of polycarbonate (PC), which is a tough transparent plastic. A typical 16mm PC plate can be seen in Figure 1.
The main application of PC multiwall plates is the glazing of architectural spaces, where both natural light and weather protection are required. As a result, the plates are exposed to wind and snow loads, which they must safely resist.
Currently, no publication describes the general performance of these plates, except manufacturer's datasheets, which are very limited to specific products and specific applications. The available publications about latticed structures relate to specific shapes, such as triangles and trapezoids, and not a general approach as presented here, which can be considered novel. The available approximated solutions for large deflections of plates are generally rather complicated, and in many cases, they involve a computational process, which is not straight-forward for field engineers. The non-linear nature of load-deflection curves is not easily represented in these solutions. Most of the research works already done do not cover the full complexity of the multiwall plates, which are shear, deformable, and orthotropic. Therefore, an engineer who needs to design a system with multiwall plates will probably face serious difficulties. This situation justifies this paper, which yields a first rough guess of these plates' performance.
To calculate the multiwall plate response to distributed load, it is necessary to know the plate's equivalent elastic properties, its dimensions (length-width), and the boundary conditions. Looking at the multiwall structure, it is obvious that the plate is orthotropic for both bending and tension, and its cross-section is transverse shear deformable. In the present study, it is assumed that all necessary equivalent elastic properties are already known. One should note that a procedure to obtain these equivalent properties of the plate is presented in Hakim and Abramovich [16].

Axes System
The axes directions are defined as shown in Figure 2, where axis x is the extrusion direction and z axis is normal to the plate surface: The origin location of the axes may be set to any convenient place.

Boundary Conditions (BCs)
For small deflections analysis, the assumption is that all in-plane stress, strains, and deflections are negligible. Therefore, the BCs, here, ignore the in-plane conditions. The most commonly used BCs are: Free (F)-no restrictions, Simply Supported (S)-no z deflection but free rotation (no bending moments), and Clamped (C)-no z deflection and no The origin location of the axes may be set to any convenient place.

Boundary Conditions (BCs)
For small deflections analysis, the assumption is that all in-plane stress, strains, and deflections are negligible. Therefore, the BCs, here, ignore the in-plane conditions. The most commonly used BCs are: Free (F)-no restrictions, Simply Supported (S)-no z deflection but free rotation (no bending moments), and Clamped (C)-no z deflection and no rotations (zero slope). S and C are the two extremes of the more complicated BC-flexible rotation support, which is rarely used.
For large deflection analysis, the in-plane BCs must be considered. The two most common BCs are: Immovable (I)-the plate edge is fixed to the support and Movable (M)the plate edges are allowed to move. It is necessary to specify both in-plane movement directions: normal to the edge and parallel to the edge. The BC used later here is SSSS-M, in which the four S stands for the four plate sides simply supported, and the M stands for the movable edges in both normal and parallel directions.
The movable (M) condition requires additional attention. When the Föppl's approximation is used, the plate in a large deflection regime is a membrane. Its deflection on movable boundary conditions should then be calculated. However, a well-known property of a membrane is that it cannot sustain in-plane compression forces, as it immediately wrinkles. However, a real plate does resist compression, as it has a bending rigidity. We, therefore, have to analyze a membrane with movable edges, which may have compression stress. Mathematically, it is possible (with the known difficulties of Föppl equations), but other practical problems would appear. Since this transversely loaded M membrane is not a common case in the literature and, perhaps, even physically not possible, no previous scientific papers that would suggest a possible solution were found. Nevertheless, an expression is suggested later in this paper.

A-B Analytical Prediction
The expressions that would predict the values for the coefficients A and B are displayed next. The variables used are: x , ν t y Tension Poisson's ratios The x, y axes origin is set as shown in Figure 2b-at the plate corner.

Small Deflection Coefficient A
The expression in (2) was derived using the Libove and Batdorf NACA Report No. 899 [17]. where The complete analysis description is presented in Appendix A.

Large Deflection Coefficient B
As presented above, the coefficient B describes the plate's large deflection response, with the plate being considered as a thin membrane. Normal membranes cannot have a movable (M) BC, so it is complicated to find previous publications displaying an expression for the membrane deflection. Wang and El-Sheikh [6] (p. 816) have analyzed an isotropic rectangular plate and have presented an approximated expression for q = Aw + Bw 3 for M BCs: Modifying this expression to an orthotropic plate and using only the B term suggests the expression: where the coefficient k includes all the numerical factors. Since (5) is based only on the first term of a multiple term series, the result is not accurate enough to correctly represent the plates response. Therefore, finite element analyses (FEA) of 80 plates with various lengths and widths were performed. The FEA software was Femap with NX-Nastran version 2021.1 with its built-in non-linear static analysis. The FEA results are given in Appendix C, while the FEA information is given in Appendix D. Typical plate arrangements and mid-point load-deflection graphs are shown in Figure 3.
The values of B in (1) were calculated with a least-squares regression from the FEA results, and the following expression is suggested for B: Note that the sum of the moduli and the sum of the length-width represent its averaged values, as the 2 division is included in K. The (a/b) is the aspect ratio of the plate. Since (5) is based only on the first term of a multiple term series, the result is not accurate enough to correctly represent the plates response. Therefore, finite element analyses (FEA) of 80 plates with various lengths and widths were performed. The FEA software was Femap with NX-Nastran version 2021.1 with its built-in non-linear static analysis. The FEA results are given in Appendix C, while the FEA information is given in Appendix D. Typical plate arrangements and mid-point load-deflection graphs are shown in Figure 3.  The values of B in (1) were calculated with a least-squares regression from the FEA results, and the following expression is suggested for B: Note that the sum of the moduli and the sum of the length-width represent its averaged values, as the 2 division is included in K. The (a/b) is the aspect ratio of the plate. The

Vacuum Chamber Test
To check the multiwall plate response to transverse distributed load, a vacuum chamber test was used, as presented in Figure 4: A 35 mm-thick wooden frame encloses a rectangular space with the required dimensions. The multiwall plate is freely placed on the frame's edges. A thin plastic sheet covers the entire device and the floor near-by. A variable-speed vacuum cleaner is connected to the internal space through a drilled hole. The vacuum created causes the plastic sheet to seal all air leakages, allowing the vacuum level to gradually increase.
An electronic vacuum sensor measures the vacuum level through another hole in the wooden frame. The sensor is connected to a controller, which displays the data. Additionally, an ultrasonic distance sensor is placed 20 cm above the plate mid-point, measures the plate deflection, and transmits it to the controller. The vacuum units are [Pa], and the distance units are [mm].
The controller has a zero button to zero the vacuum and the distance values as the test starts. During the test, the vacuum level is gradually increased, while both load and deflection are recorded. Various local buckling phenomena are also closely monitored and recorded.
The plate edges are free to move on the wood frame, creating the M movable BC. As the vacuum level increases during the test, the thin plastic sealing sheet experiences tension and presses the plate edge to the wood frame. This may change the BC from SSSS to be somewhat closer to CCCC.
Several tests were performed at the Krumbein Structures Laboratory, Faculty of Aerospace Engineering-Technion. Next, a typical test is presented.

Vacuum Chamber Test
To check the multiwall plate response to transverse distributed load, a vacuum chamber test was used, as presented in Figure 4: A 35 mm-thick wooden frame encloses a rectangular space with the required dimensions. The multiwall plate is freely placed on the frame's edges. A thin plastic sheet covers the entire device and the floor near-by. A variable-speed vacuum cleaner is connected to the internal space through a drilled hole. The vacuum created causes the plastic sheet to seal all air leakages, allowing the vacuum level to gradually increase.
An electronic vacuum sensor measures the vacuum level through another hole in the wooden frame. The sensor is connected to a controller, which displays the data. Additionally, an ultrasonic distance sensor is placed 20 cm above the plate mid-point, measures the plate deflection, and transmits it to the controller. The vacuum units are [Pa], and the distance units are [mm].
The controller has a zero button to zero the vacuum and the distance values as the test starts. During the test, the vacuum level is gradually increased, while both load and deflection are recorded. Various local buckling phenomena are also closely monitored and recorded.
The plate edges are free to move on the wood frame, creating the M movable BC. As the vacuum level increases during the test, the thin plastic sealing sheet experiences tension and presses the plate edge to the wood frame. This may change the BC from SSSS to be somewhat closer to CCCC.
Several tests were performed at the Krumbein Structures Laboratory, Faculty of Aerospace Engineering -Technion. Next, a typical test is presented.

•
Type: 10 mm-PC double wall • Nominal Area Weight 1700 g/m 2 The measured results were least-squares fitted to the expression q = Aw + Bw 3 and were compared to the theoretical curves.
Part of the tested plates data are shown in the following Figure 5, Tables 1 and 2:  The measured results were least-squares fitted to the expression q = Aw + Bw³ and were compared to the theoretical curves. Part of the tested plates data are shown in the following Figure 5, Tables 1 and 2: (a) (b)   The coefficients A and B were theoretically calculated with the expressions (2), (6), and they were compared to the measured values. The comparison is shown in Table 2: As displayed in Table 2, the calculated coefficients comply to the measured one with some differences.

FEA Results
To find the response of the plates to transversal uniform load, 80 finite element analyses were performed: 4 plate types with 4 various widths and 5 different lengths (see in Figure 3 one of the tests). The coefficients A and B for each plate were both theoretically calculated using Equations (2) and (6), respectively, and they were also found from the graphs drawn using Equation (1). A comparison of the theoretically calculated A and B coefficients and the FEA-measured coefficients is presented in Figure 6. The legend at (b) applies to all other graphs. As displayed in Table 2, the calculated coefficients comply to the measured one with some differences.

FEA Results
To find the response of the plates to transversal uniform load, 80 finite element analyses were performed: 4 plate types with 4 various widths and 5 different lengths (see in Fig. 3 one of the tests). The coefficients A and B for each plate were both theoretically calculated using Equations (2) and (6), respectively, and they were also found from the graphs drawn using Equation (1). A comparison of the theoretically calculated A and B coefficients and the FEA-measured coefficients is presented in Figure 6. The legend at (b) applies to all other graphs.  The 45 • lines represent the location of the perfect agreement between the theory and measurements. As it is shown in Figure 6, a very good agreement between the theory and analysis is found.

Large Number of Loading Tests
One of the manufacturers of PC multiwall plates is Plazit-Polygal (Plaskolite). During the years 2001-2002, they performed a large number of vacuum loading tests similar to the one described here. Plazit-Polygal has allowed the authors to use the test data for research and publication purposes. This permission is very much appreciated.
From about 250 tests, 120 tests of plates 6-16 mm thick were chosen. Each test has a set of measured load-deflection values for various width-length measurements. The analysis of every test was a linear least-squares regression that calculated the coefficients A and B in the expression q = Aw + Bw 3 . The results are listed in Appendix C. The coefficient of determination R 2 (goodness-of-fit) in all tested cases was above 0.99. These very good fits support the validity of the suggested A, B large deflection approximation (Equation (1)).
The equivalent elastic constants and moduli of the plates were measured and given with the test data in Appendix B. This information allows for the calculation of the A, B coefficients according to the theory above. The measured and calculated A, B coefficients are compared in Figures 7 and 8 The 45° lines represent the location of the perfect agreement between the theory and measurements. As it is shown in Figure 6, a very good agreement between the theory and analysis is found.

Large Number of Loading Tests
One of the manufacturers of PC multiwall plates is Plazit-Polygal (Plaskolite). During the years 2001-2002, they performed a large number of vacuum loading tests similar to the one described here. Plazit-Polygal has allowed the authors to use the test data for research and publication purposes. This permission is very much appreciated.
From about 250 tests, 120 tests of plates 6-16 mm thick were chosen. Each test has a set of measured load-deflection values for various width-length measurements. The analysis of every test was a linear least-squares regression that calculated the coefficients A and B in the expression q = Aw + Bw³. The results are listed in Appendix C. The coefficient of determination R² (goodness-of-fit) in all tested cases was above 0.99. These very good fits support the validity of the suggested A, B large deflection approximation (Equation (1)).
The equivalent elastic constants and moduli of the plates were measured and given with the test data in Appendix B. This information allows for the calculation of the A, B coefficients according to the theory above. The measured and calculated A, B coefficients are compared in Figures 7 and 8  In Figures 7 and 8, the horizontal axes are the A and B measured coefficients, while the vertical axes are the theoretically calculated coefficients. The 45° lines represent the Generally, the theory-measurement agreements are better for the coefficient A over the coefficient B, and they are also better for higher thickness over lower thickness.
Since Polygal's tests were performed more than 20 years ago, the tested samples are not available for verification anymore. Many dimensions of the plates were missing in the records, so the standard values were taken from the plate's data sheets. Nevertheless, the actual real plate dimension values almost always deviate from the standard values, as may occur in real manufacturing. These deviations can be rather significant and, for sure, influence on the plate rigidities. This is probably the main source for the inconsistencies in the reported data.
As the vacuum was measured at that time with a water manometer, which is not accurate enough, it is very possible that errors do exist in the data. These errors can be seen in Appendix C, Tables A4-A7, where the measured coefficients A should monotonically increase while the length values decrease, but in several cases, they unexpectedly decrease.
By being very comprehensive, with a large number of length/width combinations, the understanding of the plate response to transverse pressure still has some value, besides proving the q = Aw + Bw 3 response.
Note that the last plate test (10 mm Plate Faculty Lab) was performed more recently, under a better-controlled environment, and therefore, it presents a more accurate agreement.

Conclusions
Large deflection of multiwall plates, under distributed transverse pressure and SSSS-M boundary conditions, has been found to comply with the q = Aw + Bw 3 approximation rule with very good R 2 (goodness-of-fit) values.
The suggested expressions for A and B coefficients have good agreement with FEA results, while in actual tests, they appear to be more applicable for thicker plates.
The deviations in the presented experimental loading data may relate to the measurement technique. Good laboratory practice would lead to more accurate results.

. Small Deflections Coefficient A-Libove's and Reddy's Solutions
A multiwall plate has considerable transverse shear flexibility in the width direction y, while it is relatively rigid to shear in the length direction x. Nevertheless, the analysis here considers both directions, omitting the shear in the length direction only at the end of the derivation.

Appendix A.2. Libove's Solution
In the following, Poisson's ratios are designated µ x , µ y as they appear in the source, while in the main paper, they are designated ν b x , ν b y . Libove's NACA Report No. 899 [17] presents the PDEs describing this problem.
The three main equations are on page 144 of [17], between (12) and (13), without numbers. The notations are shown in [17] on page 140.
The following Equations (A1)-(A3) are simply copied from [17]: In order to have it in a more readable form, we change the shear rigidity D Qx to be S x and the D Qy to be S y . After opening the parenthesis, we get: Some simplification can be done by taking N x , N y , and N xy to be zero because no membrane forces exist in our problem.
After doing so and repeating Equations (A5) and (A6) in Equations (A8) and (A9), we get: Equations (A7)-(A9) are a linear set of three PDEs. It should be solved simultaneously for the unknowns Q x , Q y , and w for the independent load q. All the other coefficients are already known before.
According to [17], this set can be separated to be three dual variable PDEs for every unknown in terms of the load q. One of them is for w, which is a sixth order equation, as shown in [17] After replacing D y µ x with D x µ y in the above, according to [17] (8), Equation (A10) becomes: or in a shorter form: where the coefficients K i are: The following is a standard Fourier series solution of a linear PDE. The deflection boundary conditions are all w = 0, so we can assume a double sine series as the solution (Navier Solution): where W mn are coefficients to be determined. Substituting (A16) into (A14) yields: This suggests that the equation's right-hand side should also be expanded into a double sine series: Since (A19) must exist at all points x, y in the domain, so the coefficients of sin mπx a sin nπy b must be zero for every m and n. This yields: For an evenly distributed load q(x,y) = q o , the coefficients q mn are: where the coefficients K i are given in (A15). Equation (A22) presents the linear relationship between the load and the deflection in a small defection case: w max = A·q 0 (A23) Appendix A.3. Reddy's Solution J.N. Reddy [18] (p. 368), has also solved this problem. The numerical results of Reddy [18] are identical to those of Libove's NACA 899 [17].
Nevertheless, to complete this comparison, the details of Reddy's solution are repeated here with their original notations.
The elastic constant definitions used by Reddy are in Libove's terms: There are several easing assumptions that simplify the original Reddy's expressions:
No initial in-plane forces 3.
Static state-no changes in time 4.
A 11 The Navier solution with double Fourier series is: Substituting the solution and the load into the equations above yields [18] (10.2.9) (p. 368):ŝ 11 W mn +ŝ 12 X mn +ŝ 13 Y mn = Q mn s 12 W mn +ŝ 22 X mn +ŝ 23 Y mn = 0 s 13 W mn +ŝ 23 X mn +ŝ 33 Y mn = 0 (A27) where: and α m = mπ/a, β n = nπ/b. Coefficients b are now defined: Fourier coefficients are [18] (10.2.14) (p. 369): This the end of Reddy's text, the rest is a consequential result. For an evenly distributed load q(x, y) = q 0 , the coefficients Q mn are: At the plate s mid-point at location The deflection at that point is the maximal deflection w max which is: As in Libove's solution above, Reddy's solution (A33) presents the linear relationship between the load and the deflection in a small deflection case: w max = A·q 0 .

Appendix A.4. Simplified Libove's Solutions
In multiwall plates, the x-direction shear rigidity is often very high: S x → ∞ . This causes the coefficients K i to change: This change simplifies the calculation a little bit: For orthotropic plates, which are rigid for shear deformation in both directions, (A36) reduces to: while, for isotropic plates, it further reduces to: Several sources present information describing the large deflection of these plates. In order to compare the findings, it is necessary to normalize the results to a common format. Many articles present the following normalization: Which leads to: CPT: Note that the small deflection linear Classical Plate Theory (CPT) of Navier's solution for this case, with large number of terms in the summation, has: Timoshenko [3] (p. 110) states the first summation term coefficient only: N 1 = 240.38. Yankelevsky D. et al [19] (2017, Hebrew language) offer simple approximated solutions for three BCs, including the one stated here. For a single degree of freedom model SDOF, (p. 7), the plate deflection is: Walter D. Pilkey [20] (2005) p. 1001 offers an approximated solution for rectangular plate: 16a 4 q Then for a square plate: Levy Samuel [21] (1942), in NACA Report 737, solved the problem with Fourier series, creating an infinite system of non-linear algebraic equations to be solved. The truncation makes the solution an approximated one. He did that (manually!) for the first several terms for deflections up to 3.6 times the thickness. He also states that higher deflection requires more terms in order to converge accurately enough. The process is rather complex and is not easy to implement.
Using the data presented in [21] ( Table 6 at  Ishizaki Hatsuo [22] (1972) did many actual loading tests of flat glass deflections up to 10 times the thickness or breakage. The results were transferred to an Excel sheet for analysis, and the resulting expression is: The first number N 1 (small deflection coefficient) is reasonable, but the membrane coefficient N 3 is rather low, indicating a more flexible plate than in other research. ASTM E 1300 [23] (2009) "Standard Practice for Determining Load Resistance of Glass in Buildings" supplies data for flat glass under load. Using the same thicknesses and glass properties as in Ishizaki above, as well as deflections up to five times the thickness, we get the following load expression: The expression found in his work is: Eh a 4 w 3 (A49) Kaiser Rudolf [25] (1936) performed an actual test on 600 × 600 × 3.15 mm steel plate up to 2.57 times the thickness. Although only one load-deflection point is declared, some intermediate points are implicitly given. After correcting the wrong reported load data with the supplied water manometer information, the following expression was calculated: Chia Chuen-Yuan [26] (1980) (p. 64 expression 2-46, p. 65 expression 2-50) is similar to Levy [21], using the first term only (a) and eight terms (b), and shows expressions, for this case, with deflections up to two times the thickness:       Appendix D.6. Post-Processing Femap allows us to obtain the load-deflection graphs at the requested nodes. The graph data are then transferred to Excel, which produces graphs as shown in Figure 3. The Excel sheet also performs theoretical least-squares regression analyses to find the coefficients A and B for each plate. Figure A1 below shows a typical FEA screen where the mesh, BCs, loads, and deflections are presented. Figure A1. Femap typical screen.