Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Error in beam2ts() #14

Closed
ojseryd opened this issue Feb 9, 2021 · 1 comment
Closed

Error in beam2ts() #14

ojseryd opened this issue Feb 9, 2021 · 1 comment

Comments

@ojseryd
Copy link

ojseryd commented Feb 9, 2021

There is an error in the formulation of beam2ts() which lead to an error when trying to use the function.

Definition looks like this:

def beam2ts(ex,ey,ep,ed,eq=None,np=None):

But since "np" already is used as an alias for numpy all calculations done by the function lead to errors, since "np" now is either "None" or an integer.

There is also a row (calculation of 'v') which tries to use "np.np.power" which needs to be corrected to just "np.power".

Hope this can be adjusted in an upcoming release.

Below is the working code, swapping "np" to "nep", where it is used to describe the number of evaluation points:

def beam2ts(ex,ey,ep,ed,eq=None,nep=None):
"""
Compute section forces in two dimensional beam element (beam2e).

Parameters:

    ex = [x1, x2]
    ey = [y1, y2]       element node coordinates

    ep = [E,G,A,I,ks]   element properties,
                          E:  Young's modulus
                          G:  shear modulus
                          A:  cross section area
                          I:  moment of inertia

    ed = [u1, ... ,u6]  element displacements

    eq = [qx, qy]       distributed loads, local directions 

    nep                 number of evaluation points ( default=2 )
    
Returns:
      
    es = [[N1,V1,M1],   section forces, local directions, in 
          [N2,V2,M2],   n points along the beam, dim(es)= n x 3
          ..........]  

    edi = [[u1,v1,teta1],   element displacements, local directions,
           [u2,v2,teta2],   and rotation of cross section at
           .............]   in n points along the beam, dim(es)= n x 2

(Note! Rotation of the cross section is not equal to dv/dx for Timoshenko beam element)

    eci = [[x1],    local x-coordinates of the evaluation 
           [x2],    points, (x1=0 and xn=L)
           ....]

"""
EA = ep[0]*ep[2]
EI = ep[0]*ep[3]
GAK = ep[1]*ep[2]*ep[4]
alfa = EI/GAK

b = np.mat([
    [ex[1]-ex[0]],
    [ey[1]-ey[0]]
])
L = np.asscalar(np.sqrt(b.T*b))
n = np.asarray(b.T/L).reshape(2)

qx = 0.
qy = 0.
if eq != None:
    qx = eq[0]
    qy = eq[1] 
  
ne = 2

if nep != None:
    ne = nep
    
C = np.mat([
    [ 0., 0.,              0.,   1., 0., 0.],
    [ 0., 0.,              0.,   0., 0., 1.],
    [ 0., 6*alfa,          0.,   0., 1., 0.],
    [ L,  0.,              0.,   1., 0., 0.],
    [ 0., L**3,            L**2, 0., L,  1.],
    [ 0., 3*(L**2+2*alfa), 2*L,  0., 1., 0.]
])

G = np.mat([
    [ n[0], n[1], 0., 0.,   0.,   0.],
    [-n[1], n[0], 0., 0.,   0.,   0.],
    [ 0.,   0.,   1., 0.,   0.,   0.],
    [ 0.,   0.,   0., n[0], n[1], 0.],
    [ 0.,   0.,   0.,-n[1], n[0], 0.],
    [ 0.,   0.,   0., 0.,   0.,   1.]
])

M = np.ravel(C.I*(G*np.asmatrix(ed).T-np.mat([0., 0., 0., -qx*L**2/(2*EA), qy*L**4/(24*EI)-qy*L**2/(2*GAK), qy*L**3/(6*EI)]).T))
C2 = np.mat([M[0], M[3]]).T
C4 = np.mat([M[1], M[2], M[4], M[5]]).T

x = np.asmatrix(np.arange(0., L+L/(ne-1), L/(ne-1))).T
zero = np.asmatrix(np.zeros([len(x)])).T
one = np.asmatrix(np.ones([len(x)])).T

u = np.concatenate((x,one),1)*C2-qx/(2*EA)*np.power(x,2)
du = np.concatenate((one,zero),1)*C2-qx*x/EA

v = np.concatenate((np.power(x,3),np.power(x,2),x,one),1)*C4+qy/(24*EI)*np.power(x,4)-qy/(2*GAK)*np.power(x,2)
dv = np.concatenate((3*np.power(x,2),2*x,one,zero),1)*C4+qy*np.power(x,3)/(6*EI)-qy*x/GAK

teta = np.concatenate((3*(np.power(x,2)+2*alfa*one),2*x,one,zero),1)*C4+qy*np.power(x,3)/(6*EI)
dteta = np.concatenate((6*x,2*one,zero,zero),1)*C4+qy*np.power(x,2)/(2*EI)

N = EA*du
M = EI*dteta
V = GAK*(dv-teta)

es = np.concatenate((N,V,M),1)
edi = np.concatenate((u,v,teta),1)
eci = x

if np != None:
    return es,edi,eci
else:
    return es
@jonaslindemann
Copy link
Contributor

Fixed in the master branch.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants