Plane strain ring example
In this example a hollow cylinder submitted to an internal pressure $p_i$ as shown in diagram depicted below is considered. The length of the cylinder is $L_z $ m and the internal and external radious are $R_i$ and $R_e$, respectively.
Before defining the structs, the workspace is cleaned, the ONSAS directory is added to the path and scalar geometry and material parameters are defined:
close all, if ~strcmp( getenv('TESTS_RUN'), 'yes'), clear all, end
% add path
addpath( genpath( [ pwd '/../../src'] ) ) ;
% scalar parameters
% E = 1e6 ; nu = 0.3 ; p = 30e3 ; L = .75 ; Re = 0.15 ; Ri = 0.1 ;
E = 210 ; nu = 0.3 ; p = 0.01 ; L = .75 ;
global Re
global Ri
Re = 200 ; Ri = 100 ;
Linear analysis
Analytic solution
The solution displacement field is extracted from chapter 4 of (Timoshenko and Goodier, Theory of Elasticity, 3rd edition). The Navier's equation, imposing no temperature variation, no volumetric forces, and considering a radial dispalcement field leads to:
\[ \nabla (\nabla . \textbf{u}(r,\theta,z) ) = 0 \]
Due to the symmetry of the problem $\mathbf{\mathit{u_{\theta}}} = 0 $ and also $\mathbf{ \mathit{ \textbf{u} (r,\theta,z) } } = \mathbf{ \mathit{ \textbf{u}(r,z) } } $. Thus, according to the boundary conditions stated above $\mathit{u_z(r,z)=0}$ and the radial displacements field $\mathit{u_r(r)}$ only varies with $r$. Thereafter by imposing the boundary conditions stated above and substituting ($E$, $\nu$) into Lamé parameters ($\lambda=\frac{ E\nu }{(1 + 2\nu )(1 - 2\nu )}$ and $\mu=\frac{ E\nu }{(1 + 2\nu )}$) we obtain:
\[ u_r(r) = Ar + \dfrac{B}{r} \\ A = \dfrac{(1+\nu)(1-2\nu)R_i^2p_i}{E(R_e^2-R_i^2)}, \quad B = \dfrac{(1+\nu)R_i^2R_e^2p_i}{E(R_e^2-R_i^2)}\]
Numerical solution
MEB parameters
materials
The constitutive behavior of the material considered is isotropic linear elastic. Since only one material is considered, the structs defined for the materials contain only one entry:
materials = struct() ;
materials.modelName = 'elastic-linear' ;
materials.modelParams = [ E nu ] ;
elements
In this plane model, three kinds of elements are used: triangle
for the solid, edges
to add pressure loads and nodes
to set additional boundary conditions for the numerical resolution. Since three kinds of elements are used, the struct has length 3:
elements = struct() ;
elements(1).elemType = 'node' ;
elements(2).elemType = 'edge' ;
elements(2).elemCrossSecParams = L ;
elements(3).elemType = 'triangle';
elements(3).elemTypeParams = 2 ;
elements(3).elemCrossSecParams = L ;
where elemCrossSecParams
field sets the thickness of the edge and elemTypeParams
sets the plane strain triangle element.
boundaryConds
Three BCs are considered, one corresponding to a load and two for displacements. The first two BCs constrain displacements in $x$ and $y$ global directions respectively:
boundaryConds = struct() ;
boundaryConds(1).imposDispDofs = [1] ;
boundaryConds(1).imposDispVals = [0] ;
boundaryConds(2).imposDispDofs = [3] ;
boundaryConds(2).imposDispVals = [0] ;
then the third BC corresponds to the pressure. It is introduced in local
coordinates so the first entry is along the edge (tangent) and the second towards the normal vector obtained by rotating the tangent vector 90 degrees in global axis $z$:
boundaryConds(3).loadsCoordSys = 'local' ;
boundaryConds(3).loadsTimeFact = @(t) t ;
boundaryConds(3).loadsBaseVals = [ 0 p ] ;
Mesh
The mesh can be read from the msh file. However, if any changes to the mesh are desired, the .geo file can be edited and the msh file can be re-generated using GMSH.
The element properties are set using labels into GMSH follwing the MEB nomenclature. First triangle
elements have linear elastic material so entry $1$ of the materialṣ struct is assigned. Then for both node
and edge
elements any material is set. Next displacement boundary conditions are assigned to the element, since the problem is modeled into $x-y$ plane, a constrain to avoid rotation along $z$ is necessary. This is done fixing $y$ and $x$ displacements (using boundaryConds(1)
and boundaryConds(2)
as labels) on points 2 3 4 5. Finally the internal pressure is applied on the edge
elements linked with curves from one to four (Circles 1-4 in Figure). In accordance with the orientation of the curve set in GMSH, the normal vector obtained in local coordinates is $e_r$ so the internal pressure is assigned using boundaryConds(3)
. Once the mesh is created is read using:
base_msh='';
if strcmp( getenv('TESTS_RUN'),'yes') && isfolder('examples'),
base_msh=['.' filesep 'examples' filesep 'ringPlaneStrain' filesep];
end
mesh = struct();
[ mesh.nodesCoords, mesh.conecCell ] = meshFileReader( [ base_msh 'ring.msh'] ) ;
initialConds
Any non-homogeneous initial conditions are considered, thereafter an empty struct is set:
initialConds = struct();
Analysis parameters
The Newton-Raphson method is employed to solve 2 load steps. The ratio between finalTime
and deltaT
sets the number of load steps used to evaluate boundaryConds(3).loadsTimeFact
function:
analysisSettings = struct() ;
analysisSettings.methodName = 'newtonRaphson' ;
analysisSettings.stopTolIts = 30 ;
analysisSettings.stopTolDeltau = 1.0e-12 ;
analysisSettings.stopTolForces = 1.0e-12 ;
analysisSettings.finalTime = 1 ;
analysisSettings.deltaT = .5 ;
Output parameters
otherParams = struct() ;
otherParams.problemName = 'linear_PlaneStrain' ;
otherParams.plots_format = 'vtk' ;
The ONSAS software is executed for the parameters defined above and the displacement solution of each load(time) step is saved in matUs
matrix:
[ modelCurrSol, modelProperties, BCsData ] = ONSAS_init( materials, elements, boundaryConds, initialConds, mesh, analysisSettings, otherParams ) ;
%
After that the structs are used to perform the numerical time analysis
[matUs, loadFactorsMat, modelSolutions ] = ONSAS_solve( modelCurrSol, modelProperties, BCsData ) ;
Verification
The numerical and analytic solutions are compared at the final load step for the internal and external surface (since all the elements on the same surface have the same analytic solution):
% internal surface analytic solution
A = ( p * (1+nu)*(1-2*nu)*Ri^2 ) / ( E*(Re^2-Ri^2) ) ;
B = ( p * (1+nu)*Ri^2*Re^2 ) / ( E*(Re^2-Ri^2) ) ;
analyticValRi = A*Ri + B/Ri ;
% internal surface numerical solution
dofXRi = 1 ;
numericalRi = matUs( dofXRi, end ) ;
% external surface analytic solution
analyticValRe = A*Re + B/Re ;
% external surface numerical solution
dofXRe = (8-1)*6+3 ;
numericalRe = matUs( dofXRe , end ) ;
The numerical and analytical solution for the internal and external surface are plotted:
%plot parameters
lw = 2.0 ; ms = 11 ; plotfontsize = 10 ;
figure, hold on, grid on
%internal surface
plot( matUs(dofXRi,:), loadFactorsMat(:,3) , 'ro' , 'linewidth', lw,'markersize',ms )
plot( linspace(0,analyticValRi,length(loadFactorsMat(:,3) ) ) , loadFactorsMat(:,3), 'k-', 'linewidth', lw,'markersize',ms )
%internal surface
plot( matUs(dofXRe,:), loadFactorsMat(:,3) , 'ro' , 'linewidth', lw,'markersize',ms )
plot( linspace(0,analyticValRe,length(loadFactorsMat(:,3)) ) , loadFactorsMat(:,3), 'k-', 'linewidth', lw,'markersize',ms )
labx = xlabel('Displacement [m]'); laby = ylabel('\lambda(t)') ;
legend('Numeric','Analytic','location','East')
set(gca, 'linewidth', 1.2, 'fontsize', plotfontsize )
set(labx, 'FontSize', plotfontsize); set(laby, 'FontSize', plotfontsize) ;
print('output/verifLinearRingPlaneStrain.png','-dpng')
Finally the deformed configuration is illustrated:
Elastoplastic analysis
Semi-analytic solution
The solution is extracted from Hill (The mathematical theroy of plasticity, 1950). The yielding pressure $p_0$ is defined as,
\[ Y = \dfrac{2\sigma_{Y,0}}{\sqrt{3}} \\ p_0 = \dfrac{Y}{2}\left(1+\dfrac{R_i^2}{R_e^2}\right).\]
The radial displacement of the outer surface of the ring is given by,
\[ u_r(R_e) = \text{if}~~ p \leq p_0 \\ \dfrac{2 p R_e}{E \left( \dfrac{R_e^2}{R_i^2-1}\right) }( 1-\nu^2 ) \\ \text{else} \\ \dfrac{2pR_e}{E\left(\dfrac{R_e^2}{R_i^2-1}\right)}(1-\nu^2)\]
where $c$ denotes the plastic front surface in the ring and is given by the implicit function,
\[ \dfrac{p}{Y} = ln\left(\dfrac{c}{R_i}\right) + \dfrac{1}{2}\left(1-\dfrac{c^ 2}{R_e^2}\right)\]
Numerical solution
% scalar parameters
E = 210 ; nu = 0.3 ; H = 0 ; sigmaY0 = 0.24 ; L = .75 ; p = 0.01 ;
MEB parameters
materials
The constitutive behavior of the material considered is isotropic hardening. Since only one material is considered, the structs defined for the materials contain only one entry:
materials.modelName = 'isotropicHardening' ;
materials.modelParams = [ E nu H sigmaY0 ] ;
elements
The elements struct is the same as the previous model.
boundaryConds
The BC struct is the same as in the elastic-linear case. However the loadsTimeFact function can be modified to consider unloading as follows.
boundaryConds(3).loadsTimeFact = @(t) t*(t<=19) + (t-(t-19)*2)*(t>19) ;
initialConds
Any non-homogeneous initial conditions are considered, thereafter the struc is the same as in the previous example.
Mesh
The mesh can be read from the msh file. The same mesh as in the elastic-linear case is considered for this problem.
Conec = myCell2Mat( mesh.conecCell ) ;
elems = size(Conec,1)
Analysis parameters
The Newton-Raphson method is employed to solve 19 load steps. The ratio between finalTime
and deltaT
sets the number of load steps used to evaluate boundaryConds(3).loadsTimeFact
function:
analysisSettings.methodName = 'newtonRaphson' ;
analysisSettings.stopTolIts = 30 ;
analysisSettings.stopTolDeltau = 1.0e-8 ;
analysisSettings.stopTolForces = 1.0e-6 ;
analysisSettings.finalTime = 19 ;
analysisSettings.deltaT = 1 ;
Output parameters
otherParams.problemName = 'EPP_PlaneStrain' ;
otherParams.plots_format = 'vtk' ;
The ONSAS software is executed for the parameters defined above and the displacement solution of each load(time) step is saved in matUs
matrix:
[ modelCurrSol, modelProperties, BCsData ] = ONSAS_init( materials, elements, boundaryConds, initialConds, mesh, analysisSettings, otherParams ) ;
%
After that the structs are used to perform the numerical time analysis
[matUs, loadFactorsMat, modelSolutions ] = ONSAS_solve( modelCurrSol, modelProperties, BCsData ) ;
Verification
The numerical and analytic solutions are compared for the external surface (since all the elements on the same surface have the same analytic solution):
%
global Y
%
Y = 2 * sigmaY0 / sqrt(3) ;
% p0 = Y/2 * (1-a^2/b^2)
p0 = Y / 2 * (1 - Ri^2 / Re^2) ; % Yielding pressure
%
pressure_vals = loadFactorsMat(:,3) * p ;
%
cvals = zeros(length(pressure_vals),1) ;
ubAna = zeros(length(pressure_vals),1) ;
%
% Plastic front value
for i = 1:length(cvals)
p = pressure_vals(i) ;
if i == 1
% val = fsolve(@(c)c_val(c,p,Y,a,b), a) ;
val = fsolve(@(c)c_val(c,p,Y,Ri,Re), Ri) ;
else
% val = fsolve(@(c)c_val(c,p,Y,a,b), cvals(i-1)) ;
val = fsolve(@(c)c_val(c,p,Y,Ri,Re), cvals(i-1)) ;
end
cvals(i) = val ;
end
%
% Analytic radial displacement at outer surface
for i = 1:length(cvals)
p = pressure_vals(i) ;
if p < p0
% ubAna(i) = 2*p*b / ( E*( b^2/a^2-1 ) ) * (1-nu^2) ;
ubAna(i) = 2*p*Re / ( E*( Re^2/Ri^2-1 ) ) * (1-nu^2) ;
else
c = cvals(i) ;
% ubAna(i) = Y*c^2/(E*b) * (1-nu^2) ;
ubAna(i) = Y*c^2/(E*Re) * (1-nu^2) ;
end
end
Plots
% plot parameters
lw = 2.0 ; ms = 11 ; plotFontSize = 10 ;
fig = figure;
hold on, grid on
% node to plot the solution
node = 5 ;
dofX = node * 6 - 5 ;
ubNum = matUs(dofX, :) ;
%
plot(ubNum, pressure_vals, 'b-o', 'linewidth', lw,'markersize', ms)
plot(ubAna, pressure_vals, 'g-x', 'linewidth', lw,'markersize', ms)
%
legend ({'FEM', 'Analytic',}, 'location', 'east');
labx = xlabel('u_b'); laby = ylabel('p') ;
tit = title('p-u_b');
set(labx, 'fontsize', plotFontSize*.8);
set(laby, 'fontsize', plotFontSize*.8);
set(tit, 'fontsize', plotFontSize);
The numerical solution is verified for both cases:
analyticCheckTolerance = 1e-2 ;
verifBoolean = ( ( numericalRi - analyticValRi ) < analyticCheckTolerance ) && ...
( ( numericalRe - analyticValRe ) < analyticCheckTolerance ) && ...
( ( ubNum(end) - ubAna(end) ) < analyticCheckTolerance ) ;
```