Page 3 - Demo

Planar Truss Example for Comrel Add-on RCP Consult, 2021-2025 Page 3## Global flag to manage plot facilities of engineStrurelPlot:=true:## 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.with(LinearAlgebra): # https://www.maplesoft.com/support/help/maple/view.aspx?path=LinearAlgebrawith(ArrayTools): # https://www.maplesoft.com/support/help/maple/view.aspx?path=ArrayToolsif StrurelPlot then :with(Maplets[Examples]): # https://www.maplesoft.com/support/help/Maple/view.aspx?path=Maplets/Exampleswith(Maplets[Elements]):with(ColorTools):with(plottools):with(plots):plotsetup(maplet):end if:PlanarTruss := proc(INPUT)(*============================================================================ Planar Truss - Example for Maple============================================================================ 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 input - order from Stochastic ModelA1 := INPUT[1]:A2 := INPUT[2]:E1 := INPUT[3]:E2 := INPUT[4]:p := INPUT(5..10): # Maple hfarray## 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 := upperbound(IEN): # number of barsnnp := 13: # number of nodal pointsned := 2: # number of dof per nodendof := ned*nnp: # number of degrees of freedom (dof)## LM: localization matrix, nfe x dof LM := Matrix(nfe, 2*ned):for n to nfe do LM[n, 1] := IEN[n, 1]*ned-1: LM[n, 2] := IEN[n, 1]*ned: LM[n, 3] := IEN[n, 2]*ned-1: LM[n, 4] := IEN[n, 2]*ned:end do:## Deterministic rod propertiesang := arctan(200, 200): # inclination angle of the truss [deg]theta := Replicate(Array([ang, 0, -ang, 0]), 6): Remove(theta,-1): # inclination angle [deg]leng := Replicate(Array(4*[1/sqrt(2), 1]), 12): Remove(leng,-1): # bar length [m]## Stochastic rod propertiesArea := Replicate(Array([A2, A1]), 12): Remove(Area,-1): # bar cross sectional area [m2]E := Replicate(Array([E2, E1]), 12): Remove(E,-1): # Young's modulus [N/m^2]## Material propertiesk := E*Area/leng: # stiffness of each bar## Boundary conditions (supports)cc := Array([1, 2, 14]): # dof fixed/supportsdd := Array(1..ndof, i -> i): # dof freefor n from upperbound(cc) by -1 to 1 do dd := Remove(dd, cc[n]):end do:## Boundary conditions (applied forces)fc := Vector(1..ndof):nP := upperbound(p):for n to nP do fc[14 + 2*n] := -p[n]: # negative vertical forces [N]end do:for n from upperbound(cc) by -1 to 1 do fc := Remove(fc, cc[n]):end do:## Global stiffness matrixK := Matrix(ndof, ndof): # stiffness matrixT := Array(1..nfe): # transformation matrixfor n to nfe do # sin and cos of the inclination c := cos(theta[n]): s := sin(theta[n]): # coordinate transformation matrix for the bar n T[n] := Matrix([[c, s, 0, 0], [-s, c, 0, 0], [0, 0, c, s], [0, 0, -s, c]]): # local stiffness matrix for the bar n Kloc := k[n]*Matrix([[1, 0, -1, 0], [0, 0, 0, 0], [-1, 0, 1, 0], [0, 0, 0, 0]]): # global stiffness matrix assembly idx:=[entries(LM[n,1..],'nolist')]: