N2PEdgeDelamination User Manual#

Overview#

The N2PEdgeDelamination Module is designed to calculate the stresses in the free edge regions of symmetric laminated plates subjected to a uniform prescribed axial strain and temperature change. It also predicts the onset of delamination. Calculations are based on the 95028 ESDU Data Item (See [Ref-1]).

Classical laminated plate theory is widely used to calculate the elastic stresses in flat laminated plates build up from thin orthotropic layers. However, at free edges region, classical plate theory does not account for the more general three dimensional nature of the stress field which allows the free edge surface to have zero values of direct stress and shear stresses. Edge stresses arise when the forces in a layer, predominantly in the fibre direction within that layer, have to be transferred to another layer as the fibres approach the free edge (for example in laminates containing angle plies). This implies stresses in the laminate in the through-the-thickness direction which can be large enough to cause local delamination. The magnitude of the edge stresses is increased when, for given laminate thickness, a smaller number of layers is used and each contains a larger number of plies.

The module performs calculations for the edges delimiting a set elements defined as a list of N2PElements (see [Mandatory Inputs]). However, it must be noted that module does not check if the edges of the given set are free edges. In case there is an element with more than one edge (for example, a corner), the module will return the results corresponding to the more critical one.

Note

When solving RF equations, any solution greater than 50 is automatically capped at 50 for visualization purposes.

Example Usage#

1. Importing the module#

import NaxToPy as n2p
from NaxToPy.Modules.static.delamination.N2PEdgeDelamination import N2PEdgeDelamination

2. Load the model#

model = n2p.load_model(r"file path")

3. Define Elements and Load Cases#

element_list = [3026292,3026293,3026294] # List of element
n2pelem = model.get_elements(element_list)

load_case_id = [68195, 68196] # List of ID of Load Cases
n2plc = model.get_load_case(load_case_id)

4. Create and Configure N2PEdgeDelamination#

delam = N2PEdgeDelamination()
delam.Model = model  # Assign the model
delam.ElementList = n2pelem  # Assign elements
delam.LoadCases = n2plc  # Assign load cases

5. Material properties and allowables assignment (Mandatory inputs)#

mat_ID = 4  #ID of the material 
part_mat = '0' # Part ID of the material  
delam.Materials[(mat_ID,part_mat)].YoungZ = 12000  #Assign the Z elastic modulus
delam.Materials[(mat_ID,part_mat)].PoissonXZ = 0.3  # Assign the XZ Poisson ratio
delam.Materials[(mat_ID,part_mat)].PoissonYZ = 0.3  # Assign the YZ Poisson ratio 

delam.Materials[(mat_ID,part_mat)].Allowables.ShearXZ = 105  #Assign the shear XZ allowable
delam.Materials[(mat_ID,part_mat)].Allowables.ShearYZ = 105  #Assign the shear YZ allowable
delam.Materials[(mat_ID,part_mat)].Allowables.ZCompressive = 108  #Assign the compressive Z allowable
delam.Materials[(mat_ID,part_mat)].Allowables.ZTensile = 53.9  #Assign the tensile Z allowable

6. Temperature, integration distance and failure criteria(Optional inputs)#

delam.FailureCriteria = 'CI'

for LC in delam.LoadCases:
    LC.RefTemp = -30

int_distance = {}
for element in n2pelem:
    int_distance[element] = 5
delam.IntegrationDistance = int_distance

7. Run Calculations#

delam.HDF5.FilePath = r"file output"
delam.calculate()

Setup#

Prerequisites#

Python Environment: Ensure you have a compatible Python environment with the following libraries:
- NaxToPy modules
NaxToPy Framework: The module integrates with the NaxToPy library, which must be installed and configured in your environment.

Module Inputs#

The N2PEdgeDelamination class provides methods and attributes for defining the model, associating materials, and calculating interlaminar stresses and RF.

Mandatory Inputs#

Property	Type	Description
.Model	`N2PModelContent`	FEM model to be analyzed
.ListElements	`list[N2PElement]`	List of elements to analyze
.LoadCases	`list[N2PLoadCase]`	List of load cases to consider
.Material[(Material ID,Increment ID)].YoungZ	`float`	Elastic modulus in Z direction
.Material[(Material ID,Increment ID)].PoissonXZ	`float`	Poisson ratio in XZ direction
.Material[(Material ID,Increment ID)].PoissonYZ	`float`	Poisson ratio in Z direction
.Material[(Material ID,Increment ID)].Allowables.Zc	`float`	Compressive allowable in Z direction
.Material[(Material ID,Increment ID)].Allowables.Zt	`float`	Tensile allowable in Z direction
.Material[(Material ID,Increment ID)].Allowables.ShearXZ	`float`	Shear allowable in XZ direction
.Material[(Material ID,Increment ID)].Allowables.ShearYZ	`float`	Shear allowable in XZ direction

FEM Requirements#

For this module to function properly, the FEM model must include node displacements results in the output requests. Additionally, the model should accurately represent all elements and load cases where edge delamination is expected to be calculated. This module is available for Nastran and Abaqus.

Supported Type Files#

Solver	File Extension
Nastran	.op2
Nastran	.xdb
Nastran	.h5
Optistruct	.op2
Optistruct	.h3d
Abaqus	.odb

Conditional Inputs#

No conditional inputs required in this module

Optional Inputs#

Property/Method	Type	Description
.HDF5.FilePath	`str`	Path to specify the location of the created HDF5 file. If no path is provided, the file is created in memory
.HDF5.export_to_HDF5	`method`	Convert an HDF5 file that has been created in memory to an HDF5 file in disk
.N2PLoadCase[(Load Case ID)].Temperature	`float`	Temperature change with respect from an initial stress condition (Default: 0ºC)
.IntegrationDistance	`dict`	Dictionary mapping each N2PElement with its corresponding integration distance defined in [Ref-1] (Default: thickness of the laminate)
.FailureCriteria	`str`	Identification letters for the chosen failure criterion: MI (mechanical), TI (thermal) and CI (both mechanical and thermal) (Default: CI)

To export the data into a HDF5 file set the path in the property HDF5:

delam.HDF5.FilePath = "output.h5"

Warning

For FEM models with many load cases and elements to analyze, the option to export results to HDF5 must be enabled; otherwise, an exception will be raised due to memory overflow.

Output Description#

RF_c: Reserve factor in compression calculated with the expressions provided in RF calculations
RF_t: Reserve factor in tension calculated with the expressions provided in RF calculations
Mechanical interlaminar stresses: Mechanical through-the-thickness stresses (see Interlaminar mechanical and thermal loads) -\(\sigma_{zk}\) -\(\tau_{xzk}\) -\(\tau_{yzk}\)
Thermal interlaminar stresses: Thermal through-the-thickness stresses (see Interlaminar mechanical and thermal loads) -\(\sigma_{zkT}\) -\(\tau_{xzkT}\) -\(\tau_{yzkT}\)

If a path is given the results will be export in a HDF5 with NaxTo format.

Material Associaton#

Materials are automatically associated with elements.

Warning

Ensure Z elastic modulus, shear XZ and YZ modulus, and allowables are defined for all materials used.

Common Errors#

C800: Assigned model input must be a N2PModelContent instance
C801: Assigned elements must be a list of N2PElement instances
C802: Assigned elements must be a list of N2PLoadCase instances
C850, C851, C852: IntegrationDistance must be a dict mapping each N2PElement with its corresponding integration distance
C853: Failure criteria must be ‘MI’, ‘TI’ or ‘CI’
C854: Material is not an Orthotropic instance
C855: YoungZ, PoissonXZ and PoissonYZ must be defined
C856: Maximum number of iterations reached without convergence
C857: Allowables for material must be defined
C858: Load case list is a mandatory input
C859: Element list is a mandatory input
C860: Solver selected is not implemented in this module
C861: Property of the element selected must be a CompositeShell

Common Warnings#

W850: No reference temperature is set for a Load Case. By default is set to 0
W851: No integration distance is set. By default is set to the thickness of the element

Methodology#

Background Theory#

The methodology is based on the structural calculations outlined the ESDU Data Item 95028. See [Ref-1]

Previous calculations#

According to [Ref-1], to calculate the delamination stresses it is necessary to know the strain of each edge. In this case, the process to follow to obtain the strain of an edge will be the following:

Compute initial position from the nodes on the edge of the
Obtain displacements from results file.
Compute final position of the nodes by adding the displacement to the initial position

\[ \overrightarrow{P_f} = \overrightarrow{P_0} + \Delta\overrightarrow{P} \]

\[ \overrightarrow{Q_f} = \overrightarrow{Q_0} + \Delta\overrightarrow{Q} \]

Compute strain as the ratio between the initial and final length

\[ \epsilon_x = \frac{\left\| \overrightarrow{P_fQ_f} \right\| - \left\| \overrightarrow{P_0Q_0} \right\|}{\left\| \overrightarrow{P_0Q_0} \right\|} \]

Once strain is calculated on each edge, it is calculated the relative angle between global axis and the edge, which will be known as \(\beta\)

Equivalent mechanical loads#

Equivalent mechanical loads on each edge are calculated by the following process:

Obtain matrix [A], relating normal and shear in-plane forces with the in-plane deformations
Retrieve \(\epsilon_x\) from previous section
Obtain \(N_x\) from constitutive laminate equations:

\[\begin{split} \begin{bmatrix} N_x \\ N_y \\ N_{xy} \end{bmatrix} = \begin{bmatrix} A_{11} & A_{12} & A_{16} \\ A_{21} & A_{22} & A_{26} \\ A_{16} & A_{26} & A_{66} \end{bmatrix} \cdot \begin{bmatrix} \epsilon_x \\ \epsilon_y \\ \epsilon_{xy} \end{bmatrix} \end{split}\]

Considering \(N_y = N_{xy} = 0\), solve the following system for \(\epsilon_y\) and \(\gamma_xy\)

\[\begin{split} \begin{bmatrix} A_{22} & A_{26} \\ A_{26} & A_{66} \end{bmatrix} \cdot \begin{bmatrix} \epsilon_y \\ \gamma_{xy} \end{bmatrix} = \begin{bmatrix} -A_{12} \cdot \epsilon_x \\ -A_{16} \cdot \epsilon_x \end{bmatrix} \end{split}\]

Substitute in the remaining equation to obtain \(N_x\)

\[ N_x = A_{11}\cdot\epsilon_x + A_{12}\cdot\epsilon_y+A_{16}\cdot\gamma_{xy} \]

Equivalent thermal loads#

Equivalent thermal loads on each edge are calculated by the following process:

Calculate the orientation of each edge with respect to the global axes, to do this, add the orientation of the edge with respect to the global axes (\(\beta\)) to the orientation of the laminate with respect to the edge (\(\theta\))

\[ \psi = \theta + \beta \]

Calculate thermal strains

\[ e_{\theta x} = e_{\theta \alpha} \cdot cos^{2}(\psi) + e_{\theta \beta} \cdot sin^{2}(\psi) \]

\[ e_{\theta y} = e_{\theta \alpha} \cdot sin^{2}(\psi) + e_{\theta \beta} \cdot cos^{2}(\psi) \]

\[ e_{\theta xy} = (e_{\theta \alpha} -e_{\theta \beta}) \cdot 2\cdot cos(\psi) \cdot sin(\psi) \]

where:

\[ e_{\theta \alpha} = \mu_{\alpha} \cdot \Delta T ; e_{\theta \beta} = \mu_{\beta} \cdot \Delta T \]

Calculate thermal loads

\[\begin{split} \begin{bmatrix} N_x \\ N_y \\ N_{xy} \end{bmatrix} = \sum_{i=1}^{n} [\bar{Q}] \cdot \begin{bmatrix} e_{\theta x} \\ e_{\theta y} \\ e_{\theta xy} \end{bmatrix} \end{split}\]

Calculate deformation vectors

\[\begin{split} \begin{bmatrix} \epsilon_x \\ \epsilon_y \\ \gamma_{xy} \end{bmatrix} = [A]^{-1} \cdot \begin{bmatrix} N_{x} \\ N_{y} \\ N_{xy} \end{bmatrix} \end{split}\]

Interlaminar mechanical and thermal loads#

To calculate the mechanical loads between layers, we must calculate the parameters \(\lambda\) and \(\phi\) that minimize the complementary energy equation (\(\Pi_c\)). Finally, to calculate the parameters, a non-linear system of equations must be solved (the C coefficients are explained in the Ref-1):

\[ (\lambda^{4}\phi^{4}) C_1 + (\lambda^{3}\phi^{4}) C_2 + \lambda^{2} C_3 + (\lambda^{2}\phi^{2}) C_4 + \lambda C_5 + C_6 = 0 \]

\[ (\lambda^{3}\phi^{4}) C_7 + (\lambda^{2}\phi^{2}) C_8 + (\lambda \phi^{2}) C_9 + \lambda^{2} C_{10} + \lambda C_{11} + C_{12} = 0 \]

Interlaminar stresses are calculated by the following formulas (stresses marked with an uppercase bar are those computed from the classical laminate plate):

\[ \sigma_{zkA} = \frac{1}{d} \int_{0}^{d} \bar{\sigma_{yk}} \cdot \left[ \left( \frac{\lambda}{1-\lambda} \right) \phi^{2} e^{-\phi y} + \left( \frac{1}{\lambda - 1} \right) \cdot \lambda^{2}\phi^{2}e^{-\lambda \phi y} \right] \cdot \left( \frac{z^2}{2} + zB_{4k} + B_{5k} \right)\cdot dy = \frac{1}{d} \bar{\sigma_{yk}} \cdot \left( \frac{z^2}{2} + zB_{4k} + B_{5k} \right)\left[ \left( \frac{\lambda}{1-\lambda} \right) \phi (1 - e^{-\phi d}) + \left( \frac{1}{\lambda -1} \right) \cdot \lambda \phi (1 - e^{-\lambda \phi d}) \right] \]

\[ \tau_{xzkA} = \frac{1}{d} \int_{0}^{d} -\bar{\tau}_{xyk} \cdot \phi e^{-\phi y} \cdot(z + B_{2k}) \cdot dy = \frac{1}{d} \cdot \bar{\tau}_{xyk}\cdot(z + B_{2k}) \cdot (-1 + e^{-\phi d}) \]

\[ \tau_{yzkA} = \frac{1}{d} \int_{0}^{d} \bar{\sigma_{yk}} \cdot \left[ \left( \frac{\lambda}{1-\lambda} \right) \phi e^{-\phi y} + \left( \frac{1}{\lambda - 1} \right) \cdot \lambda\phi e^{-\lambda \phi y} \right] \cdot \left( z + B_{4k}\right)\cdot dy = \frac{1}{d} \bar{\sigma_{yk}} \cdot \left( z + B_{4k}\right)\cdot \left[ \left( \frac{\lambda}{1-\lambda} \right) (1 - e^{-\phi d}) + \left( \frac{1}{\lambda - 1} \right) \cdot (1 - e^{-\lambda \phi d}) \right] \]

When mechanical interlaminar loads are computed, thermal loads are calculated with the same procedure but using the thermal stresses.

RF calculations#

Depending on the type of calculation chosen (mechanical, thermal or combined), one of the following formulas will be used:

MI: Applied strain variable and temperature change constant

\[ \left( \frac{(RF_M)\tau_{xzkA} + \tau_{xzkTA}}{\tau_{xzkF}} \right)^{2} + \left(\frac{(RF_M)\tau_{yzkA} + \tau_{yzkTA}}{\tau_{yzkF}} \right)^{2} + \left( \frac{(RF_M)\sigma_{zkA} + \sigma_{zkTA}}{\sigma_{zkFt} (or \hspace{0.1cm} \sigma_{zkFc})} \right)^{2} = 1 \]

TI: Applied strain constant and variable temperature change

\[ \left( \frac{\tau_{xzkA} + (RF_T)\tau_{xzkTA}}{\tau_{xzkF}} \right)^{2} + \left( \frac{\tau_{yzkA} + (RF_T)\tau_{yzkTA}}{\tau_{yzkF}} \right)^{2} + \left( \frac{\sigma_{zkA} + (RF_T)\sigma_{zkTA}}{\sigma_{zkFt} (or \hspace{0.1cm} \sigma_{zkFc})} \right)^{2} = 1 \]

CI: Applied strain and temperature change are varied by the same factor

\[ \left( \frac{RF_C(\tau_{xzkA} + \tau_{xzkTA})}{\tau_{xzkF}} \right)^{2} + \left( \frac{RF_C(\tau_{yzkA} + \tau_{yzkTA})}{\tau_{yzkF}} \right)^{2} + \left( \frac{RF_C(\sigma_{zkA} + \sigma_{zkTA})}{\sigma_{zkFt} (or \hspace{0.1cm} \sigma_{zkFc})} \right)^{2} = 1 \]

RFs are obtained from solving the corresponding second order equation, hence leading to 2 RFs associated to each layer (one corresponding for the tension and one for the compression)

Implementation#

The main function of this module is the extraction of results using the get_result_by_LCs_Incr method from N2PModelContent. Results are processed with ndarray. The main calculation process is divided into batches and load cases. If a file path is provided, memory will be released as it iterates through the load cases.

Advanced Features#

No advanced features of the module.

References#

[Ref-1]: ESDU 95028.

This guide serves as a comprehensive reference for using the N2PEdgeDelamination Module in structural analysis tasks. For further information, refer to the source code comments and examples provided.