§ 4 Units to be coordinated


1.        Coordination issues and units to be coordinated


One disadvantage of isoparametric elements is that no matter how the number of nodes and the number of interpolation increases, only the continuity of the interpolation function itself can be guaranteed between each element. For spatial elastic mechanics problems, the strain in the deformation energy is the first-order partial derivative of the displacement), but for the variational problems involving higher-order derivatives (such as the bending problem of rods and plates, the integral of the deformation energy contains undetermined The second-order partial derivative of the deflection function), the interpolation function often requires a first-order continuity in the entire region. However, such so-called coordination problems cannot be solved by equal-parameter units. Of course, this requirement of smoothness must first be reflected at the nodes, that is, the value of the node parameter of the element should contain the relevant derivative value of the undetermined function. Such units may be collectively referred to as quasi-coordination units.

[ Type function ] Suppose  there are p nodes on an element with local numbers i = 1,2, , p . First consider the parameter value of the undetermined function u = u ( x , y , z ) and its first-order partial derivative at node i

The number of parameter values ​​for each node is r = 4 . So there are 4 p node parameter values ​​in each unit. If the interpolation function is a ternary nth degree polynomial , the number of terms or coefficients of each term is generally one, and obviously the degree n must satisfy


At this time, corresponding to the parameter values ​​of each node , the definable function is the following 4 p nth degree polynomials ( i = 1, 2 , , p ):

(i) at node i

(ii) At the remaining p - 1 nodes ji , the above sixteen functions and their first-order partial derivatives are all equal to zero.


If these type functions can be constructed, then the interpolation polynomial of u can be expressed as

If the undetermined function is ( u , v , w ) and other conditions are the same as before, its interpolation function ( ) can be expressed as

or as a vector


[ Undetermined coefficient method ] The undetermined coefficient method is to directly solve the coefficients of the interpolation polynomials from the 4p node parameter values  . Generally speaking, the number N = of complete n -degree polynomial coefficients { } will be greater than 4 p , and it is necessary to add N - 4 p conditions to the polynomial to uniquely determine these coefficients. The easiest way is to restrict the form of the polynomial as in § 3 (eg restrict the polynomial to a polynomial of degree k , that is , for x , y or z , the number of coefficients is reduced to ) or remove some The higher-order terms (that is, setting the corresponding coefficients to zero) make the number of coefficients 4 p , assuming that the restricted interpolation polynomial is rewritten as

4p equations can be obtained from 4p node parameter values

   ( i = 1 , 2 , , p )

In the formula , respectively represent the column vector composed of the polynomial terms (the number of terms is 4 p ) and the value of their partial derivatives at node i . Arrange all the node parameter values ​​of this unit into a column vector , and use U to represent the 4 p × 4 p coefficient matrix formed by the sequential arrangement of p 4 × 4 p coefficient matrices on the right side of the above equation , then the 4 p equations can be abbreviated as

From this, the coefficients of the terms can be solved

                                ( 15 )

It can be seen that, under the same restriction, the coefficients of each type of function are actually the solution in which one component takes 1 and the other components take zero.

Another way is to have a general expression for the interpolation function

Add certain constraints according to smoothness requirements or physical conditions, and assume that these constraints can be expressed as N - 4 p linear equations about

                                   ( 16 )

where Q is a ( N - 4 p ) × N coefficient matrix. For this complete n - th interpolation polynomial, the following 4p equations can also be derived from the 4p node parameter values

                                ( 17)

In the formula, the coefficient matrix of 4 p × N is represented , as long as the constraints are selected appropriately * , the general solution of the simultaneous equations ( 16 ) ( 17 ) can be obtained


Divide an N × N matrix G into two submatrices with first 4 p columns and last ( N - 4 p ) columns: G = [ ] . but

(15) or ( 17 ) represent the correspondence between the node parameter values ​​of u and the coefficients of the interpolation polynomial, and there are similar relational expressions for the undetermined functions v and w . If written in vector form, the interpolation polynomial is

                        ( 18 )

[ Generalized nodal parameters ] If the interpolation polynomial  ( 18 ) is directly used instead of ( 1 ) for element analysis , the variational equation ( 2 ) should be changed to

     In the formula, the terms { R } are equivalent to the type function , and the coefficient { a } plays the role of the node parameter value , which is called the generalized node parameter of the element . The element coefficient matrix and the equivalent external force after the ampersand in the above formula are also generalized , respectively denoted as

So there is

where A is the matrix in ( 15 ) . Similar results are obtained for interpolation polynomials .

[ Transformation of node parameter values ] The node parameter values  ​​contain partial derivatives , and their values ​​in the Cartesian coordinate system and the local coordinate system are different . For example , the node parameter values ​​in the Cartesian coordinate system

in the local coordinate system as

                    ( 20 )

In the formula, the right subscript i of the fourth-order matrix represents the value at node i . On the contrary

Note that the type function in the local coordinate system is, and the type function in the rectangular coordinate system is , from the expression

It can be seen that they have the following linear relationship

                                  ( 21 )

In fact , all parts that contain differential operations , such as the matrix B containing differential operations in the integral formula of , the directional derivative of the boundary, etc. , must be transformed accordingly under the coordinate transformation .


2.        High-order interpolation of one-dimensional elements


[ Cubic interpolation ]   The univariate cubic polynomial has four terms , and its coefficients can generally be determined by four node parameter values . Now take the two endpoints of the line segment as nodes 1 and 2 , and the node values ​​of the undetermined function and its derivatives are four node parameter values , namely

        ( i = 1 , 2 )

Take the distance coordinate ( ) as the local coordinate , so that the coordinates of the two endpoints of the line element are respectively = 0 , = 1 . Define a cubic polynomial that satisfies the following conditions as a type function :

(i)      on node i              

(ii)     on node j ( i )         


Using the symmetry of the distance coordinates from these conditions can be determined respectively


So in the local coordinate system , the interpolation polynomial can be expressed as

If you want to represent it in a rectangular coordinate system , you can use coordinate transformation

                              ( 22 )

Substitute into the above formula and expand it to get the following Hermitian cubic interpolation polynomial

[ Quatic Interpolation ]   The quintic polynomial in one variable has six terms , and its coefficients can generally be determined by the six node parameter values . Now take the node value of the undetermined function and its first and second derivatives as the node parameter value , that is

   ( i= 1,2 )

Also take the distance coordinate ( ) as the local coordinate , and define the fifth degree polynomial satisfying the following conditions as the type function :

(i)    on node i  

    (ii)   At node j (i ) its first and second derivatives are equal to zero .


From these conditions and using the symmetry of the distance coordinates, it is possible to determine

So in the local coordinate system , the interpolation polynomial can be expressed as

If it is to be represented in the rectangular coordinate system , it can also be substituted and expanded by the coordinate transformation ( 22 ) to obtain the following Hermitian quintic polynomial

The former is because the derivative has the following relationship in different coordinate systems :

[ Type function and undetermined coefficient method ]   Take the above quintic interpolation polynomial as an example . Write the general form of the quintic interpolation function

The six undetermined coefficients can be determined by the six boundary conditions

       ( i= 1 , 2 )

To be determined , it can be reduced to a set of linear equations about , whose constant terms are the values ​​of these nodal parameters .

The type function is nothing but a special interpolation function that takes a unit vector . For example, in local coordinates , taking

Then the six boundary conditions can be expressed as

Solved from this

So this special interpolation function can be written as

That is, the type function of the previous paragraph, the other type functions in the local coordinate system can also be calculated similarly .

For the one-dimensional case , it is relatively easy to use the undetermined coefficient method to find the type function . But for the two-dimensional case , especially for incomplete high-order interpolation , the definition and composition of the type function are often not so easy due to the addition of constraints. Simple . In order to avoid the direct composition of the type function , the method of generalized node parameters can be used .


3.        Higher-order interpolation of triangular elements


The local coordinate system takes the area coordinates , and makes ( ) = ( ), the definition and formation of the type function are only carried out in the local coordinate system.

[ Quadratic Interpolation ] A total of six items of   binary ( ) quadratic polynomial or ternary ( ) quadratic homogeneous . Therefore, six node parameter values ​​are required to achieve complete quadratic interpolation . In addition to giving function values ​​at the vertices of the triangle outside , its normal derivative values ​​can be given at the midpoints of the three sides , here denoting the outside normal on side = 0 ( Fig . 19.14 ).

This is the simplest quasi-coordinated plate element , which takes a parameter value at each of the six nodes , and its shape function can be defined as follows in the local coordinate system :

There are six functions with satisfying conditions :

(i)      on vertex i =1 , =0 ,( i =1 , 2,3 )  

on the midpoint of edge =0

(ii)     On the remaining relevant nodes , , the external normal derivatives of , are equal to zero .  


Now let's find the expression of the normal derivative outside each side . Since ( ) = ( ), and

In the formula, the length of the line segment with edge = 0 , that is, the distance between the two vertices j , k

From the inverse transformation matrix of the § 2 triangular element, the value of the outer normal derivative on each side can be obtained


where A is the area of ​​the triangle . For example , on the =0 side


which represents the height of the bottom edge . The result can be written uniformly as







注意,由于这些节点参数值的定义与坐标系无关,因此上述的型函数在局部的面积坐标和整体的直角坐标系中是一致的,   (i=1,2,3).

[三次插值]  三次的二元()多项式或三元()的齐次式共十项,因此需要十个节点参数值,才能做到完全的三次插值,对每个顶点i(i=1,2,3),取节点参数值为,共九个,其余一个或取在形心O的函数值,或改为一个限制条件.对后一情况,可采用同九节点等参数单元一样的限制(§3),即要求对于三元二次齐次多项式是完全的,对于()的齐次式



For the case of ten-node parameter values , the definable type functions ( i = 1 , 2 , 3 ) are as follows :

    (i)   on the centroid O ( i = 1 , 2 , 3 )               

    (ii)   on node i


    (iii)   On the remaining nodes j (i ) , its first-order partial derivatives are all equal to zero .


    The shape function of local coordinates can be obtained by using the method of undetermined coefficients



The interpolation polynomial of u can be written as

Turn to the Cartesian coordinate system , since the parameter value of the node remains unchanged , then press the formula ( 21 ) ( two-dimensional case ) to obtain its type function

( i =1 , 2 , 3 )

In the case of adding restrictions on the parameter values ​​of nine nodes , remove the centroid node and the shape function definition . As for its composition, it can follow the method of nine nodes and other parameter units, add an item to each of the above shape functions , and then request Satisfy the restriction condition ( 24 ), and thus determine a . As a result , nine types of functions can be obtained


The relationship between them and the type function of Cartesian coordinates and the expression of the interpolation polynomial are the same as in the previous case , as long as they are removed .

Note that these two elements are not coordinated . Although they are suitable for plane elastic problems , they also need to meet the criteria for convergence of the patch test to be used as plate elements . Therefore , the so-called parallel element is also added to the three-sided element with nine nodes. The restriction of triangulation ( that is, the three sides of each element in the region are parallel to three fixed directions ) can be used as a plate element . Of course, this restriction can also be relaxed to as long as parallel triangulation is performed in each sub-region of Ω . _

[ Quadic interpolation ]   A complete binary quintic polynomial or a ternary quintic homogeneous has a total of 10 coefficients , so 21 interpolation conditions need to be given to the element . First, the node parameter values ​​up to the second derivative are given at each vertex :

        ( i =1 , 2 , 3 )

There are 18 in total , and the other three conditions can be given on the three sides respectively , so that the derivative of the interpolation function remains continuous between the elements . The outer normal derivative of the interpolation polynomial on the boundary is generally a quartic polynomial , which should be determined by five independent The condition of the node is uniquely determined , and the node parameter value only provides the value at both ends , that is, there are only four independent conditions . In order to maintain continuity between the elements , restrictions can be placed on each edge . The general restrictions are as follows Two kinds :

(i)      requires a cubic polynomial at the boundary , i.e. the coefficient of its quartic term is equal to zero .  

(ii)     The value at the midpoint of the boundary is a node parameter value ( Fig . 19.16 ).  

In this way, three equations are obtained on the three sides , and together with the original 18 equations, the interpolation polynomial can be uniquely determined . The generalized node parameter method in the area coordinate system is introduced below .

Let ( ) 's fifth-order homogeneous terms be arranged as follows :


The corresponding coefficient is

And the interpolation polynomial is

Using formula ( 23 ) to find the outer normal derivative of the edge , we can get

According to the method (i), substituting into the above formula and the coefficient of the essential term is zero , then we get

in the formula

Considering the restriction of the remaining two sides = 0 , = 0 , the other two equations can be obtained in the same way . After simplification , the three equations can be written in the following form

in the formula

For the mode (ii), we can make . Substitute into the above expression to get its value at the midpoint ; the value at the midpoint of the other two boundaries can be obtained similarly . The result is the relationship between the parameter values ​​and the coefficients of the three midpoint nodes. relation can be written as


It is not difficult to list the corresponding 3 × 21 matrix Q according to the notation of ( 16 ) and ( 17 ) , and the 18 × 21 matrix corresponding to the parameter value of 18 nodes can be written as

where is a 6 × 6 matrix :


Because of its inverse matrix


The last three lines of , are easy to obtain for either method : first solve and then substitute into the latter three equations to solve . Then the linear relationship between the coefficients of the complete quintic polynomial and equivalence is obtained :

The left end is the generalized node parameter , but the { b } on the right end of the method ( ii ) is not the node parameter value , they should be changed to the value at the midpoint of the three sides . It can be seen from ( 25 ) that the last three columns of G have to be multiplied by factor

    Note that the transformation of the node parameter values ​​of the two coordinate systems is no longer ( 20 ), but to be changed to

The external normal derivative of the midpoint is independent of the coordinate system .


Four,   quadrilateral unit


[ Bicubic interpolation ( coordination plate element )]   The bicubic interpolation polynomial has a total of coefficients , and its general form can be written as

These coefficients can be uniquely determined as long as the node parameter values ​​are given at each vertex .

Considering the identity matrix in the local coordinates ( ξ , η ) and making appropriate combinations of the form function expressions of the cubic interpolation of one-dimensional units , it is not difficult to see that the form functions can be written in sequence as


It is easy to prove that they satisfy the following conditions :

(i)      equals 1 on vertex i .  

(ii)     The remaining functions and their first and mixed second derivatives are equal to zero at each node .  


So in the local coordinate system , the interpolation polynomial can be written as

Note that in two coordinate systems , the transformation between node parameter values ​​should take

And the type function correspondingly has the following relationship

                         ( i =1 , 2 , 3 , 4 )

This unit is the coordination unit .

[ Incomplete Bicubic Interpolation ] The   node parameter value is taken , and there are twelve interpolation conditions . Therefore, to uniquely determine the interpolation polynomial , it must be restricted . Generally, the following high-order items are removed from the above-mentioned complete bicubic polynomial , that is, and other four terms , while still maintaining symmetry . Note that although this interpolation function is missing four high-order terms , it still contains a complete cubic polynomial , but because the mixed derivative is not taken as nodal value , along the boundary between the two elements The normal derivative cannot guarantee its continuity. Nonetheless, it satisfies the convergence criterion of the patch test and thus can be used as the deflection mode of the plate element .

Using the undetermined coefficient method , the shape function corresponding to the local coordinate ( ) can be obtained as


In the local coordinate system , the interpolation polynomial can be written as


* The case of one or two should be changed toor.Ifspecified as a bi-degree polynomial,takeand so on.

* Notethatconjunctionwiththesubdivision of the element:someconstraints do not apply to acertainelement.The partial derivative is the node parameter value,then,a condition needs to be added.Inorder to converge,a complete quadratic polynomial should be included,anddue to symmetry,it can generally be madeto thesecond coefficient,but this constraint is for parallel coordinate axesThe isosceles right-angled triangular elementof ,it can be argued that the corresponding matrixis ​​singular.

* If a finite element solution with a certain interpolation function can obtain a solution in a constant "strain" state on any unit slice (composed of several adjacent elements), the interpolation function is said to pass the slice test, and Convergence is guaranteed. For the thin plate bending problem, if the boundary conditions corresponding to the complete quadratic polynomial are given on the boundary of any element, and the deflection is obtained by solving the result of some interpolation function, it is said to pass the slice test.

Original text