Page 3 - Demo

Planar Truss Examplefor ComrelAdd-onRCPConsult, 2021-2025Page 3(*# Global flags to manage plot facilities of engine *)StrurelPlot=True;StrurelPlotMode=3;StrurelPlotName=\_\StrurelPlotType=\Needs[\(*! Activate usage of external dll *)(*# Code based on Max Ehre work https://github.com/maxehre/polynomial_surrogatesUsed 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 distributedA1,A2: cross section of horizontal and inclined bars resp. -log-normally distributedP1-P6: negative vertical forces attacking in nodes 8-13. -Gumbel distributeduy: maximal allowed vertical deflection in node 4. constantOutput: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/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 => 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 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 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]]];(*! coordinate transformation matrix for the bar n *)