Page 3 - Demo
P. 3
Planar Truss Example for Comrel Add-on RCP Consult, 2019-2024 Page 3## Global flag to manage plot facilities of engineStrurelPlot=False## Code based on Max Ehre work - https://github.com/maxehre/polynomial_surrogates## Used with his kind permission as version from 2017 that adapted to Strurel.import numpy as np # https://numpy.orgif StrurelPlot :import matplotlib.pyplot as plt # https://matplotlib.orgimport matplotlib.cm as cmximport ctypesdef Planar_Truss(XP):'''============================================================================ Planar Truss - Example for Python============================================================================ Inputs: X = [A1, A2, E1, E2, P1, P2, P3, P4, P5, P6, u_y] E1,E2: Youngs modulus of horizontal and inclined bars resp. - log-normally distributed A1,A2: cross section of horizontal and inclined bars resp. - log-normally distributed P1-P6: negative vertical forces attacking in nodes 8-13. - Gumbel distributed u_y: maximal allowed vertical deflection in node 4. - constant Output: uout: vertical truss deflection at bottom center (N04)============================================================================ Truss description taken from Lee, S.H., Kwak, B.M. (2006). Response surface augmented moment method for efficient reliability analysis. Structural Safety 28(3), 261 - 272. https://www.sciencedirect.com/science/article/abs/pii/S0167473005000421 Based on Diego Andr%u00e9s Alvarez Mar%u00edn work - https://github.com/diegoandresalvarez/elementosfinitos/tree/master/codigo/repaso_matricial/cercha_2d N: Nodes R: Rods============================================================================ N13__R4___N12__R8___N11__R12__N10__R16__N09__R20__N08 /\\ /\\ /\\ /\\ /\\ /\\ / \\ / \\ / \\ / \\ / \\ / \\ R1 R3 R5 R7 R9 R11 R13 R15 R17 R19 R21 R23 / \\ / \\ / \\ / \\ / \\ / \\ /___R2___\\/___R6___\\/__R10___\\/__R14___\\/__R18___\\/__R22___\\ N01 N02 N03 N04 N05 N06 N07 ||\\/ uout => u_y @ N04============================================================================'''##---------------------------------------------------------------------------## Planar truss model##---------------------------------------------------------------------------## Vector inputs - order from Stochastic ModelA1 = XP[0] # cross section of horizontal barsA2 = XP[1] # cross section of inclined barsE1 = XP[2] # Youngs modulus of horizontal barsE2 = XP[3] # Youngs modulus of inclined barsP = XP[4:10] # negative vertical forces attacking in nodes 8-13. (1x6 vector)## Element, nodes and dofs association## IEN: connectivity matrix, nfe x nodes - bar 1 has nodes 1 and 13, bar 2 has nodes 1 and 2 ...IEN = np.array([[1, 13], [1, 2], [13, 2], [13, 12], [2, 12], [2, 3], [12, 3], [12, 11], [3, 11], [3, 4], [11, 4], [11, 10], [4, 10], [4, 5], [10, 5], [10, 9], [5, 9], [5, 6], [9, 6], [9, 8], [6, 8], [6, 7], [8, 7]]) - 1## Deterministic rod propertiesgrid = 4 # size of grid side in [m]ng = 12 # count of grid cellsang = np.arctan2(grid,grid) # inclination angle of the truss [rad]theta = np.tile([ang, 0, -ang, 0],int(ng/2)) # inclination angle [rad]theta = np.delete(theta,-1)leng = np.tile([grid/np.sqrt(2), grid],ng) # bar length [m]leng = np.delete(leng,-1)## FEM constantsned = 2 # number of dof per nodennp = 13 # number of nodal pointsnfe = np.shape(IEN)[0] # number of barsndof = ned*nnp # number of degrees of freedom (dof)## Dof: degrees of freedom (rows are the nodes, cols are dofs)dof = np.array(range(ndof)).reshape(nnp,ned)## Stochastic rod propertiesarea = np.tile([A2, A1],ng) # bar cross sectional area [m2]area = np.delete(area,-1)E = np.tile([E2, E1],ng) # Young's modulus [N/m^2]E = np.delete(E,-1)## Material propertiesk = E*area/leng # stiffness of each bar## Boundary conditions (supports)cc = np.array([1, 2, 14])-1 # dof fixed/supportsdd = np.array(list(set(range(ndof)) - set(cc))) # dof free## Boundary conditions (applied forces)fc = np.zeros(ndof)fc[15::2] =-np.array(P) # negative vertical forces [N]fc = np.delete(fc,cc)fc = np.expand_dims(fc, axis=1)## Global stiffness matrixK = np.zeros((ndof,ndof)) # stiffness matrix