# Fixture Design in Flexible Tooling of Aircraft Panel Based on Thin Plate Theory

Authors: Zemin Pan; Ying Liu; Zhichao Sun; Songyang Chang; Qiang Fang
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

Flexible and conformal positioning of the skin is of great significance for the efficient and low deformation assembly of aircraft panels. As the main conformal positioning part of the skin, the fixture positioning surface is complicated to design when taking into account the flexible positioning. In this study, taking the low deformation and flexible clamping of the skin under complex stress conditions as the target, a fixture positioning profile design method is proposed to predict the deformation of the skin and provide moderate support for the skin based on thin plate theory (TPT). In the face of complex stress and boundary conditions of the skin, the differential method is used to solve the stress and deformation of the skin under different support conditions, which is correctly verified by finite element simulation. The proposed method reveals the influence of the fixture positioning profile on the deflection of the aircraft skin under concentrated load (drilling force), which provides a good design theoretical basis for the fixture design for the low deformation flexible clamping of the skin.


---

## Body

## 1. Introduction

In the fuselage, wing, and other structures of aircraft, thin-walled parts are widely used to reduce the weight of the whole aircraft. The clamping deformation and assembly deviation of this kind of thin-walled parts in the assembly stage will accumulate and spread along the assembly dimension chain of products [1] and ultimately affect the quality of aircraft manufacturing [2]. So how to design a suitable fixture for thin-walled parts is of great significance in the field of aircraft assembly.Aiming at the problem of fixture design for weak stiffness parts, Cai et al. [3] came up with the “N-2-1” positioning principle in 1996. The “N-2-1” positioning principle is more suitable for the clamping of thin-walled parts than the “3-2-1” positioning principle by increasing the number of positioning elements on the main positioning surface and enhancing the overall stiffness of thin-walled parts. The increase in the number of positioning elements causes trouble in the fixture arrangement. The intelligent optimization algorithm has the characteristics of high search efficiency and convenient mathematical description, which is commonly used in fixture positioning point layout optimization. For fixture layout optimization, different fixtures develop different solutions. Minh et al. [4] presented a single-axis rod type flexible fixture system for thin-walled components in machining processes utilizing the N-2-1 location principle and proposed a geometry-based method to find the optimal location of the workpiece relative to the fixture system, which could find the optimal location within 30 s. Chen et al. [5] proposed a “N-M” positioning principle based on the “N-2-1” principle and used the genetic algorithm to search the optimal solution of fixture layout, which reduced the machining deformation of thin-walled parts by 56.5%. Abolfazl et al. [6] used the genetic algorithm to optimize the unconstrained problem of car body parts clamping error propagation and realized the optimal fixture configuration. Ma et al. [7] proposed a hybrid optimization algorithm GAOT, which generated the global optimal fixture layout, and obtained a better optimal solution than the existing optimization algorithm. Zhou et al. [8] designed the fixture layout scheme based on the hybrid particle swarm optimization algorithm with the objective of minimizing the normal deformation of metal sheet in the flexible tooling system. Wan et al. [9] applied the fast nondominated multiobjective optimization algorithm (NSGA-II) to optimize the four defined objective functions to minimize the position error of the workpiece and obtain the optimal positioner layout. Yang et al. [10] proposed a fixture location layout optimization method based on the support vector regression model and NSGA-II algorithm and realized the optimal design of support points for aircraft skin parts. Zeshan et al. [11] proposed a “N-3-2-1” positioning principle for sheet metal fixture layout optimization and combined with finite element and genetic algorithm to realize fixture layout optimization in the multipoint spot welding process. Li et al. [12] proposed a method of fixture location layout optimization design for thin-walled parts based on the Kriging agent model and flower pollination algorithm and completed the “4-2-1” fixture layout optimization design for curved thin-walled parts. Aimed at aircraft weak stiffness parts, Wang et al. [13] calculated the maximum deformation under different positioning schemes by finite element simulation and combined with the firefly algorithm to achieve the iterative optimization of positioning point layout. Li et al. [14] proposed an optimization method for the clamping scheme of aircraft thin-walled parts based on the genetic algorithm. Combined with finite element simulation, the synchronous optimization of clamping layout and clamping sequence was realized. To sum up, in the design of thin-walled parts fixture, the application of “N-2-1” positioning principle, and the positioning point layout optimization guided by the intelligent optimization algorithm, can significantly reduce the deformation of thin-walled parts in the assembly and manufacturing process and improve the final quality of products.In the processing of thin-walled parts, the workpiece and cutting tool will have elastic deformation under the action of cutting force, which will lead to the change of process parameters and affect the processing quality. Thin plate problem is a basic problem in elastic mechanics. The key to study this kind of problem is how to deal with the boundary conditions and how to get the exact analytical or numerical solution [15–17]. Based on the thin plate theory of elasticity, the machining deformation of thin-walled parts can be analyzed. The solution of thin plate problems usually includes the finite element method, meshless method, finite difference method, boundary element method, Rayleigh–Ritz method, and Galerkin method [18]. Slark et al. [19] studied the large deformation problem of simply supported thin plates by the meshless method and solved the nonlinear equations. Based on Kirchhoff plate theory, Shabana et al. [20] proposed a 12-node plate element model, in which the displacement continuity can be ensured by rigid connection between nodes, and the motion and deformation of flexible plate can be accurately described. Chen et al. [21] proposed a method of machining error compensation for thin-walled parts based on bicubic B-spline interpolation considering the influence of elastic deformation. In order to reduce the machining deformation of thin-walled titanium alloy parts, Li et al. [22] proposed a nonuniform allowance design strategy of discrete allowance volume element based on the Rayleigh Ritz method, which makes better use of the workpiece stiffness and reduces the machining error. Tang et al. [23] solved the bending problems of Kirchhoff and Winkler thin plates by using the generalized finite difference method. Numerical experiments show that the method can effectively solve the bending problems of two kinds of thin plates under different transverse loads.In addition, as a large aviation thin-walled part, the assembly sequence of aircraft panels has a crucial impact on the assembly efficiency and assembly quality of the panels. In order to improve assembly efficiency and quality, many scholars home and abroad began to study the mechanical assembly sequence as early as 4 decades ago. The optimization of the assembly sequence in the early stage was mainly realized by manual calculation, and its defects were low efficiency and single results. With the rapid development of computer computing power and the introduction of intelligent algorithms, the optimization of assembly sequence has become efficient, flexible, and rich in calculation results [24]. Yang et al. [25] developed a novel PSOBC algorithm based on the conventional PSO algorithm for complex products; this method avoids a large number of matrix calculation by establishing synchronized assembly Petri net (SAPN) to describe the precedence relationships, and the provided PSOBC not only prevents premature convergence to a high degree but also keeps a more rapid convergence rate than the standard PSO algorithm. Bahubalendruni et al. [26] proposed a novel and efficient method to obtain all valid assembly sequences and optimized assembly sequence for a given assembled product, which considers four basic predicates, namely, “liaison predicate, geometrical feasibility, mechanical feasibility, and stability” to validate each sequence, and the validity of the method is validated by different example products. Bahubalendruni et al. [27] proposed an efficient computational method to find the optimal sequence from a huge set of all assembly sequences, and the method is proven in generating optimal solutions for any given product effectively.Based on the above references, it can be found that both scholars domestic and abroad have made a lot of achievements in the research of thin-walled parts fixture design. As a kind of large thin-walled parts, how to reduce the load deformation of aircraft skin and achieve reliable clamping is the paramount problem of aircraft panel flexible tooling development. Scholars domestic and abroad have made research on the fixture layout and deformation analysis of thin-walled parts, but the research on the fixture layout mainly focuses on the point supporting flexible tooling, and the research on the typical fixture type tooling like surface supporting is too little. Few scholars established the mechanical model to guide fixture design through elastic deformation theory of aircraft skins. In this study, in order to minimize the deformation of aircraft skins under load, the research on the design of flexible fixture is carried out.


## 2. Analysis of Fixture Design

The skin of large aircraft is usually designed with variable thickness. The thickness of the skin is reduced in the area with less stress, which forms the inner surface of the skin as concave, so as to reduce the weight of the whole skin structure. Figure1 shows the skin of the front fuselage panel of a certain civil aircraft, and the blue areas are the chemical milling area of the skin. When it comes to the flexible tooling fixture, the chemical milling area of different skin may be different, which makes the fixture surface difficult to coordinate or even share the same fixture surface. In an ideal situation, the longer the length of the fixture surface is, the longer the length of the effective bonding section with the skin is, and the better the fixture stiffness can be achieved. However, due to the location of the truss axis, the length of the clamping surface is limited within a certain range, as shown in Figure 2 Therefore, it is necessary to study the quantitative relationship between the length of the surface and the effect of skin support stiffness to guide the design of the surface.Figure 1
Schematic diagram of skin chemical milling area.Figure 2
Schematic diagram of fixture with skin fitting.In this study, aiming at the design problem of the fixture, the design goal is to ensure the effect of skin support, and the solution of thin plate problem in the elasticity is used to guide the fixture surface design according to the suppression effect of the elastic deformation of the skin.


## 3. Mechanical Model Based on Thin Plate Theory

In the assembly process, the skin is affected by gravity, strap pressure, board support force, frictional force, and drilling force, so it is often necessary to establish a complex mechanical analysis model, which increase the difficulty of solving, so the establishment of a simplified mechanical model is an important part of solving engineering problems. For the case where the skin is thin and large and only the stress-strain conditions near small load area are analyzed, the classical thin plate deformation theory is adopted, where the skin radian is ignored, the skin is considered to be an ideal elastic body, and the theoretical model is based on elasticity of a rectangular thin plate. Hence, several assumptions should be given in advance [28]:(1)
Zero radian assumption: since the stress area is small when drilling, the surface radian of the skin can be ignored;(2)
Continuity assumption: it is assumed that the object is a continuous medium and there is no gap between particles;(3)
Complete elasticity assumption: the object completely obeys Hooke’s law, and the strain of the object is proportional to the stress;(4)
Isotropic assumption: the physical properties of any point in the object are the same in all directions;(5)
Uniformity assumption: the object is composed of the same type of uniform material, so that the elastic coefficient of the object does not change with the change of position coordinates.The thin plate will be bent under the action of transverse load, and the stress, strain, and displacement problems caused by it are calculated according to the bending problem of thin plate [29]. The bending problems of thin plates are generally divided into small deflection bending problems and large deflection bending problems. After the skin is positioned and clamped on the frame, the preconnection holes between the skin and truss, and the skin and corner piece need to be drilled. In this process, the fixed effect of the strap and the clamp on the skin is regarded as the fixed boundary condition of the skin, and the skin is only affected by the drilling force. The drilling force is a transverse force and perpendicular to the neutral plane. Therefore, the stress and deformation of the skin under the action of drilling force can be solved as the bending problem of thin plate.When dealing with the problem of thin plate bending, it can be considered as the problem of small deflection bending of thin plate when the stiffness of thin plate is large and the deflection is far less than the thickness under the action of external load. The Kirchhoff–Love theory [29] should be satisfied.(1)
On any normal line perpendicular to the neutral plane, the displacement of each point in the thickness of the thin plate is the same, that is, the deflectionω, which is only related to the x and y coordinates. Let the neutral plane of the thin plate be the xy plane (Figure 3), so that(1)w=wx,y.(2)
The deformation caused by the normal stressσz is negligible, so the equation can be simplified to(2)ex=1Eσx−vσyey=1Eσy−vσxγxy=21+vETxy.(3)
There is only vertical displacement and no displacement parallel to the neutral plane in the thin plate, that means(3)μz=0=υz=0=0.Figure 3
Stress analysis of thin plate microunit.The bending deformation problem of the skin is simplified to the small deflection bending problem of the thin plate, and a differential equation forω(x, y) needs to be established. As shown in Figure 3, take out a parallelepiped microunit, whose length and width are dx and dy, respectively, from the rectangular thin plate with thickness δ. The four sides of the unit are subjected to the internal force of the plate, the upper surface is subjected to the external load qdxdy, and the lower surface has no load. On the side where x is a constant, stress components sx, τxy, and τxz act; on the side where y is a constant, stress components σy, τyx (equal to τxy), and τyz act. These stresses can be synthesized into bending moments Mx, My, torque Mxy, and transverse shear forces Fsx, Fsy, respectively.Combine the stress components per unit width of the flank:(4)Mx=∫−δ2δ2σxzdz,My=∫−δ2δ2σyzdz,Mxy=Myx∫−δ2δ2τxyzdz,Combining (2) and elastic mechanics, it can be known that the generalized Hooke’s law expression for deflection ω(x, y) is(5)σx=−Ez1−v2∂2w∂x2+v∂2w∂y2,σy=−Ez1−v2∂2w∂y2+v∂2w∂x2,τxy=−Ez1−v2∂2w∂x∂y.Introducing the bending stiffness of the thin plate (the dimension isL2MT−2),(6)D=Eδ3121−v2Substituting (5) into (4), the bending moment and torque of the microunit can be obtained as(7)Mx=D−∂2w∂x2+v∂2w∂y2,My=D−∂2w∂y2+v∂2w∂x2,Mxy=Myx=−D1−v−∂2w∂x∂y.In order to keep balance of the microunit, the sum of the moments aroundx-axis and y-axis is 0, and the force in z-direction is balanced, that means the following balance conditions must be satisfied:(8)∑Mx=0,∑My=0,∑Fz=0.After omitting the high-order trace, it can be obtained that(9)∂Mx∂x+∂Myx∂y=Fsx∂My∂y+∂Mxy∂x=Fsy,and(10)∂Fsx∂x+∂Fsy∂y+q=0.Regarding Mxy=Myx, substitute (7) and (9) into equation (10) to get(11)∂4w∂x4+2∂4w∂x2∂y2+∂4w∂y4−qD=0.Arrange (11) as follows, which is the equilibrium differential equation for bending deformation of the skin with small deflection.(12)∇4w=qD,where ∇2=∂2/∂x2+∂2/∂y2 is the Laplace operator.However, the assumption for small deflection theory is that the longitudinal displacement of thin plate, which is perpendicular to the neutral plane is so small than the thickness of thin plate that the membrane force can be ignored. For a thin metal plate such as the aircraft skin, the longitudinal displacement of each point in the plane is not necessarily much less than the thickness of the plate. If so, it can be solved as a large deflection bending problem of thin plate. When establishing the equilibrium differential equation of large deflection bending problem of thin plate, it is necessary to take into account the neutral plane strain and membrane force caused by the longitudinal displacement of each point in the neutral plane.The bending differential equation of the elastic plate is(13)D∇4w=Ftx∂2w∂x2+Fty∂2w∂y2+2Ftxy∂2w∂x∂x+q.The deflectionω and the film forces (caused by lateral load q) Ftx, Fty, and Ftxy in this equation are unknown. In order to simplify the equation and reduce the number of unknowns, the stress function Φ is introduced, so that(14)Ftx=δσx=δ∂2ϕ∂y2,Fty=δσy=δ∂2ϕ∂x2,Ftxy=δτxy=−δ∂2ϕ∂x∂y.At this point, the elastic plate differential equation becomes an equation containing two unknownω and Φ, which still cannot be solved. A compatibility equation needs to be constructed for the relationship between film force and deflection. Considering the geometric relationship between the normal strain, shear strain, and three-dimensional displacement in the microunit of the thin plate, the following geometric relationship equation can be obtained:(15)εx=∂u∂x+12∂w∂x2,εx=∂u∂y+12∂w∂y2,γxy=−∂u∂x+∂u∂x+∂w∂x∂w∂y.By eliminating the lateral displacementsu and v in the equation, the following compatibility equation can be obtained:(16)∂2w∂x∂y2−∂2w∂x2∂2w∂y2=∂2εy∂x2+∂2εx∂y2−∂2γxy∂x∂y.According to (2) and (14), the compatibility equation can be simplified; then, the differential equations of elastic plate can be simultaneously combined to obtain the following differential equation system based on the large deflection theory:(17)∇4w=qD+δD∂2ϕ∂x2∂2w∂y2+∂2ϕ∂y2∂2w∂x2−2∂2ϕ∂x∂y∂2w∂x∂y,∇4ϕ=E∂2w∂x∂y2−∂2w∂x2∂2w∂y2.


## 4. Theoretical Calculation Based on the Difference Method

In the process of using small deflection theory or large deflection theory to solve the skin deformation problem, the same boundary of the skin may contain both fixed boundary and free boundary or other types of boundaries due to the panel is partiality supported by the board. Moreover, during the drilling or milling, the distribution of the load and size of the load is irregular; thus, the description of the load is also complex. To sum up, there are mathematical difficulties in solving the differential equation of skin deformation because the mathematical models are often not that accurate to include the different boundary conditions of different skin segments and the randomly distributed load conditions.Difference method is a common method to solve partial differential equations. The basic method is as follows: the domain is discretized by difference grid, the independent variable is represented by the function of the discrete variable at the point of difference grid, and the derivative or partial derivative is represented by the difference quotient of the function between the points of difference grid. The original differential equation and boundary condition are approximately expressed by difference equation in the form of algebraic equation, and then, the differential equation is replaced by the algebraic equation. It is considered that the solution of the algebraic equation is the approximate solution of the original problem.


### 4.1. Difference Calculation Based on Small Deflection Theory

For the convenience of calculation, the skin area is divided by square mesh, and the side length of square mesh is represented byΔd. As shown in Figure 4, the skin is divided into i−1×j−1 mesh regions, including i×j nodes.Figure 4
Differential mesh generation.Using the square grid as shown in Figure4, from column x = 1 to column x = i, write the difference equations for the unknown nodes in each column and list the difference equations. By using the matrix form, we can get the following results:(18)KW=Q,where K is the coefficient matrix of order (i×j)×(i×j) of the difference equations, which can be expressed as(19)K=C1B1A1B2C2D2A2ABCBAABCBA⋱ABCBAABCBAA2′D2′C2′B2′A1′B1′C1′I×J×I×J,where A,B,C,A1,B1,C1,A2,B2,C2,D2,A1′,B1′,C1′,A2′,B2′,C2′,D2′ are the block matrices, all of which are square matrices of order j×j. If the number of rows is j and the number of columns is i in this matrix, A, B, and C are the coefficient matrices corresponding to the internal i-4 column nodes, A2,B2,C2,D2,A1′,B1′,C1′,A2′,B2′,C2′,D2′ are the coefficient matrices corresponding to a column of nodes in the boundary, and A3, B3, C3, A3′, B3′, and C3′ are the coefficient matrices corresponding to the nodes of the boundary column. In the construction, the above matrix needs to be adjusted according to the actual boundary conditions.W is a i×j×1-order column matrix composed of unknown deflection wi,j at each node and is expressed as follows:(20)W=W1W2W3⋮Wii×j×1,Wi=w1,1w1,2w1,3⋮w1,jj×1.Q is the value on the right side of the difference equation and is a known condition in actual calculation, which can be expressed as(21)Q=Q1Q2Q3⋮Qii×j×1,Qi=Qi,1Qi,2Qi,3⋮Qi,jj×1,qi,j=qi,jDΔd4,where qi,j is the load value per unit area at the node (i,j).When the boundary condition is free, the virtual nodes of one line outside the boundary can be represented by two lines of nodes close to the boundary and the virtual nodes of two lines outside the boundary can be represented by three lines of nodes close to the boundary. Therefore, when constructing the coefficient matrices A, B, and C corresponding to the internali-4 column nodes, the j-4 between them acts as the conventional coefficient, which is derived from the original fourth-order difference equation. The first two lines and the second two lines of coefficient matrices A, B, and C are still missing the virtual nodes of the line outside the boundary and the two lines of virtual nodes outside the boundary. Therefore, it is necessary to rewrite the correlation coefficient in the matrix, as shown in t1–t9 and t1′–t9′ in the following formula. Since the coefficients to be rewritten in dealing with simply supported boundary and fixed boundary have been included in t1–t9 and t1′–t9′, the coefficient matrices A, B, and C can be expressed in the following general form:(22)A=t111⋱11t1′,B=t2t3t4−822−82⋱2−822−8t4′t3′t2′,c=t5t6t7t8t9−811−820−81⋱1−820−811−8t9′t8′t7′t6′t5′.After the construction of coefficient matrix K and matrix Q, the numerical solution of skin deflection under different boundary conditions can be obtained by using the MATLAB software, and then, the real solution of the original small deflection bending differential equation can be approximated. We can use the left division command in MATLAB, which has better numerical stability than calculating the inverse matrix [30]. The calculation formula is as follows:(23)W=KQ.


### 4.2. Difference Calculation Based on Large Deflection Theory

According to (17), the difference equations based on large deflection theory are composed of differential equations and compatibility equations of elastic surface. To solve this problem, we need to solve the unknown functions in the equations, deflection function Wx,y and stress function ϕx,y, according to the actual boundary condition of the mechanical model. Since the exact solutions of multiple nonlinear differential equations is difficult to obtain, based on the difference calculation method, the equilibrium differential equations of large deflection theory are transformed into the following difference equations:(24)K1W=Q,K2ϕ=G.There are two unknown matrices in the difference equations, which are deflection matrixW and stress function matrix ϕ. In the original system of differential equations, the unknown matrix ϕ is on the right side of the differential equation of elastic surface. In order to unify the coefficients, it is moved to the left side, that is, the unknown matrix K1 is included in ϕ. G is a matrix containing an unknown matrix W. K1 and K2 are the coefficient matrices to be constructed. The equations cannot be solved directly by the difference method. The approximate solution is obtained by using the idea of gradual approximation. First, the initial value of ϕ is assumed and brought into the differential equation of elastic surface to solve W by the difference method (usually, the initial value of ϕ is 0, which is solved according to the small deflection theory). Then, the solved W is brought into the compatible equation to solve the ϕ value by the difference method, and then, ϕ is brought into the differential equation of elastic surface to compare the solved W with the previous W. It is considered that the solution is completed by iterating in cycles according to this law until the W values of the previous and the following two times approach.The approximation conditions of the twoW values (Wt and Wt+1) are expressed by the Frobenius norm of the difference between the two matrices:(25)Wt−Wt+1F=∑m=1i∑n=1iwm,nt−wm,nt+12≤0.01.The differential equation of elastic surface is solved. According to the difference idea, the difference relation of each basic nodewi,j is represented by the difference diagram as shown in Figure 5.Figure 5
Difference diagram of nodal relation in large deflection theory.It can be seen that the partial differential expression of the stress functionϕx,y is needed when constructing the coefficient matrix of the difference equations. The second-order partial differential of stress function ϕx,y can be expressed by calculating the second-order degree of matrix ϕ twice with gradient function ϕ in MATLAB software. At the same time, it should be noted that in the difference equation of the basic node ωi,j, the partial differential of ϕ in its coefficient must be the second step degree of one-to-one corresponding node ϕi,j in the stress function matrix ϕ. In the coefficient matrix, the coefficients ∂2ϕ/∂x2, ∂2ϕ/∂y2, and ∂2ϕ/∂x∂y used by ωi,j in the corresponding rows of the coefficient matrix are the two-step values of ϕi,j in x, y, and xy directions, respectively. Therefore, the difference between the large deflection theory and the small deflection theory is that when the large deflection theory is used, the values of each line in the coefficient matrix are affected by ϕ and change with ϕ. In the process of solving the differential equation of elastic surface in large deflection theory, the treatment of boundary conditions is the same as that in solving the differential equation of small deflection theory. The coefficient matrix Kj of elastic surface differential equation can be expressed as follows:(26)K1=C1D1H1B2C2BDH2A3B3C3D3H3A4B4C4D4H4⋱Ai−3Bi−3Ci−3Di−3Hi−3Ai−2Bi−2Ci−2Di−2Hi−2Ai−1Bi−1Ci−1Di−1AiBiCii×j×i×j.In the matrix,Ai, Bi, Ci, Di, and Hi are the square matrices of order j, which can be given by the following equation (27).(27)Ai=tia,1b11⋱11tia,1b,Bi=ti,1dbti,1ebti,2cbti,2dti,2eti,3cti,3dti,3e⋱ti,j−2cti,j−2dti,j−2eti,j−1cti,j−1cti,j−1cbti,jcbti,jdb,Ci=ti,1fbti,1gbti,1hbti,2gbti,2fbti,2g11ti,3gti,3fti,3g1⋱1ti,j−2gti,j−2fti,j−2g11ti,j−1gti,j−1fbti,j−1fbti,jhbti,jgbti,jfb,Di=ti,1dbti,1cbti,2ebti,2dti,2cti,3eti,3dti,3c⋱ti,j−2eti,j−2dti,j−2cti,j−1eti,j−1dti,j−1cbti,jebti,jdb,Hi=tib,1b11⋱11tib,1b.Similar to the boundary treatment method based on the theory of small deflection, the corresponding coefficients of the first two and the latter two lines ofAi, Bi , Ci, Di, and Hi should be adjusted according to the change of boundary conditions, and the parameters with b subscripts are used in the matrix. Other parameters are given in Table 1.Table 1
Corresponding parameter values of coefficient matrix of large deflection theory.
ParameterCorresponding valueti,je2+1/2∂Δd2/D∂2ϕi,j/∂x∂yti,jd−8−1/2∂Δd2/D∂2ϕi,j/∂∂y2ti,je2−1/2∂Δd2/D∂2ϕi,j/∂x∂yti,jf20+2∂Δd2/D∂2ϕi,j/∂x2+2∂Δd2/D∂2ϕi,j/∂y2ti,jg−8−∂Δd2/D∂2ϕi,j/∂x2The compatibility equation is solved. The boundary condition of stress functionϕ is treated as follows: the ϕ value of the boundary node is 0, and the ϕ value of the virtual node of a row (or column) outside the boundary is the ϕ value of the symmetric node of a row (or column) inside the boundary. The construction form of coefficient matrix can refer to the construction method of coefficient matrix in the calculation of small deflection theory. The right term of the compatibility equation also deals with the deflection matrix W by using gradient function.


## 4.1. Difference Calculation Based on Small Deflection Theory

For the convenience of calculation, the skin area is divided by square mesh, and the side length of square mesh is represented byΔd. As shown in Figure 4, the skin is divided into i−1×j−1 mesh regions, including i×j nodes.Figure 4
Differential mesh generation.Using the square grid as shown in Figure4, from column x = 1 to column x = i, write the difference equations for the unknown nodes in each column and list the difference equations. By using the matrix form, we can get the following results:(18)KW=Q,where K is the coefficient matrix of order (i×j)×(i×j) of the difference equations, which can be expressed as(19)K=C1B1A1B2C2D2A2ABCBAABCBA⋱ABCBAABCBAA2′D2′C2′B2′A1′B1′C1′I×J×I×J,where A,B,C,A1,B1,C1,A2,B2,C2,D2,A1′,B1′,C1′,A2′,B2′,C2′,D2′ are the block matrices, all of which are square matrices of order j×j. If the number of rows is j and the number of columns is i in this matrix, A, B, and C are the coefficient matrices corresponding to the internal i-4 column nodes, A2,B2,C2,D2,A1′,B1′,C1′,A2′,B2′,C2′,D2′ are the coefficient matrices corresponding to a column of nodes in the boundary, and A3, B3, C3, A3′, B3′, and C3′ are the coefficient matrices corresponding to the nodes of the boundary column. In the construction, the above matrix needs to be adjusted according to the actual boundary conditions.W is a i×j×1-order column matrix composed of unknown deflection wi,j at each node and is expressed as follows:(20)W=W1W2W3⋮Wii×j×1,Wi=w1,1w1,2w1,3⋮w1,jj×1.Q is the value on the right side of the difference equation and is a known condition in actual calculation, which can be expressed as(21)Q=Q1Q2Q3⋮Qii×j×1,Qi=Qi,1Qi,2Qi,3⋮Qi,jj×1,qi,j=qi,jDΔd4,where qi,j is the load value per unit area at the node (i,j).When the boundary condition is free, the virtual nodes of one line outside the boundary can be represented by two lines of nodes close to the boundary and the virtual nodes of two lines outside the boundary can be represented by three lines of nodes close to the boundary. Therefore, when constructing the coefficient matrices A, B, and C corresponding to the internali-4 column nodes, the j-4 between them acts as the conventional coefficient, which is derived from the original fourth-order difference equation. The first two lines and the second two lines of coefficient matrices A, B, and C are still missing the virtual nodes of the line outside the boundary and the two lines of virtual nodes outside the boundary. Therefore, it is necessary to rewrite the correlation coefficient in the matrix, as shown in t1–t9 and t1′–t9′ in the following formula. Since the coefficients to be rewritten in dealing with simply supported boundary and fixed boundary have been included in t1–t9 and t1′–t9′, the coefficient matrices A, B, and C can be expressed in the following general form:(22)A=t111⋱11t1′,B=t2t3t4−822−82⋱2−822−8t4′t3′t2′,c=t5t6t7t8t9−811−820−81⋱1−820−811−8t9′t8′t7′t6′t5′.After the construction of coefficient matrix K and matrix Q, the numerical solution of skin deflection under different boundary conditions can be obtained by using the MATLAB software, and then, the real solution of the original small deflection bending differential equation can be approximated. We can use the left division command in MATLAB, which has better numerical stability than calculating the inverse matrix [30]. The calculation formula is as follows:(23)W=KQ.


## 4.2. Difference Calculation Based on Large Deflection Theory

According to (17), the difference equations based on large deflection theory are composed of differential equations and compatibility equations of elastic surface. To solve this problem, we need to solve the unknown functions in the equations, deflection function Wx,y and stress function ϕx,y, according to the actual boundary condition of the mechanical model. Since the exact solutions of multiple nonlinear differential equations is difficult to obtain, based on the difference calculation method, the equilibrium differential equations of large deflection theory are transformed into the following difference equations:(24)K1W=Q,K2ϕ=G.There are two unknown matrices in the difference equations, which are deflection matrixW and stress function matrix ϕ. In the original system of differential equations, the unknown matrix ϕ is on the right side of the differential equation of elastic surface. In order to unify the coefficients, it is moved to the left side, that is, the unknown matrix K1 is included in ϕ. G is a matrix containing an unknown matrix W. K1 and K2 are the coefficient matrices to be constructed. The equations cannot be solved directly by the difference method. The approximate solution is obtained by using the idea of gradual approximation. First, the initial value of ϕ is assumed and brought into the differential equation of elastic surface to solve W by the difference method (usually, the initial value of ϕ is 0, which is solved according to the small deflection theory). Then, the solved W is brought into the compatible equation to solve the ϕ value by the difference method, and then, ϕ is brought into the differential equation of elastic surface to compare the solved W with the previous W. It is considered that the solution is completed by iterating in cycles according to this law until the W values of the previous and the following two times approach.The approximation conditions of the twoW values (Wt and Wt+1) are expressed by the Frobenius norm of the difference between the two matrices:(25)Wt−Wt+1F=∑m=1i∑n=1iwm,nt−wm,nt+12≤0.01.The differential equation of elastic surface is solved. According to the difference idea, the difference relation of each basic nodewi,j is represented by the difference diagram as shown in Figure 5.Figure 5
Difference diagram of nodal relation in large deflection theory.It can be seen that the partial differential expression of the stress functionϕx,y is needed when constructing the coefficient matrix of the difference equations. The second-order partial differential of stress function ϕx,y can be expressed by calculating the second-order degree of matrix ϕ twice with gradient function ϕ in MATLAB software. At the same time, it should be noted that in the difference equation of the basic node ωi,j, the partial differential of ϕ in its coefficient must be the second step degree of one-to-one corresponding node ϕi,j in the stress function matrix ϕ. In the coefficient matrix, the coefficients ∂2ϕ/∂x2, ∂2ϕ/∂y2, and ∂2ϕ/∂x∂y used by ωi,j in the corresponding rows of the coefficient matrix are the two-step values of ϕi,j in x, y, and xy directions, respectively. Therefore, the difference between the large deflection theory and the small deflection theory is that when the large deflection theory is used, the values of each line in the coefficient matrix are affected by ϕ and change with ϕ. In the process of solving the differential equation of elastic surface in large deflection theory, the treatment of boundary conditions is the same as that in solving the differential equation of small deflection theory. The coefficient matrix Kj of elastic surface differential equation can be expressed as follows:(26)K1=C1D1H1B2C2BDH2A3B3C3D3H3A4B4C4D4H4⋱Ai−3Bi−3Ci−3Di−3Hi−3Ai−2Bi−2Ci−2Di−2Hi−2Ai−1Bi−1Ci−1Di−1AiBiCii×j×i×j.In the matrix,Ai, Bi, Ci, Di, and Hi are the square matrices of order j, which can be given by the following equation (27).(27)Ai=tia,1b11⋱11tia,1b,Bi=ti,1dbti,1ebti,2cbti,2dti,2eti,3cti,3dti,3e⋱ti,j−2cti,j−2dti,j−2eti,j−1cti,j−1cti,j−1cbti,jcbti,jdb,Ci=ti,1fbti,1gbti,1hbti,2gbti,2fbti,2g11ti,3gti,3fti,3g1⋱1ti,j−2gti,j−2fti,j−2g11ti,j−1gti,j−1fbti,j−1fbti,jhbti,jgbti,jfb,Di=ti,1dbti,1cbti,2ebti,2dti,2cti,3eti,3dti,3c⋱ti,j−2eti,j−2dti,j−2cti,j−1eti,j−1dti,j−1cbti,jebti,jdb,Hi=tib,1b11⋱11tib,1b.Similar to the boundary treatment method based on the theory of small deflection, the corresponding coefficients of the first two and the latter two lines ofAi, Bi , Ci, Di, and Hi should be adjusted according to the change of boundary conditions, and the parameters with b subscripts are used in the matrix. Other parameters are given in Table 1.Table 1
Corresponding parameter values of coefficient matrix of large deflection theory.
ParameterCorresponding valueti,je2+1/2∂Δd2/D∂2ϕi,j/∂x∂yti,jd−8−1/2∂Δd2/D∂2ϕi,j/∂∂y2ti,je2−1/2∂Δd2/D∂2ϕi,j/∂x∂yti,jf20+2∂Δd2/D∂2ϕi,j/∂x2+2∂Δd2/D∂2ϕi,j/∂y2ti,jg−8−∂Δd2/D∂2ϕi,j/∂x2The compatibility equation is solved. The boundary condition of stress functionϕ is treated as follows: the ϕ value of the boundary node is 0, and the ϕ value of the virtual node of a row (or column) outside the boundary is the ϕ value of the symmetric node of a row (or column) inside the boundary. The construction form of coefficient matrix can refer to the construction method of coefficient matrix in the calculation of small deflection theory. The right term of the compatibility equation also deals with the deflection matrix W by using gradient function.


## 5. Calculation Error Analysis

In this section, the “4-2-1” support layout of skin is taken as an example to evaluate the accuracy of the elastic mechanics model by using the difference calculation method based on the small deflection theory and the one based on the large deflection theory.Taking the thickness of the skin as 2 mm, two fixtures are used, and each fixture has two sections of surface to fit with the skin. The concentrated force is applied on the middle part of the skin to calculate its deflection. As shown in Figure6, the loads and constraints of the skin are shown, AB and CD are the free boundaries and AD and BC are the free-fixed composite boundaries. In the case of free-fixed composite boundary conditions, it is difficult to solve the analytical solution of the differential equation of thin plate bending deformation by using conventional methods. Therefore, using the difference method to deal with the complex boundary problems can fully reflect the convenience of the difference method to deal with the complex boundary conditions. In order to evaluate the accuracy of the theoretical calculation, ABAQUS finite element simulation is used for comparative calculation. Aluminum alloy 2024-T3 and SC8R shell elements are selected as the skin material in the simulation.Figure 6
Skin loads and constraints.First, the deformations of two kinds of the skin with dimensions of 200 × 200 mm, fixed boundary length of 80 mm, and dimensions of 400 × 400 mm, fixed boundary length of 100 mm under 65 N concentrated force, are compared. The difference method based on the large deflection theory and the difference method based on the small deflection theory are used for theoretical calculation. The difference method is 40 × 40 units. The simulation process of ABAQUS is calculated by linear analysis and nonlinear analysis with large deformation. As shown in Figures7 and 8, the curves of the maximum deflection of two kinds of skin changing with the concentrated force are shown. It can be seen that when the maximum deflection is less than 1 mm (that is, half of the skin thickness), the error of the four calculation results can be controlled within 0.10 mm. With the increase of deflection, the calculation results based on the small deflection theory and ABAQUS linear analysis are approximate, and the maximum error within 2 mm deflection is 5.76%. However, based on the small deflection theory and ABAQUS linear analysis, the maximum error within 2 mm deflection is 28.43%. This is because the strain and internal force caused by the longitudinal displacement of each point in the neutral plane are not taken into account in the calculation based on the small deflection theory. With the increase of deformation, the neutral plane stress and strain caused by longitudinal displacement increase, which makes the calculation based on small deflection theory deviate greatly from the actual situation. The calculated results based on the large deflection theory are similar to the results of ABAQUS nonlinear analysis, and the variation trend is close. The maximum error within 2 mm deflection is 12.70%. If it is considered that the finite element nonlinear analysis results with large deformation are closer to the real situation, the calculation result based on large deflection theory is more accurate in the case of large deformation of the skin.Figure 7
Curve of maximum deflection of 400 × 400 skin changing with concentrated force.Figure 8
Curve of maximum deflection of 200 × 200 skin changing with concentrated force.The deflection of ABAQUS nonlinear analysis results, calculation results based on small deflection theory, and calculation results based on large deflection theory at the red dot position shown in Figures9–11 are compared. The skin with 400 mm size and 100 mm fixed boundary length is selected to analyze the results under the action of 50 N concentrated force. Figure 9 shows the finite element simulation results, Figure 10 shows the theoretical difference calculation results of large deflection, and Figure 11 shows the theoretical difference calculation results of small deflection. Figure 12 shows the comparison of the three calculation results. It can be seen that the deflection change trend of the three calculation methods is consistent at different node positions, which indicates that the difference method based on thin plate theory not only has more accurate calculation results at the maximum deflection but also can accurately reflect the deflection at different positions of the skin. When the deflection exceeds 1 mm, the calculated results based on the large deflection theory are closer to the results of ABAQUS nonlinear analysis, which is the same as the previous conclusion. By comparison, when the actual deformation is large, the calculation results based on the small deflection theory exaggerate the deformation to a certain extent. At this time, the design and calculation according to the small deflection is safe and wasteful, which is consistent with the conclusion mentioned in the literature [12].Figure 9
Nephogram of maximum deflection deformation of finite element skin.Figure 10
Theoretical calculation results of large deflection.Figure 11
Theoretical calculation results of small deflection.Figure 12
Deflection curve of three calculation methods with node position.The calculation results based on the large deflection theory are compared under the different grid densities of 20×20, 30×30, 40×40, and 50×50, and the simulation results are shown in Figure13. In addition to the large deviation of the calculation results when using the difference grid density of 20×20, the deflection of the other three cases keeps a high consistency with the change of load, the maximum deviation is 1.3%, and the deviation is 0.03 mm, which can be ignored. However, with the increase of grid density, the computing time will be doubled. For example, under the same computing conditions, the computing time of 20×20, 30×30, 40×40, and 50×50 differential grid density is 1.5 s, 2.9 s, 8.8 s, and 20.0 s, respectively. This is due to the increase of the unknowns in the difference equations and the capacity of the coefficient matrix with the increase of the grid density. Moreover, under the same conditions, with the increase of deflection, the calculation time will also increase, which is caused by the gradual approximation calculation method of large deflection theory. Under the assumption of the initial value, the initial error of large and small deflection theory is large, so that the number of iterations increases. Therefore, it is necessary to choose the suitable grid density in different cases.Figure 13
Deflection versus load curves under different grid densities.The following conclusions can be obtained by the error evaluation of the difference method above.(1)
The results based on the small deflection theory are similar to the results of the finite element linear analysis, and the results based on the large deflection theory are also similar to the results of the finite element nonlinear analysis. It is proved that the model based on thin plate theory is accurate, and the approximate result of difference calculation is reliable. If it is considered that the results of nonlinear finite element analysis are closer to the real situation, the large deflection theory is more accurate.(2)
When the deflection is less than half of the skin thickness, the calculation results based on the large/small deflection theory are both accurate; When the deflection is more than half of the skin thickness, the calculation results based on the large deflection theory can still approach the accurate value of the simulation, and the calculation results based on the small deflection theory will exaggerate the deformation to a certain extent. At this time, if the design and calculation are carried out according to the small deflection, it is too safe and even wasteful.(3)
When the deflection is near the skin thickness, the calculation results based on the large deflection theory can approach the simulation results, but there is still 12.70% error. The possible error sources include the following points:(i)
The difference method is derived by using the Taylor series expansion of the function and omitted the higher-order terms above the second order, and its error cannot be ignored. Therefore, the more reserved terms and the smaller the mesh size is, the more accurate results can be obtained, but the amount of calculation will increase.(ii)
In the calculation process based on the large deflection theory, the treatment of the boundary conditions of the stress function is relatively rough. In practice, the fixed boundary will be subject to certain longitudinal constraints, and the stress function of the boundary is different from the approximate solution, which leads to errors.(iii)
When using the difference method, the difference of grid density will also affect the calculation results of deflection. Generally speaking, the denser the grid, the higher the calculation accuracy.


## 6. Application Examples

The application verification is carried out by taking the skin corner position of a certain type of aircraft as an example. In Figure14, a and b represent the edge size of the area and c represents the length of the fixture surface. Study the deflection after applying drilling force at the red point F1, as shown in Figure 14.Figure 14
Sketch of skin corner position.Figure15 shows the relationship between the length of the fixture surface and the maximum deflection of the skin when the two adjacent surfaces of the fixture are effective support segments (b = 280 mm) under the action of 65 N drilling force. The curve of a = 420 mm is added as the control group. It can be seen from Figure 15 that the maximum deflection of the skin decreases with the increase of the length of the board surface. Under the condition of a = 420 mm, the maximum deflection is 1.12 mm and the minimum deflection is 0.70 mm in the range of 20–110 mm. Under the condition of a = 840 mm, the maximum deflection is 1.29 mm and the minimum deflection is 0.74 mm in the range of 20–110 mm. When the length of the surface is more than 60 mm, the deflection difference of the skin is basically maintained within 0.05 mm when the distance between the two plates is different. At this time, the influence of the distance between the two plates on the deflection is very small. This shows that the skin deflection under concentrated load is mainly affected by the boundary condition of the surrounding area. When the area is large enough, the influence of the change of the far boundary on the skin deflection can be almost ignored.Figure 15
The relationship between the length of the clamping surface at the corner and the maximum deflection of the skin.Figure16 shows the difference calculation results based on the large deflection theory when a = 840 mm and the support length is 100 mm, and Figure 17 shows the finite element analysis results under the same conditions. Comparing the two figures, it can be seen that the deflection error of the theoretical calculation results and the finite element analysis results at the load point is 0.05 mm, and the error of the maximum deflection at the edge is 0.04 mm, and the overall deflection diffusion trend is the same, indicating that the theoretical calculation results are accurate.Figure 16
Theoretical difference calculation results of large deflection at corner position.Figure 17
Finite element deformation nephogram of corner position.


## 7. Conclusions

In this study, aiming at the problem of fixture design in flexible tooling, the load deformation problem of the skin supported by the fixture is transformed into the thin plate bending problem in elastic mechanics, and the mechanical models are established based on the small deflection theory and the large deflection theory, respectively. Considering the complexity of the boundary conditions of the skin under the support of the fixture, the difference method is used to solve the theoretical model. The coefficients of the difference equation are transformed into the form of coefficient matrix, and the approximate values of the deflection of each node of the skin are obtained by matrix operation. The difference results based on small deflection theory and large deflection theory are compared with the finite element simulation results. Finally, the skin corner position of a certain type of aircraft is taken as an example to verify. The support area and support effect of the fixture on the skin can be calculated and optimized by the numerical calculation method proposed in the manuscript, which guarantees the skin deformation within the allowable range when reducing the effective support section of the fixture to the skin. It has important theoretical guiding significance for the flexible design of fixture.



---

*Source: 6602155-2022-04-25.xml*