I wrote a 2D Finite Element program for Axial Loaded Plates, but the results are unexpected

Question

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):

![Text](

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][0]==in_row and nd_rnc[row][1]==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][0]=row
    nd_rnc[n][1]=col
    ndcoor[n][0]=col*(W/(nds_x-1))
    ndcoor[n][1]=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][0]=nd_a
        nds_in_ele[e][1]=nd_b
        nds_in_ele[e][2]=nd_c
    else:   #raise error
        print(nnf_msg)
    #...BUILD LOCAL STIFFNESS MATRIX
    y_bc=ndcoor[nd_b][1]-ndcoor[nd_c][1]    #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][1]-ndcoor[nd_a][1]
    y_ab=ndcoor[nd_a][1]-ndcoor[nd_b][1]
    x_cb=ndcoor[nd_c][0]-ndcoor[nd_b][0]
    x_ac=ndcoor[nd_a][0]-ndcoor[nd_c][0]
    x_ba=ndcoor[nd_b][0]-ndcoor[nd_a][0]
    x_bc=ndcoor[nd_b][0]-ndcoor[nd_c][0]
    y_ac=ndcoor[nd_a][1]-ndcoor[nd_c][1]
    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][0]=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][1]=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]

Answer 1

The primary problem is that the CST approximation has a different displacement response depending on the orientation of mesh elements relative to the applied element loading (you're only allowed to applied forces on the nodes of triangles, so distributed loads must be approximated).

You can see the effect of this by looking at only a single triangle attached to a wall with a perfectly horizontal force applied (drawn by hand and exaggerated):

Adding a second triangle, there's nothing "forcing" the triangle to compress upwards, so it will remain displaced downwards in y.

As long as you keep tessellating (repeating) this pattern to build up your beam, nothing will cancel out this asymmetry (though I think increasing mesh resolution will reduce the effect since the actual force applied to any single element will be reduced).

You can however, decide to mirror the elements to take advantage of the fact that this effect can cancel itself out:

Here, you can see that the top half will tend to try to displace the top edge downards, however the bottom wants to displace the bottom edge upwards, leaving the midline unmoved. Tesselating this 2x2 triangularized element pattern will allow you to produce a more isotropic behavior, however you should of course start to wonder why didn't you just use rectangular elements to begin with.

Somewhat luckily, you don't necessarily need this exact pattern to reduce this anisotropy; a relatively "random" orientation of elements can cancel out most of this effect, however the downside of this is that you need a lot of elements to ensure that the few elements which don't cancel out is small compared to the number of elements in the neighborhood of elements.

Answer 2

Here is how you want to test this and you need only two elements in the mesh. You want to define your left BC so it will reproduce a constant stress state as follows: assuming $u$ is the displacement in the x-direction and $v$ the displacement in the y-direction, set $u=0$ at the two nodes on this edge and $v=0$ at the bottom node on this edge. The two nodes on the right edge have an x-direction force of $P/2$. The numerical solution should give the exact solution to this simple plane stress problem as follows: if $L$, $h$, and $t$ are the length, height, and thickness of the strip, the u-displacement at the right edge is given by $u_{right}=P L/(E h t)$ and the v-displacement at the top edge is given by $-\nu u_{right}$ ($E$ and $\nu$ are the material constants).

It is required that a finite element formulation exactly represent this constant stress state in order for it to lead to an exact solution as the mesh is refined. The constant stress triangle element satisfies this requirement.

Answer 3

I think that the answer by @helloworld922 is misleading. The first image shown in the answer seems to be an effect of the Poisson effect, a contraction in one direction due to loads applied in the other direction.

If you want to obtain a state of constant stress in your simulation you need to change the boundary conditions that you are applying, namely:

• all nodes on the left should be constrained in the horizontal direction;
• only the node in the middle should be constrained in the vertical direction; and
• forces on the right are not equal because the contributions from the contiguous elements will make the middle one larger.

Using the same material properties than you, and a total load of 10000 (split as 2500, 5000 and 2500), I get the following

I used a mesh of 750 mm long instead of 800 for nicer numbers. The displacement is multiplied by 10 for visualization purposes.