Development of an engineering methodology for non-linear fracture analysis of impact-damaged pressurized spacecraft structures

Motivated by the dramatic worsening and uncertainty of orbital debris situation, the present paper is focused on the structural integrity of the spacecraft pressurized modules/pressure vessels. The objective is to develop an engineering methodology for non-linear fracture analysis of pressure wall damaged by orbital debris impact. This methodology is viewed as a key element in the survivability-driven spacecraft design procedure providing that under no circumstances will the “unzipping” occur. The analysis employs the method of singular integral equation to study the burst conditions of habitable modules, although smaller vessels containing gases at higher pressures can also be analyzed.


INTRODUCTION
he hazard from orbital debris is a growing international concern for the safety of space-based infrastructure.The series of incidents happened in last six years demonstrated that only one or two collisions can drastically change the orbital debris population.With space activity continuously running and expanding, the rate of collisions in space also increases, leading in turn to a new reality for the orbital debris environment where all functioning spacecraft are under higher risk than they were designed for.Because of the very high impact velocities and possibility to fail catastrophically when impacted, the pressurized modules and pressure vessels play a crucial role in spacecraft survivability.They are identified as the most critical components on-board spacecraft [1][2][3].Historically, considerable amounts of resources have been used to developing anti-orbital-debris shielding for such structures to avoid the pressure wall damage.However, the shielding cannot protect the spacecraft from all types of debris.Nowadays, the pressurized modules and high pressure tanks of the most heavily shielded spacecrafts are able to withstand the impact of debris up to one centimeter in diameter.The orbital debris between 1 and 10 cm in size which is too small to be tracked but large enough to cause the shielded pressure wall perforation, poses the highest risk for the spacecraft mission.As it was demonstrated experimentally, the case of both shield and pressurized wall perforation presents a potential for the pressure wall failure in an abrupt fashion [1][2][3][4].The answer to the question whether the spacecraft pressurized structure would undergo "unzipping" due to the impact of undetectable debris is crucial for the mission success or failure.Essentially, it quantifies the spacecraft survivability.Fig. 1 illustrates the survivability-driven design logic where it is assumed that impact of undetectable debris between 1 and 10 cm in size has occurred and the pressure wall is damaged.This design concept requires that when developing spacecraft, all attempts be made to prevent the accidental spacecraft breakups.The mitigation and protection measures are assessed for effectiveness through the fracture analysis (Fig. 1, T module 5).In the event that a pressure wall is predicted to "unzip", the survivability improvements can be achieved by adding more effective shielding or/and by varying the design parameters of the pressurized module.New protection measures will be evaluated by repeating the steps in the above design procedure until the "no rupture" conditions will be verified.The analysis of interaction of penetrative particles with equipment inside a spacecraft is out of scope of the current paper.The present paper is focused on the engineering methodology which allows determining the border between the simple perforation and catastrophic fracture of impact-damaged pressurized modules and pressure vessels onboard spacecraft.This methodology is viewed as a key element in the survivability-driven spacecraft design procedure providing that under no circumstances will the "unzipping" occur.Addressing this problem will not only improve the survivability of spacecraft itself but also will provide the mitigation effect since each satellite break-up causes not only the loss of space assets but the considerable addition to the orbital debris population.

MODELING OF IMPACT DAMAGE
xperimental studies have shown that the impact damage has the form of a hole surrounded by a zone of the cracklike defects (Fig. 2a, b, c).For the case of both shield and pressure wall perforation the impact damage varies from the petal hole (Fig. 2a) to the "cookie-cutter hole" (Fig. 2b).The perforation of the unshielded wall is accompanied by a zone of spall cracks adjacent to the impact hole as shown in Fig. 2c.For further analysis it suggested to model the cracked area around the penetrated hole by two radial cracks emanating from the rim of the hole perpendicularly to the applied load.The diameter of the model hole is equal to the diameter of the impact hole (D hole ), while the length of the fictitious radial cracks is bounded by a damage zone (D crack ).In cylindrical pressurized structure these two radial cracks are set to be normal to the hoop stress, i.e. along the expected fracture path (Fig. 2d).The penetration process lasts for a matter of microseconds and this process is essentially dynamic.After the appearance of an impact hole in the pre-loaded plate, the field of stress distribution around this hole does not change immediately.This transition process flows as the stress wave travels away from rim of hole.The evolution of the stress field near the hole in the perforated plate were evaluated explicitly using the Autodyn® code (Fig. 3).The obtained results are consistent E with the numerical solution of the non-steady-state problem of Kirsh [5].During the transition process the dynamic stress concentration factor K  ()=  () / increases reaching the maximum value of 3.33 and then asymptotically drops to the static value.

Solution of Singular Integral Equation
he problem of potential fracture and bursting of aerospace pressurized structures was extensively examined by the NASA Advanced Fracture Mechanics Group [1][2][3].The fracture propagation analysis was conducted analytically using the linear elastic fracture mechanics approach and numerically employing the finite element method and non-linear fracture mechanics technique.Comparison to the experimental data showed that the linear elastic fracture mechanics methods are too conservative and non-linear fracture mechanics approach is required for a more realistic treatment of the problem [1].We assumed that a single hole with two radial cracks is located in the infinite plate made of an isotropic elastic perfectly plastic material, the zones of plasticity are localized along the crack prolongations and the compressive stresses within the plastic zones  pz are equal to the tensile yield limit  y (Fig. 4a).The problem can be formulated in terms of a singular integral equation (SIE) of a form: where t and x are coordinates of the points on the crack contour L , p(x) is a crack surface traction.The unknown function g(t) is expressed using the Muskhelishvili's complex variable formulation [6] via the displacement dcontinuity across the crack contour L: Here u, v are the contour displacement components in x and y directions, respectively, is the shear modulus, E is the modulus of elasticity,  is the Poisson's ratio, is the elastic parameter for the plane stress and  T singular integral equation technique is a powerful alternative to the finite element method in the non-linear analysis of crack propagation which provides very rapid convergence of the numerical results [7][8][9][10][11].Fig. 5 summarizes the procedure of non-linear fracture analysis which employs the method of singular integral equation and includes the following basic steps: Modules 1-2: The analysis starts with specifying the design and material characteristics of the pressure wall and determining the impact hole parameters.
where H reflects the total number of links within the crack; q is the current number of link;  is the angle of link inclination and l is half of the crack link.Also, the symmetry of the problem and link angular positions (1=2=0,  3 = 4 =π) were taken into account.
Module 4: Unlike the finite element method the method of singular integral equations is free of mesh generation and only nodes are needed.The Chebyshev's nodes with normalized coordinates  and  changing from -1 to 1are generated on each link of the contour (Fig. 4c, d) where  .The open circles indicate the points ξ 1 ,.., ξ N on the crack faces where displacements are calculated.The closed circles correspond to the traction nodes η 1 , .., η N-1 .The normalized coordinates  and  change from -1 to 1.
Module 5: An efficient approach to account for the jump discontinuities of traction applied to the crack faces was proposed by Savruk [7].Following [7] the Eq. ( 1) for the case of 5-link crack is replaced by the system of singular integral equations: where M 00 ,..,M 02 , M 10 ,.., M 12 , M 20 ,..,M 22 are normalized kernels,  0 ,..,  2 are normalized derivatives of the displacement discontinuity across the crack contour  0 ()=g 0 ' (l 0 ),.., 2 ()=g 2 ' (l 2 ), and  0 ,.., 2 are normalized tractions  0 ()=p 0 (l 0 ),.., 2()=p2(l2).The last equation of the system (3) represents the condition of single-valuedness of displacements for the crack contour.Also, the symmetry of the problem and link angular positions are taken into account.Module 6: The numerical solution of the system of singular integral Eq. ( 3) is obtained by the method of mechanical quadratures [7].Functions 0(), 1(), 2() are sought in the class of functions unbounded at the ends of intervals where   u   are applied to complete the system of Eq. (3).By applying the Gauss-Chebyshev quadrature expressions the system of singular integral equations is transformed to the closed system of linear algebraic equations with 3N unknowns where N is number of the Chebyshev nodes.Module 7: The solution of closed normalized and linearized system of equations is obtained by Gauss elimination.

Calculation of length of the plastic zones
Modules 8-9: Once a solution of the linearized system of equations is obtained, the stress intensity factor (SIF) at the end of the plastic strip can be evaluated by Modules 8-9-10: The stress at the crack tips is considered to be finite.The unknown length of the plastic zones is determined from the condition that the stress intensity factor is equal to zero at the end of the plastic strip:   2 0 I K l  .The search is performed by golden section method.

Calculation of crack tip opening displacement/angle
Module 11: Once a numerical solution of the singular integral equation is obtained, the displacement can be calculated at any point on the crack faces.For the arbitrary point *  ( ) Analogously we obtain the expressions for * * 0 0 ( ) g x and * * 1 1 From Eq. ( 2) in the symmetric case we have Integrating we obtain the relation: where n is a segment number.The constants of integration C n are determined by displacement at the end of the segment: Thus the crack opening displacement (COD) for the segment Ln is defined as following Since for the plane stress (1 ) 2 4 ae G E   , the expression for * ( ) The presented technique allows determining the crack opening profile for the entire crack (Fig. 6) and calculate the opening displacement (CTOD) specifically at the crack tip: Fig. 7 illustrates the convergence of the numerical procedure.It allows calculating the crack tip opening angle as well.
Modules 12-14: The critical crack tip opening displacement is used as a fracture criterion (CTODc).Once the value of CTOD has been determined and compared with the value of CTOD c it is possible to answer the main question if there is a case of simple perforation without crack growth from the impact hole or crack propagation and subsequently catastrophic rupture.We have thus obtained the complete solution of the problem.

NUMERICAL RESULTS
his section gives the numerical examples which illustrate the application of the method of singular integral equations for the structures with cracks or crack-like damages.The Fig. 8 illustrates the evolution of the crack tip opening displacement after an impact hole was suddenly introduced in the loaded plate made of aluminum alloy 2024.Once CTOD has reached the critical value, the crack starts to propagate.The estimated speed of crack propagation in the metal (V cr ) varies in a range of 0.2 c 0 to 0.29c 0 , where c 0 is the speed of sound [12][13][14].For the calculations it was assumed that Vcr  0.27c0.

T
It is known that the ratio of the radial crack length (Lrad.cr.) to the hole diameter (Dhole) has a considerable effect on the critical stress.Fig. 9 illustrates that the singular integral equations method allows obtaining the accurate result for any specific case of (L rad.cr./D hole )-ratio.The obtained results also illustrate the fact that for L rad.cr./D hole >0.25, the hole with two radial cracks can be considered as a straight crack.In order to verify above method and illustrate its application, numerical calculations were performed and compared with results of impact and tensile tests of the 3-mm thickness specimens fabricated from alloy 2024 with ultimate tensile strength σ u =446 MPa, yield strength σ y =370 MPa, modulus of elasticity E=70000 MPa, Poisson's ratio ν=0.33 and fracture toughness K c = 53.9MPa m 1/2 .The critical CTOD was determined assuming the plane stress state and using the relation CTODc=(Kc 2 )/(σyE).To account the strain hardening effects the σy was interpreted as an average of the nominal yield stress and ultimate tensile strength.The computational analysis predicted residual strength to within 5% of the experimental data given in [4].The Fig. 10 and 11 illustrate a fair agreement of the obtained computational results with test data [15] where the specimens were perforated by 0.5 Ball projectile at ballistic velocities of 206-308 m/s and then subjected to the tensile tests.The specimens with thickness of 4.8 mm and dimension of 460×910 mm were fabricated from 7075-T6 alloy.The following input data were used for the analysis: ultimate tensile strength σu =535 MPa, yield strength σy=468 MPa, modulus of elasticity E=72000 MPa, Poisson's ratio ν=0.33 and fracture toughness K c = 63 MPa m 1/2 for transverse grain and K c = 81.6MPa m 1/2 for longitudinal grain.The Tab. 1 presents a comparison with the computational results obtained by the finite element method [1] to quantify the critical crack length in the cylindrical pressurized module experiencing 68.6 MPa hoop and 34.3 MPa longitudinal stresses respectively.The numerical analysis was performed for 2219-T87 aluminum alloy shell with the following parameters: σ u =430 MPa, σ y =343 MPa, E=73800 MPa, ν=0.33, wall thickness t s =3.17 mm, toughness at the crack initiation Kic= 68 MPa m 1/2 and toughness at the maximum load Kc max= 92 MPa m 1/2 [1].The comparisons shows that the computational results obtained by the finite element and singular integral equations methods are in a good agreement.The numerical experiments on the reinforced habitable modules of the International Space Station showed the "unzipping" of the pressure wall is unlikely.

Method
Critical crack length, mm

CONCLUSIONS
he present paper is focused on the engineering methodology which is viewed as a key element in the spacecraft design procedure providing that under no circumstances will the "unzipping" occur.A model of crack propagation in impact-damaged pressurized aerospace structure is presented.The non-linear fracture analysis is performed by the method of singular integral equations.Comparisons of the calculated results with the test data and numerical results obtained by finite element method showed good agreement.The suggested SIE-based approach is concluded to be effective way of assessing the fracture behavior of the impact damaged aerospace pressurized structures.

Figure 2 :
Figure 2: Modeling the impact holes: a) petal hole; b) "cookie-cutter hole"; c) hole with adjacent spall cracks; d) model of impact hole.

Figure 3 :
Figure 3: Snapshot of the evolution of the stress field after the hole was instantly formed in the loaded plate
for the plane strain.

Module 3 :
The piecewise traction distribution p(x) is applied to the crack surface as it is shown in Fig.4-b.It divides the contour into 5 portions (links) L0, L1, L2, L3 and L4, where each piece of the traction function is differentiable throughout each individual link.The traction-free link L 0 corresponds to the hole, links L 1 and L 3 are radial cracks and links L 2 and L 4 represent the plastic zones.The solution of the singular integral Eq. (1) must satisfy the condition of single-valuedness of displacements for the crack contour: n=0, 1, 2) expressing u 0 (), u 1 (), u 2 () in terms of the Lagrange interpolation polynomials over the Chebyshev nodes:

Figure 5 :
Figure 5: Steps of the fracture analysis.
L2 we have the following expression: of the function u2(ξ) in terms of Lagrange interpolation polynomials over the Chebyshev's nodes we obtain the expression for the function * * 2 2

Figure 8 :
Figure 8: Evolution of the crack tip opening displacement.