# Static Von-Mises Truss example

In this example the Static Von Mises Truss problem and its resolution using ONSAS are described. The aim of this example is to validate the implementations of the Newton-Raphson and Newton-Raphson-Arc-Length methods by comparing the results provided with the analytic solutions.

The structural model is formed by two truss elements with length $L$ as it is shown in the figure, with nodes $1$ and $3$ fixed, and node $2$ submitted to a nodal load $P$ and constrained to move only in the $x-z$ plane.

## Analytic solutions

The solutions for the nonlinear cases are developed in section 2.3 of (Bazzano and Pérez Zerpa, 2017). The expressions obtained for different strain measures are:

- Rotated-Engineering: $P = \dfrac{EA_o(z_2+w)\left(\sqrt{(w+z_2)^2+x_2^2}-l_o\right)}{l_o\sqrt{(w+z_2)^2+x_2^2}}$
- SVK: $P = \dfrac{EA_o (z_2+w)\left( 2 z_2 w + w^2 \right) }{ 2 l_o^3 }$

where $x_2$ and $z_2$ are the coordinates of node 2 in the reference configuration and $w$ is the vertical displacement in the $z$ direction.

## Numerical solutions

Before defining the structs, the workspace is cleared, the ONSAS directory is added to the path

```
close all, if ~strcmp( getenv('TESTS_RUN'), 'yes'), clear all, end
addpath( genpath( [ pwd '/../../src'] ) );
```

some scalar parameters are defined and computed

```
% scalar parameters
E = 210e9 ; nu = 0 ;
A = 2.5e-3 ; ang1 = 65 ; L = 2 ;
% x and z coordinates of node 2
x2 = cos( ang1*pi/180 ) * L ;
z2 = sin( ang1*pi/180 ) * L ;
```

### MEB parameters

The modelling of the structure begins with the definition of the Material-Element-BoundaryConditions (MEB) parameters.

#### materials

The `materials`

struct is initialized as empty.

`materials = {} ;`

Since for each model both bars are formed by the same material only one `materials`

struct is defined. The constitutive behavior considered in the first analysis case is the Rotated Engineering strain, then the field `hyperElasModel`

is set to:

`materials.hyperElasModel = '1DrotEngStrain' ;`

and in the field `hyperElasParams`

a vector with the parameters of the Engineering Strain model is set

`materials.hyperElasParams = [ E nu ] ;`

which in the case of this material model are the Young modulus and the Poisson ratio. The field `density`

is not set, then the default $\rho = 0$ value is considered by ONSAS.

#### elements

The `elements`

struct is initialized as empty

`elements = {} ;`

Two different types of elements are required to create the model: `node`

and `truss`

, thus, the `elements`

struct will have two entries. The type of the first entry is

`elements(1).elemType = 'node' ;`

and the second entry is

`elements(2).elemType = 'truss';`

for the geometries, the node has no geometry to assign, and the truss elements will be set as the native `circle`

cross-section, then the elemCrossSecParams field is:

`elements(2).elemCrossSecParams = { 'circle' , sqrt(A*4/pi) } ;`

#### boundaryConds

The elements are submitted to two different BoundaryConditions, then the struct `boundaryConds`

will have length two. The nodes $1$ and $3$ are fixed, without loads applied (this is the first BC), and node $2$ has a constraint in displacement and an applied load (second BC). For the displacements, the first BC corresponds to a xyz-fixed displacement,

```
boundaryConds = {} ;
boundaryConds(1).imposDispDofs = [ 1 3 5 ] ;
boundaryConds(1).imposDispVals = [ 0 0 0 ] ;
```

and the second BC corresponds to a zero displacement only in the $y$ direction.

```
boundaryConds(2).imposDispDofs = 3 ;
boundaryConds(2).imposDispVals = 0 ;
```

Regarding the loads, the second BC is set so that the final load factor is $3 \cdot 10^8$ at 1 second. The default zero density is used, then no inertial effects are considered.

```
boundaryConds(2).loadsCoordSys = 'global' ;
boundaryConds(2).loadsTimeFact = @(t) 3.0e8*t ;
boundaryConds(2).loadsBaseVals = [ 0 0 0 0 -1 0 ] ;
```

### mesh parameters

The coordinates of the nodes of the mesh are given by the matrix:

```
mesh = {} ;
mesh.nodesCoords = [ 0 0 0 ; ...
x2 0 z2 ; ...
2*x2 0 0 ] ;
```

where the columns 1,2 and 3 correspond to $x$, $y$ and $z$ coordinates, respectively, and the row $i$-th corresponds to the coordinates of node $i$.

The connectivity is introduced using the *conecCell* cell. Each entry of the cell (indexed using {}) contains a vector with the four indexes of the MEBI parameters, followed by the indexes of the nodes of the element (node connectivity). For didactical purposes each element entry is commented. First the cell is initialized:

`mesh.conecCell = cell(5,1) ;`

Then the entry of node $1$ is introduced:

`mesh.conecCell{ 1, 1 } = [ 0 1 1 1 ] ;`

the first MEB parameter (Material) is set as *zero* (since nodes dont have material). The second parameter corresponds to the Element, and a *1* is set since `node`

is the first entry of the `elements.elemType`

cell. For the BC index, we consider that node $1$ is fixed, then the first index of the `boundaryConds`

struct is used. Finally, at the end of the vector the number of the node is included (1). A similar approach is used for node $3$,

`mesh.conecCell{ 2, 1 } = [ 0 1 1 3 ] ;`

and for node $2$ only the boundary condition is changed.

`mesh.conecCell{ 3, 1 } = [ 0 1 2 2 ] ;`

Regarding the truss elements, the first material is considered, the second type of element, and no boundary conditions are applied.

```
mesh.conecCell{ 4, 1 } = [ 1 2 0 1 2 ] ;
mesh.conecCell{ 5, 1 } = [ 1 2 0 2 3 ] ;
```

#### initial Conditions

homogeneous initial conditions are considered, then an empty cell is set:

`initialConds = {} ;`

### analysisSettings

The method used in the analysis is the Newton-Raphson, then the field `methodName`

must be introduced as:

```
analysisSettings = {};
analysisSettings.methodName = 'newtonRaphson' ;
```

and the following parameters correspond to the iterative numerical analysis settings

```
analysisSettings.deltaT = 0.1 ;
analysisSettings.finalTime = 1 ;
analysisSettings.stopTolDeltau = 1e-8 ;
analysisSettings.stopTolForces = 1e-8 ;
analysisSettings.stopTolIts = 15 ;
```

### otherParams

```
otherParams = {};
otherParams.problemName = 'staticVonMisesTruss_NR_RotEng';
otherParams.plots_format = 'vtk' ;
otherParams.plots_deltaTs_separation = 2 ;
```

### Analysis case 1: Newton-Raphson with Rotated Eng Strain

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 ) ;
controlDispsNREngRot = -matUs(11,:) ;
loadFactorsNREngRot = loadFactorsMat(:,2) ;
```

### Analysis case 2: Newton-Raphson with linear elastic behavior

In this case a linear elastic behavior is assumed. Then the hyperelasmodel es overwritten

```
materials.hyperElasModel = 'linearElastic' ;
otherParams.problemName = 'staticVonMisesTruss_linearElastic';
analysisSettings.finalTime = 1.5 ;
```

and the analysis is run again

`[matUs, loadFactorsMat] = ONSAS( materials, elements, boundaryConds, initialConds, mesh, analysisSettings, otherParams ) ;`

the displacements are extracted

```
controlDispsNRLinearElastic = -matUs(11,:) ;
loadFactorsNRLinearElastic = loadFactorsMat(:,2) ;
```

and the analytic values of the load factor are computed, as well as its difference with the numerical solution

```
analyticLoadFactorsNREngRot = @(w) -2 * E*A* ...
( ( (z2+(-w)).^2 + x2^2 - L^2 ) ./ (L * ( L + sqrt((z2+(-w)).^2 + x2^2) )) ) ...
.* (z2+(-w)) ./ ( sqrt((z2+(-w)).^2 + x2^2) ) ;
difLoadEngRot = analyticLoadFactorsNREngRot( controlDispsNREngRot)' - loadFactorsNREngRot ;
```

### Analysis case 3: NR with Green Strain

In order to perform a SVK case analysis, the material is changed and the problemName is also updated

```
otherParams.problemName = 'staticVonMisesTruss_NR_Green';
materials.hyperElasModel = 'SVK' ;
analysisSettings.finalTime = 1.0 ;
lambda = E*nu/((1+nu)*(1-2*nu)) ; mu = E/(2*(1+nu)) ;
materials.hyperElasParams = [ lambda mu ] ;
```

the load history is also changed

```
boundaryConds(2).loadsTimeFact = @(t) 1.5e8*t ;
%boundaryConds(2).userLoadsFilename = 'myVMLoadFunc' ;
```

and the analysis is run

`[matUs, loadFactorsMat] = ONSAS( materials, elements, boundaryConds, initialConds, mesh, analysisSettings, otherParams ) ;`

and the displacements are extracted

```
controlDispsNRGreen = -matUs(11,:) ;
loadFactorsNRGreen = loadFactorsMat(:,2) ;
```

the analytic solution is computed

```
analyticLoadFactorsGreen = @(w) - 2 * E*A * ( ( z2 + (-w) ) .* ( 2*z2*(-w) + w.^2 ) ) ./ ( 2.0 * L^3 ) ;
difLoadGreen = analyticLoadFactorsGreen( controlDispsNRGreen )' - loadFactorsNRGreen ;
```

### Analysis case 4: NR-AL with Green Strain

```
elements(2).elemCrossSecParams{1,1} = 'rectangle' ;
elements(2).elemCrossSecParams{2,1} = [ sqrt(A) sqrt(A)] ;
```

The same loading conidition as before is used, but given by a user load function. The argument set in this case is:

In this case, the numerical method is changed for newtonRaphson arc length.

```
otherParams.problemName = 'staticVonMisesTruss_NRAL_Green' ;
analysisSettings.methodName = 'arcLength' ;
analysisSettings.finalTime = 1 ;
analysisSettings.incremArcLen = 0.15 ;
analysisSettings.iniDeltaLamb = boundaryConds(2).loadsTimeFact(.2)/100 ;
analysisSettings.posVariableLoadBC = 2 ;
```

```
[matUs, loadFactorsMat] = ONSAS( materials, elements, boundaryConds, initialConds, mesh, analysisSettings, otherParams ) ;
controlDispsNRALGreen = -matUs(11,:) ;
loadFactorsNRALGreen = loadFactorsMat(:,2) ;
analyticLoadFactorsNRALGreen = analyticLoadFactorsGreen(controlDispsNRALGreen);
difLoadGreenNRAL = analyticLoadFactorsNRALGreen' - loadFactorsNRALGreen ;
```

### Analysis case 4: NR-AL Jirasek with Green Strain

```
otherParams.problemName = 'staticVonMisesTruss_NRAL_Jirasek_Green' ;
analysisSettings.methodName = 'arcLength' ;
analysisSettings.finalTime = 1 ;
analysisSettings.incremArcLen = 0.15 ;
analysisSettings.iniDeltaLamb = boundaryConds(2).loadsTimeFact(.2)/100 ;
analysisSettings.posVariableLoadBC = 2 ;
```

Jirasek variant - Dominant dof Extracted from Jirásek & Bazant book Inelastic Analysis of Structures, 2002 Chapter 22, Numerical Methods in Plasticity Sets arcLengthFlag = 2 to secifiy Jirasek constraint method.

```
global arcLengthFlag
arcLengthFlag = 2 ;
```

The dominant dof selected for this problem correpsonds with the displacement uz of node 2.

```
global dominantDofs
dominantDofs = 11 ;
```

The scaling projection for the Jirasek method and for this selected dof is set as follows.

```
global scalingProjection
scalingProjection = -1 ;
```

```
[matUs, loadFactorsMat] = ONSAS( materials, elements, boundaryConds, initialConds, mesh, analysisSettings, otherParams ) ;
controlDispsNRAL_Jirasek_Green = -matUs(11,:) ;
loadFactorsNRAL_Jirasek_Green = loadFactorsMat(:,2) ;
analyticLoadFactorsNRAL_Jirasek_Green = analyticLoadFactorsGreen(controlDispsNRAL_Jirasek_Green);
difLoadGreenNRAL_Jirasek = analyticLoadFactorsNRAL_Jirasek_Green' - loadFactorsNRAL_Jirasek_Green ;
```

## Verification

the numerical resolution is validated for both strain measures.

```
verifBoolean = ( ( norm( difLoadEngRot ) / norm( loadFactorsNREngRot ) ) < 1e-4 ) ...
&& ( ( norm( difLoadGreen ) / norm( loadFactorsNRGreen ) ) < 1e-4 ) ...
&& ( ( norm( difLoadGreenNRAL ) / norm( loadFactorsNRALGreen ) ) < 1e-4 ) ...
&& ( ( norm( difLoadGreenNRAL_Jirasek ) / norm( loadFactorsNRAL_Jirasek_Green ) ) < 1e-4 ) ;
```

### Plots

and solutions are plotted.

```
lw = 2.0 ; ms = 11 ; plotfontsize = 18 ;
figure
plot( controlDispsNREngRot, analyticLoadFactorsNREngRot( controlDispsNREngRot) ,'b-x' , 'linewidth', lw,'markersize',ms )
hold on, grid on
plot( controlDispsNREngRot, loadFactorsNREngRot, 'k-o' , 'linewidth', lw,'markersize',ms )
plot( controlDispsNRALGreen, analyticLoadFactorsGreen( controlDispsNRALGreen ), 'g-x' , 'linewidth', lw,'markersize',ms )
plot( controlDispsNRGreen, loadFactorsNRGreen, 'r-s' , 'linewidth', lw,'markersize',ms )
plot( controlDispsNRALGreen, loadFactorsNRALGreen, 'c-^' , 'linewidth', lw,'markersize',ms )
plot( controlDispsNRAL_Jirasek_Green, loadFactorsNRAL_Jirasek_Green, 'y-*' , 'linewidth', lw,'markersize',ms )
plot( controlDispsNRLinearElastic, loadFactorsNRLinearElastic, 'm-+' , 'linewidth', lw,'markersize',ms )
labx = xlabel('Displacement w(t)'); laby = ylabel('\lambda(t)') ;
legend( 'analytic-RotEng', 'NR-RotEng','analytic-Green', 'NR-Green','NRAL-Green', 'NRAL-Jirasek-Green','LinearElastic', 'location','northoutside')
set(gca, 'linewidth', 1.0, 'fontsize', plotfontsize )
set(labx, 'FontSize', plotfontsize); set(laby, 'FontSize', plotfontsize) ;
print('output/vonMisesTrussCheck.png','-dpng')
```