WO2018196098A1 - 有限元插值函数构造方法 - Google Patents

有限元插值函数构造方法 Download PDF

Info

Publication number
WO2018196098A1
WO2018196098A1 PCT/CN2017/086830 CN2017086830W WO2018196098A1 WO 2018196098 A1 WO2018196098 A1 WO 2018196098A1 CN 2017086830 W CN2017086830 W CN 2017086830W WO 2018196098 A1 WO2018196098 A1 WO 2018196098A1
Authority
WO
WIPO (PCT)
Prior art keywords
unit
interpolation function
displacement
node
coordinate system
Prior art date
Application number
PCT/CN2017/086830
Other languages
English (en)
French (fr)
Inventor
张国祥
张帆航
Original Assignee
中南大学
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by 中南大学 filed Critical 中南大学
Priority to US15/749,021 priority Critical patent/US20190005171A1/en
Publication of WO2018196098A1 publication Critical patent/WO2018196098A1/zh

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/15—Correlation function computation including computation of convolution operations
    • G06F17/156—Correlation function computation including computation of convolution operations using a domain transform, e.g. Fourier transform, polynomial transform, number theoretic transform
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/17—Function evaluation by approximation methods, e.g. inter- or extrapolation, smoothing, least mean square method
    • G06F17/175—Function evaluation by approximation methods, e.g. inter- or extrapolation, smoothing, least mean square method of multidimensional data
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]

Definitions

  • the present invention relates to the field of analog simulation technologies, and in particular, to a method for constructing a finite element interpolation function.
  • finite element is an indispensable part of engineering analysis and design. Finite element calculation software is now widely used in various fields of structural, solid and fluid analysis engineering. In fact, finite elements are used in almost every field of engineering analysis.
  • a solution domain representing a structure or a continuum is discretized into a plurality of sub-domains (units), and interconnected by nodes on their boundaries to form a combination.
  • This part is the pre-processing part of the finite element software, that is, the unit division part, and the technology of this part is very mature.
  • the approximate field function assumed in each unit is used to slice the unknown field variables to be solved in the full solution domain.
  • the approximation function in each cell is expressed by the value of the unknown field function and its derivative at each node of the cell and its corresponding interpolation function.
  • This part is the finite element interpolation function forming part in the finite element software.
  • the finite element interpolation function is difficult to construct. There are many problems that can not construct the finite element interpolation function that satisfies the basic convergence requirements. It has always been a difficult problem in the field of finite element research.
  • the construction result of the finite element interpolation function of the same problem is not unique.
  • the selection of the finite element interpolation function has great influence on the calculation and analysis precision of the finite element software, which is directly related to the success or failure of the finite element software calculation result.
  • the finite element interpolation function is in the common boundary of adjacent units.
  • the displacement (derivative with displacement) is coordinated, that is, the displacement from the adjacent unit is the same as the displacement of the common boundary. Otherwise, the displacement conflict causes energy loss and reduces the calculation accuracy.
  • the finite element interpolation function is suitable for the surface (curve) boundary.
  • the construction of the existing interpolation function generally adopts the isoparametric coordinate method. Whether it is for a planar solid element, a three-dimensional solid element, a planar thin plate unit or a space shell unit, the calculation accuracy is low, the scope of application is limited, and coordination and/or coordination cannot be achieved. Not suitable for problems such as curve boundaries. E.g:
  • a 4-node quadrilateral unit that has been constructed based on the isoparametric coordinate method. The unit has only one complete coordination and can only be full The basic convergence requirements of the finite element calculation of the foot are low.
  • the 8-node curved quadrilateral element has been constructed.
  • the number of nodes in the unit is doubled, but only one-order complete coordination is available, which can only meet the basic convergence requirements of finite element calculation, and the calculation accuracy is low.
  • the unit is rectangular, the unit can be fully coordinated twice, but it is not suitable for the curve boundary, and the scope of use is very limited.
  • a 12-node curved quadrilateral unit constructed based on the isoparametric coordinate method.
  • the unit interpolation function has only 2 complete coordinations and the calculation accuracy is low.
  • the unit interpolation function can be fully coordinated three times, but it is not suitable for the curve boundary, and the scope of use is very limited.
  • the 20-node surface hexahedral element constructed based on the isoparametric coordinate method. Regardless of whether the unit is an arbitrary hexahedral element or a rectangular parallelepiped, the finite element interpolation function has only two complete coordinations, and the calculation accuracy is limited.
  • a 32-node surface hexahedral element constructed based on the isoparametric coordinate method. Regardless of whether the unit is an arbitrary hexahedral element or a rectangular parallelepiped, the finite element interpolation function has only two complete coordinations, and the calculation accuracy is low.
  • the current structural finite element software is based on a single isoparametric coordinate method (or area coordinate method) construction unit to solve the physical quantity (displacement, temperature, fluid and electromagnetic, etc.) interpolation function, the constructed unit interpolation function can not be high The order is complete and coordinated, even if it is complete, it is only low-order complete, and the calculation accuracy is low. For the structural problem, a high-order complete and coordinated finite element interpolation function that satisfies the basic convergence requirements of finite elements has not been constructed.
  • the object of the present invention is to disclose a finite element interpolation function construction method to improve the performance of an interpolation function.
  • the present invention discloses a finite element interpolation function construction method, including:
  • the number of coordinate elements, the number of items, and the number of times of the interpolation function equation is related to the interpolation of the target entity unit
  • the total number of node displacements is equal; at the same time, the corresponding linear transformation coordinate system in the constructed interpolation function equation is the number of times to cover the highest-order complete polynomial from low to high after the combination of the coordinate elements, and the corresponding isoparametric parameters in the constructed interpolation function equation
  • the number of times of the coordinate system portion is increased from low to high and symmetrically distributed over the highest number of times of the linear transformation coordinate system portion to cover the combination of the coordinate elements.
  • the constructed interpolation function includes any one or any combination of the following:
  • the constructed unit displacement interpolation function is:
  • u(v) a 1 +a 2 T 1 +a 3 T 2 +a 4 T 1 2 +a 5 T 1 T 2 +a 6 T 1 2 +a 7 ⁇ 2 ⁇ +a 8 ⁇ 2 ;
  • the constructed unit displacement interpolation function is:
  • the constructed unit displacement interpolation function is:
  • the constructed unit displacement interpolation function is:
  • the target solid element is a two-dimensional 4-node and each node has three high-order complete coordination arbitrary quadrilateral thin plate elements with relevant displacement components, the relevant displacement components are respectively w, ⁇ x , ⁇ y , and the constructed unit
  • the displacement interpolation function is:
  • the relevant displacement components are w, ⁇ x , ⁇ y , ⁇ xy , respectively.
  • the constructed unit displacement interpolation function is:
  • the constructed unit displacement interpolation function is:
  • T 1 , T 2 and T 3 are the coordinate axes in the linear transformation coordinate system in the unit curved surface, respectively, ⁇ , ⁇ , ⁇ are the coordinate axes in the isoparametric coordinate system; u, v, w respectively correspond to The displacements in the three local coordinate directions in the unit surface, ⁇ x and ⁇ y are the partial derivatives of w to the local coordinates x and y in the unit surface, respectively, and ⁇ xy is the second-order cross-bias derivative of w and x, y.
  • the constructed unit displacement interpolation function is:
  • the constructed unit displacement interpolation function is:
  • T 1 , T 2 and T 3 are the coordinate axes in the linear transformation coordinate system in the unit curved surface, respectively, ⁇ , ⁇ , ⁇ are the coordinate axes in the isoparametric coordinate system; u, v, w respectively correspond to The displacements in the three local coordinate directions in the unit surface, ⁇ x and ⁇ y are the partial derivatives of w to the local coordinates x and y in the unit surface, respectively, and ⁇ xy is the second-order cross-bias derivative of w and x, y.
  • the present invention also supplements the complete rigid body displacement in the displacement mode.
  • the target entity unit is a two-dimensional 4-node high-order complete coordination arbitrary quadrilateral thin plate unit, a two-dimensional 8-point high-order complete coordinated curved quadrilateral thin plate unit, a three-dimensional 4-point high-order complete coordination arbitrary quadrilateral flat thin shell unit, three-dimensional
  • the 8-point high-order fully-coordinated curved quadrilateral flat shell element, the three-dimensional 4-node high-order complete coordinated quadrilateral curved thin shell element or the three-dimensional 8-node high-order complete coordinated quadrilateral curved thin shell element the non-coordination of the cell boundary is modified by the following formula Normal angular displacement:
  • ⁇ ⁇ 23 ( ⁇ ), ⁇ ⁇ 14 ( ⁇ ), ⁇ ⁇ 34 ( ⁇ ), ⁇ ⁇ 12 ( ⁇ ) are the non-coordinated normal angular displacements of the four boundaries of the element, and the non-coordinated normal angular displacement is at the node zero.
  • the transformation formula of the global coordinate system and the linear coordinate system transformation of the present invention may be:
  • the transformation formula of the global coordinate system and the isoparametric transformation of the present invention may be:
  • the finite element is used to solve the basic unknown (the node value of the field function) by the variational principle or the weighted residual method which is equivalent to the original problem mathematical model (basic equation, boundary condition).
  • Algebraic equations or groups of ordinary differential equations can be solved by solving equations.
  • This part is assembled and solved for algebraic equations or ordinary differential equations.
  • the technology of this part is also very mature, and there are standard fixed solution modules. Thereby, the present invention has the following beneficial effects:
  • the calculation accuracy of the finite element analysis software can be greatly improved, the safety and reliability of the structural design can be improved, and the structural design can be optimized. It is more adaptable to various curved (curve) boundaries, which brings huge economic benefits for engineering, aviation and aerospace construction.
  • FIG. 1 is a flow chart of a method for constructing a finite element interpolation function disclosed in a preferred embodiment of the present invention
  • FIG. 2(a) is a schematic diagram of a global coordinate system in a planar linear transformation coordinate system
  • FIG. 2(b) is a schematic diagram of linear transformation coordinates in a planar linear transformation coordinate system
  • FIG. 3(a) is a schematic diagram of a global coordinate system in a spatial linear transformation coordinate system
  • FIG. 3(b) is a schematic diagram of linear transformation coordinates in a spatial linear transformation coordinate system
  • FIG. 4(a) is a schematic diagram of the overall coordinate system of an 8-node curved quadrilateral unit
  • FIG. 4(b) is a schematic diagram of an isoparametric coordinate system of an 8-node curved quadrilateral unit
  • FIG. 4(c) is an 8-node curved quadrilateral unit.
  • Figure 5 (a) is a schematic diagram of the overall coordinate system of a 12-node curved quadrilateral unit
  • Figure 5 (b) is a schematic diagram of an isoparametric coordinate system of a 12-node curved quadrilateral unit
  • Figure 5 (c) is a 12-node curved quadrilateral unit Schematic diagram of a linear transformation coordinate system
  • Figure 6 (a) is a schematic diagram of the overall coordinate system of a 20-node curved quadrilateral unit
  • Figure 6 (b) is a 20-node curved side four sides
  • Fig. 6(c) is a schematic diagram of the linear transformation coordinate system of the 20-node curved quadrilateral unit
  • Fig. 7(a) is a schematic diagram of the overall coordinate system of a 32-node curved quadrilateral unit
  • Fig. 7(b) is a schematic diagram of an isoparametric coordinate system of a 32-node curved quadrilateral unit
  • Fig. 7(c) is a 32-node curved quadrilateral unit Schematic diagram of a linear transformation coordinate system
  • Fig. 8(a) is a schematic diagram showing the overall coordinate system of a 4-node curved quadrilateral thin plate unit
  • Fig. 8(b) is a schematic view of an isoparametric coordinate system of a 4-node curved quadrilateral thin plate unit
  • Fig. 8(c) is a 4-node curved quadrilateral Schematic diagram of a linear transformation coordinate system of a thin plate unit
  • Figure 9 (a) is a schematic diagram of the overall coordinate system of the 8-node curved quadrilateral thin plate unit;
  • Figure 9 (b) is a schematic diagram of the isoparametric coordinate system of the 8-node curved quadrilateral thin plate unit;
  • Figure 9 (c) is an 8-node curved quadrilateral Schematic diagram of a linear transformation coordinate system of a thin plate unit;
  • Figure 10 (a) is a schematic diagram of the overall coordinate system of a spatial 4-node quadrilateral thin shell element;
  • Figure 10 (b) is a schematic diagram of an isoparametric coordinate system of a spatial 4-node quadrilateral thin shell element;
  • Figure 10 (c) is a space 4-node quadrilateral thin Schematic diagram of a linear transformation coordinate system of a shell element;
  • Figure 11 (a) is a schematic diagram of the overall coordinate system of a spatial 8-node quadrilateral thin shell element;
  • Figure 11 (b) is a schematic diagram of an isoparametric coordinate system of a spatial 8-node quadrilateral thin shell element;
  • Figure 11 (c) is a spatial 8-node quadrilateral thin Schematic diagram of a linear transformation coordinate system of a shell element;
  • Figure 12 (a) is a schematic diagram of the overall, partial and curvilinear coordinate system of a spatial 4-node quadrilateral curved thin shell element;
  • Figure 12 (b) is a schematic diagram of an isoparametric coordinate system of a spatial 4-node quadrilateral curved thin shell element;
  • Figure 12 (c) Is a schematic diagram of a linear transformation coordinate system of a space 4-node quadrilateral curved thin shell element;
  • Figure 13 (a) is a schematic diagram of the overall, partial and curvilinear coordinate system of a spatial 8-node quadrilateral curved thin shell element;
  • Figure 13 (b) is a schematic diagram of an isoparametric coordinate system of a spatial 8-node quadrilateral curved thin shell element;
  • Figure 13 (c) It is a schematic diagram of the linear transformation coordinate system of a spatial 8-node quadrilateral curved thin shell element.
  • the finite element interpolation function construction method disclosed in this embodiment includes:
  • Step S1 constructs an interpolation function by a mixed coordinate system composed of a linear transformation coordinate system and an isoparametric coordinate system.
  • the linear coordinate transformation system that is, the transformation relationship between two orthogonal coordinate systems is linear, and the isoparametric transformation system is a nonlinear coordinate transformation system.
  • the orthogonal surface coordinate transformation system existing on the curved thin shell element is equivalent to the global coordinate system of the general structural unit, and can also be transformed into a linear transformation coordinate system.
  • the unit can be changed into a rectangular (line) shape unit, so that the coordinate values of some unit nodes become simple 0 and 1, thereby reducing the difficulty of constructing the finite element interpolation function and improving the calculation.
  • the purpose of accuracy is to make the complete order of the polynomial not increase by coordinate transformation.
  • the coordinate system of the finite element interpolation function polynomial is several complete polynomials. There are several complete polynomials in the whole coordinate system, which makes the structure high. A well-ordered finite element interpolation function is possible.
  • the linear transformation coordinate system is divided into a plane linear transformation coordinate system and a spatial linear transformation coordinate system.
  • Fig. 2(a) and Fig. 2(b) The shape of the unit after coordinate transformation is shown in Fig. 2(a) and Fig. 2(b).
  • One corner of the quadrilateral element is at the origin of the coordinate, and the two corner points are located on the coordinate axis.
  • There are six undetermined coefficients o i , b i , c i , (i 1, 2) in the coordinate transformation relationship, and six linear transformation coordinate values of 0 or 1 can be transformed.
  • the shape of the unit after the coordinate change is shown in Fig. 3(a) and Fig. 3(b).
  • One corner of the hexahedral element is at the origin of the coordinate, and the three corner points are located on the coordinate axis.
  • Step S2 determining a coordinate element number, an item number, and a number of times of the interpolation function equation according to characteristics of the target entity unit, the feature including the number of known nodes and the number of node displacement components, and the number of items of the constructed interpolation function equation and the target entity unit
  • the total number of displacements of the interpolated correlation nodes is equal; at the same time, the corresponding linear transformation coordinate system in the constructed interpolation function equation is the highest order sub-polynomial with the number of times to cover the combination of each coordinate element from low to high, corresponding to the constructed interpolation function equation
  • the number of times of the isoparametric coordinate system portion is increased from low to high and symmetrically distributed over the highest number of times of the linear transformation coordinate system portion to cover the combination of the coordinate elements.
  • the present embodiment constructs a finite element interpolation function based on a linear transformation coordinate and an isoparametric coordinate coordinate system.
  • the complete term of the finite element interpolation polynomial uses linear transformation coordinates, while the extra terms use isoparametric coordinates.
  • the linear transformation coordinates are represented by isoparametric coordinates, and then substituted into the finite element interpolation polynomial.
  • the number of polynomials represented by the obtained isoparametric coordinates is not higher than the complete order of the finite element interpolation polynomial, thus ensuring that the finite element interpolation function has C 0 at the element boundary.
  • the constructed element displacement interpolation function is :
  • u(v) a 1 +a 2 T 1 +a 3 T 2 +a 4 T 1 2 +a 5 T 1 T 2 +a 6 T 1 2 +a 7 ⁇ 2 ⁇ +a 8 ⁇ 2 .
  • T 1 , T 2 and T 3 are the coordinate axes in the linear transformation coordinate system of the unit curved surface, respectively, ⁇ , ⁇ , ⁇ are the coordinate axes in the isoparametric coordinate system; u, v, w respectively correspond to the unit The displacements in the three local coordinate directions in the surface, ⁇ x and ⁇ y are the partial derivatives of the local coordinates x and y in the unit curved surface, respectively, and ⁇ xy is the second-order cross-bias derivative of w and x, y;
  • the total number of interpolation-related node displacements corresponding to u(v), u(v, w) is equal to the number of known nodes, and the total number of interpolation-related node displacements corresponding to w is the product of the number of known nodes and the relevant displacement component, and subsequent Let me repeat.
  • the unit displacement interpolation function can establish a system of equations according to the displacement of the unit nodes, and determine the undetermined coefficients a 1 to a 8 in a simultaneous solution, which has a 2nd order completeness, and the 8-node quadrilateral unit can only have 2 order completeness, which is more than the traditional unit displacement.
  • the completeness of the interpolation function is 1 step higher.
  • the completeness of the element displacement interpolation function is improved once, and its convergence performance and anti-distortion performance will be greatly improved.
  • the common boundary of two adjacent elements is a straight line and the nodes in the boundary are equally divided, the common boundary is coordinated, regardless of the other boundary shapes of the unit. Therefore, when dividing the unit, the common boundary of the two adjacent units is kept as a straight line. And the nodes in the boundary are divided into boundaries, and the free outer boundary of the unit has no coordination requirement, which can be a curve. At this time, the two-dimensional 8-node element displacement interpolation function has high-order complete coordination and is suitable for the curve boundary, and does not increase. The difficulty of unit division.
  • This two-dimensional 8-node quadrilateral curved edge unit can degenerate a 6-node triangular curved edge unit.
  • the constructed unit displacement interpolation function is :
  • the unit displacement interpolation function can establish a system of equations according to the displacement of the unit nodes, and the simultaneous solution determines the undetermined coefficients a 1 to a 12 , which have 3 order completeness, which is 2 orders higher than the completeness of the traditional unit displacement interpolation function, and its convergence performance and The anti-distortion performance is very good.
  • the common boundary of two adjacent elements is a straight line and the nodes in the boundary are equally divided, the common boundary is coordinated, regardless of the other boundary shapes of the unit. Therefore, when dividing the unit, the common boundary of the two adjacent units is kept as a straight line. And the nodes in the boundary are divided into boundaries, and the free boundary of the unit has no coordination requirement, which can be a curve.
  • the two-dimensional 12-node unit displacement interpolation function has high-order complete coordination and can be adapted to arbitrary curve boundaries, and there is no The difficulty of dividing the unit.
  • This two-dimensional 12-node quadrilateral curved edge unit can degenerate a 9-node triangular curved edge unit.
  • the unit displacement interpolation function can establish a system of equations according to the displacement of the unit nodes, and the simultaneous solution determines the undetermined coefficients a 1 to a 20 , which have 2 order completeness, which is 1 order higher than the completeness of the traditional unit displacement interpolation function.
  • the curved hexahedral element can be automatically degenerated into a curved pentahedron unit and a curved side tetrahedral unit; the six quadrilateral side surfaces can be automatically degenerated into a triangular side surface.
  • the common angular line of the unit is required to be a straight line and the middle node is divided into its common angular line.
  • the unit displacement interpolation function is coordinated on the common edge line, regardless of the shape of other parts of the unit.
  • the common quadrilateral side surface of the unit is a plane, the four sides are straight lines, and the nodes in the side are divided into their boundaries.
  • the unit displacement interpolation function is coordinated on the common quadrilateral side plane, regardless of the shape of other planes of the unit.
  • the triangular side surface common to the unit is a plane
  • the side of the triangle can be a curve
  • the unit displacement interpolation function is coordinated on the side plane of the common triangle, regardless of the shape of other faces of the unit.
  • the sides of the triangle can be curved conditions, making it possible to construct degenerate tetrahedral elements suitable for any curved boundary.
  • the plane hexahedral element is used inside the structure, and the outer surface boundary of the degenerate tetrahedral unit, the pentahedral unit and the hexahedral unit exposed surface can be used to ensure the coordination of the unit displacement interpolation function.
  • the constructed unit displacement interpolation function is :
  • the unit displacement interpolation function can establish a system of equations according to the displacement of the unit nodes, and the simultaneous solution determines the undetermined coefficients a 1 to a 32 , which have 3 order completeness, which is 2 orders higher than the completeness of the traditional unit displacement interpolation function.
  • the curved hexahedral element can also be automatically degenerated into a curved pentahedral unit and a curved side tetrahedral unit; the six quadrilateral side surfaces can be automatically degenerated into triangular side surfaces to fit the curved structural boundary.
  • the common angular line of the unit is required to be a straight line and the middle node is divided into its common angular line.
  • the unit displacement interpolation function is coordinated on its common angular line, regardless of the shape of other parts of the unit.
  • the common quadrilateral side surface of the unit is a plane, the four sides are straight lines, and the nodes in the side are divided into their sides.
  • the unit displacement interpolation function is coordinated on the common quadrilateral side plane, regardless of the shape of other planes of the unit.
  • the triangular side surface common to the unit is a plane
  • the side of the triangle can be a curve
  • the unit displacement interpolation function is coordinated on the side plane of the common triangle, regardless of the shape of other faces of the unit.
  • the sides of the triangle can be curved conditions, making it possible to construct degenerate tetrahedral elements suitable for any curved boundary.
  • the plane hexahedral element is used inside the structure, and the outer surface boundary of the degenerate tetrahedral unit, the pentahedral unit and the hexahedral unit exposed surface can be used to ensure the coordination of the unit displacement interpolation function.
  • the target entity unit is a two-dimensional 4-node and each node has three correlation displacement components (the three-node parameter is w, ⁇
  • the displacement displacement function of the constructed unit when x , ⁇ y , w, ⁇ x , ⁇ y are related, in this case, the total number of displacements is the product of the number of known nodes and the number of node displacement components, the high-order complete coordination of arbitrary quadrilateral thin plate elements for:
  • the unit displacement interpolation function can establish a system of equations according to the displacement of the unit nodes, and the simultaneous solution determines the undetermined coefficients a 1 to a 12 , which have a 3 order completeness, and the complete order is higher, and the unit displacement interpolation function is in the adjacent unit.
  • the deflection and tangential corners on the common boundary can be coordinated.
  • the conventional quadrilateral thin plate unit has only a rectangular unit and a triangular unit, and is not suitable for any polygonal line boundary, and the application range is very limited.
  • the unit displacement interpolation function coordinates the deflection and tangential rotation angle on the element boundary, but the normal rotation angle is not coordinated, that is, the C 1 order non-coordination problem.
  • This problem has been considered as an unsolvable problem. Therefore, in the regular isoparametric coordinate system, the present embodiment uses only a very special correction function to correct only the non-coordinating normal rotation angle of one boundary of the unit, and does not affect the displacement and rotation angle of other boundaries of the unit. It does not affect the normal corner value of other boundaries of the unit. Therefore, only the unit has a non-coordinated normal corner common boundary, and the non-coordinated normal corner of the unit free boundary does not have to be corrected.
  • the present embodiment can correct the non-coordinated normal angular displacement of the cell boundary as follows:
  • ⁇ ⁇ 23 ( ⁇ ), ⁇ ⁇ 14 ( ⁇ ), ⁇ ⁇ 34 ( ⁇ ), ⁇ ⁇ 12 ( ⁇ ) are the non-coordinated normal angular displacements of the four boundaries of the element, and the non-coordinated normal angular displacement is at the node zero.
  • the non-coordinated normal angular displacement of any boundary of the unit can be selectively corrected.
  • the target entity unit is a two-dimensional 4-node and each node has four correlation displacement components (four-node parameters are w, ⁇ x , ⁇ y , ⁇ xy , w, ⁇ x , ⁇ y , ⁇ xy are related.
  • the total number of displacements is the product of the number of known nodes and the number of node displacement components, where ⁇ xy is w versus x, y
  • the constructed element displacement interpolation function is:
  • the unit displacement interpolation function can establish a system of equations according to the displacement of the unit nodes, and determine the undetermined coefficients a 1 to a 16 by the simultaneous solution. Considering the four-node parameters w, ⁇ x , ⁇ y , ⁇ xy , the interpolation function of the unit displacement can be guaranteed.
  • the normal corners on the cell boundary are coordinated without special modification; and the 3rd order completeness is obtained, and the complete order is higher.
  • the displacement interpolation function of the element is not only the deflection and the tangential rotation angle on the common boundary of the adjacent elements.
  • the normal corners are also coordinated, but the thickness of the plate unit and the material properties are required to be the same. Otherwise, new inconsistencies will be caused.
  • the constructed element displacement interpolation function is:
  • the unit displacement interpolation function can establish a system of equations according to the displacement of the unit nodes, and the simultaneous solution determines the undetermined coefficients a 1 to a 24 , which have 5 order completeness, and the complete order is high, which is suitable for any curved boundary.
  • the unit displacement interpolation function can be a curve to fit the curve boundary.
  • the unit displacement interpolation function The deflection and the tangential corner are coordinated on the common boundary of the unit, but the normal corner cannot be coordinated and must be corrected.
  • the outer boundary of the element can be a curve.
  • the present embodiment corrects the non-coordinated normal angular displacement of the boundary of the unit:
  • ⁇ ⁇ 23 ( ⁇ ), ⁇ ⁇ 14 ( ⁇ ), ⁇ ⁇ 34 ( ⁇ ), ⁇ ⁇ 12 ( ⁇ ) is the non-coordinated normal angular displacement of the cell boundary, and the non-coordinated normal angular displacement is zero at the node.
  • the non-coordinated normal angular displacement of any boundary of the unit can be selectively corrected.
  • the target solid element is a three-dimensional 4-node high-order complete coordination arbitrary quadrilateral flat thin shell element (where w, ⁇ x , ⁇ y is related), that is, for a 4-node planar arbitrary quadrilateral thin shell unit, as shown in Fig. 10 (a ), Fig. 10(b) and Fig. 10(c), by means of coordinate transformation, the displacement interpolation function of any quadrilateral planar thin plate element can be converted into a spatial thin shell element displacement interpolation function, wherein the node displacement vector is between the two coordinate systems.
  • the conversion relationship is:
  • the key is the construction of the displacement interpolation function of the arbitrary quadrilateral thin shell element. Based on the linear transformation coordinate method and the isoparametric coordinate mixing method of the present invention, the unit displacement interpolation function can be assumed as:
  • the interpolation function can establish a system of equations according to the displacement of the unit nodes, and determine the undetermined coefficients by simultaneous solution; have 3 order completeness, complete order, and the deflection of the thin shell element displacement interpolation function on the common boundary of adjacent units It can be coordinated with the tangential corner.
  • the displacement and tangential rotation angles of the thin-shell element displacement interpolation function are coordinated at the element boundary, but the normal rotation angle is not coordinated, that is, the C 1 order non-coordination problem. Therefore, in the regular isoparametric coordinate system of the present embodiment, only a non-coordinating normal rotation angle of a boundary of the unit is corrected by a very special correction function, which does not affect the deflection of the unit boundary and the tangential rotation angle value. It is a normal corner value that does not affect the other boundaries of the unit. Therefore, only the non-coordinating normal corners of the unit common boundary can be corrected, and the non-coordinated normal corner of the unit free boundary does not have to be corrected.
  • the present embodiment corrects the non-coordinated normal angular displacement of the boundary of the thin shell element by the following formula:
  • ⁇ ⁇ 23 ( ⁇ ), ⁇ ⁇ 14 ( ⁇ ), ⁇ ⁇ 34 ( ⁇ ), ⁇ ⁇ 12 ( ⁇ ) is the non-coordinated normal angular displacement of the cell boundary, and the non-coordinated normal angular displacement is zero at the node.
  • the non-coordinated normal angular displacement of any boundary of the unit can be selectively corrected.
  • the target solid element is a three-dimensional 8-node high-order complete coordination curve quadrilateral flat thin shell element (where w, ⁇ x , ⁇ y is related), that is, for an 8-node planar arbitrary quadrilateral thin shell unit, as shown in Fig. 11 (a As shown in Fig. 11(b) and Fig. 11(c), the displacement interpolation function of any quadrilateral planar thin plate element can be converted into a spatial thin shell element displacement interpolation function by coordinate transformation.
  • the node displacement vector is transformed between two coordinate systems:
  • the key is the construction of the displacement interpolation function of the arbitrary quadrilateral thin shell element. Based on the linear transformation coordinate method and the isoparametric coordinate mixing method of the present invention, the unit displacement interpolation function can be assumed as:
  • the above interpolation function can establish a system of equations according to the displacement of the unit nodes, and determine the undetermined coefficients by a simultaneous solution; it has a 5th order completeness, and the complete order is high, which is suitable for the curved boundary.
  • the unit displacement interpolation function The deflection and the tangential corner are coordinated on the common boundary of the unit, but the normal corner cannot be coordinated and must be corrected.
  • the deflection and tangential rotation angle of the element displacement interpolation function are coordinated at the element boundary, but the normal rotation angle is uncoordinated, that is, the C 1 order non-coordination problem. Therefore, in the regular isoparametric coordinate system of the present embodiment, only a non-coordinating normal rotation angle of a boundary of the unit is corrected by a very special correction function, which does not affect the deflection of the unit boundary and the tangential rotation angle value. It is a normal corner value that does not affect the other boundaries of the unit. Therefore, only the non-coordinating normal corners of the unit common boundary can be corrected, and the non-coordinated normal corner of the unit free boundary does not have to be corrected.
  • the present embodiment corrects the non-coordinated normal angular displacement of the boundary of the unit:
  • ⁇ ⁇ 23 ( ⁇ ), ⁇ ⁇ 14 ( ⁇ ), ⁇ ⁇ 34 ( ⁇ ), ⁇ ⁇ 12 ( ⁇ ) is the non-coordinated normal angular displacement of the element boundary, and the non-coordinated normal angular displacement is zero at the node.
  • the non-coordinated normal angular displacement of any boundary of the unit can be selectively corrected.
  • the target solid element is a three-dimensional 4-node high-order fully coordinated quadrilateral curved thin shell element (where w, ⁇ x , ⁇ y are related), as shown in Fig. 12(a), Fig. 12(b) and Fig. 12(c) ), it can be assumed that the overall coordinates of any point inside the curved thin shell element are:
  • N i ( ⁇ , ⁇ ) is a conventional shape function.
  • the rigid body displacement ⁇ ' Ri ⁇ of the node in the global coordinate system is given by the rigid body motion of the unit microbody, and the motion of the microcentre centroid includes rotation about three coordinate axes and translation along the coordinate axis.
  • the rigid body motion of the curved shell element is:
  • the rigid body displacement of the rigid body motion obtained by the dynamic method is:
  • x 0 , y 0 , z 0 are the global coordinates of the centroid of the unit micro-body.
  • the rigid body displacement of the joint generated by the rigid body motion is:
  • L i is a transformation matrix
  • the rigid body motion of the centroid of the curved thin shell element is:
  • the unit node displacement vector converts the relationship between the global coordinates and the curve coordinates as:
  • the total displacement field of the unit after the displacement of the rigid body is:
  • I is a unit matrix of 20 x 20.
  • the rigid body displacement does not produce strain, so the strain matrix in the orthogonal main curve coordinate system is:
  • B is the strain matrix of the elements in the orthogonal main curve coordinate system.
  • the rigid body displacement does not generate nodal forces, for static condensation, and the finite element equation is established according to the principle of virtual work:
  • the node displacement solution equation of the global coordinate system can be obtained.
  • ⁇ , ⁇ , ⁇ are the arc lengths of the curve coordinates.
  • ⁇ , ⁇ , ⁇ are the arc lengths of the curve coordinates.
  • the coordinate mixing method is similar to the planar thin plate problem. Linear coordinate transformation and isoparametric coordinate transformation are performed on the arc length coordinates of the curve. It can be assumed that the element displacement interpolation function is:
  • the equations are established according to the displacement of the unit nodes, and the simultaneous solution determines the undetermined coefficients.
  • the unit displacement interpolation function has 3 order completeness, and the complete order is higher.
  • the unit displacement interpolation function can coordinate the deflection and tangential rotation angle on the common boundary of adjacent units.
  • the unit displacement interpolation function coordinates the deflection and tangential rotation angle on the element boundary, but the normal rotation angle is not coordinated, that is, the C 1 order non-coordination problem. Therefore, in the regular isoparametric coordinate system, this embodiment uses a very special correction function to correct only the non-coordinating normal rotation angle of one boundary of the unit, and does not affect the displacement and rotation angle of other boundaries of the unit. It does not affect the normal corner value of other boundaries of the unit. Therefore, only the unit has a non-coordinated normal corner common boundary, and the non-coordinated normal corner of the unit free boundary does not have to be corrected.
  • the present embodiment corrects the non-coordinated normal angular displacement of the boundary of the unit:
  • ⁇ ⁇ 23 ( ⁇ ), ⁇ ⁇ 14 ( ⁇ ), ⁇ ⁇ 34 ( ⁇ ), ⁇ ⁇ 12 ( ⁇ ) are the non-coordinated normal angular displacements of the four boundaries of the element, and the non-coordinated normal angular displacement is at the node zero.
  • the non-coordinated normal angular displacement of any boundary of the unit can be selectively corrected.
  • the target solid element is a three-dimensional 8-node high-order fully coordinated quadrilateral curved thin-shell element (where w, ⁇ x , ⁇ y are related), that is, for an 8-node curved thin-shell element, as shown in Fig. 13(a) and 13(b) and Figure 13(c), assume that the overall coordinates of any point inside the curved shell element are:
  • N i ( ⁇ , ⁇ ) is a conventional shape function.
  • the rigid body displacement ⁇ ' Ri ⁇ of the node in the global coordinate system is given by the rigid body motion of the unit microbody, and the motion of the microcentre centroid includes rotation about three coordinate axes and translation along the coordinate axis.
  • the rigid body motion of the shell element is:
  • the rigid body displacement of the rigid body motion obtained by the dynamic method is:
  • x 0 , y 0 , z 0 are the global coordinates of the centroid of the unit micro-body.
  • the rigid body displacement of the joint generated by the rigid body motion is:
  • L i is a transformation matrix
  • the rigid body motion of the centroid of the curved thin shell element is:
  • the unit node displacement vector converts the relationship between the global coordinates and the curve coordinates as:
  • the total displacement field of the unit after the displacement of the rigid body is:
  • I is a unit matrix of 40 x 40.
  • the rigid body displacement does not produce strain, so the strain matrix in the orthogonal main curve coordinate system is:
  • B is the strain matrix of the elements in the orthogonal main curve coordinate system.
  • the rigid body displacement does not generate nodal forces, for static condensation, and the finite element equation is established according to the principle of virtual work:
  • the node displacement solution equation of the global coordinate system can be obtained.
  • ⁇ , ⁇ , ⁇ are the arc lengths of the curve coordinates, for an 8-node arbitrary quadrilateral single Element, based on the original linear transformation coordinate method and the isoparametric coordinate method, similar to the planar thin plate problem, linear coordinate transformation and isoparametric coordinate transformation on the arc length coordinate of the curve, the assumption that the element displacement interpolation function is:
  • the equations can be established according to the displacement of the unit nodes, and the undetermined coefficients are determined by the simultaneous solution.
  • the unit displacement interpolation function has 5 order completeness, and the complete order is high, which is suitable for any curved boundary.
  • the unit displacement interpolation function The deflection and the tangential corner can be coordinated on the common boundary of the unit, but the normal corner cannot be coordinated and must be corrected.
  • ⁇ ⁇ 23 ( ⁇ ), ⁇ ⁇ 14 ( ⁇ ), ⁇ ⁇ 34 ( ⁇ ), ⁇ ⁇ 12 ( ⁇ ) is the non-coordinated normal angular displacement of the cell boundary, and the non-coordinated normal angular displacement is zero at the node.
  • the non-coordinated normal angular displacement of any boundary of the unit can be selectively corrected.
  • the present embodiment is based on the linear transformation coordinate and the isoparametric coordinate hybrid coordinate method, instead of constructing the interpolation function of the physical quantity based on a single coordinate system, so that the calculation precision of the finite element analysis software can be greatly improved, and the structural design is safe and reliable.

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Algebra (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • Operations Research (AREA)
  • Computing Systems (AREA)
  • Complex Calculations (AREA)
  • Pharmaceuticals Containing Other Organic And Inorganic Compounds (AREA)

Abstract

一种有限元插值函数构造方法，涉及模拟仿真技术领域，方法包括：以线性变换坐标系与等参坐标系所组成的混合坐标系构造插值函数（S1）；根据目标实体单元的特征确定插值函数方程式的坐标元数、项数和次数，所述特征包括已知节点数和节点位移分量数，且所构造插值函数方程式的项数与目标实体单元的插值相关节点位移总数相等；同时，所构造插值函数方程式中对应线性变换坐标系部分为各项次数以覆盖各坐标元组合后从低到高递增的最高次完整多项式，所构造插值函数方程式中对应等参坐标系部分的各项次数在所述线性变换坐标系部分的最高次数上以覆盖各坐标元组合后从低到高递增并呈对称性分布（S2）。该方法提高了插值函数的性能。

Description

有限元插值函数构造方法 技术领域
本发明涉及模拟仿真技术领域，尤其涉及一种有限元插值函数构造方法。
背景技术
目前有限元是工程分析和设计必不可少的重要组成部分，有限元计算软件现已广泛应用于结构、固体和流体分析工程的各个领域。事实上，有限元几乎在工程分析的每个领域都得到应用。
在工程或物理问题的数学模型(基本变量、基本方程、求解域和边界条件等)确定以后，有限元法作为对其进行分析的数值计算方法可归纳如下三部分：
(1)将一个表示结构或连续体的求解域离散为若干个子域(单元)，并通过它们边界上的结点相互联结成为组合体。该部分为有限元软件的前处理部分，即单元划分部分，该部分的技术已很成熟。
(2)用每个单元内所假设的近似函数来分片地表示全求解域内待求的未知场变量。而每个单元内的近似函数由未知场函数及其导数在单元各个结点上的数值和与其对应的插值函数来表达。该部分即为有限元软件中的有限元插值函数形成部分，有限元插值函数构造难度大，有很多问题还不能构造出满足基本收敛要求的有限元插值函数，一直是有限元研究领域的难题。
同一问题的有限元插值函数构造结果不是唯一，有限元插值函数的选用对有限元软件的计算分析精度影响很大，直接关系到有限元软件计算结果的成败。构造高精度有限元插值函数有三个关键性条件：一是有限元插值函数所用多项式的完备(完整)阶数(次数)越高计算精度越高；二是有限元插值函数在相邻单元公共边界上的位移(含位移的导数)协调，即从相邻单元插值同一公共边界的位移要一致，否则，位移冲突引起能量损失，降低计算精度；三是有限元插值函数要适合曲面(曲线)边界。目前还无法构造能同时满足以上条件的有限元插值函数，构造高阶完备协调的有限元插值函数非常困难。
现有的插值函数的构造一般采用等参坐标法，无论是对于平面实体单元、三维实体单元、平面薄板单元或空间壳单元，都存在计算精度低，适用范围有限且不能兼顾协调性和/或不适应曲线边界等问题。例如：
1)、基于等参坐标法已构造出的4节点四边形单元。该单元只有1次完备协调，只能满 足有限元计算的基本收敛要求，计算精度低。
2)、基于等参坐标法构已造出的8节点曲边四边形单元。该单元节点数增加一倍，但还只具有1阶完备协调，只能满足有限元计算的基本收敛要求，计算精度低。当该单元为矩形时，单元能2次完备协调，但又不适合曲线边界，使用范围非常有限。
3)、基于等参坐标法已构造出的12节点曲边四边形单元。该单元插值函数只有2次完备协调，计算精度低。当该单元为矩形时，单元插值函数能3次完备协调，但不适合曲线边界，使用范围非常有限。
4)、基于等参坐标法构造出的8节点任意六面体单元，当该单元为任意六面体单元时，适合于折线形边界，但单元插值函数只有1阶完备性，只能满足有限元计算的基本收敛要求，计算精度低。当该单元为长方体时，单元插值函数能2次完备协调，但又不适合于折线形边界，使用范围非常有限。
5)、基于等参坐标法构造出的20节点曲面六面体单元，不管该单元为任意六面体单元还是长方体，有限单元插值函数都只有2次完备协调，计算精度提高有限。
6)、基于等参坐标法构造出的32节点曲面六面体单元，不管该单元为任意六面体单元还是长方体，有限单元插值函数都只有2次完备协调，计算精度低。
7)、基于等参坐标法只构造出的4节点(w相关的三节点参数w，θx，θy和四节点参数w，θx，θy，θxy)三次完备矩形薄板单元位移插值函数，尽管单元位移插值函数的完备阶数较高，但单元边界法向转角位移不协调，也不适合任意折线边界，适用范围非常有限。
8)、基于等参坐标法，还不能构造出二次完备的4节点任意四边形薄板单元，更无法解决单元协调问题。
9)、基于等参坐标法，目前只能将平面4节点非协调矩形和三角形薄板单元位移插值函数，经坐标变换法用于空间薄壳结构，适用范围非常有限，且不协调。
10)、基于等参坐标法构已造出的三维8节点协调低阶完备四边形超参数曲面壳单元位移插值函数和三维8节点低阶完备协调曲面四边形相对自由度壳单元位移插值函数。这两种单元位移插值函数协调，厚薄壳结构通用，但只具有1阶完备性，计算精度低，当壳的厚度趋向薄曲壳时，存在剪切“锁死”和薄膜“锁死”等问题。
综上，目前结构有限元软件都是基于单一的等参坐标法(或面积坐标法)构造单元求解物理量(位移、温度、流体和电磁等)的插值函数，所构造的单元插值函数不能既高阶完备又协调，即使完备也只是低阶完备，计算精度较低。对于结构问题还没有构造出满足有限元基本收敛要求的高阶完备又协调有限单元插值函数。
发明内容
本发明目的在于公开一种有限元插值函数构造方法，以提高插值函数的性能。
为实现上述目的，本发明公开了一种有限元插值函数构造方法，包括：
以线性变换坐标系与等参坐标系所组成的混合坐标系构造插值函数；
根据目标实体单元的特征确定插值函数方程式的坐标元数、项数和次数，所述特征包括已知节点数和节点位移分量数，且所构造插值函数方程式的项数与目标实体单元的插值相关节点位移总数相等；同时，所构造插值函数方程式中对应线性变换坐标系部分为各项次数以覆盖各坐标元组合后从低到高递增的最高次完整多项式，所构造插值函数方程式中对应等参坐标系部分的各项次数在所述线性变换坐标系部分的最高次数上以覆盖各坐标元组合后从低到高递增并呈对称性分布。
本发明中，可选的，所构造的插值函数包括以下的任意一项或任意组合：
1)、当所述目标实体单元为二维8节点高阶完备协调四边形曲边单元时，所构造单元位移插值函数为：
u(v)＝a1+a2T1+a3T2+a4T1 2+a5T1T2+a6T1 2+a7ξ2η+a8ξη2；
2)、当所述目标实体单元为二维12节点高阶完备协调四边形曲边单元时，所构造单元位移插值函数为：
Figure PCTCN2017086830-appb-000001
3)、当所述目标实体单元为三维20节点高阶完备协调曲面六面体单元时，所构造单元位移插值函数为：
Figure PCTCN2017086830-appb-000002
4)、当所述目标实体单元为三维32节点高阶完备协调曲面六面体单元时，所构造单元位移插值函数为：
Figure PCTCN2017086830-appb-000003
5)、当所述目标实体单元为二维4节点且各节点有3个相关位移分量的高阶完备协调任意四边形薄板单元时，相关位移分量分别为w、θx、θy，所构造单元位移插值函数为：
Figure PCTCN2017086830-appb-000004
6)、当所述目标实体单元为二维4节点且各节点有4个相关位移分量的高阶完备协调任意四边形薄板单元时，相关位移分量分别为w、θx、θy、θxy，所构造单元位移插值函数为：
Figure PCTCN2017086830-appb-000005
7)、当所述目标实体单元为二维8节点且各节点有3个相关位移分量(w，θx，θy相关)的高阶完备协调曲边四边形薄板单元时，所构造单元位移插值函数为：
Figure PCTCN2017086830-appb-000006
8)、当所述目标实体单元为三维4节点高阶完备协调任意四边形平板薄壳单元(其中w，θx，θy相关)时，所构造单元位移插值函数为：
Figure PCTCN2017086830-appb-000007
9)、当所述目标实体单元为三维8节点高阶完备协调曲线四边形平板薄壳单元(其中w，θx，θy相关)时，所构造单元位移插值函数为：
Figure PCTCN2017086830-appb-000008
其中，上述各方程式中，T1、T2、T3分别为单元曲面内线性变换坐标系中坐标轴，ζ、η、ζ分别为等参坐标系中坐标轴；u、v、w分别对应单元曲面内三个局部坐标方向上的位移，θx、θy分别为w对单元曲面内局部坐标x、y的偏导数，θxy为w对x、y二阶交叉偏导数。
本发明中，对于工程上普遍采用的空间薄壳，采用合适的正交曲线坐标及相应的几何方程，根据上述单元的原理，像平面问题一样在空间正交曲线坐标系中直接构造高阶完备协调曲面薄壳单元，计算单元刚度矩阵，再进行空间坐标转换；具体包括：
10)、当所述目标实体单元为三维4节点高阶完备协调四边形曲面薄壳单元(其中w，θx，θy 相关)时，所构造单元位移插值函数为：
Figure PCTCN2017086830-appb-000009
11)、当所述目标实体单元为三维8节点高阶完备协调四边形曲面薄壳单元(其中w，θx，θy相关)时，所构造单元位移插值函数为：
Figure PCTCN2017086830-appb-000010
其中，上述各方程式中，T1、T2、T3分别为单元曲面内线性变换坐标系中坐标轴，ξ、η、ζ分别为等参坐标系中坐标轴；u、v、w分别对应单元曲面内三个局部坐标方向上的位移，θx、θy分别为w对单元曲面内局部坐标x、y的偏导数，θxy为w对x、y二阶交叉偏导数。
优选的，对于曲线坐标系上的壳单元，本发明中还在位移模式中补充完整的刚体位移。
进一步的，当目标实体单元为二维4节点高阶完备协调任意四边形薄板单元、二维8点高阶完备协调曲边四边形薄板单元、三维4点高阶完备协调任意四边形平板薄壳单元、三维8点高阶完备协调曲边四边形平板薄壳单元、三维4节点高阶完备协调四边形曲面薄壳单元或三维8节点高阶完备协调四边形曲面薄壳单元时，按下式修正单元边界的非协调法向转角位移：
Figure PCTCN2017086830-appb-000011
其中，Δθη23(η)，Δθη14(η)，Δθξ34(ξ)，Δθξ12(ξ)为单元四个边界的非协调法向转角位移，且该非协调法向转角位移在节点处为零。
可选的，本发明整体坐标系与线性坐标系变换的变换公式可为：
二维情况下，
Figure PCTCN2017086830-appb-000012
坐标变换关系中有6个待定系数oi，bi，ci，(i＝1，2)由变换6个为0或1的线性变换坐标值的方程组确定。
三维情况下，
Figure PCTCN2017086830-appb-000013
坐标变化关系有12个待定系数oi，bi，ci，di，(i＝1，2，3)由变换12个为0或1的线性变换坐标值的方程组确定。
可选的，本发明整体坐标系与等参系变换的变换公式可为：
Figure PCTCN2017086830-appb-000014
基于上述所构造的插值函数，然后有限元通过和原问题数学模型(基本方程、边界条件)等效的变分原理或加权余量法，建立求解基本未知量(场函数的结点值)的代数方程组或常微分方程组，求解方程组即可得问题解答。该部分即为代数方程组或常微分方程组组装和求解，该部分的技术也很成熟，已有标准定型的求解模块。藉此，本发明具有以下有益效果：
基于线性变换坐标与等参坐标混合坐标法，而不是基于单一坐标系统构造求解物理量的插值函数，可使有限元分析软件的计算精度大幅度提高，提高结构设计的安全可靠性，优化结构设计，更能适应各种曲面(曲线)边界，从而为工程、航空和航天等建设带来巨大的经济效益。
下面将参照附图，对本发明作进一步详细的说明。
附图说明
构成本申请的一部分的附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：
图1是本发明优选实施例公开的有限元插值函数构造方法流程图；
图2(a)是平面线性变换坐标系统中的整体坐标系示意图；图2(b)是平面线性变换坐标系统中线性变换坐标的示意图；
图3(a)是空间线性变换坐标系统中的整体坐标系示意图；图3(b)是空间线性变换坐标系统中线性变换坐标的示意图；
图4(a)是8节点曲边四边形单元的整体坐标系示意图；图4(b)是8节点曲边四边形单元的等参坐标系示意图；图4(c)是8节点曲边四边形单元的线性变换坐标系示意图；
图5(a)是12节点曲边四边形单元的整体坐标系示意图；图5(b)是12节点曲边四边形单元的等参坐标系示意图；图5(c)是12节点曲边四边形单元的线性变换坐标系示意图；
图6(a)是20节点曲边四边形单元的整体坐标系示意图；图6(b)是20节点曲边四边 形单元的等参坐标系示意图；图6(c)是20节点曲边四边形单元的线性变换坐标系示意图；
图7(a)是32节点曲边四边形单元的整体坐标系示意图；图7(b)是32节点曲边四边形单元的等参坐标系示意图；图7(c)是32节点曲边四边形单元的线性变换坐标系示意图；
图8(a)是4节点曲边四边形薄板单元的整体坐标系示意图；图8(b)是4节点曲边四边形薄板单元的等参坐标系示意图；图8(c)是4节点曲边四边形薄板单元的线性变换坐标系示意图；
图9(a)是8节点曲边四边形薄板单元的整体坐标系示意图；图9(b)是8节点曲边四边形薄板单元的等参坐标系示意图；图9(c)是8节点曲边四边形薄板单元的线性变换坐标系示意图；
图10(a)是空间4节点四边形薄壳单元的整体坐标系示意图；图10(b)是空间4节点四边形薄壳单元的等参坐标系示意图；图10(c)是空间4节点四边形薄壳单元的线性变换坐标系示意图；
图11(a)是空间8节点四边形薄壳单元的整体坐标系示意图；图11(b)是空间8节点四边形薄壳单元的等参坐标系示意图；图11(c)是空间8节点四边形薄壳单元的线性变换坐标系示意图；
图12(a)是空间4节点四边形曲面薄壳单元的整体、局部和曲线坐标系示意图；图12(b)是空间4节点四边形曲面薄壳单元的等参坐标系示意图；图12(c)是空间4节点四边形曲面薄壳单元的线性变换坐标系示意图；
图13(a)是空间8节点四边形曲面薄壳单元的整体、局部和曲线坐标系示意图；图13(b)是空间8节点四边形曲面薄壳单元的等参坐标系示意图；图13(c)是空间8节点四边形曲面薄壳单元的线性变换坐标系示意图。
具体实施方式
以下结合附图对本发明的实施例进行详细说明，但是本发明可以由权利要求限定和覆盖的多种不同方式实施。
实施例1
如图1所示，本实施例公开的有限元插值函数构造方法包括：
步骤S1、以线性变换坐标系与等参坐标系所组成的混合坐标系构造插值函数。
线性坐标变换系统即两直角坐标系统的变换关系为线性，等参坐标变换系统为非线性坐标变换系统。已有在曲面薄壳单元上的正交曲面坐标变换系统相当于一般结构单元的整体坐标系统，也可以变换为线性变换坐标系统。
线性变换坐标变换系统的贡献：一是能将单元变为有一直角面(线)形状单元，使部分单元节点坐标值变为简单的0和1，从而达到降低有限元插值函数构造难度和提高计算精度目的；二是使多项式的完备阶数经坐标变换不升高，有限元插值函数多项式在变换后的坐标系统是几次完备多项式，在整体坐标系统就有几次完备多项式，从而使得构造高阶完备的有限元插值函数成为可能。
线性变换坐标系统分平面线性变换坐标系统和空间线性变换坐标系统。
(a)、针对平面线性变换坐标系统，设坐标变换关系为：
Figure PCTCN2017086830-appb-000015
单元经过坐标变换后的形状如图2(a)及图2(b)所示，四边形单元一个角点位于坐标原点，两个角点位于坐标轴上。坐标变换关系中有6个待定系数oi，bi，ci，(i＝1，2)，可变换6个为0或1的线性变换坐标值。
(b)、针对空间线性变换坐标系统，设坐标变化关系为：
Figure PCTCN2017086830-appb-000016
单元经过坐标变化后的形状如图3(a)及图3(b)所示，六面体单元一个角点位于坐标原点，三个角点位于坐标轴上。坐标变化关系有12个待定系数oi，bi，ci，di，(i＝1，2，3)，可变换12个为0或1的线性变换坐标值。
另一方面，该步骤中，整体坐标系与等参系变换的变换公式为：
Figure PCTCN2017086830-appb-000017
步骤S2、根据目标实体单元的特征确定插值函数方程式的坐标元数、项数和次数，该特征包括已知节点数和节点位移分量数，且所构造插值函数方程式的项数与目标实体单元的插值相关节点位移总数相等；同时，所构造插值函数方程式中对应线性变换坐标系部分为各项次数以覆盖各坐标元组合后从低到高递增的最高次完整多项式，所构造插值函数方程式中对应等参坐标系部分的各项次数在所述线性变换坐标系部分的最高次数上以覆盖各坐标元组合后从低到高递增并呈对称性分布。
特殊情况下，如果对应线性变换坐标系部分的各项总数刚好等于所构造插值函数方程式的总项数，所构造插值函数方程式中将不出现等参坐标系部分。
藉此，本实施例基于线性变换坐标与等参坐标混合坐标系统构造有限元插值函数。在有限单元插值多项式的选择中，有限单元插值多项式的完备项使用线性变换坐标，而多余项使用等参坐标。将线性变换坐标用等参坐标表示，再代入有限单元插值多项式，所得等参坐标表示的多项式次数不高于有限单元插值多项式的完备阶数，从而保证有限元插值函数在单元边界上具有C0阶协调，否则，如果有限单元插值多项式存在非完备项，且不用等参坐标替换，将不能保证有限元插值函数在单元边界上C0阶协调，更不可能保证有限元插值函数在单元边界上C1阶协调。从而，非常巧妙地解决构造高阶完备有限元插值函数而导致不协调的难题。
下面针对各目标实体单元的构造函数分数如下：
1)、如图4(a)、图4(b)及图4(c)所示，当目标实体单元为二维8节点高阶完备协调四边形曲边单元时，所构造单元位移插值函数为：
u(v)＝a1+a2T1+a3T2+a4T1 2+a5T1T2+a6T1 2+a7ξ2η+a8ξη2。
下述各方程式中，T1、T2、T3分别为单元曲面内线性变换坐标系中坐标轴，ξ、η、ζ分别为等参坐标系中坐标轴；u、v、w分别对应单元曲面内三个局部坐标方向上的位移，θx、θy分别为w对单元曲面内局部坐标x、y的偏导数，θxy为w对x、y二阶交叉偏导数；且下述各情况中，对应u(v)、u(v，w)的插值相关节点位移总数与已知节点数相等，对应w的插值相关节点位移总数为已知节点数与相关位移分量的乘积，后续不再赘述。
该单元位移插值函数可根据单元节点位移建立方程组，联立求解确定待定系数a1～a8，其具有2阶完备性，8节点四边形单元也只能具有2阶完备性，比传统单元位移插值函数的完备性高1阶。单元位移插值函数的完备性提高1次，其收敛性能和抗畸变性能会有非常大的提高。
当两相邻单元的公共边界为直线且边界中节点均分边界时，该公共边界即协调，与单元其它边界形状无关，因此，在划分单元时，只要保持两相邻单元的公共边界为直线且边界中节点均分边界，而单元的自由外边界没有协调性要求可为曲线，此时的二维8节点单元位移插值函数即具有高阶完备协调性，且适合于曲线边界，并没有增加单元划分的难度。此二维8节点四边形曲边单元可退化6节点三角形曲边单元。
2)、如图5(a)、图5(b)及图5(c)所示，当目标实体单元为二维12节点高阶完备协调四边形曲边单元时，所构造单元位移插值函数为：
Figure PCTCN2017086830-appb-000018
该单元位移插值函数可根据单元节点位移建立方程组，联立求解确定待定系数a1～a12，其具有3阶完备性，比传统单元位移插值函数的完备性高2阶，其收敛性能和抗畸变性能非常好。当两相邻单元的公共边界为直线且边界中节点均分边界时，该公共边界即协调，与单元其它边界形状无关，因此，在划分单元时，只要保持两相邻单元的公共边界为直线且边界中节点均分边界，而单元的自由边界没有协调性要求可为曲线，此时的二维12节点单元位移插值函数即具有高阶完备协调性，并能适合于任意曲线边界，并没有曾加单元划分的难度。此二维12节点四边形曲边单元可退化9节点三角形曲边单元。
3)、如图6(a)、图6(b)及图6(c)所示，当目标实体单元为三维20节点高阶完备协调曲面六面体单元时，所构造单元位移插值函数为：
Figure PCTCN2017086830-appb-000019
该单元位移插值函数可根据单元节点位移建立方程组，联立求解确定待定系数a1～a20，其具有2阶完备性，比传统单元位移插值函数的完备性高1阶。该曲边六面体单元可以自动退化为曲边五面体单元、曲边四面体单元；六个四边形侧曲面可以自动退化为三角形侧曲面。
有关单元位移插值函数在单元边界处的协调问题分为三种情况：
(a)六面体的公共棱角线协调
与平面问题类似，要求单元公共棱角线为直线且中节点均分其公共棱角线，单元位移插值函数在公共棱角线上即协调，与单元其它部位的形状无关。
(b)六面体的四边形侧曲面协调
要求单元公共的四边形侧曲面为平面、四条边均为直线，且边中的节点均分其边界，此时单元位移插值函数在公共四边形侧平面上即协调，与单元其它面的形状无关。
(c)六面体的三角形侧曲面协调
只要求单元公共的三角形侧曲面为平面，三角形的边可为曲线，单元位移插值函数在公共三角形侧平面上即协调，与单元其它面的形状无关。三角形的边可为曲线条件，使得构造适合于任意曲边边界的退化四面体单元成为可能。
在单元划分时，结构内部用平面六面体单元，可用退化后的四面体单元、五面体单元和六面体单元外露曲面模似结构的外部曲面边界，从而就可保证了单元位移插值函数的协调，同时也能很好地模似结构的曲面边界。与常规的单元划分方法一致，并没有增加单元划分的难度。
4)、如图7(a)、图7(b)及图7(c)所示，当所述目标实体单元为三维32节点高阶完备协调曲面六面体单元时，所构造单元位移插值函数为：
Figure PCTCN2017086830-appb-000020
该单元位移插值函数可根据单元节点位移建立方程组，联立求解确定待定系数a1～a32，其具有3阶完备性，比传统单元位移插值函数的完备性高2阶。该曲边六面体单元同样可以自动退化为曲边五面体单元、曲边四面体单元；六个四边形侧曲面可以自动退化为三角形侧曲面，以适合曲面结构边界。
有关单元位移插值函数在单元边界处的协调问题分为三种情况：
(a)六面体的公共棱角线协调
与平面问题类似，要求单元公共棱角线为直线且中节点均分其公共棱角线，单元位移插值函数在其公共棱角线上既协调，与单元其它部位的形状无关。
(b)六面体的四边形侧曲面协调
要求单元公共的四边形侧曲面为平面、四条边均为直线，且边中的节点均分其边，此时单元位移插值函数在公共四边形侧平面上既协调，与单元其它面的形状无关。
(c)六面体的三角形侧曲面协调
只要求单元公共的三角形侧曲面为平面，三角形的边可为曲线，单元位移插值函数在公共三角形侧平面上即协调，与单元其它面的形状无关。三角形的边可为曲线的条件，使得构造适合于任意曲边边界的退化四面体单元成为可能。
在单元划分时，结构内部用平面六面体单元，可用退化后的四面体单元、五面体单元和六面体单元外露曲面模似结构的外部曲面边界，从而就可保证了单元位移插值函数的协调，同时也能很好地模似结构的曲面边界。与常规的单元划分方法一致，并没有增加单元划分的难度。
5)、如图8(a)、图8(b)及图8(c)所示，当目标实体单元为二维4节点且各节点有3个相关位移分量(三节点参数为w，θx，θy，w，θx，θy相关，该情况下，位移总数为已知节点数和节点位移分量数的乘积)的高阶完备协调任意四边形薄板单元时，所构造单元位移插值函数为：
Figure PCTCN2017086830-appb-000021
该单元位移插值函数可根据单元节点位移建立方程组，联立求解确定待定系数a1～a12，其具有3阶完备性，完备的阶数较高，该单元位移插值函数在相邻单元的公共边界上挠度和切向转角能协调。而传统的四边形薄板单元只有矩形单元和三角形单元，不适合于任意折线形边界，应用范围非常有限。
该单元位移插值函数的在单元边界上挠度和切向转角协调，但法向转角不协调，即C1阶非协调性问题，该问题一直被认为是不能解决的难题。为此，本实施例在规则的等参坐标系里，用一种非常特殊的修正函数仅对单元一条边界的非协调法向转角进行修正，不影响单元其它边界的位移和转角值，关键是不影响单元其它边界的法向转角值，因此，可只对单元有非协调法向转角公共边界的进行修正，单元自由边界的非协调法向转角不必修正。
优选地，本实施例可按下式修正单元边界的非协调法向转角位移：
Figure PCTCN2017086830-appb-000022
其中，Δθη23(η)，Δθη14(η)，Δθξ34(ξ)，Δθξ12(ξ)为单元四个边界的非协调法向转角位移，且该非协调法向转角位移在节点处为零。可选择性修正单元任一边界的非协调法向转角位移。
6)、如图8(a)、图8(b)及图8(c)所示，当目标实体单元为二维4节点且各节点有4个相关位移分量(四节点参数为w，θx，θy，θxy，w，θx，θy，θxy相关，该情况下，位移总数为已知节点数和节点位移分量数的乘积，其中，θxy为w对x、y二阶交叉偏导数)的高阶完备协调任意四边形薄板单元时，所构造单元位移插值函数为：
Figure PCTCN2017086830-appb-000023
该单元位移插值函数可根据单元节点位移建立方程组，联立求解确定待定系数a1～a16，其考虑四节点参数w，θx，θy，θxy，可保证单元位移插值函数的在单元边界上法向转角协调，而不必进行特殊修正；且具有3阶完备性，完备的阶数较高，该单元位移插值函数在相邻单元的公共边界上，不仅挠度和切向转角协调，而且法向转角也协调，但要求板单元厚度和材料性质一致，否则，会引起新的不协调问题。
7)、如图9(a)、图9(b)及图9(c)所示，当目标实体单元为二维8节点且各节点有3个相关位移分量的高阶完备协调曲边四边形薄板单元时，所构造单元位移插值函数为：
Figure PCTCN2017086830-appb-000024
该单元位移插值函数可根据单元节点位移建立方程组，联立求解确定待定系数a1～a24，其具有5阶完备性，完备的阶数较高，适合于任意曲线形边界。在划分单元时，只要保持两相邻单元的公共边界为直线且边界中节点均分边界，而单元的自由外边界没有协调性要求可为曲线，以适合曲线边界，此时，单元位移插值函数在单元公共边界上挠度和切向转角才协调，但法向转角不能协调须修正。
当两相邻单元的公共边界为直线且边界中节点均分边界时，该单元位移插值函数的在单元边界上挠度和切向转角协调，但法向转角不协调，即C1阶非协调性问题，该问题一直被认为是不能解决的难题。为此，本实施例在规则的等参坐标系里，用一种非常特殊的修正函数可仅对单元一条边界的非协调法向转角进行修正，不影响单元其它边界的挠度和切向转角值，关键是不影响单元其它边界的法向转角值，可只对有非协调法向转角单元的公共边界进行修正，单元自由边界的非协调法向转角不必修正，因此，单元外边界可为曲线，以适合曲线边界。优选地，本实施例按下式修正单元边界的非协调法向转角位移：
Figure PCTCN2017086830-appb-000025
其中，Δθη23(η)，Δθη14(η)，Δθξ34(ξ)，Δθξ12(ξ)为单元边界的非协调法向转角位移，该非协调法向转角位移在节点处为零。可选择性修正单元任一边界的非协调法向转角位移。
8)、当目标实体单元为三维4节点高阶完备协调任意四边形平板薄壳单元(其中w，θx，θy相关)时，即对于4节点平面任意四边形薄壳单元，如图10(a)、图10(b)及图10(c)所示，通过坐标变换，可将任意四边形平面薄板单元位移插值函数转换为空间薄壳单元位移插值函数，其中节点位移向量在两坐标系之间转换关系为：
a′i＝La；ai＝LTa′i ai＝[ui vi wi θxi θyi θzi]；a′i＝[u′i v′i w′i θ′xi θ′yi θ′zi]
Figure PCTCN2017086830-appb-000026
其中，a′i，ai为节点位移；T，L，λ为转换矩阵。
单元刚度矩阵和荷载列向量的转换关系为：
K′e＝TKeT；Q′e＝TQ
其它转换关系按常规法进行实施。
关键是任意四边形薄壳单元位移插值函数的构造，基于本发明的线性变换坐标法与等参坐标混合法，可假定单元位移插值函数为：
Figure PCTCN2017086830-appb-000027
该插值函数可根据单元节点位移建立方程组，联立求解确定待定系数；具有3阶完备性，完备的阶数较高，且该薄壳单元位移插值函数在相邻单元的公共边界上的挠度和切向转角能协调。
该薄壳单元位移插值函数的在单元边界上挠度和切向转角协调，但法向转角不协调，即C1阶非协调性问题。为此，本实施例在规则的等参坐标系里，用一种非常特殊的修正函数仅对单元一条边界的非协调法向转角进行修正，不影响单元边界的挠度和切向转角值，关键是不影响单元其它边界的法向转角值，因此，可只对单元公共边界的非协调法向转角进行修正，单元自由边界的非协调法向转角不必修正。
优选地，本实施例按下式修正薄壳单元边界的非协调法向转角位移：
Figure PCTCN2017086830-appb-000028
其中，Δθη23(η)，Δθη14(η)，Δθξ34(ξ)，Δθξ12(ξ)为单元边界的非协调法向转角位移，该非协调法向转角位移在节点处为零。可选择性修正单元任一边界的非协调法向转角位移。
9)、当目标实体单元为三维8节点高阶完备协调曲线四边形平板薄壳单元(其中w，θx，θy相关)时，即对于8节点平面任意四边形薄壳单元，如图11(a)、图11(b)及图11(c)所示，通过坐标变换，可将任意四边形平面薄板单元位移插值函数转换为空间薄壳单元位移插值函数。其中节点位移向量在两坐标系之间转换关系为：
a′i＝Lai；ai＝LTa′i
ai＝[ui vi wi θxi θyi θzi]；a′i＝[u′i v′i w′i θ′xi θ′yi θ′zi]
Figure PCTCN2017086830-appb-000029
其中，a′i，ai为节点位移；T，L，λ为转换矩阵。
单元刚度矩阵和荷载列向量的转换关系为：
K′e＝TKeT；Q′e＝TQ
其它转换关系按常规法进行实施。
关键是任意四边形薄壳单元位移插值函数的构造，基于本发明的线性变换坐标法与等参坐标混合法，可假定单元位移插值函数为：
Figure PCTCN2017086830-appb-000030
上述插值函数可根据单元节点位移建立方程组，联立求解确定待定系数；其具有5阶完备性，完备的阶数较高，适合于曲线形边界。在划分单元时，只要保持两相邻单元的公共边界为直线且边界中节点均分边界，而单元的自由外边界没有协调性要求可为曲线，以适合曲线边界，此时，单元位移插值函数在单元公共边界上挠度和切向转角才协调，但法向转角不能协调须修正。
当两相邻单元的公共边界为直线且边界中节点均分边界时，该单元位移插值函数的在单元边界上挠度和切向转角协调，但法向转角不协调，即C1阶非协调性问题。为此，本实施例在规则的等参坐标系里，用一种非常特殊的修正函数仅对单元一条边界的非协调法向转角进行修正，不影响单元边界的挠度和切向转角值，关键是不影响单元其它边界的法向转角值，因此，可只对单元公共边界的非协调法向转角进行修正，单元自由边界的非协调法向转角不必修正。
优选地，本实施例按下式修正单元边界的非协调法向转角位移：
Figure PCTCN2017086830-appb-000031
其中Δθη23(η)，Δθη14(η)，Δθξ34(ξ)，Δθξ12(ξ)为单元边界的非协调法向转角位移，且该非 协调法向转角位移在节点处为零。可选择性修正单元任一边界的非协调法向转角位移。
10)、当目标实体单元为三维4节点高阶完备协调四边形曲面薄壳单元(其中w，θx，θy相关)时，如图12(a)、图12(b)及图12(c)所示，可假定曲面薄壳单元内部任意一点的总体坐标为：
Figure PCTCN2017086830-appb-000032
其中，Ni(ξ，η)为常规的形函数。
在整体坐标系中节点的刚体位移{δ′Ri}由单元微体的刚体运动给出，微体形心的运动包括绕三个坐标轴的转动和沿坐标轴的平动。曲面薄壳单元形心的刚体运动为：
{V′R}＝{u′0 v′0 w′0 θ′x0 θ′y0 θ′z0}T
在整体坐标系里，按动力学方法得到这个刚体运动产生的节点刚体位移为：
{δ′Ri}＝{u′Ri v′Ri w′Ri θ′xRi θ′yRi θ′zRi}T＝[TRi]{V′R}
Figure PCTCN2017086830-appb-000033
其中，转换矩阵TR的子矩阵：
[TR]＝{TR1 TR2 TR3 TR4}T
Figure PCTCN2017086830-appb-000034
式中，x0，y0，z0为单元微体形心整体坐标。
在曲线坐标系里，刚体运动产生的节点刚体位移为：
{δRi}＝{uRi vRi wRi θxRi θyRi θzRi}T＝[Li]T[TRiIL0]T{VR}
Figure PCTCN2017086830-appb-000035
其中，Li为转换矩阵。
Figure PCTCN2017086830-appb-000036
Figure PCTCN2017086830-appb-000037
式中λx′αi＝cos(x′，α)等是x、y、z轴在正交主曲线坐标系α，β，δ的各方向余弦。
在曲线坐标系里，曲面薄壳单元形心的刚体运动为：
{VR}＝{u0 v0 w0 θx0 θy0 θz0}T
单元节点位移向量在整体坐标和曲线坐标之间转换关系为：
Figure PCTCN2017086830-appb-000038
在曲线坐标系里，补充刚体位移后的单元位移场全量为：
Figure PCTCN2017086830-appb-000039
式中I是20×20的单位矩阵。
刚体位移不产生应变，所以，在正交主曲线坐标系中应变矩阵为：
Figure PCTCN2017086830-appb-000040
式中，B是正交主曲线坐标系中单元的应变矩阵。
刚体位移也不产生节点力，作静力凝聚，按虚功原理建立有限元方程：
Figure PCTCN2017086830-appb-000041
其中，
Figure PCTCN2017086830-appb-000042
是在正交主曲线坐标系里单元的节点载荷，
Figure PCTCN2017086830-appb-000043
是在正交主曲线坐标系里由单元原有的刚度矩阵
Figure PCTCN2017086830-appb-000044
经位移扩展后得到的，即
Figure PCTCN2017086830-appb-000045
其中，
Figure PCTCN2017086830-appb-000046
Figure PCTCN2017086830-appb-000047
Figure PCTCN2017086830-appb-000048
在整体坐标系刚度矩阵和荷载向量的转换关系为：
Figure PCTCN2017086830-appb-000049
集成总体坐标系的各个单元刚度矩阵和荷载向量，就可以得到整体坐标系的节点位移求解方程。
在正交主曲线坐标系(α，β，δ)里，α，β，δ为曲线坐标的弧长，对于4节点任意四边形单元，如图11所示，基于独创的线性变换坐标法与等参坐标混合法，与平面薄板问题相似，对曲线弧长坐标进行线性坐标变换与等参坐标变换，可假定单元位移插值函数为：
Figure PCTCN2017086830-appb-000050
根据单元节点位移建立方程组，联立求解确定待定系数。该单元位移插值函数具有3阶完备性，完备的阶数较高，该单元位移插值函数在相邻单元的公共边界上挠度和切向转角能协调。该单元位移插值函数的在单元边界上挠度和切向转角协调，但法向转角不协调，即C1阶非协调性问题。为此，在规则的等参坐标系里，本实施例用一种非常特殊的修正函数仅对单元一条边界的非协调法向转角进行修正，不影响单元其它边界的位移和转角值，关键是不影响单元其它边界的法向转角值，因此，可只对单元有非协调法向转角公共边界的进行修正，单元自由边界的非协调法向转角不必修正。
优选地，本实施例按下式修正单元边界的非协调法向转角位移：
Figure PCTCN2017086830-appb-000051
其中，Δθη23(η)，Δθη14(η)，Δθξ34(ξ)，Δθξ12(ξ)为单元四个边界的非协调法向转角位移，且该非协调法向转角位移在节点处为零。可选择性修正单元任一边界的非协调法向转角位移。
11)、当目标实体单元为三维8节点高阶完备协调四边形曲面薄壳单元(其中w，θx，θy相关)时，即对于8节点曲面薄壳单元，如图13(a)、图13(b)及图13(c)所示，假定曲壳单元内部任意一点的总体坐标为：
Figure PCTCN2017086830-appb-000052
其中，Ni(ξ，η)为常规的形函数。
在整体坐标系中节点的刚体位移{δ′Ri}由单元微体的刚体运动给出，微体形心的运动包括绕三个坐标轴的转动和沿坐标轴的平动。壳单元形心的刚体运动为：
{V′R}＝{u′0 v′0 w′0 θ′x0 θ′y0 θ′z0}T
在整体坐标系里，按动力学方法得到这个刚体运动产生的节点刚体位移为：
{δ′Ri}＝{u′Ri v′Ri w′Ri θ′xRi θ′yRi θ′zRi}T＝[TRi]{V′R}
Figure PCTCN2017086830-appb-000053
其中，转换矩阵TR的子矩阵：
[TR]＝{TR1 TR2 TR3 TR4 TR5 TR6 TR7 TR8}T
Figure PCTCN2017086830-appb-000054
式中x0，y0，z0为单元微体形心整体坐标。
在曲线坐标系里，刚体运动产生的节点刚体位移为：
{δRi}＝{uRi vRi wRi θxRi θyRi θzRi}T＝[Li]T[TRiIL0]T{VR}
Figure PCTCN2017086830-appb-000055
其中，Li为转换矩阵。
Figure PCTCN2017086830-appb-000056
Figure PCTCN2017086830-appb-000057
式中λx′αi＝cos(x′，α)等是x、y、z轴在正交主曲线坐标系α，β，δ的各方向余弦。
在曲线坐标系里，曲面薄壳单元形心的刚体运动为：
{VR}＝{u0 v0 w0 θx0 θy0 θz0}T
单元节点位移向量在整体坐标和曲线坐标之间转换关系为：
Figure PCTCN2017086830-appb-000058
在曲线坐标系里，补充刚体位移后的单元位移场全量为：
Figure PCTCN2017086830-appb-000059
式中I是40×40的单位矩阵。
刚体位移不产生应变，所以，在正交主曲线坐标系中应变矩阵为：
Figure PCTCN2017086830-appb-000060
式中B是正交主曲线坐标系中单元的应变矩阵。
刚体位移也不产生节点力，作静力凝聚，按虚功原理建立有限元方程：
Figure PCTCN2017086830-appb-000061
其中，
Figure PCTCN2017086830-appb-000062
是在正交主曲线坐标系里单元的节点载荷，
Figure PCTCN2017086830-appb-000063
是在正交主曲线坐标系里由单元原有的刚度矩阵
Figure PCTCN2017086830-appb-000064
经位移扩展后得到的，即
Figure PCTCN2017086830-appb-000065
其中，
Figure PCTCN2017086830-appb-000066
Figure PCTCN2017086830-appb-000067
Figure PCTCN2017086830-appb-000068
在整体坐标系刚度矩阵和荷载向量的转换关系为：
Figure PCTCN2017086830-appb-000069
集成总体坐标系的各个单元刚度矩阵和荷载向量，就可以得到整体坐标系的节点位移求解方程。
在正交主曲线坐标系(α，β，δ)里，α，β，δ为曲线坐标的弧长，对于8节点任意四边形单 元，基于独创的线性变换坐标法与等参坐标混合法，与平面薄板问题相似，对曲线弧长坐标进行线性坐标变换与等参坐标变换，可假定单元位移插值函数为：
Figure PCTCN2017086830-appb-000070
可根据单元节点位移建立方程组，联立求解确定待定系数。该单元位移插值函数具有5阶完备性，完备的阶数较高，适合于任意曲线形边界。在划分单元时，只要保持两相邻单元的公共边界为直线且边界中节点均分边界，而单元的自由外边界没有协调性要求可为曲线，以适合曲线边界，此时，单元位移插值函数在单元公共边界上挠度和切向转角能协调，但法向转角不能协调须修正。
当两相邻单元的公共边界为直线且边界中节点均分边界时，该单元位移插值函数的在单元边界上挠度和切向转角协调，但法向转角不协调，即C1阶非协调性问题，一个世纪以来，一直被认为是不能解决的难题。为此，在规则的等参坐标系里，首次提出用一种非常特殊的修正函数可仅对单元一条边界的非协调法向转角进行修正，不影响单元其它边界的挠度和切向转角值，关键是不影响单元其它边界的法向转角值，可只对有非协调法向转角单元的公共边界进行修正，单元自由边界的非协调法向转角不必修正，因此，单元外边界可为曲线，以适合曲线边界。
经理论分析首次提出按下式修正单元边界的非协调法向转角位移：
Figure PCTCN2017086830-appb-000071
其中，Δθη23(η)，Δθη14(η)，Δθξ34(ξ)，Δθξ12(ξ)为单元边界的非协调法向转角位移，该非协调法向转角位移在节点处为零。可选择性修正单元任一边界的非协调法向转角位移。
综上，本实施例基于线性变换坐标与等参坐标混合坐标法，而不是基于单一坐标系统构造求解物理量的插值函数，可使有限元分析软件的计算精度大幅度提高，提高结构设计的安全可靠性，优化结构设计，更能适应各种曲面(曲线)边界，从而为工程、航空和航天等建设带来巨大的经济效益。
以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员 来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (7)

  1. 一种有限元插值函数构造方法，其特征在于，包括：
    以线性变换坐标系与等参坐标系所组成的混合坐标系构造插值函数；
    根据目标实体单元的特征确定插值函数方程式的坐标元数、项数和次数，所述特征包括已知节点数和节点位移分量数，且所构造插值函数方程式的项数与目标实体单元的插值相关节点位移总数相等；同时，所构造插值函数方程式中对应线性变换坐标系部分为各项次数以覆盖各坐标元组合后从低到高递增的最高次完整多项式，所构造插值函数方程式中对应等参坐标系部分的各项次数在所述线性变换坐标系部分的最高次数上以覆盖各坐标元组合后从低到高递增并呈对称性分布。
  2. 根据权利要求1所述的有限元插值函数构造方法，其特征在于，所构造的插值函数包括以下的任意一项或任意组合：
    1)、当所述目标实体单元为二维8节点高阶完备协调四边形曲边单元时，所构造单元位移插值函数为：
    u(v)＝a1+a2T1+a3T2+a4T1 2+a5T1T2+a6T1 2+a7ξ2η+a8ξη2；
    2)、当所述目标实体单元为二维12节点高阶完备协调四边形曲边单元时，所构造单元位移插值函数为：
    Figure PCTCN2017086830-appb-100001
    3)、当所述目标实体单元为三维20节点高阶完备协调曲面六面体单元时，所构造单元位移插值函数为：
    Figure PCTCN2017086830-appb-100002
    4)、当所述目标实体单元为三维32节点高阶完备协调曲面六面体单元时，所构造单元位移插值函数为：
    Figure PCTCN2017086830-appb-100003
    5)、当所述目标实体单元为二维4节点且各节点有3个相关位移分量的高阶完备协调任意四边形薄板单元时，相关位移分量分别为w、θx、θy，所构造单元位移插值函数为：
    Figure PCTCN2017086830-appb-100004
    6)、当所述目标实体单元为二维4节点且各节点有4个相关位移分量的高阶完备协调任意四边形薄板单元时，相关位移分量分别为w、θx、θy、θxy，所构造单元位移插值函数为：
    Figure PCTCN2017086830-appb-100005
    7)、当所述目标实体单元为二维8节点且各节点有3个相关位移分量的高阶完备协调曲边四边形薄板单元时，相关位移分量分别为w、θx、θy，所构造单元位移插值函数为：
    Figure PCTCN2017086830-appb-100006
    其中，上述各方程式中，T1、T2、T3分别为单元曲面内线性变换坐标系中坐标轴，ξ、η、ζ分别为等参坐标系中坐标轴；u、v、w分别对应单元曲面内三个局部坐标方向上的位移，θx、θy分别为w对单元曲面内局部坐标x、y的偏导数，θxy为w对x、y二阶交叉偏导数。
  3. 根据权利要求1所述的有限元插值函数构造方法，其特征在于，对于工程上普遍采用的空间薄壳，采用合适的正交曲线坐标及相应的几何方程，根据上述单元的原理，像平面问题一样在空间正交曲线坐标系中直接构造高阶完备协调曲面薄壳单元，计算单元刚度矩阵，再进行空间坐标转换；具体包括：
    1)、当所述目标实体单元为三维4节点且各节点w位移分量有3个相关位移分量的高阶完备协调任意四边形平板薄壳单元时，相关位移分量分别为w、θx、θy，所构造单元位移插值函数为：
    Figure PCTCN2017086830-appb-100007
    2)、当所述目标实体单元为三维8节点且各节点w位移分量有3个相关位移分量的高阶完备协调曲线四边形平板薄壳单元时，相关位移分量分别为w、θx、θy，所构造单元位移插值函数为：
    Figure PCTCN2017086830-appb-100008
    3)、当所述目标实体单元为三维4节点且各节点w位移分量有3个相关位移分量的高阶完备协调四边形曲面薄壳单元时，相关位移分量分别为w、θx、θy，所构造单元位移插值函数为：
    Figure PCTCN2017086830-appb-100009
    4)、当所述目标实体单元为三维8节点且各节点w位移分量有3个相关位移分量的高阶完备协调四边形曲面薄壳单元时，相关位移分量分别为w、θx、θy，所构造单元位移插值函数为：
    Figure PCTCN2017086830-appb-100010
    其中，上述各方程式中，T1、T2、T3分别为单元曲面内线性变换坐标系中坐标轴，ξ、η、ζ分别为等参坐标系中坐标轴；u、v、w分别对应单元曲面内三个局部坐标方向上的位移，θx、θy分别为w对单元曲面内局部坐标x、y的偏导数，θxy为w对x、y二阶交叉偏导数。
  4. 根据权利要求3所述的有限元插值函数构造方法，其特征在于，还包括：
    对于曲线坐标系上的壳单元，在位移模式中补充完整的刚体位移。
  5. 根据权利要求1至4任一所述的有限元插值函数构造方法，其特征在于，还包括：
    当所述目标实体单元为二维4节点高阶完备协调任意四边形薄板单元、二维8点高阶完备协调曲边四边形薄板单元、三维4点高阶完备协调任意四边形平板薄壳单元、三维8点高阶完备协调曲边四边形平板薄壳单元、三维4节点高阶完备协调四边形曲面薄壳单元或三维8节点高阶完备协调四边形曲面薄壳单元，且单元各节点只有3个相关位移分量时，按下式修正单元边界的非协调法向转角位移：
    Figure PCTCN2017086830-appb-100011
    其中，Δθη23(η)，Δθη14(η)，Δθξ34(ξ)，Δθξ12(ξ)为单元四个边界的非协调法向转角位移，且该非协调法向转角位移在节点处为零。
  6. 根据权利要求5所述的有限元插值函数构造方法，其特征在于，整体坐标系与线性坐标系变换的变换公式为：
    二维情况下，
    Figure PCTCN2017086830-appb-100012
    坐标变换关系中有6个待定系数oi，bi，ci，(i＝1，2)由变换6个为0或1的线性变换坐标值的方程组确定；
    三维情况下，
    Figure PCTCN2017086830-appb-100013
    坐标变化关系有12个待定系数oi，bi，ci，di，(i＝1，2，3)由变换12个为0或1的线性变换坐标值的方程组确定。
  7. 根据权利要求5所述的有限元插值函数构造方法，其特征在于，整体坐标系与等参系变换的变换公式为：
    Figure PCTCN2017086830-appb-100014
PCT/CN2017/086830 2017-04-26 2017-06-01 有限元插值函数构造方法 WO2018196098A1 (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
US15/749,021 US20190005171A1 (en) 2017-04-26 2017-06-01 Method for constructing finite element interpolation function

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
CN201710307000.X 2017-04-26
CN201710307000.XA CN107193780A (zh) 2017-04-26 2017-04-26 有限元插值函数构造方法

Publications (1)

Publication Number Publication Date
WO2018196098A1 true WO2018196098A1 (zh) 2018-11-01

Family

ID=59873371

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/CN2017/086830 WO2018196098A1 (zh) 2017-04-26 2017-06-01 有限元插值函数构造方法

Country Status (3)

Country Link
US (1) US20190005171A1 (zh)
CN (1) CN107193780A (zh)
WO (1) WO2018196098A1 (zh)

Families Citing this family (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US10394977B2 (en) * 2014-06-06 2019-08-27 Robert E. Spears Method and apparatus for shape-based energy analysis of solids
US10910325B2 (en) * 2017-05-29 2021-02-02 Intel Corporation Integrated circuit packages with conductive element having cavities housing electrically connected embedded components
CN108416079A (zh) * 2018-01-15 2018-08-17 中南大学 优化有限元软件计算精度的非协调插值函数构造方法、系统及存储介质
CN111292422A (zh) * 2020-03-02 2020-06-16 上海索辰信息科技有限公司 一种基于曲面法线改进网格生成质量的方法
CN112613211B (zh) * 2020-12-22 2022-10-21 郑州大学 平面结构中任意三角形单元的变形分解方法
CN114913218B (zh) * 2022-04-29 2024-04-30 中国人民解放军63891部队 一种基于极小模网络的烟幕红外遮蔽面积测量方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102682152A (zh) * 2011-03-01 2012-09-19 利弗莫尔软件技术公司 使用扩展的无网格有限元法对结构性能的数值仿真
JP2014153999A (ja) * 2013-02-12 2014-08-25 Toyota Communication Systems Co Ltd 有限要素法解析方法及び有限要素法解析装置
CN104317985A (zh) * 2014-09-19 2015-01-28 大连理工大学 一种基于界带有限元和拉格朗日坐标的流体仿真方法

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5497451A (en) * 1992-01-22 1996-03-05 Holmes; David Computerized method for decomposing a geometric model of surface or volume into finite elements
JP2744888B2 (ja) * 1993-09-10 1998-04-28 インターナショナル・ビジネス・マシーンズ・コーポレイション ３−ｄオブジェクトを領域に区画する方法及びシステム
US6330704B1 (en) * 1999-02-09 2001-12-11 Coyote Systems, Inc. Method and apparatus for three dimensional interconnect analysis
US6876956B1 (en) * 1999-08-31 2005-04-05 California Institute Of Technology Method and system for thin-shell finite-element analysis
US6816820B1 (en) * 1999-09-24 2004-11-09 Moldflow Ireland, Ltd. Method and apparatus for modeling injection of a fluid in a mold cavity
US10296683B2 (en) * 2013-11-06 2019-05-21 Thornton Tomasetti, Inc. Computer implemented apparatus and method for finite element modeling using hybrid absorbing element

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102682152A (zh) * 2011-03-01 2012-09-19 利弗莫尔软件技术公司 使用扩展的无网格有限元法对结构性能的数值仿真
JP2014153999A (ja) * 2013-02-12 2014-08-25 Toyota Communication Systems Co Ltd 有限要素法解析方法及び有限要素法解析装置
CN104317985A (zh) * 2014-09-19 2015-01-28 大连理工大学 一种基于界带有限元和拉格朗日坐标的流体仿真方法

Also Published As

Publication number Publication date
CN107193780A (zh) 2017-09-22
US20190005171A1 (en) 2019-01-03

Similar Documents

Publication Publication Date Title
WO2018196098A1 (zh) 有限元插值函数构造方法
Tornabene et al. Static analysis of laminated composite curved shells and panels of revolution with a posteriori shear and normal stress recovery using generalized differential quadrature method
Dhatt et al. Finite element method
Lee et al. A new strain smoothing method for triangular and tetrahedral finite elements
Vizotto Computational generation of free-form shells in architectural design and civil engineering
Wildberger Universal hyperbolic geometry I: trigonometry
Eliseev et al. Finite deformation of thin shells in the context of analytical mechanics of material surfaces
Xie et al. A novel triangular boundary crack front element for 3D crack problems based on 8-node serendipity element
CN102096736A (zh) 一种基于渐近变分法的复合材料层合板仿真及优化方法
Wang et al. A 4-node quasi-conforming Reissner–Mindlin shell element by using Timoshenko's beam function
Kiseleva et al. Comparison of scalar and vector FEM forms in the case of an elliptic cylinder
Li et al. A new 8-node quadrilateral spline finite element
Yuan et al. The inverse mapping and distortion measures for 8-node hexahedral isoparametric elements
WO2019137236A1 (zh) 优化有限元软件计算精度的非协调插值函数构造方法、系统及存储介质
Kwak et al. Free vibration analysis of laminated rectangular plates with varying thickness using Legendre-radial point interpolation method
Dey et al. An efficient plate element for the vibration of composite plates
Peng et al. Advances in the Base Force Element Method
Abenda et al. Real Soliton Lattices of the Kadomtsev Petviashvili II Equation and Desingularization of Spectral Curves: The Gr TP (2, 4) Case
Chen et al. A 3D pyramid spline element
Sanz et al. Three-dimensional BEM for piezoelectric fracture analysis
Borodzik et al. Plane algebraic curves of arbitrary genus via Heegaard Floer homology
Walkington A C^1 Tetrahedral Finite Element without Edge Degrees of Freedom
Wankhade Geometric nonlinear analysis of skew plates using finite element method
Schoop A simple nonlinear flat element for large displacement structures
Dorninger A nonlinear layered shell finite element with improved transverse shear behavior

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 17908083

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 17908083

Country of ref document: EP

Kind code of ref document: A1

32PN Ep: public notification in the ep bulletin as address of the adressee cannot be established

Free format text: NOTING OF LOSS OF RIGHTS PURSUANT TO RULE 112(1) EPC (EPO FORM 1205N DATED 17/12/2019)

122 Ep: pct application non-entry in european phase

Ref document number: 17908083

Country of ref document: EP

Kind code of ref document: A1