Page 3 - Demo

Planar TrussExamplefor ComrelAdd-onRCPConsult, 2021-2025Page 3## Global variables to manage plot facilities of engine$StrurelPlot=true$StrurelPlotName='PlanarTruss_';$StrurelPlotType='.png';$StrurelPlotMode=2require'numo/narray'# https://github.com/rubynumo/numonarrayrequire'numo/linalg'# https://github.com/rubynumo/numolinalgincludeNumo 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.defPlanarTruss(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 distributedA1,A2: cross section of horizontal and inclined bars resp. lognormally distributedp1p6: negative vertical forces attacking in nodes 813. Gumbel distributedOutput: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/S0167473005000421Based on Diego Andr%u00e9s Alvarez Mar%u00edn work -https://github.com/diegoandresalvarez/elementosfinitos/tree/master/codigo/repaso_matricial/cercha_2dN: NodesR: 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 813. (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)forn in0..nfe1@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=@Edelete(1)## Material propertiesk=@E*areab/leng## Global stiffness matrix@K=DFloat.zeros(ndof,ndof)# stiffness matrix@T=[]# transformation matrixforn in0..nfe-1# sin and cos of the inclinationco=Mathcos(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[1ndof].delete(c)# dof free## Boundary conditions (applied forces)f=DFloat.zeros(ndof)forn in0@Plength1f[15+2*n]=-@P[n]# negative vertical forces [N]endf=f.delete(c)## Solve the system of equations@ac=DFloatzeros(3)