Merge pull request #299 from ONSAS/jorge

changes in roorda and arch examples

.github/workflows/tests.yml

@@ -9,7 +9,7 @@ jobs:
       - name: run mechanical-work tests
         uses: joergbrech/moxunit-action@v1.1
         with:
-          src: examples/staticVonMisesTruss examples/linearPlaneStrainExample examples/uniformCurvatureCantilever examples/uniaxialExtension src src/elements src/vtk src/mesh
+          src: examples/staticVonMisesTruss examples/linearPlaneStrain examples/uniformCurvatureCantilever examples/uniaxialExtension src src/elements src/vtk src/mesh
           tests: ./test/runTestProblems_moxunit_disp.m
   # tests-disp-coverage:
   #   runs-on: [ubuntu-latest]

examples/circularArch/onsasExample_circularArch.m

@@ -0,0 +1,99 @@
+
+close all, clear all ; addpath( genpath( [ pwd '/../../src'] ) );
+
+% structural parameters
+E  = 210e9 ;		% Young's modulus (Pa)
+nu = 0.3 ;		% Poisson ratio
+a  = .05; 		% Rectangular section width (mm)
+b  = .01; 		% Rectangular section height (mm)
+A  = a*b; 		% Rectangular section Area (mm2)
+R  = 1 ;
+fiApoyo = 30 ;
+
+nEM = 10; 			% numero elementos finitos mitad de arco
+
+%Material Definitions
+materials(1).hyperElasModel  = '1DrotEngStrain' ;
+materials(1).hyperElasParams = [ E nu ] ;
+
+%Element Definitions
+elements(1).elemType = 'node' ;
+elements(2).elemType = 'frame';
+elements(2).elemTypeGeometry = [2 a b ] ;
+
+%Boundary Conditions Definitions
+boundaryConds(1).imposDispDofs = [ 1 3 5 ] ;	% pinned node
+boundaryConds(1).imposDispVals = [ 0 0 0 ] ;
+boundaryConds(2).imposDispDofs =  3 ;			% node restrained in X-Y plane
+boundaryConds(2).imposDispVals =  0 ;
+boundaryConds(2).loadsCoordSys = 'global' ;
+boundaryConds(2).loadsTimeFact = @(t) t ;
+
+boundaryConds(2).loadsBaseVals = 1e4*[ .001 0 0 0 -1 0 ] ;	% Fx = -12.8e3 (N), My = 12.8e3*ecc (Nmm)
+																% Note that LBA Load is: Pcr = 12.8e3 (N)
+
+%Nodal Coordinates Definitions. X(mm)	Y(mm)	Z(mm)
+
+anguloTotalArco = (180-fiApoyo) - fiApoyo
+deltaTheta = anguloTotalArco / ( 2*nEM )
+
+for i=1:(2*nEM+1)
+	anguloi = (180-fiApoyo) - (i-1)*deltaTheta ;
+	mesh.nodesCoords(i,:) = [  R*cos( pi / 180 * anguloi )  0 R*sin( pi / 180 * anguloi )   ] ; 	% Coords of Nodes in vertical member, incl corner node #(m+1)
+end
+
+largo = 0 ;
+for j=1:2*nEM
+  largo = largo + norm( mesh.nodesCoords(j+1,:) - mesh.nodesCoords(j,:) );
+end
+perimArcoTeo = R*anguloTotalArco*pi/180 ;
+
+%Conectivity Definitions
+mesh.conecCell = cell(5,1) ;
+         					   %M E B I / Node
+mesh.conecCell{ 1 } = [ 0 1 1 0  1   ] ; 	% Node at coord (0,0,0)
+mesh.conecCell{ 2 } = [ 0 1 1 0  2*nEM+1 ] ;	% Node at coord (L,0,L)
+mesh.conecCell{ 3 } = [ 0 1 2 0  nEM+1   ] ;	%
+
+for j=1:2*nEM					   %M E B I / Nodes
+	mesh.conecCell{ 4+j-1 } = [ 1 2 0 0  j j+1 ] ; % frame finite elements
+end
+
+initialConds = [];	% no initial conditions for static analysis
+
+%Static Analysis Parameter Definitions
+analysisSettings.methodName    = 'arcLength' ;
+analysisSettings.deltaT        =   0.02  ;
+analysisSettings.finalTime      =   1   ;
+analysisSettings.stopTolDeltau =   1e-8 ;
+analysisSettings.stopTolForces =   1e-8 ;
+analysisSettings.stopTolIts    =   30   ;
+
+analysisSettings.incremArcLen = .02 ;
+analysisSettings.iniDeltaLamb = boundaryConds(2).loadsTimeFact(1)/100 ;
+analysisSettings.posVariableLoadBC = 2 ;
+
+otherParams.problemName = 'circularArch_AL';
+otherParams.plotsFormat = 'vtk' ;
+
+%Analysis case 1: Solution of Stable Branch with ArcLength
+[matUs, loadFactorsMat] = ONSAS( materials, elements, boundaryConds, initialConds, mesh, analysisSettings, otherParams ) ;
+controlDispsALstab =  -matUs(6*nEM+5,:) ; % rotation wrt Y-axis at node m+1
+loadFactorsALstab  =  loadFactorsMat(:,2) ;
+
+return
+%Plot Load Displacement Curves
+lw = 1.0 ; ms = 5 ; plotfontsize = 12 ;
+figure
+plot( controlDispsALstab, loadFactorsALstab, 'k-' , 'linewidth', lw,'markersize',ms )
+hold on
+plot( controlDispsALunstab, loadFactorsALunstab, 'r-' , 'linewidth', lw,'markersize',ms )
+hold off
+
+labx = xlabel('rotation @Corner node (rad)');   laby = ylabel('\lambda(t)') ;
+legend( 'stable branch', 'unstable branch','location','southeast') ;
+set(gca, 'linewidth', 1.0, 'fontsize', plotfontsize ) ;
+set(labx, 'FontSize', plotfontsize); set(laby, 'FontSize', plotfontsize) ;
+title('Roorda Frame / Load-Displacement curves / With imperfectons') ;
+grid on ;
+print('output/RoordaFrame.png','-dpng')

examples/RoordaFrame_Beams/onsasExample_RoordaFrame.m → examples/roordaFrame/onsasExample_RoordaFrame.m

@@ -1,19 +1,26 @@
 % Static Arclength Analysis of Roorda Frame
 % 2021.10.06
-
 close all, clear all ; addpath( genpath( [ pwd '/../../src'] ) );
+
 % structural parameters
-E = 210e3 ;		% Young's modulus (MPa)
-nu = 0.3 ;		% Poisson ratio
-a = 50; 		% Rectangular section width (mm)
-b = 10; 		% Rectangular section height (mm)
+E = 210e3 ;	% Young's modulus (MPa)
+nu = 0.3  ;	% Poisson ratio
+a = 50 ; 		% Rectangular section width (mm)
+b = 10 ; 		% Rectangular section height (mm)
 A = a*b; 		% Rectangular section Area (mm2)
-L = 1000 ; 		% Frame members length (mm)
-ecc = 1; 		% Point Load eccentricity (mm)
+L = 1000 ;	% Frame members length (mm)
 m = 4; 			% number of finite elements per member
 
+vececc = [ 1 .25 ] ;
+vecpAL = [ 2 .5 ] ;
+
+figure, grid on, hold on
+
+for indecc = 1:length(vececc)
+  ecc = vececc(indecc);
+
 %Material Definitions
-materials(1).hyperElasModel  = 'hiperelastic' ;
+materials(1).hyperElasModel  = '1DrotEngStrain' ;
 materials(1).hyperElasParams = [ E nu ] ;
 
 %Element Definitions
@@ -29,7 +36,7 @@
 boundaryConds(2).loadsCoordSys = 'global' ;
 boundaryConds(2).loadsTimeFact = @(t) t ;
 
-boundaryConds(2).loadsBaseVals = 12.8e3*[ 0 0 0 ecc -1 0 ] ;	% Fx = -12.8e3 (N), My = 12.8e3*ecc (Nmm)
+boundaryConds(2).loadsBaseVals = 12.8e3*[ 0 0 0 ecc*1 -1 0 ] ;	% Fx = -12.8e3 (N), My = 12.8e3*ecc (Nmm)
 																% Note that LBA Load is: Pcr = 12.8e3 (N)
 
 %Nodal Coordinates Definitions. X(mm)	Y(mm)	Z(mm)
@@ -43,7 +50,7 @@
 	mesh.nodesCoords(i,:) = [   (i-m-1)*dl 0  L  ] ;  	% Coords of Nodes in horizontal member, excl corner node #(m+1)
 end
 
-%Conectivity Definitions										   
+%Conectivity Definitions
 mesh.conecCell = cell(5,1) ;
 					   %M E B I / Node
 mesh.conecCell{ 1 } = [ 0 1 1 0  1   ] ; 	% Node at coord (0,0,0)
@@ -65,7 +72,7 @@
 analysisSettings.stopTolIts    =   30   ;
 
 analysisSettings.finalTime     = 1.1 ;
-analysisSettings.incremArcLen = 2 ;
+analysisSettings.incremArcLen = vecpAL(indecc) ;
 analysisSettings.iniDeltaLamb = boundaryConds(2).loadsTimeFact(1)/100 ;
 analysisSettings.posVariableLoadBC = 2 ;
 
@@ -89,16 +96,17 @@
 
 %Plot Load Displacement Curves
 lw = 1.0 ; ms = 5 ; plotfontsize = 12 ;
-figure
-plot( controlDispsALstab, loadFactorsALstab, 'k-' , 'linewidth', lw,'markersize',ms )
-hold on
-plot( controlDispsALunstab, loadFactorsALunstab, 'r-' , 'linewidth', lw,'markersize',ms )
-hold off
+
+plot( controlDispsALstab, loadFactorsALstab, 'k-x' , 'linewidth', lw,'markersize',ms )
+plot( controlDispsALunstab, loadFactorsALunstab, 'r-o' , 'linewidth', lw,'markersize',ms )
 
 labx = xlabel('rotation @Corner node (rad)');   laby = ylabel('\lambda(t)') ;
 legend( 'stable branch', 'unstable branch','location','southeast') ;
 set(gca, 'linewidth', 1.0, 'fontsize', plotfontsize ) ;
 set(labx, 'FontSize', plotfontsize); set(laby, 'FontSize', plotfontsize) ;
 title('Roorda Frame / Load-Displacement curves / With imperfectons') ;
 grid on ;
-print('output/RoordaFrame.png','-dpng')
+
+end
+
+print('output/RoordaFrame.png','-dpng')

src/assembler.m

@@ -116,7 +116,7 @@
 
 			[ Finte, Ke ] = linearStiffMatBeam3D(elemNodesxyzRefCoords, elemTypeGeometry, density, hyperElasParams, elemDisps ) ;
 
-		else
+		elseif strcmp( hyperElasModel, '1DrotEngStrain')
 
       [ fs, ks, stressElem ] = elementBeamForces( elemNodesxyzRefCoords, elemTypeGeometry, [ 1 hyperElasParams ], u2ElemDisps( Ut       , dofselem ) , ...
                                                u2ElemDisps( Udott    , dofselem ) , ...
@@ -126,7 +126,8 @@
       if density > 0
         Fmase = fs{3} ;  Ce    = ks{2} ;   Mmase = ks{3} ;
       end
-
+    else
+      error('wrong hyperElasModel for frame element.')
 		end
 
   % ---------  triangle solid element -----------------------------

test/runTestProblems_local.m

@@ -6,7 +6,7 @@
 
 keyfiles = { 'staticVonMisesTruss/onsasExample_staticVonMisesTruss.m';
             'uniformCurvatureCantilever/onsasExample_uniformCurvatureCantilever.m' ;
-            'linearPlaneStrainExample/onsasExample_linearPlaneStrain.m' ;
+            'linearPlaneStrain/onsasExample_linearPlaneStrain.m' ;
             'uniaxialExtension/onsasExample_uniaxialExtension.m'  } ;
 
 current  = 1 ;   verifBoolean = 1 ;  testDir = pwd ;