Study on Failure Criteria and the Numerical Simulation Method of a Coal-Based Carbon Foam under Multiaxial Loading

: Coal-based carbon foam (CCF) has broad application prospects in aerospace, composite material tooling and other ﬁelds. However, the lack of failure criteria limits its promotion. In previous studies, the failure criteria of similar materials were proposed, but there are some limitations. This paper proposes improved failure criteria based on macro-mechanical tests. Furthermore, uniaxial and multiaxial loading tests were carried out to obtain accurate failure criteria of CCF. Finally, 3-points bending tests of the CCF sandwich structure were conducted and their ﬁnite element models (FEMs) were established. The CCF test results show that the mechanical properties of CCF are transversely isotropic. The failure criteria in this paper can accurately predict the stress when the CCF fails. The error band boundary formula caused by the dispersion of the material were also given. The maximum load P max calculated by the failure surface (3684 N) was only 4.7% larger than the mean value measured by the test (3518 N), and all of the P max measured by the test (3933 N, 3640 N, 3657 N, 3269 N, 3091 N) were between the maximum value (4297 N) and minimum value (3085 N) calculated by the error band boundary formula, which means that the failure criteria have good precision.


Introduction
Coal-based carbon foam (CCF) is a high-temperature-resistant foam with good comprehensive performance [1][2][3]. Its comprehensive performance makes it suitable for the hypersonic aircraft's thermal protection system (TPS). The Parker Solar Probe has already used CCF as its TPS material [4,5], and many projects and scholars [6][7][8][9][10][11][12] have also widely studied other styles of carbon foam TPS. The thermal expansion coefficient of CCF is similar to Invar, which can be used to manufacture high-precision composite material forming tooling. In contrast, the cost of CCF tooling is about 50% lower than Invar [13,14]. CCF also has a broad application prospect in other structures [2], which makes CCF great value to research.
Failure analysis is an essential part of precise CCF structural design, which requires the clarification of the mechanical properties and failure criteria of CCF. Sihn et al. [15] measured the compression and shear properties of CCFs provided by different manufacturers and found that they all exhibit transverse isotropy: the mechanical properties in the foaming direction are significantly different from those perpendicular to the foaming direction. K. Li et al. [16] used the space frame matrix and tetrakaidecahedron element methods to establish a mesoscopic mechanical model of the three-dimensional open-cell CCF. They found that the CCF's elastic properties depended on the material's relative density, the shape and size of the ligament section, and Young's modulus and Poisson's The difficulty in establishing the failure criteria of CCF is that CCF is a kind of transverse isotropic cellular foam. Gibson et al. [24] studied the macro failure criteria of both isotropic and anisotropic reticulated vitreous carbon foam (RVC) made of resin. This study measured the stresses in each direction of the isotropic RVC under various multiaxial loads and calculated mean stress (p) and von Mises stress (q) at the time of failure to obtain the isotropic RVC failure surface in p − q space. Equivalent mean stress (p) and equivalent von Mises stress (q) were proposed according to the definitions of p and q, respectively. It is discovered that the failure surface of anisotropic RVC inp −q space has a similar shape function to that of isotropic RVC in p − q space. Due to the different mesostructures of RVC and CCF, the failure criteria of RVC cannot be applied to CCF. In addition, the proportion of shear stress inq was not thoroughly discussed in their study, which requires modification in the definition ofq.
Researches on macro failure criteria of carbon foam have been hard to find since Gibson's work. However, there are many studies on the macro failure criteria of other anisotropic foams. M. Alkhader and M. Vural [25,26] studied the macro failure criteria of ideal open-cell anisotropic foams (IOAF) under biaxial and multiaxial loading. Their criteria can accurately predict the test results of Voronoi foam. The IOAF contains no defects, so its tensile strength S t and compressive strength S c are equal. Most plastic foams are less sensitive to defects, so the ratio of compressive strength S c to tensile strength S t of plastic foam is relatively low (S c /S t = 1.09 in Voronoi foam), which can be regarded as an approximate IOAF. M. Vural et al. [27,28] studied the macro failure criteria of transversely isotropic PVC foams under multiaxial loading. They have carried out many experiments, and their results have well verified the accuracy of the criteria. Their PVC foam's value of S c /S t is 1.22~1. 24, which can be regarded as IOAF. Other scholars [29][30][31][32] also studied the criteria of plastic foams (such as metal foams, resin foams, etc.) and achieved valuable results. However, CCF cannot be considered as IOAF, whose value of S c /S t is 2.95~3.47 in our study, because CCF is a brittle foam sensitive to defects. Therefore, the criteria of CCF need further study.
According to the achievement of Gibson, the failure criteria of CCF should probably be similar to isotropic cellular foam. Nowadays, isotropic cellular foam failure criteria models are extensively and deeply studied by scholars, among which "Crushable Foam" stands out and is commonly used for foam structure failure analysis. This model is developed by Deshpande and Fleck [33]. In their study, Deshpande and Fleck [33,34] investigated the failure criteria of metallic and PVC foams under uniaxial and hydrostatic compression loading. The failure surface in the p − q space is a semi-ellipse with a central axis on the horizontal axis. Therefore, it is reasonable to assume that the failure surface of CCF in thê p −q space is also a semi-ellipse.
This study aims to establish failure criteria of CCF based on basic macro-mechanical properties that can be used for CCF structure analysis. The novelty in this paper lies in using a modified equivalent von Mises stress method to transform the failure criteria of "Crushable Foam" into CCF, filling the gap of CCF in related fields. Finally, the Abaqus finite element model, including the CCF failure criteria, was established by a user-defined subroutine, and three-point bending tests of the CCF sandwich structure were conducted to verify the criteria.

Material
The production process of CCF used in this paper was introduced in our previous study [35]. Its mesoscopic structure is presented in Figure 1, which shows that the cells are approximately circular in the section perpendicular to the foaming direction and have higher ellipticity in the section along the foaming path. Therefore, it can be assumed that the CCF we studied is a transverse isotropic foam, and the test results later also confirmed this assumption.
can accurately predict the test results of Voronoi foam. The IOAF contains no defects, so its tensile strength St and compressive strength Sc are equal. Most plastic foams are less sensitive to defects, so the ratio of compressive strength Sc to tensile strength St of plastic foam is relatively low (Sc/St = 1.09 in Voronoi foam), which can be regarded as an approximate IOAF. M. Vural et al. [27,28] studied the macro failure criteria of transversely isotropic PVC foams under multiaxial loading. They have carried out many experiments, and their results have well verified the accuracy of the criteria. Their PVC foam's value of Sc/St is 1.22~1.24, which can be regarded as IOAF. Other scholars [29][30][31][32] also studied the criteria of plastic foams (such as metal foams, resin foams, etc.) and achieved valuable results. However, CCF cannot be considered as IOAF, whose value of Sc/St is 2.95~3.47 in our study, because CCF is a brittle foam sensitive to defects. Therefore, the criteria of CCF need further study.
According to the achievement of Gibson, the failure criteria of CCF should probably be similar to isotropic cellular foam. Nowadays, isotropic cellular foam failure criteria models are extensively and deeply studied by scholars, among which "Crushable Foam" stands out and is commonly used for foam structure failure analysis. This model is developed by Deshpande and Fleck [33]. In their study, Deshpande and Fleck [33,34] investigated the failure criteria of metallic and PVC foams under uniaxial and hydrostatic compression loading. The failure surface in the p−q space is a semi-ellipse with a central axis on the horizontal axis. Therefore, it is reasonable to assume that the failure surface of CCF in the −p q space is also a semi-ellipse.
This study aims to establish failure criteria of CCF based on basic macro-mechanical properties that can be used for CCF structure analysis. The novelty in this paper lies in using a modified equivalent von Mises stress method to transform the failure criteria of "Crushable Foam" into CCF, filling the gap of CCF in related fields. Finally, the Abaqus finite element model, including the CCF failure criteria, was established by a user-defined subroutine, and three-point bending tests of the CCF sandwich structure were conducted to verify the criteria.

Material
The production process of CCF used in this paper was introduced in our previous study [35]. Its mesoscopic structure is presented in Figure 1, which shows that the cells are approximately circular in the section perpendicular to the foaming direction and have higher ellipticity in the section along the foaming path. Therefore, it can be assumed that the CCF we studied is a transverse isotropic foam, and the test results later also confirmed this assumption.

Failure Criteria for CCF
The stress-strain relationship of CCF under rectangular coordinate system before failure is where σ is the stress vector, ε is the strain vector, C is the stiffness matrix, σ i and τ j are the normal stress and shear stress, respectively, with i = x, y, z and j = xy, yz, xz. The z direction is defined as the direction of foaming, which is also the one with the strongest mechanical properties. The stiffness matrix C can be expressed as where E i , G j and ν j are the Young modulus, shear modulus and Poisson's ratio, respectively, with i = x, y, z and j = xy, yz, xz, where E x = E y , G xz = G yz and ν xz = ν yz . Gibson et al. [24] used the mean stress p and von Mises stress q used to describe the failure criteria of isotropic RVC, which are expressed as and the equivalent mean stressp and equivalent von Mises stressq based on p and q were further proposed to describe the failure criteria of anisotropic RVC, which are expressed as [24] where S ic is the compressive strength and S j is the shear strength. Tests were conducted to obtain the failure criteria of anisotropic RVC inp −q space, including uniaxial compression and tension loading, biaxial compression loading, and axisymmetric loading with all the possible combinations of axial and radial stresses. The advantage of this method is that failure criteria are characterized by the macroscopic mechanical properties of foam, which is convenient for engineers to measure and analyze the failure of foam structure in finite element software without too much research on the microscopic structure of the foam. However, the contribution of shear stress to foam failure is understudied. The shear stress item in Equation (6) is k xy 2 + k xz 2 + k zy 2 , and its coefficient is defined as strength factor k s in this research, whose value is 1 in Equation (6). This k s reflects the contribution ratio of normal stress and shear stress to foam failure, the smaller k s is, the less likely shear stress will lead to material failure, and vice versa. Assuming an extreme situation, k s = 0, Equation (6) will have no shear stress item; only the normal stress causes the foam failure, which is impossible. On the contrary, if the value of k s is vast, then the shear stress will dominate the failure, while the influence of normal stress can be ignored, which is also against common sense. Therefore, the value of k s must be carefully defined, and its value method will be discussed later in this study. As mentioned above, tests with shear stress are not included in the previous study [24], which means that it can only be proved that Equation (6) is valid in the case of multiaxial loading with only normal stress, and cannot be extended to loading with shear stress.
According to the above reasons, the modified equivalent von Mises stressq m is expressed aŝ where S it is the tension strength. In addition to the change of k s , the expression of Equation (9) in tension is also slightly different from Equation (7), because the difference between tensile strength and compressive strength in CCF huge, while it is small in RVC. When the CCF is under uniaxial compression loading, the normal stress perpendicular to the loading direction and all shear stresses are 0. In this case,p andq m always satisfŷ q m = 3p, according to Equations (5) and (8). In particular,p = 1/3 andq m = 1 when the foam fails, no matter which axis is under loading. Similarly,p = −1/3 andq m = 1 when the foam fails under uniaxial tension loading. When the CCF is subjected to compressive shear loading, all stresses except the shear stress on the loaded surface is 0, which makeŝ p = 0 andq m = √ k s /2 when the shear failure occurs. To sum up, the failure surface of CCF inp −q m space must go through the point (−1/3, 1) (1/3, 1) and (0, √ k s /2). The failure surface of the isotropic cellular foam in the p − q space is a semi-ellipse with major axis on horizontal axis [25]. Based on this, it is assumed that the failure surface of CCF inp −q m space is also such a semi-ellipse with points (−1/3, 1) (1/3, 1) and (0, √ k s /2) on it; the expression of failure surface inp −q m space iŝ where a is the shape parameter of the failure surface, whose value needs to be obtained by fitting Equation (16) with at least one additional point on the failure surface. This additional point can only be obtained by multiaxial normal stress loading, because all the failure points of uniaxial loading were used for the acquisition of Equation (16) and the value of k s is temporarily unknown. Once the value of a is fitted, the value of k s is determined.

Experiment for Failure Criteria
The test matrix and specimens are shown in Table 2 and Figure 2, respectively. Uniaxial tests were conducted to obtain the basic macro mechanical properties of CCF in different directions, including tension strength S it , compressive strength S ic , shear strength S j , Young modulus E i , shear modulus G j and Poisson's ratio ν j , with the z-axis is the foaming direction. All the tests are carried out on the hydraulic fatigue machine. In order to obtain the quasi-static mechanical properties of CCF, the loading speed was set as 1 mm/min, The Aerospace 2023, 10, 721 6 of 17 relationship between the mechanical properties of CCF and the strain rate can be studied in the future.
Uniaxial compression ASTM C365 [37] ASTM E132 [38] 30 × 30 × 20  It is noted that the standard ASTM C365 [29] proposes that it can determine whether to attach tabs to the compressed surface of the specimen according to material properties. For brittle foams like CCF, it is necessary to attach tabs to obtain real strength, because there must be a large number of broken cell wall structures distributed on the outer surface of the cut CCF specimen. The load-bearing capacity of these cell walls is significantly weaker than that of the intact cell wall. Therefore, if the specimen is directly compressed without the tabs, its damage will occur on the surface of the specimen (as shown in Figure  3a), and the measured strength is relatively lower than the actual value. If tabs are pasted on the foam, the adhesive penetrates into the surface cells and mixes with the carbon skeleton to form an adhesive layer. Since the strength of the adhesive is much greater than the strength of the foam itself, the damage will occur inside the CCF (as shown in Figure 3b), and the measured strength reflects the actual performance of CCFs with intact mesostructure. In this paper, 1 mm thick steel tabs were attached to the surface of the CCF specimen in compression tests. It is noted that the standard ASTM C365 [29] proposes that it can determine whether to attach tabs to the compressed surface of the specimen according to material properties. For brittle foams like CCF, it is necessary to attach tabs to obtain real strength, because there must be a large number of broken cell wall structures distributed on the outer surface of the cut CCF specimen. The load-bearing capacity of these cell walls is significantly weaker than that of the intact cell wall. Therefore, if the specimen is directly compressed without the tabs, its damage will occur on the surface of the specimen (as shown in Figure 3a), and the measured strength is relatively lower than the actual value. If tabs are pasted on the foam, the adhesive penetrates into the surface cells and mixes with the carbon skeleton to form an adhesive layer. Since the strength of the adhesive is much greater than the strength of the foam itself, the damage will occur inside the CCF (as shown in Figure 3b), and the measured strength reflects the actual performance of CCFs with intact mesostructure.
In this paper, 1 mm thick steel tabs were attached to the surface of the CCF specimen in compression tests.
face of the cut CCF specimen. The load-bearing capacity of these cell walls is significantly weaker than that of the intact cell wall. Therefore, if the specimen is directly compressed without the tabs, its damage will occur on the surface of the specimen (as shown in Figure  3a), and the measured strength is relatively lower than the actual value. If tabs are pasted on the foam, the adhesive penetrates into the surface cells and mixes with the carbon skeleton to form an adhesive layer. Since the strength of the adhesive is much greater than the strength of the foam itself, the damage will occur inside the CCF (as shown in Figure 3b), and the measured strength reflects the actual performance of CCFs with intact mesostructure. In this paper, 1 mm thick steel tabs were attached to the surface of the CCF specimen in compression tests.  The passive confining pressure test is a simple and effective multiaxial loading test, Kolluri et al. [40] used this method to study the failure criteria of aluminum foam, and Gibson et al. [24] used a similar method to study RVC. This test was conducted to obtain The passive confining pressure test is a simple and effective multiaxial loading test, Kolluri et al. [40] used this method to study the failure criteria of aluminum foam, and Gibson et al. [24] used a similar method to study RVC. This test was conducted to obtain the failure point with multiaxial normal stress loading, as mentioned above. During the test, the specimen is placed in a steel groove. There is an interference fit between the specimen and the groove, and the magnitude of interference is 0.05 mm. Before the specimen is put into the groove, the groove needs to be heated to 80 • C, and the inner wall is coated with lubricating oil. The testing machine applies a compressive load through a compact, as shown in Figure 4. The wall thickness of the groove is 15 mm and is made from steel, so the stiffness of the groove is much larger than that of the CCF specimen. Therefore, it can be considered that the displacement perpendicular to the loading direction of the specimen is limited to 0 during the compression process. When the loading direction is the z direction of CCF, it has ε x = ε y , and the relationship among σ x , σ y , and σ z before failure is expressed as where σ z can be obtained from the test. When the loading direction is the x direction of CCF, it has ε z = ε y , the relationship among σ x , σ y , and σ z before failure is expressed as where σ x can be obtained from the test. Because of transverse isotropy, an expression similar to Equation (13) can be obtained when the y direction is loaded.
Aerospace 2023, 10, x FOR PEER REVIEW 8 of 18 the failure point with multiaxial normal stress loading, as mentioned above. During the test, the specimen is placed in a steel groove. There is an interference fit between the specimen and the groove, and the magnitude of interference is 0.05 mm. Before the specimen is put into the groove, the groove needs to be heated to 80 °C, and the inner wall is coated with lubricating oil. The testing machine applies a compressive load through a compact, as shown in Figure 4. The wall thickness of the groove is 15 mm and is made from steel, so the stiffness of the groove is much larger than that of the CCF specimen. Therefore, it can be considered that the displacement perpendicular to the loading direction of the specimen is limited to 0 during the compression process. When the loading direction is the z direction of CCF, it has  = xy , and the relationship among  x ,  y , and  z before failure is expressed as where  z can be obtained from the test. When the loading direction is the x direction of CCF, it has  = z y , the relationship among  x ,  y , and  z before failure is expressed as where  x can be obtained from the test. Because of transverse isotropy, an expression similar to Equation (13) can be obtained when the y direction is loaded.

Loading direction
Steel groove Loading block Figure 4. The loading method of passive multiaxial compression tests. Figure 5 shows the stress-strain curves of the uniaxial loading tests. It can be seen that CCF exhibits good linear elasticity before failure, so the modulus and strength of CCF correspond to the slope and maximum stress of the linear segment in stress-strain curves,   Figure 5 shows the stress-strain curves of the uniaxial loading tests. It can be seen that CCF exhibits good linear elasticity before failure, so the modulus and strength of CCF correspond to the slope and maximum stress of the linear segment in stress-strain curves, respectively. The modulus and strength of CCF are shown in Tables 3 and 4. The tensile modulus and tensile strength in the x direction of CCF are 3.25% and 7.8% different from those in the y direction, respectively, and 33.6% and 24.3% different from those in the z direction. The compressive modulus and compressive strength in the x direction of CCF are 4.90% and 2.22% different from those in the y direction, respectively, while the tensile modulus and tensile strength in the z direction are 28.3% and 35.8% different. It is clear that the difference in mechanical properties between the x and y directions is much smaller than the dispersion. Therefore, the CCF studied in this paper can be considered a transversely isotropic material.  The modulus and strength of CCF are shown in Tables 3 and 4. The tensile modulus and tensile strength in the x direction of CCF are 3.25% and 7.8% different from those in the y direction, respectively, and 33.6% and 24.3% different from those in the z direction. The compressive modulus and compressive strength in the x direction of CCF are 4.90% and 2.22% different from those in the y direction, respectively, while the tensile modulus and tensile strength in the z direction are 28.3% and 35.8% different. It is clear that the difference in mechanical properties between the x and y directions is much smaller than the dispersion. Therefore, the CCF studied in this paper can be considered a transversely isotropic material. Poisson's ratio was obtained by simultaneously measuring the lateral and longitudinal displacements of the specimen in uniaxial compression tests with ν xy = 0.23 and ν zx = 0.56. The parameter G yz , S yz , and ν zy were not measured in this study due to the lack of specimen, and their value can be considered equal to G xz , S xz , and ν zx respectively, because of the transversely isotropy of CCF. The results of passive confining pressure tests are listed in Table 5. The stresses perpendicular to the loading direction are calculated by Equations (12) and (13). So far, the stresses in all directions of CCF failure in uniaxial and multiaxial loading tests have been obtained. The failure points inp −q m space of uniaxial tension, uniaxial compression, and passive confining pressure tests, listed in Table 6, were calculated according to Equations (5) and (8)- (10). We mark these points in thep −q m coordinate system, as shown in Figure 6. The failure surface iŝ q 2 m = a p 2 − 1/9 + 1 = −1.468p 2 + 1.163, a = −1.468 (14) The upper bound of the error band iŝ

Results and Discussion
The lower bound of the error band iŝ Then, we got k s = 2.326, which was calculated by Equation (11). So far, theq m of the shear test can be calculated. As shown in Table 6, these failure points are also marked in thep −q m coordinates. It can be found that these points are distributed near the failure surface and are all in the error band.   CCF is a brittle foam, which causes its strength to be susceptible to material defects and has a large dispersion of strength. In the −ˆm pq space, the failure points under the same loading method are distributed on a straight line passing through the origin of the coordinates. This is because the tests used in this paper belong to the proportional loading of the stress in each direction, and the relationship between p and ˆm q always maintains a linear relationship. The error we provide is very significant, and it shows that 95% of the failure points of CCF fall within the range of the error band. When designing the CCF structure in engineering, it should at least ensure that the (ˆ, m pq) point will not exceed Equation (16). CCF is a brittle foam, which causes its strength to be susceptible to material defects and has a large dispersion of strength. In thep −q m space, the failure points under the same loading method are distributed on a straight line passing through the origin of the coordinates. This is because the tests used in this paper belong to the proportional loading of the stress in each direction, and the relationship betweenp andq m always maintains a linear relationship. The error we provide is very significant, and it shows that 95% of the failure points of CCF fall within the range of the error band. When designing the CCF structure in engineering, it should at least ensure that the (p,q m ) point will not exceed Equation (16).

3-Points Bending Tests and Finite Element Model Simulation
Foam materials are often used as the core of sandwich structures under bending load in engineering. Therefore, we conducted three-point bending tests of the CCF sandwich structure and established a finite element model simulating the test process. The effectiveness of the failure criteria of CCF in analyzing the failure of the CCF structure in practice can be verified by comparing the test results and simulation results. The specimens are shown in Figure 7, designed according to ASTM C393 [41], whose numbers are BD-1~5. The thickness of CCF is 15 mm. The thickness direction of the core material is the same as the foaming direction of the foam.

3-Points Bending Tests and Finite Element Model Simulation
Foam materials are often used as the core of sandwich structures under bending load in engineering. Therefore, we conducted three-point bending tests of the CCF sandwich structure and established a finite element model simulating the test process. The effectiveness of the failure criteria of CCF in analyzing the failure of the CCF structure in practice can be verified by comparing the test results and simulation results. The specimens are shown in Figure 7, designed according to ASTM C393 [41], whose numbers are BD-1~5. The thickness of CCF is 15 mm. The thickness direction of the core material is the same as the foaming direction of the foam.  The face sheet of the CCF sandwich structure is made of carbon fiber/bismaleimide resin composite material (CFRC), the grade is T800/IS2101, whose properties are listed in Table 7. The thickness of each layer is 0.18 mm, and the layup is [0/45/90/-−45]s. The maximum load Pmax and failure mode were obtained in the tests. Finite element models (FEMs) based on the three-point bending tests were built to verify the accuracy of the failure criteria of CCF. The CCF constitutive model with failure criteria is established through the VUMAT subroutine, whose process is shown in Figure  8.  (5) and (8)- (11). When the stress state point ˆ( , ) m pq reaches or exceeds the failure surface, the subroutine judges that the The face sheet of the CCF sandwich structure is made of carbon fiber/bismaleimide resin composite material (CFRC), the grade is T800/IS2101, whose properties are listed in Table 7. The thickness of each layer is 0.18 mm, and the layup is [0/45/90/−45] s . The maximum load P max and failure mode were obtained in the tests. Finite element models (FEMs) based on the three-point bending tests were built to verify the accuracy of the failure criteria of CCF. The CCF constitutive model with failure criteria is established through the VUMAT subroutine, whose process is shown in Figure 8. The CCF failure criteria needs 19 parameters: Young's modulus (E x , E y , E z ), shear modulus (G xy , G xz , G yz ), Poisson's ratio (ν xy , ν zx , ν zy ), tension strength (S xt , S yt , S zt ), compression strength (S xc , S yc , S zc ), shear strength (S xy , S xz , S yz ) and shape parameter a. All the parameters above were obtained and can be defined in the subroutine. The Young's modulus, shear modulus and Poisson's ratio form the stiffness matrix C, according to Equation (2). The tension strength, compression strength, shear strength and shape parameter form the failure surface according to Equations (5) and (8)- (11). When the stress state point (p,q m ) reaches or exceeds the failure surface, the subroutine judges that the CCF has failed. The primary failure mode of the CCF sandwich structure under bending load is a shear failure, which makes the crack continue to expand instead of close after failure. Therefore, there is no compressive stress on the failure section of CCF, so it can be considered that the failure elements have no bearing capacity and the stiffness matrix can be set to 0. The subroutine was used in the Abaqus explicit dynamic analysis.
vent an hourglass issue in large deformation. The modulus of steel and CFRC were defined by the "isotropy" and "engineering constant" modules, respectively, and their strength were not defined cause the damage only occurred in CCF. The mesh convergence is related to the mesh density. This FEM has 16,000 elements in CFRC and 20,000 in CCF, whose result is 0.8% larger than that of FEM with 12,000 elements in CFRC and 20,000 in CCF, which means the mesh convergence can be guaranteed. The interface between the panels and the core was not damaged, so the elements at the junction of the panels and the core were set as common nodes. The contact surfaces of block-base and block-panel were all defined as a "hard contact" module. A displacement constraint was imposed on the center base to move downward at a speed of 1 mm/min. The FEM is shown in Figure 9. Eight-node brick element (C3D8R) was used in the whole model, while the stiffness-based hourglass control was used in the CCF core to prevent an hourglass issue in large deformation. The modulus of steel and CFRC were defined by the "isotropy" and "engineering constant" modules, respectively, and their strength were not defined cause the damage only occurred in CCF. The mesh convergence is related to the mesh density. This FEM has 16,000 elements in CFRC and 20,000 in CCF, whose result is 0.8% larger than that of FEM with 12,000 elements in CFRC and 20,000 in CCF, which means the mesh convergence can be guaranteed. The interface between the panels and the core was not damaged, so the elements at the junction of the panels and the core were set as common nodes. The contact surfaces of block-base and block-panel were all defined as a "hard contact" module. A displacement constraint was imposed on the center base to move downward at a speed of 1 mm/min.
A total of four, finite element models were calculated, numbered FEM-1~4. The differences among FEM-1~3 is the damage rule in the subroutine. FEM-1 used the failure surface (Equation (14)) as its damage rule, so that it would calculate the average value of P max in the CCF sandwich structure under three-point bending. FEM-2 and FEM-3 used the upper bound of the error band (Equation (15)) and the lower bound of the error band (Equation (16)) as their damage rules, respectively, so that they would obtain the maximum and minimum values of the error band of P max with a 95% confidence rate. FEM-4 replaceŝ q m withq, proposed by Gibson [24], to compare the difference in accuracy between our method and previous methods. It is worth noting that Equation (14) is obtained through the test of only normal stress loading, which has no concern with the value of k s , so the equations of the failure surface inp −q space andp −q m space were the same. Therefore, the damage rule used in FEM-4 was also Equation (14). The differences among FEM-1~4 are listed in Table 8. A total of four, finite element models were calculated, numbered FEM-1~4. The differences among FEM-1~3 is the damage rule in the subroutine. FEM-1 used the failure surface (Equation (14)) as its damage rule, so that it would calculate the average value of Pmax in the CCF sandwich structure under three-point bending. FEM-2 and FEM-3 used the upper bound of the error band (Equation (15)) and the lower bound of the error band (Equation (16)) as their damage rules, respectively, so that they would obtain the maximum and minimum values of the error band of Pmax with a 95% confidence rate. FEM-4 replaces ˆm q with q , proposed by Gibson [24], to compare the difference in accuracy between our method and previous methods. It is worth noting that Equation (14) is obtained through the test of only normal stress loading, which has no concern with the value of s k , so the equations of the failure surface in −p q space and −ˆm pq space were the same.
Therefore, the damage rule used in FEM-4 was also Equation (14). The differences among FEM-1~4 are listed in Table 8.

Validation Results
Load-displacement curves and Pmax of three-point bending tests and simulations are shown in Figure 10. FEM-1~4 overlap because the difference between these finite element models is the failure criteria used to calculate Pmax, and other parameters in the model are the same. It can be seen that all the Pmax values measured by the test (3933 N, 3640 N, 3657 N, 3269 N, 3091 N) were between the maximum value (4297 N) and minimum value (3085 N), which, respectively, were calculated by FEM-2 and FEM-3. The mean value of Pmax calculated by FEM-1 (3684 N) is only 4.7% larger than the mean value of Pmax measured by the test (3518 N), which shows that the failure criteria proposed in this paper have

Validation Results
Load-displacement curves and P max of three-point bending tests and simulations are shown in Figure 10. FEM-1~4 overlap because the difference between these finite element models is the failure criteria used to calculate P max , and other parameters in the model are the same. It can be seen that all the P max values measured by the test (3933 N, 3640 N, 3657 N, 3269 N, 3091 N) were between the maximum value (4297 N) and minimum value (3085 N), which, respectively, were calculated by FEM-2 and FEM-3. The mean value of P max calculated by FEM-1 (3684 N) is only 4.7% larger than the mean value of P max measured by the test (3518 N), which shows that the failure criteria proposed in this paper have sufficient precision. The P max calculated by the Gibson method (FEM-4) is 4808 N, which is 30.5% larger than the mean value of P max measured by the test.
The mechanism of FEM-4 producing larger errors is as follows: by substituting the stresses at the initial failure node in FEM-1 (as shown in Table 9) into Equations (5) and (6), we then obtainedp = 0.0439 andq = 0.6359 in the Gibson method, respectively. By substituting thep andq into the failure surface (Equation (14)), we then obtained q 2 = 0.4043 < 1.1602 = −1.468p 2 + 1.163. This means that when our failure criteria determine that the CCF is damaged, Gibson's failure criteria also considers that the CCF is yet to be damaged. This is due to the large difference betweenq andq m that the k s in q is defaulted to be 1 by Gibson, and the k s inq m is 2.326 according to our study. The k s represents the contribution rate of shear stress to failure. Thus, in the case of shear stress dominated failure, the failure load calculated byq must be greater than the shear load calculated byq m . This shows that the modified equivalent von Mises stressq m that we proposed can make the failure criteria have better accuracy. yet to be damaged. This is due to the large difference between q and ˆm q that the s k in q is defaulted to be 1 by Gibson, and the s k in ˆm q is 2.326 according to our study. The s k represents the contribution rate of shear stress to failure. Thus, in the case of shear stress dominated failure, the failure load calculated by q must be greater than the shear load calculated by ˆm q . This shows that the modified equivalent von Mises stress ˆm q that we proposed can make the failure criteria have better accuracy. Figure 10. Load-displacement curves and Pmax of 3-points bending tests and simulations. Table 9. Stresses at the initial failure node in FEM-1. The failure mode of the three-point bending test and its finite element model are shown in Figure 11. It can be seen that the failure mode of the CCF sandwich structure in the three-point bending test is the shear damage of the core material, and longitudinal cracks are generated in the core material and propagate to the interface between the sheet and the core on both sides. The FEM can simulate the initial damage and the propagation process. However, the longitudinal cracks in FEM-1 are along the thickness direction, while the longitudinal cracks in the test results have a certain angle with the thickness direction. This shows that, although our FEM can accurately calculate the load at the initial damage, it failed to simulate the damage propagation. This may be because the VUMAT subroutine that we used cannot simulate the mechanical behavior of CCF well after damage, and further research is needed in the future.
The shear stress xz  dominated the failure of CCF, so Figure 12 shows the contours of xz  before and after the initial damage of FEM-1. It can be seen from Figure 12a that the initial damage occurred on the lower surface under tension. It can be seen from Table   Figure 10. Load-displacement curves and P max of 3-points bending tests and simulations. The failure mode of the three-point bending test and its finite element model are shown in Figure 11. It can be seen that the failure mode of the CCF sandwich structure in the three-point bending test is the shear damage of the core material, and longitudinal cracks are generated in the core material and propagate to the interface between the sheet and the core on both sides. The FEM can simulate the initial damage and the propagation process. However, the longitudinal cracks in FEM-1 are along the thickness direction, while the longitudinal cracks in the test results have a certain angle with the thickness direction. This shows that, although our FEM can accurately calculate the load at the initial damage, it failed to simulate the damage propagation. This may be because the VUMAT subroutine that we used cannot simulate the mechanical behavior of CCF well after damage, and further research is needed in the future.
The shear stress τ xz dominated the failure of CCF, so Figure 12 shows the contours of τ xz before and after the initial damage of FEM-1. It can be seen from Figure 12a that the initial damage occurred on the lower surface under tension. It can be seen from Table 9 that τ xz is smaller than S xz (2.61 MPa) measured in the test. This is because σ x /S xt = 0.575 and the contribution of σ x to damage cannot be ignored. Figure 12b shows the stress contour of the damaged unit after initial damage. It can be found that the position of the unit with the largest τ xz moves upward, so the damaged element propagated upward. At the same time, the maximum τ xz increases (2.309 MPa), which may be caused by the stress concentration and the decrease of the shear plane in the CCF. 9 that xz  is smaller than xz S (2.61 MPa) measured in the test. This is because / =0.575 x xt S  and the contribution of x  to damage cannot be ignored. Figure 12b shows the stress contour of the damaged unit after initial damage. It can be found that the position of the unit with the largest xz  moves upward, so the damaged element propagated upward. At the same time, the maximum xz  increases (2.309 MPa), which may be caused by the stress concentration and the decrease of the shear plane in the CCF.

Conclusions
Failure criteria of core-based carbon foam are proposed in this paper. It is applied to a finite element model and can be summarized as follows: • The failure criteria can be obtained by the following steps: (a) Obtain the basic mechanical properties of CCF in all directions through uniaxial loading tests, including strength, modulus, and Poisson's ratio, and then we obtain p and ˆm q independent of shear stress; (b) Get the equation of failure surface in −ˆm pq space by fitting the results of multiaxial loading tests without shear load; (c) Calculate the s k according 9 that xz  is smaller than xz S (2.61 MPa) measured in the test. This is because / =0.575 x xt S  and the contribution of x  to damage cannot be ignored. Figure 12b shows the stress contour of the damaged unit after initial damage. It can be found that the position of the unit with the largest xz  moves upward, so the damaged element propagated upward. At the same time, the maximum xz  increases (2.309 MPa), which may be caused by the stress concentration and the decrease of the shear plane in the CCF.

Conclusions
Failure criteria of core-based carbon foam are proposed in this paper. It is applied to a finite element model and can be summarized as follows: • The failure criteria can be obtained by the following steps: (a) Obtain the basic mechanical properties of CCF in all directions through uniaxial loading tests, including strength, modulus, and Poisson's ratio, and then we obtain p and ˆm q independent of shear stress; (b) Get the equation of failure surface in −ˆm pq space by fitting the results of multiaxial loading tests without shear load; (c) Calculate the s k according

Conclusions
Failure criteria of core-based carbon foam are proposed in this paper. It is applied to a finite element model and can be summarized as follows:

•
The failure criteria can be obtained by the following steps: (a) Obtain the basic mechanical properties of CCF in all directions through uniaxial loading tests, including strength, modulus, and Poisson's ratio, and then we obtainp andq m independent of shear stress; (b) Get the equation of failure surface inp −q m space by fitting the results of multiaxial loading tests without shear load; (c) Calculate the k s according to the equation of failure surface, and then we obtain the complete equations ofp andq m ; • A large number of material-level tests were carried out to measure the basic mechanical properties of CCF, and the equation of the failure surface and its upper and lower bounds of error bands were obtained; • CCF constitutive model with failure criteria was established through the VUMAT subroutine. It was used to predict the maximum failure load and failure mode of the CCF sandwich structure under a three-point bending load. The analysis results are in good agreement with experimental results and have higher accuracy than previous methods.
In conclusion, our CCF failure criteria have sufficient precision and can be obtained from simple macro-mechanical tests, and they are also available for engineering design analysis.