p_shell

Contents Functions

p_shell

Purpose

Element property function for shells and plates (flat shells)

Syntax

il = p_shell('default');
il = p_shell('database ProId name'); 
il = p_shell('dbval ProId name');
il = p_shell('dbval -unit TM ProId name');
il = p_shell('dbval -punit TM ProId name');
il = p_shell('SetDrill 0',il);

Description

This help starts by describing the main commands : p_shell Database and Dbval. Supported p_shell subtypes and their formats are then described.

Database,Dbval] ...

p_shell contains a number of defaults obtained with the database and dbval commands which respectively return a structure or an element property row. You can select a particular entry of the database with using a name matching the database entries.

You can also automatically compute the properties of standard shells with

kirchhoff e	Kirchhoff shell of thickness e.
mindlin e	Mindlin shell of thickness e.
laminate MatIdi Ti Thetai	Specification of a laminate property by giving the different ply MatId, thickness and angle. By default the z values are counted from -thick/2, you can specify another value with a z0.

You can append a string option of the form -f i to select the appropriate shell formulation. For example, you will obtain the element property row with ProId 100 associated with a .1 thick Kirchhoff shell (with formulation 5) or the corresponding Mindlin plate use

 il = p_shell('database 100 MindLin .1')
 il = p_shell('dbval 100 kirchhoff .1 -f5')
 il = p_shell('dbval 100 laminate z0=-2e-3 110 3e-3 30 110 3e-3 -30')
 il = fe_mat('convert SITM',il);
 il = p_shell(il,'dbval -unit TM 2 MindLin .1') % set in TM, provide data in SI
 il = p_shell(il,'dbval -punit TM 2 MindLin 100') % set in TM, provide data in TM

For laminates, you specify for each ply the MatId, thickness and angle.

Shell format description and subtypes

Element properties are described by the row of an element property matrix or a data structure with an .il field containing this row (see section 7.4). Element property functions such as p_shell support graphical editing of properties and a database of standard properties.

For a tutorial on material/element property handling see section 4.5.1. For a reference on formats used to describe element properties see section 7.4.

p_shell currently only supports two subtypes

1 : standard isotropic

  [ProID type   f d O   h   k   MID2 RatI12_T3 MID3 NSM Z1 Z2 MID4]

type		identifier obtained with fe_mat('p_shell','SI',1).
f		0 use default of element. For other formulations the specific help for each element (quad4, tria3, ...), each formulation specifies integration rule.
d	-1	no drilling stiffness. The element DOFs are the standard translations and rotations at all nodes (DOFs .01 to .06). The drill DOF (rotation .06 for a plate in the xy plane) has no stiffness and is thus eliminated by fe_mk if it corresponds to a global DOF direction. The default is d=1 (d is set to 1 for a declared value of zero).
	d	arbitrary drilling stiffness with value proportional to d is added. This stiffness is often needed in shell problems but may lead to numerical conditioning problems if the stiffness value is very different from other physical stiffness values. Start with a value of 1. Use il=p_shell('SetDrill d',il) to set to d the drilling stiffness of all p_shell subtype 1 rows of the property matrix il.
h		plate thickness.
k	k	shear correction factor (default 5/6, default used if k is zero). This correction is not used for formulations based on triangles since tria3 is a thin plate element.
RatI12_T3		Ratio of bending moment of inertia to nominal T3/I12 (default 1).
NSM		Non structural mass per unit area.
MID2		material property for bending. Defauts to element MatId if equal to 0.
MID3		material property for transverse shear.
z1,z2		(unused) offset for fiber computations.
MID4		material property for membrane/bending coupling.

Shell strain is defined by the membrane, curvature and transverse shear (display with p_shell('ConstShell')).

{

є_xx

є_yy

2 є_xy

κ_xx

κ_yy

2 κ_xy

γ_xz

γ_yz

}=[

N,x	0	0	0	0
0	N,y	0	0	0
N,y	N,x	0	0	0
0	0	0	0	N,x
0	0	0	−N,y	0
0	0	0	−N,x	N,y
0	0	N,x	0	−N
0	0	N,y	N	0

] {

} (9.23)

2 : composite

  [ProID type   Z0 NSM SB FT TREF GE LAM MatId1 T1 Theta1 SOUT1 ...]

ProID	Section property identification number.
type	Identifier obtained with fe_mat('p_shell','SI',2).
Z0	Distance from reference plate to bottom surface.
NSM	Non structural mass per unit area.
SB	Allowable shear stress of the bonding material.
FT	Failure theory.
TREF	Reference temperature.
Eta	Hysteretic loss factor.
LAM	Laminate type.
MatIdi	MatId for ply i, see m_elastic 1, m_elastic 5, ...
Ti	Thickness of ply i.
Thetai	Orientation of ply i.
SOUTi	Stress output request for ply i.

Note that this subtype is based on the format used by NASTRAN for PCOMP and the formulation used for each topology is discussed in each element (see quad4, tria3). You can use the DbvalLaminate commands to generate standard entries.

{

} = [

A	B	0
B	D	0
0	0	F

] {

} (9.24)

setTheta

When dealing with laminated plates, the classical approach uses a material orientation constant per element. OpenFEM also supports more advanced strategies with orientation defined at nodes but this is still poorly documented.

The material orientation is the reference for plies. Any angle defined in a laminate command is an additional rotation. In the example below, the element orientation is rotated 30 degrees, and the ply another 30. The fibers are thus oriented 60 degrees in the xy plane. Stresses are however given in the material orientation thus with a 30 degree rotation. Per ply output is not currently implemented.

The element-wise material angle is stored for each element. In column 7 for tria3, 8 for quad4, ... The setTheta command is a utility to ease the setting of these angles. By default, the orientation is done at element center. To use the mean orientation at nodes use command option -strategy 2.

model=ofdemos('composite');
model.il = p_shell('dbval 110 laminate 100 1 30'); % single ply

% Define material angle based on direction at element
MAP=feutil('getnormalElt MAP -dir1',model);
bas=basis('rotate',[],'rz=30;',1);
MAP.normal=MAP.normal*reshape(bas(7:15),3,3)';
model=p_shell('setTheta',model,MAP);

% Obtain a MAP of material orientations
MAP=feutil('getnormalElt MAP -dir1',model);
feplot(model);fecom('showmap',MAP)

% Set elementwise material angles using directions given at nodes. 
% Here a global direction
MAP=struct('normal',ones(size(model.Node,1),1)*bas(7:9), ...
    'ID',model.Node(:,1),'opt',2);
model=p_shell('setTheta',model,MAP);

% Using an analytic expression to define components of 
% material orientation vector at nodes
data=struct('sel','groupall','dir',{{'x-0','y+.01',0}},'DOF',[.01;.02;.03]);
model=p_shell('setTheta',model,data);
MAP=feutil('getnormalElt MAP -dir1',model);
feplot(model);fecom('showmap',MAP)

model=p_shell('setTheta',model,0) is used to reset the material orientation to zero.

Technically, shells use the of_mk('BuildNDN') rule 23 which generates a basis at each integration point. The first vector v1x,v1y,v1z is built in the direction of r lines and v2x,v2y,v2z is tangent to the surface and orthogonal to v1. When a InfoAtNode map provides v1x,v1y,v1z, this vector is projected (NEED TO VERIFY) onto the surface and v2 taken to be orthogonal.