Thermal buckling of thin injection-molded FRP plates with fiber orientation varying over the thickness

The different thermo-elastic properties of glass fibers and polymer matrices can generate residual thermal stresses in injection-molded fiber-reinforced plastic (FRP) objects. During cooling from mold to room temperature, these stresses can be relaxed by large deformations resulting from an instability of the unwarped configuration (i.e., buckling). This article investigates the thermal buckling of thin FRP disks via an analytical formulation based on the Foppl-von Karman theory. Expanding on our previous work, cylindrical orthotropy with material parameters varying over the disk thickness is assumed in order to account for thickness dependency of the glass fiber orientation distribution. A disk parameter generalizing the thermal anisotropy ratio for homogeneous orthotropic disks is introduced and its relation with the occurrence and periodicity of buckling is discussed. This is done for a skin-coreskin model, for which the core-to-total thickness ratio is defined. For fiber orientation distributions typical of injection-molded disks, it is found that there exists a value of the thickness ratio for which no buckling occurs. It is also demonstrated that the periodicity of the first buckling mode is described by the generalized thermal anisotropy ratio, thus extending the results obtained for a homogeneous fiber orientation distribution. Improvements in the accuracy of the predictions for experimental data available in the literature when using the skin-core-skin model are shown. Finally, we study the relation between buckling temperature and disk thickness and propose an expression for the dependence of the normalized buckling temperature on the thermal anisotropy ratio. Results of FEM simulations are used to validate the proposed expression, proving its applicability and accuracy.


Introduction
Fiber-reinforced plastics (FRPs) are a class of materials consisting of a thermoplastic matrix reinforced by short glass fibers. They are often processed by injecting the molten polymer and the glass fibers into a mold (injection molding). The specific process parameters determine the orientation of the fibers in the final product, which in turn translates to its thermo-elastic properties (see for example Tseng and Osswald (1994)). As the injection-molded part cools down, it tends to shrink while building up stresses (Adhikari et al. (2016)). The problem of plates with residual thermal stresses has been widely investigated (see Hetnarski (2014) for a comprehensive review). Residual stresses can be relaxed by sudden large deformation, a phenomenon referred to as buckling. In general, it is possible to identify a parameter describing the magnitude of the internal stresses. Buckling occurs when this parameter exceeds a certain critical value, the buckling threshold.
In a previous work (Gualdi et al. (2021)), we studied the thermal buckling of thin plates with uniform fiber orientation distribution through the wall thickness. A framework for the buckling analysis based on the von Kármán theory with constitutive linearity in the thermo-elastic material behavior (Hooke-Duhamel model) was formulated. The framework was specified for the particular case of a free disk with cylindrical orthotropy and uniform thermo-elastic properties through its thickness. A general technique to determine post-buckling deflections based on a perturbation approach and the Fourier-Galerkin method was proposed. It was found that the periodicity of the first buckling mode for orthotropic disks is fully determined by the anisotropy ratio of the elastic moduli (or equivalently by the ratio of the coefficients of thermal expansion), in agreement with experimental evidence available in the literature Koyama (1994, 1996a,b)). By averaging the thermo-elastic properties of the material over the disk thickness, the model can be applied to systems where the fiber orientation distribution is a function of the vertical coordinate. However, a more general formulation accounting for thickness-dependent properties is required to properly investigate their effects on the buckling behavior.
The orientation of the glass fibers in injection-molded products with thin walls has been extensively investigated. Regardless of the mold geometry, the fibers tend to align in the flow direction or perpendicular to it (transverse direction) during the injection process depending on their distance from the mold walls. For flat parts, the orientation can be investigated by taking cross-sections perpendicular to their mid-planes. It is reported in the lit-erature that the fibers are generally parallel to the walls with predominant alignment along the transverse direction close to the center and along the flow direction close to the surface. This is shown for reference geometries in Oumer and Mamat (2012) (rectangular cavity), Lionetto et al. (2021) (tensile bar), Hamanaka et al. (2017) (rectangular plate), and Kikuchi and Koyama (1996a) (disk), and for a complex geometry in Baranowski et al. (2019) (bearing part for automotive applications). Variations in the predominant orientation are often observed close to the lateral edges of the mold. This effect is highlighted for fiber-reinforced PA66 disks in Kikuchi and Koyama (1996a), where it is shown to become stronger for thicker disks and remain secondary for thinner geometries. In this work, we will assume that the in-plane variations in the fiber orientation distribution can be neglected, i.e., we will consider geometries with a large lateral-dimension-over-thickness ratio. When this assumption is satisfied, the fiber orientation distribution results in a 2D orthotropic material model with parameters varying through the thickness (vertical direction z). The material properties in the in-plane directions of orthotropy at a given z are chosen such that the macroscopic material behavior resulting from the complex in-plane distribution is reproduced. An equivalent description is that of a continuous sequence of infinitely-thin orthotropic layers z = constant. In the following, we will refer to this representation as orthotropic stack or stack model, in contrast to orthotropic single layer or single-layer model when the material parameters are constant through the thickness. The latter yields the formulation of our previous work as recalled in the paragraph above. Even though the material properties in the stack model are allowed to vary arbitrarily though the thickness, they are usually approximated with piecewise constant functions. The skin-core-skin three-layer model is one of the most used representations. More accurate descriptions can generally be obtained by adding a shell layer between the core and the skin, resulting in the five-layer skin-shell-core-shell-skin model. Techniques to translate continuous fiber orientation distributions to such models can be found for example in Lionetto et al. (2021) and Hamanaka et al. (2017). These stack models have been introduced here as approximations of injection-molded geometries, but they can be also used to describe exactly composite plates made of orthotropic layers. In fact, the general formulation provided in our previous work can be applied to such problems too.
The buckling behavior of products with material properties varying through their thickness has been widely studied for composite materials. Often, it is investigated how the buckling parameter depends on the layer properties, in order to limit the amount of warpage by specific design of the stack build-up. It is generally found via analytical considerations based on linear elasticity that the ratios between the thicknesses of the various layers are crucial in the determination of the onset of buckling. This is shown for example in Lopatin and Morozov (2008), where a symmetrical sandwich panel with orthotropic core subjected to in-plane compression is considered.
When the stack is made of orthotropic layers, also the ratio between the elastic moduli in the directions of orthotropy becomes a relevant parameter. Such result can be found again in Lopatin and Morozov (2008) but also in Wankhade and Niyogi (2020), where the buckling of composite plates with orthotropic laminae is addressed via numerical methods. While layer thickness and level of orthotropy can be easily tuned in laminated compos-ites, in injection-molded products they are the result of the flow path and process conditions and thus harder to control. Nevertheless, they are expected to play analogous roles in determining the buckling characteristics of the product. The relation between layer thickness, level of orthotropy, and observed warpage upon cooling is discussed for injection-molded disks in Kikuchi and Koyama (1996b). Based on experiments and numerical simulations, the authors correlate the critical buckling temperature to Rα , the ratio between the coefficients of thermal expansion in the radial and tangential directions averaged over the disk thickness. This parameter is in fact a combination of the layer thickness ratios and levels of anisotropy, and the two contributions cannot be analyzed separately with the proposed approach. A modified linear correlation between Rα and the normalized buckling temperature is also proposed.
The present work investigates the relation between thickness-dependent fiber orientation distributions and buckling behavior in thin injection-molded disks via equivalent stack models. The effects of layer thickness and orthotropy level are investigated independently and summarized by a disk parameter δ, which is a generalized thermal anisotropy ratio and plays an analogous role to that of Rα in Kikuchi and Koyama (1996b). The formulation and solution technique presented in our previous work (Gualdi et al. (2021)) is used. When comparing the results to those of the single-layer formulation, an improvement in the accuracy of the predictions is observed when using a skin-core-skin formulation. Analytical considerations are used to improve the modified linear correlation between averaged thermal anisotropy and buckling temperature proposed in Kikuchi and Koyama (1996b). The range of validity of the new expression is extended and the fit with the available numerical results is improved.
The paper is organized as follows. The Föppl-von Kármán equations for a stack of orthotropic circular layers are derived in Section 2. The solution technique is also recalled in the same section. In Section 3, the role of the generalized thermal anisotropy ratio δ is discussed for a skin-core-skin geometry. Results are compared to those of the single-layer model from our previous work. Finally, the correlation between normalized buckling temperature and thermal anisotropy ratio is studied and an improved fitting function is proposed.

The Föppl-von Kármán equations for continuous orthotropic disks
Let us consider a thin disk of thickness H and radius R. Let the center of the coordinate system coincide with the center of the midplane of the disk.
The radial, tangential, and normal (or vertical) directions are indicated with r, ϕ, and z, respectively. We assume that the material behavior is described by the Hooke-Duhamel law (linear thermo-elasticity) and that the thermoelastic properties are symmetric with respect to the midplane of the plate.
Thermal stresses are taken proportional to Θ = Θ (z) = Θ (−z), the (symmetric) temperature difference due to cooling down of the disk, during which the temperature decreases from its higher initial value T i of the hot mold to its final lower value T f (the colder room temperature). Since large deflections w can occur during cooling, the mechanics of deformation is best described by the Föppl-von Kármán theory. A general coordinate-free formulation of the Föppl-von Kármán equations for thin anisotropic and inhomogeneous plates was derived in Gualdi et al. (2021). A solution technique based on perturbation schemes and Galerkin-Fourier approximations was also presented.
The equations were then specified for the particular case of a cylindrically orthotropic disk with constant thermo-elastic properties through the thickness.
Here, we want to describe disks with thickness-dependent fiber orientation distributions. We will still assume that each cross-section z = constant is cylindrically orthotropic, but we will allow the material properties to be functions of z. Under these assumptions, the Föppl-von Kármán equations take the same structure as the particular homogeneous case in Gualdi et al. (2021): (1) The expressions of the differential operators∆ 2 1 ,∆ 2 2 , and [·, ·] are given in Appendix A and differ from the ones in Gualdi et al. (2021) only in the material-dependent coefficients. The same holds for the four boundary conditions of vanishing normal and shear resultants, bending moment, and generalized shear force. Their explicit expressions are given in Appendix A. The unknown deflection w and Airy stress function χ (both functions of ρ := r/R and ϕ) are dimensionless. They can be obtained from their dimensional forms w dim and χ dim via the following relations: The previous relations and the differential operators in (1) involve geometrical averages and ratios of the thermo-elastic properties in the two directions of orthotropy. Their definitions are analogous to those in Gualdi et al. (2021) and are here just recalled: For each value of z, ω = γ = 1 represents isotropic behavior, ω, γ < 1 corresponds to a predominant radial fiber orientation, and ω, γ > 1 to a tangential one. As the present model is formulated in terms of the composite thermoelastic properties, it can also be applied when anisotropic properties are the result of the presence of different reinforcements (e.g., platelets) and/or processing conditions (e.g., flow orientation).
Since the structure of the system and its boundary conditions are the same as for the more specific case described in Gualdi et al. (2021), the same solution technique discussed there can be here applied directly. Therefore, w and χ are sought for by expanding them in powers of the small parameter ε as follows: where the parameter ε is defined as Here, Θ c is the minimum temperature difference required for buckling (still unknown) andΘ is the average of Θ(z) over the thickness. The dimensionless buckling parameter µ and its critical value µ c are defined by a more general version of Equation (40) from Gualdi et al. (2021): The coefficients a 0 , b 0 , c 0 , k 1,0 , and k 2,0 are functions of the thermo-elastic parameters of the material only and their expressions are given in Appendix A. For a temperature decrease,Θ > 0. However, the sign of its normalization constant depends on the considered fiber orientation distribution and thus both positive and negative values of µ c are of interest. We also define the dimensionless parameters λ and δ as As it will be demonstrated later, they will determine the buckling mode and its occurrence just like ω and γ do in the uniform case.
After solving the system (1) order by order, the approximate solution is fully determined up to the second order and reads The complete expressions for the solutions at each order are given in Appendix B.

Results
Many commercial FEM software packages allow to specify thicknessdependent fiber orientation distribution tensors through a layer approach: the geometry is described as a stack of discrete layers, each characterized by a 2D material model. In the case of thin injection-molded parts, a commonly used discretization is the symmetric skin-core-skin geometry. This representation has proven successful to capture the mechanical behavior of parts with continuous through-the-thickness fiber orientation distributions while keeping the amount of required input parameters limited. In the following, we will thus focus on such three-layer geometry. The current model however is more general and can be applied to geometries with any number of layers.
Let h be the thickness of the core layer and (H − h)/2 that of each skin layer (H is thus the total disk thickness). We define a thickness ratio α as: The values of the thermo-elastic properties in each layer are constant and will be denoted with the subscripts S or C when referring to the skin or core layer, respectively. First, we will draw some general conclusions regarding the occurrence of buckling; then we will compare the predictions to those of the equivalent single-layer model; finally, the dependence of the buckling temperature on thickness and thermal anisotropy ratio will be discussed.

Occurrence and periodicity of buckling
The occurrence of buckling is determined by the pre-buckled stresses, which are described by the zeroth-order Airy stress function: as follows from Equation (B.2) with µ = µ(Θ) according to Equation (9), while the parameters λ and δ were introduced in Equation (10). When the material properties are such that χ (0) ≡ 0, the disk behaves macroscopically as if it were isotropic and no stress is built-up during cooling. In our previous work, the condition for no buckling was formulated in terms of ω (ratio of elastic moduli) and γ (ratio of the coefficients of thermal expansion) for uniform orthotropic disks. Their values follow the same trends with respect to fiber orientation: ω, γ < 1 for a radial alignment, ω, γ > 1 for a tangential alignment, ω = γ = 1 for an isotropic distribution. Therefore, the two parameters are interchangeable in discussing the occurrence of buckling.
When multiple layers are considered, λ (generalization of ω) and δ (generalization of γ) become independent, meaning that λ = 1 does not imply δ = 1 and vice versa. It can be seen from Equation (14) that the stress field does vanish for δ → 1 ± , irrespective of the value of λ. Therefore, δ is identified here as the relevant parameter to describe the occurrence of buckling. Since its definition involves all the thermo-elastic constants and their geometrical averages, the overall isotropic behavior is thus the result of all these factors combined. Consequently, if one were to control the degree of orientation in the layers or their thickness, the macroscopic isotropic state would be reached and buckling would be avoided. This is numerically shown in Figure 1, where all the thermo-elastic properties are kept constant (values provided in Table   1) and the thickness ratio α is varied from 0 to 1. A predominant tangential fiber orientation is assumed in the core and a radial one in the skin layer.
Three cases are analyzed: in case I, the orientation is stronger in the core, and conversely stronger in the skin in case III; case II describes a balanced configuration, with the same degree of orientation in all layers. The buckling temperature Θ c is plotted as a function of α. The value α iso of the thickness ratio for which no buckling occurs (isotropic thickness ratio) can be identified as the asymptote in the buckling temperature. On the second axis, the corresponding value of δ is displayed. As expected, δ → 1 when the isotropic thickness ratio is approached.
The equivalence of (ω, γ) and δ extends to the determination of the buckling mode too. In the uniform case, axisymmetric buckling (coffee-cup) was observed for ω, γ < 1, whereas non-axisymmetric buckling of period π (saddle) occurred for ω, γ > 1. Analogous results hold for the skin-core-skin formulation: buckling is axisymmetric for α < α iso (when δ < 1) and nonaxisymmetric of period π for α > α iso (when δ > 1). These results hold irrespective of the value of λ.
The analogy between the uniform and multi-layer geometries and the results discussed above are summarized in Table 2.

Comparison with single-layer model
In Kikuchi and Koyama (1994), the type and period of buckling was linked to the anisotropy in coefficients of thermal expansion Rα = k r /k ϕ .
Experiments were performed by injection-molding disks of PA66 with various reinforcements and the macroscopic orthotropic thermo-elastic properties were measured. In our previous work, it was demonstrated that the proposed model for uniform orthotropic disks is able to correctly capture the buckling mode and excellent quantitative agreement with the results of linear thermoelastic FEM simulations was shown.
Additional simulations were performed in Kikuchi and Koyama (1996a) to investigate the influence of the disk thickness on the buckling temperature.
The fiber orientation distribution was first derived by solving the flow field in the mold and then used as input for the simulations. The values of the thermal parameters averaged over the thickness of the disks were provided, together with the elastic constants corresponding to the highest degree of orientation. The calculations showed that the buckling temperature increases with the disk thickness. The thinnest disks (1.5 to 3 mm) buckled axisymmetrically (coffee cup), whereas the thicker disks (4 and 5 mm) did not show any warpage, as their buckling threshold was higher than the applied temperature difference. These results were already analyzed in our previous work with the model for uniform disks (single-layer formulation) by using the measured thickness-averaged coefficients of linear thermal expansion. For the elastic moduli and Poisson's ratios, the values corresponding to the highest degree of orientation were used. Since the actual material properties vary through the thickness and show less orientation than in the ideal case of perfect alignment, the buckling temperature was underestimated and the deflection was overestimated. The fiber orientation distribution through the disk thickness was studied for the thinnest (1.5 mm) and thickest (5 mm) disks. Since no buckling was predicted for the thickest disk, we focus on the thinnest one. As explained in the next paragraphs, we will translate the fiber orientation distribution to an orientation tensor for a three-layer geometry. The tensor will then be used to estimate the thickness-dependent elastic parameters needed as inputs for the present model.
The fiber alignment is described in Kikuchi and Koyama (1996a) by the layer averaged orientation angle φ , defined as the volume average per layer of the fiber angle with respect to the radial direction. When φ = 0, all the fibers in the layer are radially oriented, whereas when φ = π/2 they are all in the tangential direction. In Figure 2, φ is plotted against the normalized thicknessz = z/(H/2) for the 1.5mm-thick disk. As φ = π/4 corresponds to an average fiber orientation of 45 • , approximately 15.4% of the disk shows a predominantly tangential fiber alignment and the remaining 84.6% a radial one. This implies that in the corresponding skin-core-skin model one should use α = 0.154. To better describe the fiber orientation distribution in each layer, additional information about the distribution of the angles around their averages should be considered. However, in the following we will assume that all the fibers in a layer share the same orientation angle φ . This is done in order to maintain the focus on the validation of the proposed model. As it is demonstrated in Advani and Tucker (1987), the symmetric orientation tensor in cylindrical coordinates for a planar uniform fiber orientation distribution with angle θ with respect to the radial direction has elements α ij given by α 11 = cos 2 θ; α 22 = sin 2 θ; α 12 = 0.
Let φ S and φ C be the averages of φ over the skin and core layers, respectively (see Figure 2). By assuming that the orientation of all fibers is given by either of the two averages, the orientation tensor in each layer can be computed directly from Equation (15) by taking either θ = φ S or θ = φ C .
Once the orientation tensor has been estimated, the thermo-elastic constants can be calculated. This is done numerically in Digimat-MF by providing the matrix and fiber thermo-elastic properties, the filler content, and the orientation tensor. Linear thermo-elastic behavior is assumed. Mori-Tanaka ho-mogenization is used in the software to calculate the resulting composite behavior. The inputs used in these calculations are summarized in Table 3 and are extracted from Kikuchi and Koyama (1996a) and material datasheets, unless otherwise specified. The outcome of this procedure is shown in Table   4. The estimated constants are used as inputs for the three-layer model with α = 0.154. In agreement with previous results, the model correctly predicts axisymmetric buckling. The buckling threshold and magnitude of warpage at the disk edge are reported in Table 5 and compared to the FEM results and our previous single-layer analysis. Despite the approximation introduced in the estimation of the thermo-elastic parameters, the skin-core-skin geometry significantly improves the quantitative agreement of the predicted buckling temperature difference with the FEM results. The magnitude of the deflection is underpredicted, which is believed to be caused by the conservative calculation of the coefficients of thermal expansion for the composite. A better match can be obtained by improving the calculation of the composite properties, for which more information about the fiber orientation distribution is required.

Buckling temperature, thickness, and anisotropy
The dependence of the buckling temperature difference (∆T ) b on the disk thickness h (Θ c and H in our notation, respectively) and on the average thermal anisotropy ratio Rα is discussed and quantified in Kikuchi and Koyama (1996b) by means of FEM simulations for disks of varying thickness and fiber orientation distribution.
To investigate the thickness dependence, the fiber orientation distribution was kept fixed while h was varied. It was observed that the ratio (∆T ) b /h 2 is constant. This property is typical of composite laminates, where the laminar thickness can be changed freely without affecting the thermo-elastic properties. Conversely, injection-molded disks of different thickness will also show a different fiber orientation distribution and hence material properties.
In the following, we will thus focus on composite laminates.
In our formulation, the buckling temperature for a laminate of total thickness h L is given by Equation (9) (after some rearrangements): All the terms on the right including µ c depend on the material parameters and relative layer thickness only: the integrals involved in the definitions are already normalized by the appropriate power of H and the ra-tioD/H 3 is independent of H too. This confirms the numerical results of Kikuchi and Koyama (1996b). The previous equation also highlights the dependence on the disk radius R, as Θ c · R 2 is independent of R.
The second part of Kikuchi and Koyama (1996b) focuses on the relation between (∆T ) b /h 2 L and the average thermal anisotropy ratio Rα L . By performing FEM simulations of disks with various thicknesses and fiber orientation distributions, it was shown that linear correlations exist between (∆T ) b /h 2 L and Rα L in the intervals 0.26 ≤ Rα L ≤ 0.75 and 1.78 ≤ Rα L ≤ 2.20, separately. This is shown in Figure 3 together with the Rsquared values for each fit. The two linear functions where combined in the following general expression: The values of the coefficients a, b, and c in each Rα L interval are given in the original publication. The exponential term was introduced to reproduce the divergent behavior of the buckling temperature around Rα L = 1. We will now show that our formulation allows to derive analytically a different relation which improves the fit to the experimental data and includes the asymptotic behavior near Rα L = 1.
Since the attention is put on the average thermal anisotropy ratio, the disk can be analyzed with our previous single-layer formulation, where the ratio of the coefficients of thermal expansion was indicated by γ 2 ≡ Rα L .
Evaluation of Equation (16) for a single layer yields: where C is a function of the ratio of the elastic moduli ω, the geometric average of the Poisson's ratiosν and of the coefficients of thermal expansioñ k, and of the disk radius R. In principle, the value of C depends on the considered data point of Figure 3, as each simulation refers to a different fiber orientation distribution. This is evident when comparing data points for Rα L < 1 and Rα L > 1: (∆T ) b /h 2 L should be always positive, implying that in the first case C > 0 whereas in the second C < 0. The lack of detailed information about the thermo-elastic parameters in each simulation does not allow to better quantify the C dependence on them (which is also partially hidden in the numerical calculation of µ c ). However, let us assume that the dependence is weak in the two regions Rα L < 1 and Rα L > 1 separately, so that the following approximation holds: The constants C 1 and C 2 are positive and they are determined by fitting the experimental data of Figure 3, yielding C 1 = 7.88 and C 2 = 14.21. The fitted function is plotted on the same figure together with its statistical coefficient of determination R-squared. When Rα L < 1, the fit to the FEM results is excellent, with an R-squared value higher than that of the linear fit. It can thus be concluded that the dependence of C on the thermo-elastic constants is weak for Rα L < 1 and the proposed model of Equation (19) is accurate.
When Rα L > 1, the R-squared value drops from 0.997 to 0.609, indicating that the dependence of C on the thermo-elastic parameters might be stronger than for Rα L < 1. As mentioned earlier, the statement cannot be further validated analytically as µ c is the solution of a generalized eigenvalue problem and scarce information about the fiber orientation distribution for each simulation is provided. Nevertheless, the proposed fit extends the linear one of Equation (17) to the whole interval Rα L > 1, even though more data is needed to better determine its accuracy. The divergent behavior for Rα L = 1 is already incorporated in the choice of the fitting function in Equation (19) and does not need to be accounted for with the introduction of corrective terms.

Conclusions
The thermal buckling of thin injection-molded FRP plates has been investigated via an analytical approach based on the Föppl-von Kármán theory.
The presence of short glass fibers is modeled by using cylindrical orthotropy with thermo-elastic material parameters varying through the disk thickness.
The formulation was further specified for a skin-core-skin geometry. It was shown that for injection-molded disks there exists a value α iso of the core-tototal thickness ratio α for which the disk can be considered macroscopically isotropic, i.e., no buckling occurs. The period of the first buckling mode is determined by the predominant fiber orientation as described by the generalized thermal anisotropy ratio δ: when δ < 1 (predominant radial fiber orientation), buckling is axisymmetric (coffee cup); when δ > 1 (predominant tangential fiber orientation), buckling is non-axisymmetric of period π (saddle); when δ = 1, also α = α iso and no buckling occurs. This result generalizes that obtained for uniform disks, where the same description holds by replacing δ with the elastic and thermal anisotropy ratios ω and γ, respectively.
The skin-core-skin formulation was used for model validation on a case previously analyzed with a uniform fiber orientation distribution. Improvements in the accuracy of the predictions were obtained, showing the importance of including a more detailed description of the fiber orientation distribution through the disk thickness. Based on analytical considerations, the relations between buckling temperature, disk thickness, and thermal anisotropy were explored. An expression correlating the normalized buckling temperature (∆T ) b /h 2 L to the thermal anisotropy ratio Rα L was proposed, improving the modified linear fit previously proposed in the literature.  Table 1.

Acknowledgements
This work was conducted at the DSM Material Science Center research facility in Geleen, the Netherlands. The support of DSM, and especially of Dr. Lucien Douven for the thermo-elastic calculations with Digimat-MF, is highly appreciated. This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.  Kikuchi and Koyama (1996a). Missing data close to the disk center and surface has been replaced by a constant extrapolation. Dashed lines correspond to the averages of φ over the core and skin layers (α = 0.154). Case Layer  Table 1: Thermoelastic properties for the cases of Figure 1.

A. Appendix A
All the results in this Appendix are obtained analogous to Gualdi et al. (2021).

A.1. Differential operators in polar coordinates
In the dimensionless Föppl-von Kármán equations (1) for orthotropic disks we havẽ with λ as defined in Equation (10).

B.2. First order
where w (1) is sought for by using a Legendre-Fourier expansion in combination with a Galerkin approach: