TLDR: I used Python to write a 2D Finite Element program using 'Constant Strain Triangles' and my beam keeps pointing slightly upwards instead of straight sideways (like the force). I'm new to FEA and have very little linear algebra experience so I don't have the insight to know if I did something fundamentally wrong.
So, for now, this program is meant to simulate the strain and displacement of nodes in a thin plate (or beam) that's in tension due to a distributed external force, i.e. configurations that look something like this (force in image obviously isn't distributed but you get the idea):
I used the constant strain triangles method because triangular elements will be convenient for the next part of the project when the plates aren't simple rectangles. My main resource has been the lecture and example here (which is pretty much the same as the info here).
I ran the program and the displacements in the x-direction for each node seemed reasonable, but each node seems wants to drift 'upwards' instead of bowing inwards due to the Poisson Effect: (Pardon my homemade graphic). As you can see, the beam w/force applied is tilted upwards, which I found very strange. It does the same thing for different height/width beams, and if I add more nodes. (see Edit)
I'm new to FEA in general (haven't even used a commercial package) and I have extremely limited experience with linear algebra. What did I do that led to this?
- I know CST is apparently not as accurate as other methods, but would that cause this problem?
- I read (after the fact) that the triangular elements should be as equilateral-as-possible, so do I have a problem since my elements are right triangles?
- Does traction have something to do with this? The term keeps popping up in lectures I read but I don't fully understand what it means.
- What else should I look into?
Thanks in advance to anyone who looks into this. I tried to read-up on the problem before posting but I came up empty, so I figured I'd try to post here. Any help is appreciated!
Edit: I successfully fixed this problem by adjusting my meshing algorithm so that the repeating pattern was mirrored, as suggested in the checked answer. Also, it seems that elements with legs closer in length work better. The output of my program is shown below: the plate bows inward symmetrically about the neutral axis now. I don't have a graphic of a long beam like I originally mentioned, but I tried it out and it also worked. Thanks to everyone who had suggestions!
Original code (Python):
import graphics as gr import numpy as np import math import matplotlib.pyplot as plt #constants P=10000.0 #Load (Newtons) W=800 #Width of Beam (mm) H=50 #Height of Beam (mm) Z=0.05 #Thickness of Beam (mm) E_beam=10**5 #Beam Elastic Modulus pr_beam=0.45 #Poissons Ratio of the beam nds_x=4 #number of nodes extending in the horizontal direction nds_y=3 #number of nodes extending in the vertical direction nnds=nds_x*nds_y #total number of nodes ndof=nnds*2 #total number of degrees of freedom in the whole system nele=2*(nds_x-1)*(nds_y-1) #total number of elements eper=2*(nds_x-1) #elements per element row ndcoor=np.zeros((nnds,2)) #Table which stores the INITIAL coordinates (in terms of mm) for each node nd_rnc=np.zeros((nnds,2)) #Table which stores the 'row and column' coordinates for each node nds_in_ele=np.zeros((nele, 3)) #the nodes which comprise each element glbStiff=np.zeros((ndof,ndof)) #global stiffness matrix (GSM) lst_wallnds= #List of nodes (indices) which are coincident with the rigid wall on the left lst_wallnds.clear() lst_walldofs= #dofs indices of nodes coincident with the rigid wall lst_walldofs.clear() lst_endnds= #nodes on the free edge of the beam lst_endnds.clear() nnf_msg='Node not found!' #Function 'node_by_rnc' returns the index of the node which has the same row and column as the ones input (in_row, in_col) def node_by_rnc(in_row, in_col, start_mrow): #'start_mrow' == where the func starts searching (for efficiency) run=True row=start_mrow while run==True: if row>nnds-1: run=False elif nd_rnc[row]==in_row and nd_rnc[row]==in_col: run=False else: row=row+1 if row>nnds-1: return nnf_msg #returns error message else: return row #Function 'add_to_glbStiff' takes a local stiffness matrix and adds the value of each 'cell' to the corrosponding cell in the GSM def add_to_glbStiff(in_mtrx, nd_a, nd_b, nd_c): global glbStiff #First column in local stiffness matrix (LSM) is the x-DOF of Node A, second is the y-DOF of Node A, third is the x-DOF of Node B, etc. (same system for rows; the matrix is symmetric) dofs=[2*nd_a, 2*nd_a+1, 2*nd_b, 2*nd_b+1, 2*nd_c, 2*nd_c+1] #x-DOF for a node == 2*[index of the node]; y-DOF for node == 2*[node index]+1 for r in range(0,6): #LSMs are always 6x6 for c in range(0,6): gr=dofs[r] #gr == row in global stiffness matrix gc=dofs[c] #gc == column in global stiffness matrix glbStiff[gr][gc]=glbStiff[gr][gc]+in_mtrx[r][c] #Add the value of the LSM 'cell' to what's already in the corrosponding GSM cell for n in range(0,nnds): #puts node coordinates and rnc indices into matrix row=n//nds_x col=n%nds_x nd_rnc[n]=row nd_rnc[n]=col ndcoor[n]=col*(W/(nds_x-1)) ndcoor[n]=row*(H/(nds_y-1)) if col==0: lst_wallnds.append(n) elif col==nds_x-1: lst_endnds.append(n) for e in range(0,nele): #FOR EVERY ELEMENT IN THE SYSTEM... #...DETERMINE NODES WHICH COMPRISE THE ELEMENT erow=e//eper #erow == the row which element 'e' is on eor=e%eper #element number on row (i.e. eor==0 means the element is attached to rigid wall) if eor%2==0: #downwards-facing triangle nd_a_col=eor/2 nd_b_col=eor/2 nd_c_col=(eor/2)+1 nd_a=node_by_rnc(erow, nd_a_col, nds_x*erow) nd_b=node_by_rnc(erow+1, nd_b_col, nds_x*erow) nd_c=node_by_rnc(erow, nd_c_col, nds_x*erow) else: #upwards-facing triangle nd_a_col=(eor//2)+1 nd_b_col=(eor//2)+1 nd_c_col=eor//2 nd_a=node_by_rnc(erow+1, nd_a_col, nds_x*(erow+1)) nd_b=node_by_rnc(erow, nd_b_col, nds_x*erow) nd_c=node_by_rnc(erow+1, nd_c_col, nds_x*(erow+1)) if nd_a!=nnf_msg and nd_b!=nnf_msg and nd_c!=nnf_msg: #assign matrix element values if no error nds_in_ele[e]=nd_a nds_in_ele[e]=nd_b nds_in_ele[e]=nd_c else: #raise error print(nnf_msg) #...BUILD LOCAL STIFFNESS MATRIX y_bc=ndcoor[nd_b]-ndcoor[nd_c] #used "a, b, c" instead of "1, 2, 3" like the the example PDF; ex: 'y_bc' == 'y_23' == y_2 - y_3 y_ca=ndcoor[nd_c]-ndcoor[nd_a] y_ab=ndcoor[nd_a]-ndcoor[nd_b] x_cb=ndcoor[nd_c]-ndcoor[nd_b] x_ac=ndcoor[nd_a]-ndcoor[nd_c] x_ba=ndcoor[nd_b]-ndcoor[nd_a] x_bc=ndcoor[nd_b]-ndcoor[nd_c] y_ac=ndcoor[nd_a]-ndcoor[nd_c] detJ=x_ac*y_bc-y_ac*x_bc Ae=0.5*abs(detJ) D=(E_beam/(1.0-(pr_beam**2.0)))*np.array([[1.0, pr_beam, 0.0],[pr_beam, 1.0, 0.0],[0.0, 0.0, (1-pr_beam)/2.0]]) B=(1.0/detJ)*np.array([[y_bc, 0.0, y_ca, 0.0, y_ab, 0.0],[0.0, x_cb, 0.0, x_ac, 0.0, x_ba],[x_cb, y_bc, x_ac, y_ca, x_ba, y_ab]]) BT=np.transpose(B) locStiff=Z*Ae*np.matmul(np.matmul(BT,D),B) #...ADD TO GLOBAL STIFFNESS MATRIX add_to_glbStiff(locStiff, nd_a, nd_b, nd_c) #Deleting contrained DOFs from the GSM nwnds=len(lst_wallnds) #number of wall nodes for w in range(0,nwnds): #Populates list of all DOFs which have 0 displacement (the corrosponding rows and columns get completely erased from GSM) lst_walldofs.append(2*lst_wallnds[w]) lst_walldofs.append(2*lst_wallnds[w]+1) glbStiff=np.delete(np.delete(glbStiff, lst_walldofs, 0), lst_walldofs, 1) #delete the rows and columns corrosponding to the DOFs that are fixed #Keeping track of what rows (and columns) in the 'new' GSM corrospond to which DOF indices lst_frdofs=np.zeros(ndof) #lst_frdofs = List of "Free" DOFS i.e. DOFs NOT coincident with the wall for d in range(0,ndof): lst_frdofs[d]=d #Before deleting fixed DOFs: [the global index for each DOF] == [the corrosponding row/column in the GSM]... lst_frdofs=np.delete(lst_frdofs,lst_walldofs) #...after deleting the fixed DOF rows/columns: 'lst_frdofs' stores the global index for each DOF in the row corrosponding the the row in the GSM #Specifying the Load lpn=P/nds_y #Load per Node (on free end) mtrx_P=np.zeros(ndof) #The vector which stores the input force values for each DOF for en in range(0, len(lst_endnds)): mtrx_P[2*lst_endnds[en]]=lpn #Applies a force of 'lpn' to each node on the free end in the X-direction mtrx_P=np.delete(mtrx_P, lst_walldofs) #Deletes the rows corrosponding to the DOFs that were deleted from the GSM #Solve for q for each DOF mtrx_q=np.linalg.solve(glbStiff, mtrx_P) #Determining the final locations of each node nd_disp=np.zeros((nnds,2)) #Tabulating how much each node moved in the x and y directions for g in range(0,len(lst_frdofs)): gdof=lst_frdofs[g] if gdof%2==0: #even global DOF index -> displacement in the x-direction nd=int(gdof/2) #nd == node which the DOF (gdof) belongs to nd_disp[nd]=mtrx_q[g] #add the displacement to the table/matrix else: #odd global DOF index -> displacement in the y-direction nd=int((gdof-1)/2) nd_disp[nd]=mtrx_q[g] fnl_ndcoor=np.add(ndcoor, nd_disp) #[Final coordinates (in terms of mm) for each node] = [Original coordinates for that node] + [the displacement of the node]