Basic Course on finite element analysis (Zeng) Note One-two-dimensional rod element Finite element program (Python-based)

Zeng Teacher's "Basic course of finite element Analysis" The third chapter has the derivation of two-dimensional rod element, and combines an example to analyze and solve the program based on MATLAB. But I feel the MATLAB code in the book is a bit wordy, and some of the implementation methods are more troublesome, such as already know the starting point and the end point coordinates, still need to give the unit local coordinates and the overall coordinate angle, this is not necessary. So I reconstructed the program in Python, not to translate the MATLAB code in the book into Python (in fact it can be done completely, and soon!). ）。 For example, I use object-oriented thinking, the bar unit is written into a class, so the idea is clearer.

#!/usr/bin/python#-*-coding:utf-8-*-ImportMathImportNumPy as NP sqrt, cos, sin, pi=math.sqrt, Math.Cos, Math.sin, Math.PI"Pre-processing"Nodenumber= 4KK= Np.zeros (2*nodenumber,nodenumber))"""boundary condition, u represents the displacement vector of the node, if the displacement of a certain degree of freedom is unknown, then the place is filled with ' u_unknown ', F is the load vector of the node, if the load of a certain degree of freedom is unknown, the place is filled with ' F_unknown '"""U= Np.array ([0, 0,'U_unknown', 0,'U_unknown','U_unknown', 0, 0], dtype=object) F= Np.array (['F_unknown','F_unknown', 2e4,'F_unknown', 0, -2.5e4,'F_unknown','F_unknown'], dtype=object)classbar2d:"""defines the two-dimensional rod element class, which contains basic information about the member: E modulus of elasticity, A-bar element area, the node number of the starting point of the I unit, the node number of the end of the J element X1 The coordinates of the Y1 starting point, the coordinates of the X2 y2 end point, the position of the DOF element in the overall stiffness matrix The length of the element, the sine of the cosine direction of the cos sin element, the K element stiffness matrix"""    def __init__(self, E, A, x1, y1, x2, Y2, I, j): Self. E=E Self. A=A self.i=I self.j=J#defines an array of degrees of freedom of a single-rigid matrix converted to the total matrix of degrees of freedomSelf. DOF = Np.array ([2*i-2, 2*i-1, 2*j-2, 2*j-1],dtype=np.int16) self. L= sqrt ((x1-x2) **2 + (y1-y2) **2) Self.cos= (x2-x1)/Self . L Self.sin= (y2-y1)/Self . L L, C, s=Self . L, Self.cos, Self.sin self.k= (e*a/l) *np.array ([[C*c, C*s,-c*c,-c*S], [C*s, S*s,-c*s,-s*s], [-c*c,-c*s, C*c, c*s], [-c*s,-s*s, c*s, s*S]])"define a function to solve the element stress"    defStress (U):#obtaining the 1*4 displacement matrix of two nodes of the unit from the Displacement matrix UU =U (self. DOF) E, L, C, S=Self . E, self. L, SELF.C, Self.s T= Np.array ([-C,-S, c, S]) Self.bar_stress= e/l *Np.dot (T, u)returnself.bar_stress"defining a total matrix integration function"defbar2d2node_assembly (KK, bar): forN1inchXrange (4):         forN2inchXrange (4): Kk[bar. DOF[N1], bar. DOF[N2]]+=bar.k[n1, N2]returnKK'Solving node Displacements'defnode_disaplacement (KK, U, F):#get the reduced total matrixDel_row_col = Np.where (U = =0) Kk_delrow=Np.delete (KK, Del_row_col, 0) Kk_delcol= Np.delete (Kk_delrow, del_row_col,1) KK=Kk_delcol#get nodal force matrix of node displacement positionf = f[np.where (U = ='U_unknown')] U=Np.linalg.solve (KK, F) u[np.where (U=='U_unknown')] =ureturnU'To solve the nodal force, the function must be called after the node displacement U has been obtained.'defNode_force (KK, U): F=Np.dot (KK, U)returnF

Still taking examples from the book

The discretization of the structure is consistent with the book.

In the program bar2d this class represents the rod unit, when instantiated, need to provide information:

Elastic modulus E, area A, start I and end J numbers, and their respective coordinates.

So for this example, here are some information:

E, a, x1, y1, x2, y2, x3, Y3, x4, Y4 = 2.95e11, 0.0001, 0, 0, 0.4, 0, 0.4, 0.3, 0, 0.3= bar2d (E, A, x1, y1, x2, y 2, 1, 2= Bar2d (E, A, X3, y3, X2, Y2, 3, 2= Bar2d (E, A, x1, y1, X3, Y3, 1, 3= bar2d (E, A, X4, Y4 , X3, Y3, 4, 3= [Bar1, Bar2, Bar3, BAR4] for in bars:    bar2d2node_assembly (KK, BAR)

The corresponding function is then called to solve the displacement and load:

# Solving Displacements U = node_disaplacement (KK, U, F)#  solution node force F = Node_force (KK, u)

The final results are as follows:

It can be seen that the final result is consistent with the results obtained in the book, because the unit units are in M, so the displacement Unit is M.

--------------------------------------------------------------------------

Writing this code, reminds me of the time I went to Huihui teacher's "Finite element method" this course. Remember that the basic part of the previous chapters was taught by Mr. Jin, and the following chapters were for a small summer semester taught by Professor Sun Lizhi of the University of California, Irvine, who taught in English all the time. It's the only time I've ever listened to all the English lessons, though a little bit hard, but most of them can be understood. Foreign teacher's teaching ideas and domestic teachers still have some differences, foreign teachers generally from the basic attribute method, I remember that at that time sun teacher used a lot of space to talk about the weighted residual method and the Galerkin method. Although now forget almost, but it is an unforgettable experience, hot summer, Zhongshan Hospital, sweating Sun professor.

The teacher in the middle of the job, is a three-dimensional rod element structure of the example, let us write a program to solve the example, and then use the general finite element program analysis, compare the results. The programming language is unlimited, but all of us have chosen Matlab automatically. Originally I was forced lattice, I also want to use FORTRAN, but weigh a bit or forget, learn too troublesome. Finally, we chose Matlab. In retrospect, the code that was written was too bad. Using a day to get started with Matlab, and then gig out, in order to declare the variable with the tag of the use of N multi-line eval function, unbearable.

Basic Course on finite element analysis (Zeng) Note One-two-dimensional rod element Finite element program (Python-based)