Page 3 - Demo

Planar Truss Example for Comrel Add-on RCP Consult, 2021-2025 Page 3## Global flag to manage plot facilities of engine$StrurelPlot=falserequire 'numo/narray' # https://github.com/ruby-numo/numo-narrayrequire 'numo/linalg' # https://github.com/ruby-numo/numo-linalginclude Numo if $StrurelPlotrequire 'matplotlib/pyplot' # https://github.com/mrkn/matplotlib.rbend## 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.def PlanarTruss(xp)=begin============================================================================ Planar Truss Example============================================================================ Inputs: X = [A1, A2, E1, E2, p1, p2, p3, p4, p5, p6] 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 Output: 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=============================================================================end##---------------------------------------------------------------------------## Planar truss model##---------------------------------------------------------------------------## Vector inputs - order from Stochastic Model@A1=xp[0] # cross section of horizontal bars@A2=xp[1] # cross section of inclined bars@E1=xp[2] # Youngs modulus of horizontal bars@E2=xp[3] # Youngs modulus of inclined bars@P =xp[4..9] # 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=[[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]]## FEM constantsnfe=@IEN.length; # number of barsnnp=13; # number of nodal pointsned=2; # number of dof per nodendof=nnp*ned; # number of degrees of freedom (dof)## LM: localization matrix, nfe x dof @LM=Int32.zeros(nfe,2*ned)for n in 0..nfe-1 @LM[n,0]=@IEN[n][0]*ned-1; @LM[n,1]=@IEN[n][0]*ned @LM[n,2]=@IEN[n][1]*ned-1; @LM[n,3]=@IEN[n][1]*nedend## Deterministic rod propertiessizeg=2nsec=6ang=Math.atan2(sizeg,sizeg) # inclination angle of the truss [deg]theta=DFloat[*[ang,0,-ang,0]*nsec] # inclination angle [deg]theta=theta.delete(-1)leng=DFloat[*[sizeg*2/Math.sqrt(sizeg),sizeg*2]*nsec*2] # bar length [m]leng=leng.delete(-1)## Stochastic rod propertiesareab=DFloat[*[@A2,@A1]*nsec*2] # bar cross sectional area [m2]areab=areab.delete(-1)@E=DFloat[*[@E2,@E1]*nsec*2] # Young's modulus [N/m^2]@E=@E.delete(-1)## Material propertiesk=@E*areab/leng## Global stiffness matrix@K=DFloat.zeros(ndof,ndof) # stiffness matrix@T=[] # transformation matrixfor n in 0..nfe-1 # sin and cos of the inclination co=Math.cos(theta[n]) si=Math.sin(theta[n]) # coordinate transformation matrix for the bar n @T[n]=DFloat[*[[co, si, 0, 0],[-si, co, 0, 0],[0, 0, co, si],[0, 0, -si, co]]] # local stiffness matrix for the bar n @Kloc=k[n]*DFloat[*[[1, 0, -1, 0],[0, 0, 0, 0],[-1, 0, 1, 0],[0, 0, 0, 0]]] # global stiffness matrix assembly @K[@LM[n,0..]-1,@LM[n,0..]-1]+=((@T[n].transpose).dot(@Kloc)).dot(@T[n])end## Boundary conditions (supports)c=Int32[1,2,14]-1 # dof fixed/supportsd=Int32[1..ndof].delete(c) # dof free## Boundary conditions (applied forces)f=DFloat.zeros(ndof)for n in 0..@P.length-1 f[15+2*n]=-@P[n] # negative vertical forces [N]endf=f.delete(c)## Solve the system of equations@ac=DFloat.zeros(3)@Kcc=@K[c,c]