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=True;Needs[\(*! Activate usage of external dll *)(*# 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. *)PlanarTruss[XP_]:=Module[{},(*============================================================================ Planar Truss - Example for Mathematica============================================================================ Inputs: X = [A1, A2, E1, E2, P1, P2, P3, P4, P5, P6, uy] 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 uy: 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 => uy @ N04============================================================================*)(*#--------------------------------------------------------------------------- # Planar truss model #---------------------------------------------------------------------------*)(*# Vector inputs - order from Stochastic Model *)A1=XP[[1]]; (*! cross section of horizontal bars *)A2=XP[[2]]; (*! cross section of inclined bars *)E1=XP[[3]]; (*! Youngs modulus of horizontal bars *)E2=XP[[4]]; (*! Youngs modulus of inclined bars *)P= XP[[5;;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={{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 constants *)nfe=Length[IEN]; (*! number of bars *)nnp=13; (*! number of nodal points *)ned=2; (*! number of dof per node *)ndof=nnp*ned; (*! number of degrees of freedom (dof) *)(*# LM: localization matrix, nfe x dof *)LM=Table[0,{nfe},{2*ned}];For[n=1,n<=nfe,n++, 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;];(*# Deterministic rod properties *)grid=4;ng=12;ang=ArcTan[grid,grid]; (*! inclination angle of the truss [deg] *)theta=Flatten[ConstantArray[{ang,0,-ang,0},ng/2]]; (*! inclination angle [deg] *)theta=Delete[theta,-1];leng=Flatten[ConstantArray[{4/Sqrt[2],4},ng]]; (*! bar length [m] *)leng=Delete[leng,-1];(*# Stochastic rod properties *)Unprotect[E];area=Flatten[ConstantArray[{A2,A1},ng]]; (*! bar cross sectional area [m2] *)area=Delete[area,-1];E=Flatten[ConstantArray[{E2,E1},ng]]; (*! Young's modulus [N/m^2] *)E=Delete[E,-1];(*# Material properties *)k=E*area/leng; (*! stiffness of each bar *)(*# Boundary conditions (supports) *)cc={1,2,14}; (*! dof fixed/supports *)dd=Range[ndof]; (*! dof free *)For[n=Length[cc],n>0,n--,dd=Delete[dd,cc[[n]]];];(*# Boundary conditions (applied forces) *)fc=Array[0&,ndof];For[n=1,n<=Length[P],n++,fc[[14+2*n]]=-P[[n]];]; (*! negative vertical forces [N] *)For[n=Length[cc],n>0,n--,fc=Delete[fc,cc[[n]]];];(*# Global stiffness matrix *)Unprotect[K];K=Table[0,{ndof},{ndof}]; (*! stiffness matrix *)T=Array[Table[0,{4},{4}]&,nfe]; (*! transformation matrix *)For[n=1,n<=nfe,n++, (*! sin and cos of the inclination *) c=Cos[theta[[n]]]; s=Sin[theta[[n]]];