Programming for solving plane rigid frame based on M ATL AB

. Based on the idea of the matrix displacement method, this paper designs a program which can be used to solve the internal force of the continuous beam and rigid frame with MATLAB. It mainly demonstrates how to design a program to realize the matrix displacement method with MATLAB. In addition, some techniques are included in order to realize the correspondence between the manual calculation and the computer calculation, such as “Using lambda to locate”, “Crossing out rows and columns” and visual design. Therefore, based on the structural mechanics, combined with the principle of matrix displacement method, this paper shows the whole process from inputting the information of the rigid frame to solving the internal force of the rigid frame to outputting the bending moment diagram using MATLAB as the programming tool.


Introduction
Since the publication of Mathematical Principles of Nature Philosophy, Newton's law of motion and law of universal gravitation have been systematically established and elaborated. After that, structural mechanics has gradually developed under the framework system of classical mechanics. In the development history of structural mechanics, to solve the statically indeterminate structure, the force method was established by Maxwell in 1864. Though with its easy rationale, the choice of basic system requires the engineers' extensive experience and human intervention. At the beginning of the 20th century, the displacement method [1][2][3][4], which is more advantageous than the force method, has been developed to solve high-order statically indeterminate structures. However, with the increase of basic unknowns, the manual method of the displacement method is still unbearable. In order to solve this problem, in 1930 Hardy Cross invented the progressive analysis dethronement-Moment Distribution Method to simplify the calculation, which brought great convenience to the engineering calculation at that time. Then, with the development of the computer, the progressive calculation was replaced by the matrix displacement method, which is more suitable for programming. The matrix displacement method greatly enhances the efficiency of solving the high-order statically indeterminate structure and liberates the designer from the tedious calculation to a great extent. At the same time, in 1943, the mathematician Courant put forward the prototype of the finite element method. Subsequently, with the development of Aerospace Engineering, the finite element method gradually formed a complete system. Besides, because the matrix displacement method is the elementary form of the finite element method [5][6][7][8]. The matrix displacement method is also called the member finite element method.
Using MATLAB to program a matrix displacement method program can give full play to its advantages of matrix operation and can be applied to the force solution of a bridge or a building structure. Although there is complete finite element calculation software on the market, the internal finite element calculation steps are "black boxes" for users between the input and output. Therefore, for some personalized problems which are not in the scope of the general program, it still needs to be studied [9].
Therefore, based on the structural mechanics, combined with the principle of matrix displacement method [10], this paper shows the whole process from inputting the information of the rigid frame to solving the internal force of a rigid frame to outputting the bending moment diagram using MATLAB as the programming tool.

Information input
First of all, some basic variables need to be involved in the operation of the program. In the input information function -"TypeInInformation" function, it contains the basic information.
Among them, input the number of members that make up the rigid frame, and then input the information of each member in turn, including "length"-length of current member, "EI"-bending rigidity of current member, "EA"-tensile rigidity of current member, "alpha"-rotation angle of current member (define alpha in radian system; horizontal right as 0; clockwise as positive), "loadtype"-type of load of current member (load type "0"-No load; "1"-Uniform load (full load); "2"-Concentrated load (the load is considered to be applied on the 1/2 of the member ); "loadvalue"-value of load (if uniform load-Unit: M/N; if concentrated load-Unit: N).
After inputting the basic information, the function calculates the "EAlength" (EA/length), "Fp0"( )equivalent node load vector. (Since one member has two nodes and each node has three force components (N, V, M), "Fp0" is a vector with six rows and one column). "Fp0" is calculated according to the "loadtype" and "loadvalue". "1" represents uniform load, with left node N = 0, V = − . In addition, introduce "P0"( ) -a new variable that inverse sign of "Fp0"( ) for clear expression purpose.

Element stiffness matrix in local coordinate system
As mentioned above, since a member has two nodes and each node has three force components (N, V, M), " " is a vector with six rows and one column. These six force components correspond to six displacement components, that is, the displacement at the nodes of a member is also a vector of six rows and one column. According to the displacement method, it can be deduced that there is a relationship between displacement and force: = . is the element stiffness matrix in the local coordinate system. Obviously, it's a matrix of six rows and six columns.
While in the program, the cell matrix "k0=cell(n,1)" (corresponding , n rows, and 1 column, each row is a matrix, n is the member number), "P0=cell(n,1)" (corresponding ) is defined to represent the element stiffness matrix and equivalent node load of each member in the local coordinate system. Similarly, all variables (including "length", "EI", "EA", "alpha", "loadtype", "loadvalue", "EAlength") that need to be inputed in the "TypeInInformation" function are defined as "zeros (n, 1)" (n rows and 1 column vector, i.e. n members, each member stores one information above).
After the definition is completed, loop the "TypeInInformation" function n times, which means, input the above basic information of n members.
According to structural mechanics, is determined by the information above. So after the input is completed, loop n times, make each (k0) =

Element stiffness matrix in global coordinate system
However, as a rigid frame, each member needs to be integrated together, which requires a global coordinate system, so these stiffness matrix need to be transformed into the global coordinate system. In order to distinguish, the element stiffness matrix in the local coordinate system is represented by (k0 in the program) and the element stiffness matrix in the global coordinate system is represented by (ke in the program). The rules applied to the definition of (P0) and (Pe) are the same. According to structural mechanics, element coordinate transformation matrix T (T is already known by mathematical derivation) is used to transform stiffness matrix and equivalent node load from local( / ) to global( / ), so that = (1) = − = ∆ (2) the element stiffness matrix in global coordinate system is formed.
While in the program, the cell matrix "ke=cell (n, 1)"(corresponding to , n rows and 1 column, each row is a matrix, n is the number of members), "P0 = cell (n, 1)" (corresponding to -, n rows and 1 column, each row is a vector), "T = cell (n, 1)" (element coordinate transformation matrix of each member, T is related to the "alpha") are defined. The angle is defined as "alpha=zeros (n, 1)". After these preparation, the element stiffness matrix in the global coordinate system can be obtained according to the above formula.

Global stiffness matrix in global coordinate system K
In the traditional displacement method, the influence of each node displacement( ∆ ) on all node force(F) is considered respectively, and then the node force of each node are stacked, and then the relationship( = ∆) is formed. However, the "Element Integration Method" is used to consider the individual contribution of one element (or one member) to "F" at one time, and then stack them. That is to say, each element(or member) only controls six nodal force at its two nodes, which corresponds, the element stiffness matrix in global coordinate system (6 × 6 matrix). And ultimately we need a global stiffness matrix K( 3(n + 1) × 3(n + 1) matrix). So we need to put into K each time (attention, we need to follow certain rules, that is, use "lambda" to locate), and record it as (also a 3(n + 1) × 3(n + 1) matrix, and except for the place for , which is a 6 × 6 matrix, the rest of is 0). Then K=∑ , and that is called the "element integration method" to form K.
Besides, the purpose of location vector "lambda" is that because the previous "loop n times" refers to "n" members, but there are a total of "n+1" nodes, that correspond, "n+1" equivalent node force vectors. Therefore, we need to use " lambda " to map the contribution of each member to .

The formation of ∆
Based on the displacement method of structural mechanics, in order to solve statically indeterminate structures, the basic equations of displacement method are needed Here, K and P are the global stiffness matrix and the equivalent node load. Theoretically, we only need to solve ∆ according to the following formula ∆= −1 (5) However, since some bearings may constrain one or more directions of displacement, the displacement of these constrained directions must be "0". While in the manual calculation, because it is simple and can be judged directly by experience that some displacements are "0", the displacements which are "0" will not be encoded (this can also be shown in Chapter 3: Example) in the first place, so that "a problem" can be avoided.
"A problem" is: in the computer program, which node displacement is " 0 " will not be judged by experience in advance. In fact, it does not conform to the characteristics of the program: mechanical, which is a major difference between the human brain and the computer. Then, all nodes should be encoded directly from beginning to end in programming, which brings a problem: the number of equations ( ∆= ) is likely to be greater than the number of unknowns (∆), then the equation has no solution. So here we need to add an additional program process: "Crossing out rows and columns".

A solution to too many equations -"Crossing out rows and columns"
The rule of "Crossing out rows and columns" is to delete the rows and columns in K matrix when this row or column's displacement is "0", and delete the row in P vector when this row's displacement is "0" (P has only 1 column). In this way, " the number of equations" = "the number of unknowns".
To implement it in the program, define the constraint condition "ConstraintCondition= zeros(n+1,1)" (note that all of the following discussed are node conditions, n members have n + 1 nodes), and the counter "counter=zeros(n+1,1)" (the function of the counter is similar to "lambda", which records how many rows and columns are crossed out in the last loop, locates them in K and P, and cross them out).
After the definition, enter the loop "for j=1:n+1", loop n+1 times. Enter the " value " of constraint condition of each node, where " 0 " indicates no constraint; " 1 " indicates xdirection constraint; "2" indicates y-direction constraint; "3" indicates x&y direction constraint; "4" indicates three direction constraints. Take the case of " 1 -> xdirection constraint" for example, execute in one loop if j == 1 % if it is the first node counter (j, 1) = 0; % Because no row is crossed out before, the counter remains 0. end if j > 1 % if it is not the first node counter (j, 1) = counter (j-1, 1); % The counter maintains the last loop's value. end K(3*j-3+1-counter(j, 1),:) = []; % Cross out the corresponding rows in K. Here, "3*j-3" represents how many rows are in front of the j-th node, "+1" represents the first row of the j-th node, "-counter (j, 1)" represents to subtract the number of rows that have been crossed out according to previous counter records. In MATLAB syntax, "K(, : ) = [] " means to delete the specified row. Therefore, the function of this line of code is to find out the row to be deleted (Step 1: due to the x-direction constraint, the displacement at "3*j-3+1" is "0", so the corresponding row and column in K need to be deleted. Step2: "-counter(j,1)" subtracts the number of previously deleted rows and columns to get the real ones to be deleted), and then delete them. K (:, 3*j-3+1-counter(j, 1)) = []; % Same as deleting rows, this line of code is to locate and delete columns in K. P (3*j-3+1-counter(j, 1),:) = []; % Similarly, this line of code is to locate and delete the corresponding rows in P. if j==1 counter (j, 1) = 1; % The number of rows and columns that have been crossed out above is needed to record at the end of the loop, and because the case takes "1 -> x-direction constraint" as an example, the value of counter of the first node is 1. end if j >1 counter (j, 1) = counter (j-1, 1)+1; % If it is not the first node, then counter value is added by 1 based on the original value. end end The above method is based on the case of "1 -> x-direction constraint". If the constraint is different, we only need to modify it slightly. For example, "2 -> y-direction constraint" only needs to change the line " K (3*j-3+1-counter(j,1),:) = []; " to " K(3*j-3+2counter(j,1),:) = []; " , which represents the corresponding row whose displacement in Y direction is also 0. The same rule is true for the column. For the "3 -> X&Y direction constraint", only another one row and one column need to be crossed out (that is, two rows and two columns in total), and the counter adds 2 more each time. It's also the same rule for "4 -> three direction constraints": three rows and three columns need to be crossed out, and the counter adds 3 more each time.
Judge all j nodes according to the above crossing out rows and columns rules successively. After execution, the remaining ∆= is not the previous ∆= , but ∆= that some of the rows and columns are deleted. For the sake of distinction, write as ∆ = .

Solve and Recover ∆
With the equation ∆ = , the number of equations is equal to the number of unknowns, which can be solved mathematically. Utilize ∆ = −1 (6) Value of ∆ can be calculated directly. However, some rows with 0 items are deleted in ∆ . To get the real node displacement ∆, we need to add the deleted rows with 0 items.

The formation of &
After getting the full ∆, comparing its form with and , we can find that the definition of ∆ is in terms of "nodes", so its form is a 3 × (n + 1) vector. However, the previous and are considered and defined in terms of "members". So it requires some changes that help converted ∆ to ∆ (will be described in detail in chapter 1.4.1). Then, according to the formula = ∆ − (7) = (8) The force of each member in the global coordinate system and the force of each member in the local coordinate system can be calculated.

Visual Design
At this point, the theoretical calculation is finished when (F0) has been solved and the purpose of matrix displacement method is achieved.
The bending moment is an important factor affecting the safety performance of the structure. So after getting the data, the bending moment diagram can be drawn automatically with programming so that it can be shown more intuitive to the reader. In the program, the last function, "DrawM" function, is designed to draw the bending moment diagram according to the bending moment information from . It is divided into two steps: Positioning and Connecting. First, locate the control point of each member "LocationOfMember" (it involves the length, angle and other information of the member), then locate the control point of the bending moment diagram "LocationOfM" (it involves the location, bending moment value, load type and other information of the member). Specifically, the number of control points need to be judged by the load type. "0" (No load) is the direct connection between the two control points. "1" (Uniform load) is three control points are connected with parabola, and "2" (Concentrated load) is three control points are connected with a polyline. After the position of control points of the member and the position of control points of the bending moment are pointed, the bending moment diagram can be drawn by using the "plot" command to connect those points.

Example
According to structural mechanics, the steps of calculating the plane rigid frame with the matrix displacement method are as follows: (1) Sort out the original data and code the members and frames locally and globally (2) Forming element stiffness matrix in local coordinate system (3) Forming element stiffness matrix in global coordinate system (4) Forming integral stiffness matrix K with element integration method (5) Calculating the element equivalent node load vector in the local coordinate system, convert it to the element equivalent node load vector in the global coordinate system, and form the global equivalent node load vector P with the element integration method (6) Solving the equation ∆= , and get the displacement vector of nodes ∆ (7) Solving the internal force vector of the nodes Here is an example which solved from the perspective of manual calculation and computer calculation.
Try to solve the internal force of the frame as shown in figure1. It is assumed that each member is of rectangular section, beam 2 × ℎ 2 =0.5m × 1.26 , column 1 × ℎ 1 =0.5m× 1 .  figure 2. Members CA, AB, and BD are 1, 2 and 3 respectively. Note that since nodes C and D are fixed bearings, which displacement are "0", so they are not encoded in the first place in the manual calculation, which directly avoids the problem of "equation number > unknown number" when solving the equation ∆= . Therefore, node A is encoded as (1,2,3) and node B is encoded as (4,5,6). For the computer calculation, step (1), enter the input information function to input the basic information of members in the rigid frame, as shown in figure 3.
After inputting the basic information, enter the process previously described for processing. As is stated in chapter 1.3, it is impossible to judge whether there is a node displacement of 0 (or not to judge it conform the mechanical characteristic of computer calculation and is more convenient for a program) during the computer calculation, so the formed K is a 12 × 12 matrix (as is highlighted in figure 3), which needs to be crossed out some rows and columns. And step (5) is also completed simultaneously with step (2)  It can be seen that the error between the computer calculation result and the manual calculation result is within the allowable range, the program is accurate.
The moment diagram drawn by 'DrawM' is shown in figure4 below, which is consistent with the manually calculated moment diagram.

Conclusion
At this point, the introduction of the program is completed. This program is suitable for continuous beam and rigid frame. Its main function is to demonstrate how to correspond to the manual calculation of the matrix displacement method with the programmed computer calculation. In order to achieve that, in addition to the basic programming, some skills are included, such as "Using lambda to locate", "Crossing out rows and columns" and visual design.
In addition, there are some small details in the program that need to be paid attention to, such as symbol naming: variables with 0, such as "P0", represent variables without a "T" change in the local coordinate system. Variables with "e", such as "Pe", represent variables without a "T" change in the global coordinate system. Besides, variables used in the program contain a set of information, which is usually a cell matrix, representing the information set of all members(n) or all nodes(n+1). Accordingly, it is necessary to loop the variables n times (number of members) or n+1 times (number of nodes) to get the full information. Detailed variable name corresponding relationship is shown on table 1. The purpose of the program is to consider the initial correspondence between the manual calculation and the computer calculation, to solve the internal force in the program. But further modification and expansion can be made based on that. First, the influence of bearing displacement and temperature change, elastic bearing and shear deformation can be considered to get more accurate results. Secondly, the program of deformation under load can be designed according to the formula " = − " of material mechanics. Finally, although the program has a "DrawM" function to draw the moment diagram, it can be more precise and well developed based on the MATLAB GUI. The intelligent development of the visual graphical user interface and the automatic drawing function of the internal force diagram can further increase its universality and intelligence.