From b5c2821602ef324afba2fc0e849ac9e4329a8da1 Mon Sep 17 00:00:00 2001 From: mforets Date: Fri, 14 Aug 2020 19:03:01 -0300 Subject: [PATCH] agregar ejemplo Neo-Hookean --- .../onsasExample_uniaxialCompression_NH.m | 101 ++++++++++++++++++ sources/assembler.m | 66 ++++++------ sources/cosseratNH.m | 26 +++++ sources/elementTetraSolid.m | 50 +++++---- 4 files changed, 184 insertions(+), 59 deletions(-) create mode 100644 examples/onsasExample_uniaxialCompression_NH.m create mode 100644 sources/cosseratNH.m diff --git a/examples/onsasExample_uniaxialCompression_NH.m b/examples/onsasExample_uniaxialCompression_NH.m new file mode 100644 index 000000000..f9cba30a0 --- /dev/null +++ b/examples/onsasExample_uniaxialCompression_NH.m @@ -0,0 +1,101 @@ +%% Example uniaxialSolid +% Linear elastic solid submitted to uniaxial loading. +% Geometry given by $L_x$, $L_y$ and $L_z$, tension $p$ applied on +% face $x=L_x$. +% We use Neo-Hookean constitutive law. The solid is in a compression regime. +%% + +% uncomment to delete variables and close windows +clear all, close all + +%% General data +dirOnsas = [ pwd '/..' ] ; +problemName = 'uniaxialCompression_Manual' ; + +%% Structural properties + +% compression applied and x, y, z dimensions +p = -3 ; Lx = 1 ; Ly = 1 ; Lz = 1 ; + +% an 8-node mesh is considered with its connectivity matrix +Nodes = [ 0 0 0 ; ... + 0 0 Lz ; ... + 0 Ly Lz ; ... + 0 Ly 0 ; ... + Lx 0 0 ; ... + Lx 0 Lz ; ... + Lx Ly Lz ; ... + Lx Ly 0 ] ; + +Conec = {[ 0 1 1 0 0 5 8 6 ]; ... % loaded face + [ 0 1 1 0 0 6 8 7 ]; ... % loaded face + [ 0 1 0 0 1 4 1 2 ]; ... % x=0 supp face + [ 0 1 0 0 1 4 2 3 ]; ... % x=0 supp face + [ 0 1 0 0 2 6 2 1 ]; ... % y=0 supp face + [ 0 1 0 0 2 6 1 5 ]; ... % y=0 supp face + [ 0 1 0 0 3 1 4 5 ]; ... % z=0 supp face + [ 0 1 0 0 3 4 8 5 ]; ... % z=0 supp face + [ 1 2 0 0 0 1 4 2 6 ]; ... % tetrahedron + [ 1 2 0 0 0 6 2 3 4 ]; ... % tetrahedron + [ 1 2 0 0 0 4 3 6 7 ]; ... % tetrahedron + [ 1 2 0 0 0 4 1 5 6 ]; ... % tetrahedron + [ 1 2 0 0 0 4 6 5 8 ]; ... % tetrahedron + [ 1 2 0 0 0 4 7 6 8 ] ... % tetrahedron + } ; + + +% ====================================================================== +% --- MELCS parameters --- + +materialsParams = cell(1,1) ; % M +elementsParams = cell(1,1) ; % E +loadsParams = cell(1,1) ; % L +crossSecsParams = cell(1,1) ; % C +springsParams = cell(1,1) ; % S + +% --- Material parameters --- +E = 1 ; nu = 0.3 ; +materialsParams{1} = [ 0 3 E nu ] ; + +% --- Element parameters --- +elementsParams{1,1} = [ 5 ] ; +elementsParams{2,1} = [ 4 1 ] ; # the second index indicates that we take the complex-step computation expression + +% --- Load parameters --- +loadsParams{1,1} = [ 1 1 p 0 0 0 0 0 ] ; + +% --- CrossSection parameters --- + +% ---------------------------------------------------------------------- +% --- springsAndSupports parameters --- +springsParams{1, 1} = [ inf 0 0 0 0 0 ] ; +springsParams{2, 1} = [ 0 0 inf 0 0 0 ] ; +springsParams{3, 1} = [ 0 0 0 0 inf 0 ] ; + +% ====================================================================== + + +%% --- Analysis parameters --- +stopTolIts = 30 ; +stopTolDeltau = 1.0e-12 ; +stopTolForces = 1.0e-12 ; +targetLoadFactr = 1 ; +nLoadSteps = 10 ; + +numericalMethodParams = [ 1 stopTolDeltau stopTolForces stopTolIts ... + targetLoadFactr nLoadSteps ] ; + +controlDofs = [ 7 1 1 ] ; + +%% Output parameters +plotParamsVector = [ 3 ] ; +printflag = 2 ; + +reportBoolean = 0; + +% --- Analytic sol --- +analyticSolFlag = 0 ; % not available + +%% run ONSAS +acdir = pwd ; cd(dirOnsas); ONSAS, cd(acdir) ; +% -------------------------------------------------------- diff --git a/sources/assembler.m b/sources/assembler.m index 9ca88cd92..89884227d 100644 --- a/sources/assembler.m +++ b/sources/assembler.m @@ -1,4 +1,4 @@ -% Copyright (C) 2019, Jorge M. Perez Zerpa, J. Bruno Bazzano, Jean-Marc Battini, Joaquin Viera, Mauricio Vanzulli +% Copyright (C) 2019, Jorge M. Perez Zerpa, J. Bruno Bazzano, Jean-Marc Battini, Joaquin Viera, Mauricio Vanzulli % % This file is part of ONSAS. % @@ -33,15 +33,15 @@ % ----------------------------------------------- if booleanCppAssembler CppAssembly - + % ----------------------------------------------- % --- octave/matlab assembler --- % ----------------------------------------------- else % ----------------------------------------------- nElems = size(Conec, 1) ; nNodes = length( Ut ) / 6 ; - - + + % ==================================================================== % --- 1 declarations --- % ==================================================================== @@ -53,7 +53,7 @@ Fint = zeros( nNodes*6 , 1 ) ; Fmas = zeros( nNodes*6 , 1 ) ; Fvis = zeros( nNodes*6 , 1 ) ; - + % ------- tangent matrix ------------------------------------- case 2 %~ if octaveBoolean @@ -66,7 +66,7 @@ valsC = zeros( nElems*24*24, 1 ) ; valsM = zeros( nElems*24*24, 1 ) ; - counterInds = 0 ; % counter non-zero indexes + counterInds = 0 ; % counter non-zero indexes % ------- matrix with stress per element ---------------------------- case 3 @@ -85,50 +85,50 @@ % extract element properties elemMaterialParams = materialsParamsMat( Conec( elem, 5) , : ) ; elemrho = elemMaterialParams( 1 ) ; - elemConstitutiveParams = elemMaterialParams( 2:end ) ; + elemConstitutiveParams = elemMaterialParams( 2:end ) ; elemElementParams = elementsParamsMat( Conec( elem, 6),: ) ; typeElem = elemElementParams(1) ; [numNodes, dofsStep] = elementTypeInfo ( typeElem ) ; - + % obtains nodes and dofs of element nodeselem = Conec( elem, 1:numNodes )' ; dofselem = nodes2dofs( nodeselem , 6 ) ; dofselemRed = dofselem ( 1 : dofsStep : end ) ; elemDisps = u2ElemDisps( Ut , dofselemRed ) ; - - elemCoords = coordsElemsMat( elem, 1:dofsStep:( numNodes*6 ) ) ; - + + elemCoords = coordsElemsMat( elem, 1:dofsStep:( numNodes*6 ) ) ; + if typeElem == 2 || typeElem == 3 elemCrossSecParams = crossSecsParamsMat ( Conec( elem, 4+4 ) , : ) ; end stress = [] ; - + switch typeElem % ----------- truss element ------------------------------------ case 2 [ fs, ks, stress ] = elementTrussInternForce( elemCoords, elemDisps, elemConstitutiveParams, elemCrossSecParams, paramOut ) ; - + Finte = fs{1} ; Ke = ks{1} ; - + if solutionMethod > 2 booleanConsistentMassMat = elemElementParams(2) ; - + dotdotdispsElem = u2ElemDisps( Udotdott , dofselem ) ; [ Fmase, Mmase ] = elementTrussMassForce( elemCoords, elemrho, elemCrossSecParams(1), elemElementParams(2), paramOut, dotdotdispsElem ) ; % %~ Fmase = fs{3} ; Ce = ks{2} ; Mmase = ks{3} ; %~ Fmase = fs{3} ; Ce = ks{2} ; Mmase = ks{3} ; end - - + + % ----------- frame element ------------------------------------ case 3 @@ -136,17 +136,17 @@ u2ElemDisps( Udott , dofselem ) , ... u2ElemDisps( Udotdott , dofselem ), elemrho ) ; Finte = fs{1} ; Ke = ks{1} ; - + if solutionMethod > 2 Fmase = fs{3} ; Ce = ks{2} ; Mmase = ks{3} ; end - - + + % --------- tetrahedron solid element ----------------------------- case 4 - + [ Finte, Ke, stress ] = elementTetraSolid( elemCoords, elemDisps, ... - elemConstitutiveParams, paramOut ) ; + elemConstitutiveParams, paramOut, elemElementParams(2)) ; end % case tipo elemento % ------------------------------------------- @@ -164,11 +164,11 @@ end case 2 - + for indRow = 1:length( dofselemRed ) - - entriesSparseStorVecs = counterInds + (1:length( dofselemRed) ) ; - + + entriesSparseStorVecs = counterInds + (1:length( dofselemRed) ) ; + indsIK ( entriesSparseStorVecs ) = dofselemRed( indRow ) ; indsJK ( entriesSparseStorVecs ) = dofselemRed ; valsK ( entriesSparseStorVecs ) = Ke( indRow, : )' ; @@ -179,10 +179,10 @@ valsC( entriesSparseStorVecs ) = Ce ( indRow, : )' ; end end - + counterInds = counterInds + length( dofselemRed ) ; end - + case 3 StressVec( elem, (1:length(stress) ) ) = stress ; @@ -206,7 +206,7 @@ dampingMat(1:2:end) = nodalDispDamping ; dampingMat(2:2:end) = nodalDispDamping * 0.01 ; end - + switch paramOut case 1, @@ -231,7 +231,7 @@ Assembled{1} = K ; - if solutionMethod > 2 + if solutionMethod > 2 valsM = valsM (1:counterInds) ; valsC = valsC (1:counterInds) ; M = sparse( indsIK, indsJK, valsM , size(KS,1), size(KS,1) ) ; @@ -240,7 +240,7 @@ M = sparse(size( K ) ) ; C = sparse(size( K ) ) ; end - + Assembled{2} = C ; Assembled{3} = M ; @@ -248,7 +248,7 @@ Assembled{1} = StressVec ; end - + end % if booleanCppAssembler % ---------------------------------------- @@ -261,7 +261,7 @@ % % ============================================================================== -function nodesmat = conv ( conec, coordsElemsMat ) +function nodesmat = conv ( conec, coordsElemsMat ) nodesmat = [] ; nodesread = [] ; diff --git a/sources/cosseratNH.m b/sources/cosseratNH.m new file mode 100644 index 000000000..5c57eb5cc --- /dev/null +++ b/sources/cosseratNH.m @@ -0,0 +1,26 @@ +function [S, ConsMat] = cosseratNH( consParams, Egreen, consMatFlag) + +young = consParams(1) ; +nu = consParams(2) ; + +lambda = young * nu / ( (1 + nu) * (1 - 2*nu) ) ; +shear = young / ( 2 * (1 + nu) ) ; + +C = 2*Egreen + eye(3); % Egreen = 1/2 (C - I) +invC = inv(C); +detC = det(C); % TODO use analyDet ? +J = sqrt(detC); +S = shear * (eye(3) - invC) + lambda * log(J) * invC; + +if consMatFlag == 0 % only stress computed + ConsMat = [] ; + +elseif consMatFlag == 1 % complex-step computation expression + + ConsMat = zeros(6,6); + ConsMat = complexStepConsMat( 'cosseratNH', consParams, Egreen ) ; + +else + error("the analytical expression for the Neo-Hookean constitutive law is not available") + +end diff --git a/sources/elementTetraSolid.m b/sources/elementTetraSolid.m index 617b163c4..26efefc99 100644 --- a/sources/elementTetraSolid.m +++ b/sources/elementTetraSolid.m @@ -1,4 +1,4 @@ -% Copyright (C) 2019, Jorge M. Perez Zerpa, J. Bruno Bazzano, Jean-Marc Battini, Joaquin Viera, Mauricio Vanzulli +% Copyright (C) 2019, Jorge M. Perez Zerpa, J. Bruno Bazzano, Jean-Marc Battini, Joaquin Viera, Mauricio Vanzulli % % This file is part of ONSAS. % @@ -19,7 +19,7 @@ % function [ Finte, KTe, stress ] = elementTetraSolid( ... - elemCoords, elemDisps, elemConstitutiveParams, paramOut ) + elemCoords, elemDisps, elemConstitutiveParams, paramOut, consMatFlag ) Finte = zeros(12,1) ; booleanKTAnalytic = 1 ; @@ -36,7 +36,7 @@ eleCoordSpa = tetCoordMat + eleDispsMat ; xi = 0.25 ; wi = 1/6 ; - + % matriz de derivadas de fun forma respecto a coordenadas isoparametricas deriv = shapeFuns( xi, xi , xi , 1 ) ; @@ -47,13 +47,13 @@ if vol<0, vol, error('Element with negative volume, check connectivity.'), end - + %~ E = params(1) ; %~ nu = params(2) ; %~ mu = E / (2.0 * ( 1.0 + nu ) ) ; %~ Bulk = E / (3.0 * ( 1.0 - 2.0 * nu ) ) ; - + %~ eyetres = eye(3) ; %~ eyevoig = zeros(6,1) ; %~ eyevoig(1:3) = 1.0 ; @@ -61,43 +61,41 @@ %~ ConsMat = zeros(6,6) ; %~ ConsMat (1:3,1:3) = Bulk + 2* mu * ( eyetres - 1.0/3.0 ) ; %~ ConsMat (4:6,4:6) = mu * eyetres ; - + %~ Kml = BMat' * ConsMat * BMat * tetVol ; - + %~ KTe([1:2:end], [1:2:end]) = Kml ; %~ strain = BMat * Ue ; %~ stress = ConsMat * strain ; %~ Fint = stress' * BMat * tetVol ; - + %~ Finte(1:2:end) = Fint ; funder = inv(jacobianmat)' * deriv ; - + H = eleDispsMat * funder' ; F = H + eye(3) ; Egreen = 0.5 * ( H + transpose( H ) + transpose( H ) * H ) ; - global consMatFlag - if elemConstitutiveParams(1) == 2 % Saint-Venant-Kirchhoff compressible solid - + [ S, ConsMat ] = cosseratSVK( elemConstitutiveParams(2:3), Egreen, consMatFlag ) ; elseif elemConstitutiveParams(1) == 3 % Neo-Hookean Compressible [ S, ConsMat ] = cosseratNH ( elemConstitutiveParams(2:3), Egreen, consMatFlag ) ; end - + matBgrande = BgrandeMats ( funder , F ) ; - + Svoigt = mat2voigt( S, 1 ) ; - + Finte = transpose(matBgrande) * Svoigt * vol ; - + strain = zeros(6,1); stress = Svoigt ; @@ -116,7 +114,7 @@ Kgl( (i-1)*3+3 , (j-1)*3+3 ) = matauxgeom(i,j); end end - + KTe = Kml + Kgl ; @@ -135,16 +133,16 @@ matBgrande = zeros(6, 12) ; matBgrande(1:3,:) = [ diag( deriv(:,1) )*F' diag( deriv(:,2) )*F' diag( deriv(:,3) )*F' diag( deriv(:,4) )*F' ]; - + for k=1:4 %~ for i=1:3 % fila %~ for j=1:3 % columna %~ matBgrande ( i , (k-1)*3 + j ) = deriv(i,k) * F(j,i) ; %~ end - %~ end + %~ end + - %~ for j=1:3 %~ matBgrande ( 4:6 , (k-1)*3 + j ) = [ deriv(2,k)*F(j,3)+deriv(3,k)*F(j,2) ; ... %~ deriv(1,k)*F(j,3)+deriv(3,k)*F(j,1) ; ... @@ -166,7 +164,7 @@ function BMat = BMats ( deriv ) BMat = zeros(6,12) ; - + for k = 1:4 for i = 1:3 @@ -181,7 +179,7 @@ BMat ( 6 , (k-1)*3 + 1 ) = deriv(2,k) ; BMat ( 6 , (k-1)*3 + 2 ) = deriv(1,k) ; - + end @@ -193,13 +191,13 @@ else derivOrder = varargin(1){1} ; end - + A = zeros(4,4) ; A(:,1) = 1.0 ; A(:,2:4) = tetcoordmat' ; - + invA = inv(A) ; - + vol = det(A) / 6.0 ; - + fun = invA(2:4,:) ;