Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Implement Neo-Hookean + add example uniaxial compression solid #128

Merged
merged 1 commit into from
Aug 14, 2020
Merged

Implement Neo-Hookean + add example uniaxial compression solid #128

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
101 changes: 101 additions & 0 deletions examples/onsasExample_uniaxialCompression_NH.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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) ;
% --------------------------------------------------------
66 changes: 33 additions & 33 deletions sources/assembler.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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.
%
Expand Down Expand Up @@ -33,15 +33,15 @@
% -----------------------------------------------
if booleanCppAssembler
CppAssembly

% -----------------------------------------------
% --- octave/matlab assembler ---
% -----------------------------------------------
else
% -----------------------------------------------
nElems = size(Conec, 1) ; nNodes = length( Ut ) / 6 ;


% ====================================================================
% --- 1 declarations ---
% ====================================================================
Expand All @@ -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
Expand All @@ -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
Expand All @@ -85,68 +85,68 @@
% 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

[ fs, ks, stress ] = elementBeamForces( elemCoords, elemCrossSecParams, elemConstitutiveParams, solutionMethod, u2ElemDisps( Ut , dofselem ) , ...
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
% -------------------------------------------
Expand All @@ -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, : )' ;
Expand All @@ -179,10 +179,10 @@
valsC( entriesSparseStorVecs ) = Ce ( indRow, : )' ;
end
end

counterInds = counterInds + length( dofselemRed ) ;
end

case 3
StressVec( elem, (1:length(stress) ) ) = stress ;

Expand All @@ -206,7 +206,7 @@
dampingMat(1:2:end) = nodalDispDamping ;
dampingMat(2:2:end) = nodalDispDamping * 0.01 ;
end

switch paramOut

case 1,
Expand All @@ -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) ) ;
Expand All @@ -240,15 +240,15 @@
M = sparse(size( K ) ) ;
C = sparse(size( K ) ) ;
end

Assembled{2} = C ;
Assembled{3} = M ;

case 3
Assembled{1} = StressVec ;

end

end % if booleanCppAssembler
% ----------------------------------------

Expand All @@ -261,7 +261,7 @@
%
% ==============================================================================

function nodesmat = conv ( conec, coordsElemsMat )
function nodesmat = conv ( conec, coordsElemsMat )
nodesmat = [] ;
nodesread = [] ;

Expand Down
26 changes: 26 additions & 0 deletions sources/cosseratNH.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
Loading