Merge pull request #312 from ONSAS/jorge

agrego ejemplos von mises y springmass

.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/linearPlaneStrain examples/uniformCurvatureCantilever examples/uniaxialExtension examples/nonlinearPendulum src src/elements src/vtk src/mesh
+          src: examples/staticVonMisesTruss examples/springMass examples/linearPlaneStrain examples/uniformCurvatureCantilever examples/uniaxialExtension examples/nonlinearPendulum src src/elements src/vtk src/mesh
           tests: ./test/runTestProblems_moxunit_disp.m
   # tests-disp-coverage:
   #   runs-on: [ubuntu-latest]

examples/dynamicVonMises/centralDiffDynVonMises.m

@@ -0,0 +1,56 @@
+
+function [u, normalForce, t] = centralDiffDynVonMises(rho, Lx, L0, Lc, Ic, Ac, E, kc, m, c, g, tf, dt)
+
+  mb = L0*Ac*rho ;
+  Lz = sqrt( L0^2 - Lx^2 ); %m
+
+%% Defino Vector de fuerzas Internas: fint(u)  - u = [u1 , u2]^T
+Fint = @(u) E*Ac*L0*(u(1)^2+u(2)^2-2*Lx*u(1)+2*Lz*u(2))/2/L0^4*[-Lx+u(1);Lz+u(2)]+[kc*u(1);0];
+
+%% Defino Vector de fuerzas Externas: gravedad
+ft = @(t) [0;-(m+mb)/2*g]; %N
+
+%% Defino Matriz de Masa Concentrada
+M = [mb 0 ; 0 (mb+m)/2] ;
+
+%% Defino Matriz de Amortiguamiento
+C = [c/10 0 ; 0 c];
+
+%% Defino Condiciones Iniciales
+t0 = 0;
+u0 = [0;0];
+v0 = [0;0];
+ac0 = M\(ft(t0)-C*v0-Fint(u0)); % de ec de movimiento Mu.. + Fint(u) = ft
+
+% Inicializacion Difrerencia Centrada
+
+a0 = 1/dt^2; a1=1/2/dt; a2=2*a0; a3=1/a2;
+
+nTimes = tf/dt
+
+Meff = a0*M+a1*C;
+M2 = a0*M-a1*C;
+
+% Comienza Marcha en el Tiempo usando Diferencia Centrada
+
+k = 2 ; % index at which equilibrium is done. in this case t=0
+
+epsg = @(u) (u(1)^2+2*Lz*u(2)-2*Lx*u(1)+u(2)^2)/2/L0^2;
+
+u      = zeros(2, nTimes+2) ;
+acc    = zeros(2, nTimes+2) ;
+vel    = zeros(2, nTimes+2) ;
+normalForce = zeros(1, nTimes+2) ;
+t      = -dt:dt:tf ;
+
+u(:,k-1) = u0 - dt*v0 + a3*ac0; % solution u(-dt) at time -dt
+u(:,k  ) = u0; % u(0) %
+
+while k<=(nTimes+1)
+    feff      = ft(t(k)) - Fint(u(:,k)) +a2*M*u(:,k) - M2*u(:,k-1) ;
+    u(:,k+1) = Meff\feff ; % sets solution at t+dt
+    acc(:,k) = a0*(u(:,k+1)-2*u(:,k)+u(:,k-1)); % computers acc and vel at t
+    vel(:,k) = a1*(u(:,k+1)-u(:,k-1));
+    normalForce(k+1) = E * Ac * epsg( u(:,k+1 ) ) ;
+    k=k+1;
+end

examples/dynamicVonMises/onsasExample_dynamicVonMises.m

@@ -0,0 +1,124 @@
+%md
+%mdBefore defining the structs, the workspace is cleaned, the ONSAS directory is added to the path and scalar auxiliar parameters are defined.
+close all, clear all ; addpath( genpath( [ pwd '/../../src'] ) );
+% scalar parameters
+
+%%% Parametros Estructura
+rho = 7850; % kg/m3 (acero)
+Lx  = .374/2;
+L0  = .205 ;
+
+Lc  = .240; %m
+Ic  = .0254*.0032^3/12; %m4
+Ac  = .0254*.0032; %m2
+E   = 200000e6 %Pa (acero)
+kc  = 3*E*Ic/Lc^3; %N/m
+mConc   = 1.4; %kg Pandeo incipiente en 1.4
+c   = 0; %kg/s (amortiguamiento por friccion juntas y arrastre pesa)
+g   = 9.81; %m/s2
+
+tf  = 1.0;
+dt  = .000025; % sec
+
+[u, normalForce, times ] = centralDiffDynVonMises(rho, Lx, L0, Lc, Ic, Ac, E, kc, mConc, c, g, tf, dt);
+
+mb = L0*Ac*rho; %kg
+Lz = sqrt( L0^2 - Lx^2 ); %m
+
+rhoBarraMasa = mConc*.5 / (Lz*Ac);
+
+M = [mb 0 ; 0 (mb+mConc)/2]
+
+
+nu = .3 ;
+materials(1).hyperElasModel  = '1DrotEngStrain' ;
+materials(1).hyperElasParams = [ E nu ] ;
+materials(1).density =  rho ;
+
+materials(2).hyperElasModel  = '1DrotEngStrain' ;
+materials(2).hyperElasParams = [ 1e10*E nu ] ;
+materials(2).density = rhoBarraMasa ;
+
+materials(3).hyperElasModel  = '1DrotEngStrain' ;
+materials(3).hyperElasParams = [ kc*L0/Ac nu ] ;
+materials(3).density = rho ;
+
+
+elements(1).elemType = 'node' ;
+
+elements(2).elemType = 'truss';
+elements(2).elemTypeGeometry = [2 sqrt(Ac) sqrt(Ac) ] ;
+elements(2).elemTypeParams = 0 ;
+
+boundaryConds(1).imposDispDofs = [ 3 5 ] ;
+boundaryConds(1).imposDispVals = [ 0 0 ] ;
+
+boundaryConds(2).imposDispDofs =  [ 1 3 ] ;
+boundaryConds(2).imposDispVals =  [ 0 0 ] ;
+boundaryConds(2).loadsCoordSys = 'global'                  ;
+boundaryConds(2).loadsTimeFact = @(t) 1                    ;
+boundaryConds(2).loadsBaseVals = [ 0 0 0 0 -(mConc+mb)/2*g 0 ] ;
+
+boundaryConds(3).imposDispDofs =  [ 1 3 ] ;
+boundaryConds(3).imposDispVals =  [ 0 0 ] ;
+
+boundaryConds(4).imposDispDofs =  [ 1 3 5 ] ;
+boundaryConds(4).imposDispVals =  [ 0 0 0 ] ;
+
+initialConds                = struct() ;
+
+mesh.nodesCoords = [   0  0   0   ; ...
+                      Lx  0  Lz   ; ...
+                      Lx  0  Lz-Lz; ...
+                      -L0  0  0    ] ;
+
+mesh.conecCell = { } ;
+mesh.conecCell{ 1, 1 } = [ 0 1 1 0   1   ] ;
+mesh.conecCell{ 2, 1 } = [ 0 1 2 0   2   ] ;
+mesh.conecCell{ 3, 1 } = [ 0 1 3 0   3   ] ;
+mesh.conecCell{ 4, 1 } = [ 0 1 4 0   4   ] ;
+
+mesh.conecCell{ 4+1, 1 } = [ 1 2 0 0   1 2 ] ;
+mesh.conecCell{ 4+2, 1 } = [ 2 2 0 0   2 3 ] ;
+mesh.conecCell{ 4+3, 1 } = [ 3 2 0 0   1 4 ] ;
+
+analysisSettings.methodName    = 'newmark' ;
+%md and the following parameters correspond to the iterative numerical analysis settings
+analysisSettings.deltaT        =   0.0005  ;
+analysisSettings.finalTime      =   1    ;
+analysisSettings.stopTolDeltau =   1e-8 ;
+analysisSettings.stopTolForces =   1e-8 ;
+analysisSettings.stopTolIts    =   10   ;
+analysisSettings.alphaNM      =   0.25   ;
+analysisSettings.deltaNM      =   0.5   ;
+
+%md
+%md### otherParams
+otherParams.problemName = 'dynamicVonMisesTruss';
+otherParams.plotsFormat = 'vtk' ;
+%md
+%md### Analysis case 1: NR with Rotated Eng Strain
+%md In the first case ONSAS is run and the solution at the dof of interest is stored.
+[matUs, loadFactorsMat] = ONSAS( materials, elements, boundaryConds, initialConds, mesh, analysisSettings, otherParams ) ;
+
+uONSAS = [ matUs(1,:) ; matUs(6+5,:) ] ;
+
+timesONSAS    = 0:analysisSettings.deltaT:analysisSettings.finalTime ;
+
+
+subplot(3,1,1)
+plot(times(1:10:end),1000*u(1,1:10:end))
+hold on
+plot(timesONSAS,1000*uONSAS(1,:),'r' )
+xlabel('t [s]'); ylabel('u_1 [mm]');
+axis( [0 2 1e3*min(u(1,:))*1.1 1e3*max(u(1,:))*1.1] );
+subplot(3,1,2)
+hold on
+plot(times(1:10:end),1000*u(2,1:10:end))
+plot(timesONSAS,1000*uONSAS(2,:),'r' )
+xlabel('t [s]'); ylabel('u_2 [mm]');
+axis([0 2 1e3*min(u(2,:))*1.1 1e3*max(u(2,:))*1.1]);
+subplot(3,1,3)
+plot(times(1:10:end),normalForce(1:10:end))
+xlabel('t [s]'); ylabel('Directa [N]')
+axis([0 2 min(normalForce)*1.1 max(normalForce)*1.1]);

examples/springMass/myLoadSpringMass.m

@@ -0,0 +1,8 @@
+function f = myLoadSpringMass( t)
+
+%Force data
+omegaBar = 2   ;
+p0       = 0.1 ;
+
+f=zeros(12,1);
+f(7) = p0 *sin( omegaBar*t) ; 

examples/springMass/onsasExample_springMass.m

@@ -0,0 +1,127 @@
+% ------------------------------------
+% springmass example
+% Notation and analytical based on chapter 3 from
+% Ray W. Clough and Joseph Penzien, Dynamics of Structures, Third Edition, 2003
+% ------------------------------------
+
+close all, clear all ; addpath( genpath( [ pwd '/../../src'] ) );
+
+otherParams.problemName = 'springMass' ;
+otherParams.plotsFormat = 'vtk' ;
+
+% spring mass system
+k        = 39.47 ;
+p0       = 40    ;
+c        = 0.1   ;
+m        = 1     ;
+omegaBar = 4*sqrt(k/m) ;
+p0       = 40    ;
+u0       = 0.0   ; % initial displacement
+
+% parameters for truss model
+l   = 1   ;
+A   = 0.1 ;
+rho = m * 2 / ( A * l ) ;
+E   = k * l /   A       ;
+
+omegaN = sqrt( k / m );
+xi     = c / m  / ( 2 * omegaN ) ;
+nodalDamping = c ;
+
+freq   = omegaN / (2*pi)      ;
+TN     = 2*pi / omegaN        ;
+dtCrit = TN / pi              ;
+
+% numerical method params
+stopTolDeltau = 1e-10           ;
+stopTolForces = 1e-10           ;
+stopTolIts    = 30              ;
+alphaHHT      = 0;
+% ------------------------------------
+
+
+materials(1).hyperElasModel  = '1DrotEngStrain' ;
+materials(1).hyperElasParams = [ E 0 ] ;
+materials(1).density  = rho ;
+
+elements(1).elemType = 'node' ;
+elements(2).elemType = 'truss';
+elements(2).elemTypeGeometry = [2 sqrt(A) sqrt(A) ] ;
+elements(2).elemTypeParams = 0 ;
+
+boundaryConds(1).imposDispDofs =  [ 1 3 5 ] ;
+boundaryConds(1).imposDispVals =  [ 0 0 0 ] ;
+
+boundaryConds(2).imposDispDofs =  [ 3 5 ] ;
+boundaryConds(2).imposDispVals =  [ 0 0 ] ;
+boundaryConds(2).loadsCoordSys = 'global'                  ;
+boundaryConds(2).loadsTimeFact = @(t) p0*sin( omegaBar*t )                    ;
+boundaryConds(2).loadsBaseVals = [ 1 0 0 0 0 0 ] ;
+
+% initial conditions
+initialConds.nonHomogeneousInitialCondU0    = [ 2 1 u0    ] ;
+
+mesh.nodesCoords = [  0  0  0 ; ...
+                      l  0  0 ] ;
+mesh.conecCell = { } ;
+mesh.conecCell{ 1, 1 } = [ 0 1 1 0   1   ] ;
+mesh.conecCell{ 2, 1 } = [ 0 1 2 0   2   ] ;
+mesh.conecCell{ 3, 1 } = [ 1 2 0 0   1 2   ] ;
+
+
+analysisSettings.methodName    = 'newmark' ;
+%md and the following parameters correspond to the iterative numerical analysis settings
+analysisSettings.deltaT        =   0.005  ;
+analysisSettings.finalTime      =   1*2*pi/omegaN   ;
+analysisSettings.stopTolDeltau =   1e-8 ;
+analysisSettings.stopTolForces =   1e-8 ;
+analysisSettings.stopTolIts    =   10   ;
+analysisSettings.alphaNM      =   0.25   ;
+analysisSettings.deltaNM      =   0.5   ;
+
+[matUsNewmark, loadFactorsMat] = ONSAS( materials, elements, boundaryConds, initialConds, mesh, analysisSettings, otherParams ) ;
+
+if (c == 0) && (p0 == 0) % free undamped
+  analyticSolFlag = 1 ;
+  analyticFunc = @(t)   (   u0 * cos( omegaN * t )  ) ;
+  analyticCheckTolerance = 2e-1 ;
+else
+  beta   = omegaBar / omegaN ;
+  omegaD = omegaN * sqrt( 1-xi^2 ) ;
+
+  G1 = (p0/k) * ( -2 * xi * beta   ) / ( ( 1 - beta^2 )^2 + ( 2 * xi * beta )^2 ) ;
+  G2 = (p0/k) * (  1      - beta^2 ) / ( ( 1 - beta^2 )^2 + ( 2 * xi * beta )^2 ) ;
+  if u0 < l
+    A  = u0 - G1 ;
+    B  =  (xi*omegaN*A - omegaBar*G2 ) / (omegaD);
+  else
+    error('this analytical solution is not valid for this u0 and l0');
+  end
+  analyticSolFlag = 1 ;
+  analyticFunc = @(t) ...
+     ( A * cos( omegaD   * t ) + B  * sin( omegaD   * t ) ) .* exp( -xi * omegaN * t ) ...
+    + G1 * cos( omegaBar * t ) + G2 * sin( omegaBar * t ) ;
+    analyticCheckTolerance = 5e-2 ;
+end
+
+times = 0:analysisSettings.deltaT:(analysisSettings.finalTime+analysisSettings.deltaT) ;
+
+valsAnaly = analyticFunc(times) ;
+valsNewmark = matUsNewmark(6+1,:) ;
+
+analysisSettings.methodName    = 'alphaHHT' ;
+analysisSettings.alphaHHT      =   0   ;
+
+[matUsHHT, loadFactorsMat] = ONSAS( materials, elements, boundaryConds, initialConds, mesh, analysisSettings, otherParams ) ;
+
+valsHHT = matUsHHT(6+1,:) ;
+
+verifBooleanNewmark =  ( ( norm( valsAnaly - valsNewmark    ) / norm( valsAnaly   ) ) <  analyticCheckTolerance )
+verifBooleanHHT     =  ( ( norm( valsAnaly - valsNewmark     ) / norm( valsAnaly  ) ) <  analyticCheckTolerance )
+verifBoolean = verifBooleanHHT && verifBooleanNewmark ;
+
+figure
+hold on, grid on
+plot(valsAnaly,'b-x')
+plot(valsNewmark,'r-o')
+plot(valsHHT,'g-s')

src/assembler.m

@@ -63,14 +63,13 @@
 %  --- 2 loop assembly ---
 % ====================================================================
 
-density = materials.density ;
-
 for elem = 1:nElems
   mebiVec = Conec( elem, 1:4) ;
 
   % extract element properties
-  hyperElasModel      = materials(mebiVec(1)).hyperElasModel ;
-  hyperElasParams     = materials(mebiVec(1)).hyperElasParams ;
+  hyperElasModel   = materials(mebiVec(1)).hyperElasModel  ;
+  hyperElasParams  = materials(mebiVec(1)).hyperElasParams ;
+  density          = materials(mebiVec(1)).density         ;
 
   elemType         = elements(mebiVec(2)).elemType         ;
   elemTypeParams   = elements(mebiVec(2)).elemTypeParams   ;

test/runTestProblems_moxunit_disp.m

@@ -28,3 +28,7 @@
 function test_5
   onsasExample_nonlinearPendulum
   assertEqual( verifBoolean, true );
+
+function test_6
+    onsasExample_springMass
+    assertEqual( verifBoolean, true );