有限单元法及程序设计（Finite element method and program design）

1、有限单元法及程序设计（Finite element method and program design）Finite element method and program designIntroduction1. mechanical analysis method: analytic method, numerical methodThe finite element method, the actual structure shape and the load is complex, mostly by analytical method and numerical method is d

2、ifficult to develop, with the progress of computer, a new numerical analysis method which is developed.2 basic steps:(1) structural discretization: the structure is divided into finite units by lines or planes from the set.(2) element analysis: the relation between nodal displacement and nodal force

3、 (element stiffness matrix) of the derived element.(3) overall analysis: the whole structure of each element is analyzed, and the relation between the displacement of the structural point and the nodal force is derived.3 programming steps:Ask questions and work out a solutionStructural mathematical

4、modelDraw a program flow chartProgrammingCompiler debuggerTrial verification program4., according to the national standard (GB-1526-89) procedures, processes, icons, quasi symbols and regulations:EMBED PBrushThe graph represents the starting and ending points of the program flow chart;Diagrams show

5、the input and output of data information;A data set to be completed before a series of operations is performed by a graph;The chart shows the criteria;Graphs represent various processing functions, such as mathematical operations;The chart shows the path and direction of the process.Finite element m

6、ethod and program design of bar structuresThe finite element method of plane member elementBasic concepts of finite element methodThe basic idea: first after (the first unit analysis, then the overall analysis)Basic concept: whole number: node endpoint number by natural number 1, 2, 3,. (in the over

7、all coordinate system xOy)Local number: each unit is marked with I, J (in the local coordinates of the unit, EMBED, Equation.3, in the same direction as the overall coordinate system).3.Fe=, Ke, Delta, e, where: Ke = EMBED, Equation.3, element stiffness matrix, elements are stiffness coefficientsDel

8、ta e=, EMBED, Equation.3 element, rod section, displacement arrayRod end force array of Fe= EMBED Equation.3 unitK, EMBED, Equation.3, =P (1-7)K=, EMBED, Equation.3, integral stiffness matrix, EMBED, Equation.3 = EMBED, Equation.3, EMBED, Equation.3 displacement arrayP=, EMBED, Equation.3 node load

9、array3. finite element displacement analysis of continuous beams needs to be consideredStiffness integration method:(1-3) the K expansion order and the expanding element are 0, and the element contribution matrix is obtainedUnit 1: K = EMBED Equation.3 Unit II: K = = EMBED Equation.3The element cont

10、ribution matrix is superimposed to form the whole stiffness matrixK=, K + K = EMBED Equation.3The introduction of support conditions at both endsWithout considering the constraints, the overall stiffness matrix is obtained, and the main diagonal element K, EMBED, Equation.3 is changed to 1, line I,

11、column J, and the rest elements are changed to 0,The corresponding load element is also changed to 0.Treatment of non nodal loadsThe analysis is carried out by equivalent nodal load:Each node (including two nodes) is constrained to stop the rotation of the node, and the constraint moments are the su

12、m of the fixed end moments of the related units which are handed over to the node, and the clockwise is positiveRemove an additional constraint (P EMBED, Equation.3, the same as the constraint moment, in the opposite direction)The two part of the rod end bending moment is superimposedElement stiffne

13、ss matrix in second local coordinate systemsGeneral unitThe elastic modulus, the sectional moment of inertia and the sectional area of the element EQ, oac (L, e) are E, I and A, respectively. The I and j ends of the unit have three bars, end forces, EMBED, Equation.3 (i.e. axial forces, shear and be

14、nding moments) and their corresponding three bar end displacements, EMBED, Equation.3, as shown in figure 1-7. Figure EMBED Equation.3 as the element local coordinate, I point at the origin of EMBED, Equation.3 axis and the rod axis coincide, provided by I to j for EMBED positive Equation.3 axis, Eq

15、uation.3 axis by EMBED clockwise 90? EMBED Equation.3 direction. The positive direction of force and displacement is shown in figure 1-7.EMBED PBrushIn this unit, the element end force array and the rod end displacement array are respectivelyEMBED Equation.3 element rod end force arrayEMBED Equation

16、.3 rod end displacement arrayIn order to derive the relationship between the end force of a general element and the displacement at the end of a rod, we consider the following two cases respectively.Firstly, the relationship between axial force EMBED Equation.3 and axial displacement EMBED Equation.

17、3 of two rod ends is analyzed. According to Hookes law, yesEMBED Equation.3 (a)Secondly, the relationship between rod end bending moment EMBED Equation.3 and rod end shear EMBED Equation.3 and rod end angle EMBED Equation.3 and lateral displacement EMBED Equation.3 are considered. According to the r

18、otation displacement equation of the displacement method of structural mechanics, and according to the symbols and the sign and sign given in this section, theEMBED Equation.3 (b)Combine (a), (b) two, and form the matrix, as followsEMBED, Equation.3, EMBED, Equation.3 (1-17)The upper form can be abb

19、reviated as EMBED Equation.3 (1-18)The element stiffness matrix isEMBED Equation.3 (1-19)Properties of element stiffness matrixEach element represents the rod end force caused by the displacement of the rod end of the unit. The physical meaning of any element EMBED, Equation.3 (R, s, 1 to 6) is the

20、value of the force at the end of the rod at the end of the rod at the end of the s when the displacement component is equal to 1 RIs the symmetric matrix, whose element is EMBED Equation.3Singular matrix, whose element determinant equals zero, that is, EMBED Equation.3Characterized by rapid differen

21、tiationAxial force element: a unit that considers only the end displacement of the axial member and the force at the end of the rodCoordinate transformation of element stiffness matrixSet up the element stiffness matrix and element stiffness equation is the local coordinate system in EMBED Equation.

22、3, the general structure of members, each unit is divided by the analysis of the local coordinate system is obviously different. Therefore, when we study the condition of structural equilibrium and the condition of continuous deformation, we must select a unified coordinate system xOy, which is call

23、ed the whole coordinate system. At the same time, the element stiffness matrix established in the local coordinate system must be transformed into the Dan Yuangang matrix of the whole coordinate system.Fig. 1-8a) and Fig. 1-8b) represent the rod end force components of the E, EQ, oac, and xOy in the

24、 local coordinate system EMBED, Equation.3 and the whole coordinate system, respectively.In order to derive the relationship between the rod end forces Xi, Yi, Mi and the EMBED Equation.3 in the local coordinate system, Xi and Yi are projected onto the EMBED Equation.3 axis respectivelyEMBED, Equati

25、on.3, a)In the formula, an alpha represents the angle between the X axis and the EMBED Equation.3 axis, which is clockwise.EMBED, PBrush, figure 1-8In the two coordinate system, couple components unchanged, i.e.EMBED, Equation.3, b)In the same way, the end force of the j end of the EQ oac (E) is ava

26、ilableEMBED, Equation.3, c)Combine (a), b), c) and express in matrix form, availableEMBED Equation.3 (1-24)This formula is the transformation formula of the end force of the unit rod in the two coordinate systemsEMBED Equation.3 (1-25)Equation: the element end force array in the EMBED Equation.3 loc

27、al coordinate systemElement end force array in EMBED Equation.3 integral coordinate systemEMBED Equation.3 unit coordinate transformation matrix (1-26)T is an orthogonal matrix whose inverse matrix is equal to its transpose matrixThat is, EMBED Equation.3Unit unknown quantity codingIn order to facil

28、itate the programming calculation, it is necessary to count the displacement components of the nodes according to certain rules. The node displacement of the structure has the difference between the free node displacement and the bearing node displacement (also called the bearing node displacement).

29、 The displacement of the free node is unknown. There are two kinds of ways to solve the displacement of unknown nodes: the pre process method and the post-processing method.EMBED PBrushWhen the structure is analyzed by post-processing, the displacements of all points are unknown, and the nodal displ

30、acement array is (see Figure 1-9)EMBED Equation.3Pix, Piy and Pi respectively represent the theta role in node i (i=1, 2, 3, 4) on the horizontal force, vertical force and moment. It is stipulated that the positive direction of the nodal force Pix and Piy is the same as that of the overall coordinat

31、e system X and y, and that Pi theta is positive at clockwise direction; the positive direction of the nodal displacement is consistent with the positive direction of the nodal force.After the stiffness equation of each element is obtained, the displacement equation of the whole structure can be esta

32、blished according to the equilibrium condition of the node and the continuity condition of the displacementEMBED Equation.3 (1-37)Or EMBED Equation.3 (1-38)amongEMBED Equation.3 (1-39)Is the overall stiffness matrix of the structure, or the original stiffness matrix of the structure.In the establish

33、ment of the overall stiffness equation (1-38), assuming that all point displacement is unknown, a whole structure of ignorance, and under the action of external force, in addition to the elastic deformation, may also occur at this time, the rigid body displacement, node displacement is not uniquely

34、determined. This formula (1-39) is a singular matrix, and can not be reversed. Therefore, (1-38) it is not possible to find the displacement of the node.In the rigid frame shown in Fig. 2-9a, the nodes 1 and 4 are the fixed points, so the nodal displacement is known and the support condition is zero

35、. The support condition is introduced into the overall stiffness equationEMBED Equation.3 (1-40)Can be divided into two groups of equations, a set ofEMBED Equation.3 (1-41)It can find the unknown node displacement, delta 2, Delta 3. Make a group ofEMBED Equation.3 (1-42)Reaction equation. Using the

36、formula (1-41) to calculate the nodal displacement 2 and 3, and substituting the upper formula, the unknown support reaction force can be calculated.For the general bar structure, these steps can be analyzed. No matter how many nodes structure displacement, after adjusting the order of the total, it

37、 can be divided into two groups: one group includes all of the unknown node displacement component to a F representation; another group to support node displacement component in Delta R. Accordingly, all nodes are divided into two groups, and corresponding to the delta F is known as the nodal force

38、array, represented by PF; and corresponding to the delta R is the bearing node force array, expressed in PR. HenceEMBED Equation.3In conjunction with the above analysis method, the elements of the overall stiffness matrix K0 are rearranged, and K0, Delta 0=P0 can be writtenEMBED Equation.3 (1-43)Exp

39、anded upperEMBED, Equation.3 (1-44), style (1-45)Formula (1-44) is the modified equation of overall stiffness, which differs from the formula (1-38) in that the supporting conditions are introduced.Post process method: the whole stiffness matrix is formed by the element stiffness matrix and the stif

40、fness equation is introduced. Then, the supporting position is introduced and the position displacement of the node is solvedPre process: only the unknown free node displacement component number is obtained, and the node displacement in the array does not contain the known constraint node displaceme

41、nt componentThe rigid frame section with a composite junction is shown in Figure 1-10For three units, the numbers are 1, 2, 3, and so onThe arrows between the poles represent the direct direction of the local coordinate system,The frame structure number is 15. Consider the following unitsNodal displ

42、acement component number. The first treatment method is adoptedAs follows:The independent displacement components are numbered by natural numbers only,Displacement number. If some displacement components are connectedConditions or direct bar axial rigidity conditions are restricted to one anotherEqu

43、al，The number is the same.At the support, some displacement components are zero due to rigid constraints, and the number of displacements is zero.Therefore, figure 1-10 is numbered as followsThe number of rod unit number node module unit number beginning end displacement components AB 12000123 BC 2

44、23123456 DC in 54000457 in the computer program, the two ends of the element node number using two-dimensional array JE (I, e), called the two ends of the element node number array.JE (1, e) = the node number at the beginning of the cellJE (2, e) = the node number at the end of the unitIn this caseJ

45、E (1, 1) =1, JE (2, 1) =2JE (1, 2) =2, JE (2, 2) =3JE (1, 3) =5, JE (2, 3) =4Displacement number of displacement component of any node can be represented by two-dimensional array JN (I, J). It is called the node displacement number arrayDisplacement number of JN (1, J) = node j along X directionDisp

46、lacement number of JN (2, J) = node j along Y directionDisplacement number of JN (2, J) = node j angular displacementIn this case, for the third node,JN (1, 3) = 4JN (2, 3) = 5JN (3, 3) = 6The EQ oac unit (E,) at the beginning and end of displacement, a row (beginning in front), the digital element

47、orientation array, which can be used to form the coordination condition of displacement and displacement between corresponding nodes conveniently. Its expansion type isEMBED Equation.3In the formula, the D elements M1 to MD are the displacement numbers corresponding to the displacement components at

48、 both ends of element EQ, oac (E). In this case,EMBED Equation.3The overall stiffness matrix of plane structuresAfter the element analysis and the element stiffness matrix, the whole analysis is needed. The first approach as an example, the discrete element to form the original structure, so as to s

49、atisfy the equilibrium conditions of displacement of the structure node continuity conditions and force, so as to obtain the stiffness equation of modified structure is given earlier type (1-44)KFF F=PFType: KFF called the whole stiffness matrix of the modified structure; F, PF were free node displa

50、cement and free nodal load vector. When it is known not to free node displacement caused by misunderstanding, type (1-44) is often referred to as the overall stiffness equation of KFF, referred to as the overall stiffness matrix, F and PF were referred to as node displacement and load vector array.

51、In order to write conveniently, the following table is often omitted.A stiffness matrix integrated overall stiffness matrix, usually by direct stiffness method, the calculation process is divided into two steps, first calculate the contribution of each element matrix, and then put them together, it

52、is concluded that the whole stiffness matrix. However, in the actual calculation, the inconvenience, reason is the calculation needs first contribution matrix Ke of all units are saved, and the overall stiffness matrix of order K Ke, which takes up a lot of storage capacity, so the actual operation,

53、 the edge location, advancesMethod. The principle has not changed, and the result is the same.EMBED PBrushIn Figure 1-10, the cell matrix of unit 1 isEMBED Equation.3The element displacement component number is marked on the upper and the right of the element stiffness matrix in the formula. Because

54、 the overall stiffness matrix elements are arranged according to displacement component number, according to the treatment method, the element stiffness matrix corresponding to the component number of zero elements into the overall stiffness matrix, non-zero numbers indicate the remaining elements i

55、n the middle of the overall stiffness of the moment of the row and column number.thereforeEMBED, Equation.3, EMBED, Equation.3, EMBED, Equation.3The stiffness matrix of the unit is.EMBED Equation.3The location of each unit in K isEMBED Equation.3The stiffness matrix of the element isEMBED Equation.3

56、The position of each element in K isEMBED, Equation.3, EMBED, Equation.3, EMBED, Equation.3According to the above positioning method, the stiffness matrix of the three elements and the corresponding stiffness matrix of the whole element are obtainedEMBED Equation.3When electromechanical computing is

57、 adopted,Set K to zero, this is K=0 EMBED Equation.3;The K EQ EQ EQ S (1) related elements, in accordance with the seats, up to K;K EQ S (II) related elements, continue to follow the seats, up to K;K EQ S (3) related elements, continue to follow the seats, up to K, the overall stiffness matrix finally completed.Find the overall stiffness matrix of the rigid frame shown in Figure 1-12 K.Each rod has the same size.A=0.5m2 I