notebook.community

Edit and run

1D Structural element (Rod in PtFEM.jl) with a mesh consisting of 3 finite elements (Line in PtFEM.jl) and 4 nodes





In [1]:



    


using PtFEM, Plots, DataFrames



    











In [2]:



    


l = 1.0                      # Length [m]
q = 5.0                      # Distributed load [N/m]
EA = 10.0                    # [Ns/m];
N = 4                        # Number of nodes
els = 3                      # Number of elements
ell = [1/els for i in 1:els] # Length of each element;

Specify the type of structural element and the type of finite element used:





In [3]:



    


struct_el = :Rod
fin_el = :Line



    


















    
Out[3]:








:Line

















In [4]:



    


?Rod



    


















    






search: Rod prod prod! produce cumprod cumprod! broadcast broadcast!











    
Out[4]:







Rod
Concrete 1D structural element with only axial stresses.


Constructor
Rod(nels, np_types, nip, fin_el)



Arguments
* nels::Int64             : Number of fin_els (stored in field nxe)
* np_types::Int64         : Number of different property types
* nip::Int64              : Number of integration points
* fin_el::FiniteElement   : Line(nod, nodof)



?StructuralElement  : Help on structural elements
?FiniteElement      : Help on finite element types
?Line               : Help on a Line finite element

In this example the mesh of the structural element will consist of 3 element, each only made of a single material and the finite element type is Line:





In [5]:



    


np_types = 1       # A single property type
EA = 1.0e5         # Axial strain stiffness
nip = 1            # Numerical integration will use only a single integration point (explained in future notebooks);



    











In [6]:



    


?Line



    


















    






search: Line LineNumberNode linearindices vline hline vline! hline! inline











    
Out[6]:







Line (Interval)
1D type finite element


Constructor
Line(nod, nodof)
Line(nodof)



Arguments
* nod::Int64       : Number of nodes for finite element, defaults to 2
* nodof::Int64     : Number of degrees of freedom per node



?FiniteElement      : Help on finite element types



















In [7]:



    


nod = 2        # Each line finite element has 2 nodes
nodof = 1      # Each node has a single degree of freedom (deflection in the x direction);



    











In [8]:



    


data = Dict(
:struc_el => Rod(els, np_types, nip, Line(nod, nodof)),
  :properties => [EA;],
  :x_coords => 0.0:l/els:l,
  :support => [(N, [0])],
  :loaded_nodes => [(i, repeat([0.0], inner=nodof)) for i in 1:N],
  :eq_nodal_forces_and_moments => [(i, repeat([0.0], inner=nodof*nod)) for i in 1:els]
);
In above input dictionary we've also specified the support at node 4. In this case the support is clamped, i.e. the deflection in the x direction is fixed. Dictionary entries :loaded_nodes and :eq_nodal_forces_and_moments are initial set to 0.0. These are updated below where applicable.

In this example there are only distributed loads. Otherwise set data[:loaded_nodes] to the external forces and moment directly applied to the nodes, e.g.:

data[:loaded_nodes][2] = (2, [5.0])
Determine the distributed loads contribution to nodal forces:




In [9]:



    


for node in 1:els
  data[:eq_nodal_forces_and_moments][node][2][1] = 1/2*q*l/els
  data[:eq_nodal_forces_and_moments][node][2][2] = 1/2*q*l/els
end  

# Add the equivalent distributed forces and moment to loaded_nodes entry

for node in 1:N-1
  data[:loaded_nodes][node][2][1] += data[:eq_nodal_forces_and_moments][node][2][1]
  data[:loaded_nodes][node+1][2][1] += data[:eq_nodal_forces_and_moments][node][2][2]
end



    











In [10]:



    


F = data[:loaded_nodes]



    


















    
Out[10]:








4-element Array{Tuple{Int64,Array{Float64,1}},1}:
 (1,[0.833333])
 (2,[1.66667]) 
 (3,[1.66667]) 
 (4,[0.833333])

Vector F contains the nodal loads except the support force in node 4 (which is unknown in general)

In this case it is easy to see the support force will be -5 [N], but in most cases this is not that trivial and we would like the program to compute such forces and/or moments.

In EEM: f = {F[1], F[2], F[3], F[4] + Fs[4]} # where Fs[4] is the unknown support force needed in node 4 to prevent a deflection.





In [11]:



    


m = FE4_1(data);



    


















    






There are 3 equations and the skyline storage is 5.


















In [12]:



    


m.displacements



    


















    
Out[12]:








4×1 Array{Float64,2}:
 2.5e-5    
 2.22222e-5
 1.38889e-5
 0.0       

















In [13]:



    


m.actions



    


















    
Out[13]:








3×2 Array{Float64,2}:
 0.833333  -0.833333
 2.5       -2.5     
 4.16667   -4.16667 

















In [14]:



    


dis_df = DataFrame(
  x_translation = m.displacements[:, 1],
)



    


















    
Out[14]:







x_translation
12.499999999999999e-5
22.2222222222222217e-5
31.3888888888888884e-5
40.0

















In [15]:



    


fm_df = DataFrame(
    normal_force_1 = m.actions[:, 1],
    normal_force_2 = m.actions[:, 2],
    normal_force_1_corrected = m.actions[:, 1],
    normal_force_2_corrected = m.actions[:, 2]
);

Correct element forces and moments for equivalent nodal forces and moments introduced for loading between nodes





In [16]:



    


if :eq_nodal_forces_and_moments in keys(data)
    eqfm = data[:eq_nodal_forces_and_moments]
    k = data[:struc_el].fin_el.nod * data[:struc_el].fin_el.nodof
    for t in eqfm
      for i in 1:k
        fm_df[t[1], i+2] = round(fm_df[t[1], i] - t[2][i], 5)
      end
    end
end
fm_df



    


















    
Out[16]:







normal_force_1normal_force_2normal_force_1_correctednormal_force_2_corrected
10.8333333333333321-0.8333333333333321-0.0-1.66667
22.5-2.51.66667-3.33333
34.166666666666664-4.1666666666666643.33333-5.0

















In [17]:



    


### Plot the results



    











In [18]:



    


gr()



    


















    
Out[18]:








Plots.GRBackend()

















In [19]:



    


x = 0.0:l/els:l;
u = dis_df[:x_translation]
N = vcat(fm_df[:normal_force_1_corrected][1], fm_df[:normal_force_2_corrected]);



    











In [20]:



    


p = Vector{Plots.Plot{Plots.GRBackend}}(2)
titles = ["EEM fig 1.1 u(x)", "EEM fig 1.1 N(x)"];



    











In [21]:



    


p[1]=plot(N, x, yflip=true, xflip=true, xlab="Normal force [N]", ylab="x [m]", color=:red,
    fill=true, fillalpha=0.1, leg=false, title=titles[2])
for i in 1:els
    plot!(p[1], [fm_df[:normal_force_2][i], fm_df[:normal_force_2][i]], [(i-1)*l/els, i*l/els], color=:blue)
end



    











In [22]:



    


p[2] = plot(u, x, xlim=(0.0, 0.00004), yflip=true, xlab="Displacement [m]", ylab="x [m]",
fillto=0.0, fillalpha=0.1, leg=false, title=titles[1]);



    











In [23]:



    


defl(x) = q/(2*EA).*(l^2 .- x.^2)
x1 = 0.0:0.01:l
plot!(p[2], defl(x1), x1, xlim=(0.0, 0.00004), color=:red);



    











In [24]:



    


plot(p..., layout=(1,2))



    


















    
Out[24]:









  
    
  



  
    
  



  
    
  




























-5


-4


-3


-2


-1


0


0.0


0.2


0.4


0.6


0.8


1.0


EEM fig 1.1 N(x)


Normal force [N]


x [m]








  
    
  


























0.00000


0.00001


0.00002


0.00003


0.00004


0.0


0.2


0.4


0.6


0.8


1.0


EEM fig 1.1 u(x)


Displacement [m]


x [m]






















In [ ]:

Content source: goedman/RobGoedmansNotebooks.jl

Similar notebooks:

notebook.community | gallery | about