Semi-analytical Models for Cones and Cylinder (compmech.conecyl)#

The ConeCyl class embodies all the methods and functions necessary to perform:

  • linear static analysis

  • linear buckling analysis

  • non-linear static analysis

in conical and cylindrical structures.

The implementation focused on laminate composite shells with a constant laminate constitutive relation. This means that the semi-analtical models were derived doing an integration using a constant \([F]\) for the whole domain. Recalling that \([F]\) correlates the strains and the distributed stresses by the relation:

\[\{N\} = [F] \{\varepsilon\}\]

The ConeCyl object#

class compmech.conecyl.ConeCyl[source]#
Attributes:
E11
F
F_reuse
Fc
H
K
L
LA
MLA
Nxxtop
P
P_inc
T
T_inc
Ts
Xs
alphadeg
alpharad
analysis
bc
betadeg
betarad
c0
cosa
cs
eigvals
eigvecs
excluded_dofs
excluded_dofs_ck
force_orthotropic_laminate
forces
forces_inc
funcnum
h
increments
inf
is_cylinder
k0
k0uk
k0uu
kG
kG0
kG0_Fc
kG0_P
kG0_T
kGuu
kL
kLuu
kTuk
kTuu
kphitBot
kphitTop
kphixBot
kphixTop
kuBot
kuTop
kvBot
kvTop
kwBot
kwTop
lam
laminaprop
laminaprops
m0
m1
m2
model
n0
n2
name
ni_method
ni_num_cores
nt
nu
num0
num_eigvalues
num_eigvalues_print
nx
out_num_cores
outputs
pdC
pdLA
pdT
phit
phix
plyt
plyts
r1
r2
s
sina
size
stack
tLAdeg
tLArad
thetaTdeg
thetaTrad
u
uTM
v
w
with_k0L
with_kLL
xiLA
zero

Methods

SPLA(PLs[, NLgeom, plot])

Runs the Single Perturbation Load Approach (SPLA)

add_SPL(PL[, pt, thetadeg, increment])

Add a Single Perturbation Load \(\{{F_{PL}}_i\}\)

add_force(x, thetadeg, fx, ftheta, fz[, ...])

Add a punctual force

apply_shim(thetadeg, width, thick[, ncpts])

Distributes the axial load in order to simulate a shim

calc_fext([inc, kuk, silent])

Calculates the external force vector \(\{F_{ext}\}\)

calc_fint(c[, inc, m, return_u, silent])

Calculates the internal force vector \(\{F_{int}\}\)

calc_full_c(cu[, inc])

Returns the full set of Ritz constants

calc_kT(c[, inc, silent])

Calculates the tangent stiffness matrix

eigen([c, tol, kL, kG, combined_load_case])

Performs a non-linear eigenvalue analysis at a given state

exclude_dofs_matrix(k[, return_kkk, ...])

Makes the partition of the dofs for prescribed displacements

fit_Nxxtop(ts, us[, update_Nxxtop])

Adjusts the axial load distribution for a desired top edge displacement

from_DB(name)

Load cone/cylinder data from the local database

get_size()

Calculates the size of the stiffness matrices

lb([c, tol, combined_load_case])

Performs a linear buckling analysis

plot(c[, invert_x, plot_type, vec, ...])

Contour plot for a Ritz constants vector.

plotAbaqus(frame, fieldOutputKey, vec, ...)

Print a field output for a cylinder/cone model from Abaqus

save()

Save the ConeCyl object using pickle

static([NLgeom, silent])

Static analysis for cones and cylinders

strain(c[, xs, ts, gridx, gridt, inc])

Calculates the strain field

stress(c[, xs, ts, gridx, gridt, inc])

Calculates the stress field

uvw(c[, xs, ts, gridx, gridt, inc])

Calculates the displacement field

calc_k0
SPLA(PLs, NLgeom=True, plot=False)[source]#

Runs the Single Perturbation Load Approach (SPLA)

A set of non-linear results will be

Parameters:
PLs: list

The perturbation loads used to build the knock-down curve. It must be a list of float values.

NLgeombool, optional

Flag passed to the static() method that tells whether a geometrically non-linear analysis is to be performed.

Returns:
curveslist

The sequence of curves, one curve for each perturbation load given in the input parameter PLs. Each curve in the list is a dict object with the keys:

Key

Description

'wall_time_s'

The wall time for the non-linear analysis

'name'

The name of the curve. Ex: 'PL = 1. N'

'cs'

A list with a vector of Ritz constants for each load increment needed

'increments'

A list with the values of increments needed

'wPLs'

A list with the normal displacement at the perturbation load application point for each load increment

'uTMs'

A list containing the axial displacement for each load increment

'Fcs'

A list containing the axial reaction force for each load increment

Notes

The curves are stores in the ConeCyl parameter outputs['SPLA_curves'].

add_SPL(PL, pt=0.5, thetadeg=0.0, increment=False)[source]#

Add a Single Perturbation Load \(\{{F_{PL}}_i\}\)

Adds a perturbation load to the ConeCyl object, the perturbation load is a particular case of the punctual load with only a normal component.

Parameters:
PLfloat

The perturbation load value.

ptfloat, optional

The normalized position along the \(x\) axis in which the new SPL will be included.

thetadegfloat, optional

The angular position of the SPL in degrees.

incrementbool, optional

If this perturbation load should be incrementally applied in a non-linear analysis.

Notes

Each single perturbation load is added to the forces parameter of the ConeCyl object, which may be changed by the analyst at any time.

add_force(x, thetadeg, fx, ftheta, fz, increment=False)[source]#

Add a punctual force

Adds a force vector \(\{f_x, f_\theta, f_z\}^T\) to the forces parameter of the ConeCyl object.

Parameters:
xfloat

The \(x\) position.

thetadegfloat

The \(\theta\) position in degrees.

fxfloat

The \(x\) component of the force vector.

fthetafloat

The \(\theta\) component of the force vector.

fzfloat

The \(z\) component of the force vector.

incrementbool, optional

If this punctual force should be incrementally applied in a non-linear analysis.

apply_shim(thetadeg, width, thick, ncpts=10000)[source]#

Distributes the axial load in order to simulate a shim

The axial load distribution \({N_{xx}}_{top}\) will be adjusted such that the resulting displacement \(u\) at \(x=0\) (top edge) will look similar to a case where a shim is applied.

Parameters:
thetadegfloat

Position in degrees of the center of the shim.

widthfloat

Circumferential width of the shim.

thickfloat

Thickness of the shim.

ncptsint, optional

Number of control points used in the least-squares routine.

Returns:
tsnp.ndarray

Positions \(\theta\) of the control points.

usnp.ndarray

Displacements \(u\) of the control points.

Notes

This function changes the Nxxtop parameter of the current ConeCyl object. Returning ts and us is useful for post processing purposes only.

Examples

>>> ts, us = cc.apply_shim(0., 25.4, 0.1)
calc_fext(inc=1.0, kuk=None, silent=False)[source]#

Calculates the external force vector \(\{F_{ext}\}\)

Recall that:

\[\{F_{ext}\}=\{{F_{ext}}_0\} + \{{F_{ext}}_\lambda\}\]

such that the terms in \(\{{F_{ext}}_0\}\) are constant and the terms in \(\{{F_{ext}}_\lambda\}\) will be scaled by the parameter inc.

Parameters:
incfloat, optional

Since this function is called during the non-linear analysis, inc will multiply the terms \(\{{F_{ext}}_\lambda\}\).

kuknp.ndarray, optional

Obsolete, created for displacement controlled analyses, but the implementation has not been finished, see exclude_dofs_matrix().

silentbool, optional

A boolean to tell whether the msg messages should be printed.

Returns:
fextnp.ndarray

The external force vector

calc_fint(c, inc=1.0, m=1, return_u=True, silent=False)[source]#

Calculates the internal force vector \(\{F_{int}\}\)

The following attributes will affect the numerical integration:

Attribute

Description

ni_num_cores

int, number of cores used for the numerical integration

ni_method
str, integration method:

  • 'trapz2d' for 2-D Trapezoidal’s rule

  • 'simps2d' for 2-D Simpsons’ rule

nx

int, number of integration points along the \(x\) coordinate

nt

int, number of integration points along the \(\theta\) coordinate
Parameters:
cnp.ndarray

The Ritz constants that will be used to compute the internal forces.

incfloat, optional

Load increment, necessary to calculate the full set of Ritz constants using calc_full_c().

minteger, optional

A multiplier to be applied to nx and nt, if one whishes to use more integration points.

return_ubool, optional

If the internal force vector corresponsing to the unknown set of Ritz constants should be returned.

silentbool, optional

A boolean to tell whether the msg messages should be printed.

Returns:
fintnp.ndarray

The internal force vector.

calc_full_c(cu, inc=1.0)[source]#

Returns the full set of Ritz constants

When prescribed displacements take place the matrices and the Ritz constants are partitioned like:

k = | kkk    kku |
    | kuk    kuu |

and the corresponding Ritz constants:

c = | ck |
    | cu |

This function adds the set of known Ritz constants (ck) to the set of unknown (cu) based on the prescribed displacements.

Parameters:
cunp.ndarray

The set of unknown Ritz constants

incfloat, optional

Load increment, necessary to calculate the full set of Ritz constants.

Returns:
cnp.ndarray

The full set of Ritz constants.

calc_kT(c, inc=1.0, silent=False)[source]#

Calculates the tangent stiffness matrix

The following attributes will affect the numerical integration:

Attribute

Description

ni_num_cores

int, number of cores used for the numerical integration

ni_method
str, integration method:

  • 'trapz2d' for 2-D Trapezoidal’s rule

  • 'simps2d' for 2-D Simpsons’ rule

nx

int, number of integration points along the \(x\) coordinate

nt

int, number of integration points along the \(\theta\) coordinate
Parameters:
cnp.ndarray

The Ritz constants vector of the current state.

incfloat, optional

Load increment, necessary to calculate the full set of Ritz constants using calc_full_c().

silentbool, optional

A boolean to tell whether the msg messages should be printed.

Returns:
kTuusparse matrix

The tangent stiffness matrix corresponding to the unknown degrees of freedom.

eigen(c=None, tol=0, kL=None, kG=None, combined_load_case=None)[source]#

Performs a non-linear eigenvalue analysis at a given state

The following attributes of the ConeCyl object will affect the non-linear eigenvalue analysis:

Attribute

Description

num_eigenvalues

Number of eigenvalues to be extracted

num_eigvalues_print

Number of eigenvalues to print after the analysis is completed

Additionally, the non-linear analysis parameters described in static() will affect the integration of the non-linear matrices kL and kG if they are not given as input parameters.

Parameters:
combined_load_caseint or None, optional

It tells whether the linear buckling analysis must be computed considering combined load cases, each value will tell the algorithm to rearrange the linear matrices in a different way. The valid values are 1, or 2, where:

  • 1 : find the critical axial load for a fixed torsion load

  • 2 : find the critical axial load for a fixed pressure load

  • 3 : find the critical torsion load for a fixed axial load

Notes

The extracted eigenvalues are stored in the eigvals parameter of the ConeCyl object and the \(i^{th}\) eigenvector in the eigvecs[i-1, :] parameter.

exclude_dofs_matrix(k, return_kkk=False, return_kku=False, return_kuk=False)[source]#

Makes the partition of the dofs for prescribed displacements

Makes the following partition of a given matrix:

k = | kkk    kku |
    | kuk    kuu |
Parameters:
kscipy.sparse.coo_matrix

Matrix to be partitioned.

return_kkkbool, optional

If the region \(kkk\) must be returned.

return_kkubool, optional

If the region \(kku\) must be returned.

return_kukbool, optional

If the region \(kuk\) must be returned.

Returns:
outdict

A dict object containing the keys for the corresponding sub-matrices kkk, kku, kuk, kuu. The sub-matrix out['kuu'] is a scipy.sparse.csr_matrix, while the others are 2-D np.ndarray objects.

fit_Nxxtop(ts, us, update_Nxxtop=True)[source]#

Adjusts the axial load distribution for a desired top edge displacement

Parameters:
tsnp.ndarray

Corrdinates \(\theta\) of each control point.

usnp.ndarray

Desired displacement \(u\) for each control point.

update_Nxxtopbool, optional

Tells whether self.Nxxtop should be updated.

Returns:
Nxxtopnp.ndarray

The coefficients for the \({N_{xx}}_{top}\) function.

from_DB(name)[source]#

Load cone/cylinder data from the local database

Parameters:
namestr

A key contained in the ccs dictionary of module compmech.conecyl.conecylDB.

get_size()[source]#

Calculates the size of the stiffness matrices

The size of the stiffness matrices can be interpreted as the number of rows or columns, recalling that this will be the size of the Ritz constants’ vector \(\{c\}\), the internal force vector \(\{F_{int}\}\) and the external force vector \(\{F_{ext}\}\).

Returns:
sizeint

The size of the stiffness matrices.

lb(c=None, tol=0, combined_load_case=None)[source]#

Performs a linear buckling analysis

The following parameters of the ConeCyl object will affect the linear buckling analysis:

Attribute

Description

num_eigenvalues

Number of eigenvalues to be extracted

num_eigvalues_print

Number of eigenvalues to print after the analysis is completed
Parameters:
combined_load_caseint, optional

It tells whether the linear buckling analysis must be computed considering combined load cases, each value will tell the algorithm to rearrange the linear matrices in a different way. The valid values are 1, or 2, where:

  • 1 : find the critical axial load for a fixed torsion load

  • 2 : find the critical axial load for a fixed pressure load

  • 3 : find the critical torsion load for a fixed axial load

Notes

The extracted eigenvalues are stored in the eigvals parameter of the ConeCyl object and the \(i^{th}\) eigenvector in the eigvecs[i-1, :] parameter.

plot(c, invert_x=False, plot_type=1, vec='w', deform_u=False, deform_u_sf=100.0, filename='', ax=None, figsize=(3.5, 2.0), save=True, add_title=True, title='', colorbar=False, cbar_nticks=2, cbar_format=None, cbar_title='', cbar_fontsize=10, aspect='equal', clean=True, dpi=400, texts=[], xs=None, ts=None, gridx=300, gridt=300, num_levels=400, inc=1.0, vecmin=None, vecmax=None)[source]#

Contour plot for a Ritz constants vector.

Parameters:
cnp.ndarray

The Ritz constants that will be used to compute the field contour.

vecstr, optional

Can be one of the components:

  • Displacement: 'u', 'v', 'w', 'phix', 'phit',

    'magnitude'

  • Strain: 'exx', 'ett', 'gxt', 'kxx', 'ktt', 'kxt', 'gtz', 'gxz'

  • Stress: 'Nxx', 'Ntt', 'Nxt', 'Mxx', 'Mtt', 'Mxt', 'Qt', 'Qx'

deform_ubool, optional

If True the contour plot will look deformed.

deform_u_sffloat, optional

The scaling factor used to deform the contour.

invert_xbool, optional

Inverts the \(x\) axis of the plot. It may be used to match the coordinate system of the finite element models created using the desicos.abaqus module.

plot_typeint, optional

For cylinders only 4 and 5 are valid. For cones all the following types can be used:

  • 1: concave up (with invert_x=False) (default)

  • 2: concave down (with invert_x=False)

  • 3: stretched closed

  • 4: stretched opened (\(r \times \theta\) vs. \(H\))

  • 5: stretched opened (\(\theta\) vs. \(H\))

savebool, optional

Flag telling whether the contour should be saved to an image file.

dpiint, optional

Resolution of the saved file in dots per inch.

filenamestr, optional

The file name for the generated image file. If no value is given, the \(name\) parameter of the ConeCyl object will be used.

axAxesSubplot, optional

When ax is given, the contour plot will be created inside it.

figsizetuple, optional

The figure size given by (width, height).

add_titlebool, optional

If a title should be added to the figure.

titlestr, optional

If any string is given add_title will be ignored and the given title added to the contour plot.

colorbarbool, optional

If a colorbar should be added to the contour plot.

cbar_nticksint, optional

Number of ticks added to the colorbar.

cbar_format[ None | format string | Formatter object ], optional

See the matplotlib.pyplot.colorbar documentation.

cbar_fontsizeint, optional

Fontsize of the colorbar labels.

cbar_titlestr, optional

Colorbar title. If cbar_title == '' no title is added.

aspectstr, optional

String that will be passed to the AxesSubplot.set_aspect() method.

cleanbool, optional

Clean axes ticks, grids, spines etc.

xsnp.ndarray, optional

The \(x\) positions where to calculate the displacement field. Default is None and method _default_field is used.

tsnp.ndarray, optional

The theta positions where to calculate the displacement field. Default is None and method _default_field is used.

gridxint, optional

Number of points along the \(x\) axis where to calculate the displacement field.

gridtint, optional

Number of points along the \(theta\) where to calculate the displacement field.

num_levelsint, optional

Number of contour levels (higher values make the contour smoother).

incfloat, optional

Load increment, necessary to calculate the full set of Ritz constants using calc_full_c().

vecminfloat, optional

Minimum value for the contour scale (useful to compare with other results). If not specified it will be taken from the calculated field.

vecmaxfloat, optional

Maximum value for the contour scale.

Returns:
axmatplotlib.axes.Axes

The Matplotlib object that can be used to modify the current plot if needed.

plotAbaqus(frame, fieldOutputKey, vec, nodes, numel_cir, elem_type='S4R', ignore=[], ax=None, figsize=(3.3, 3.3), save=True, aspect='equal', clean=True, plot_type=1, outpath='', filename='', npzname='', pyname='', num_levels=400)[source]#

Print a field output for a cylinder/cone model from Abaqus

This function is intended to be used with models created by the DESICOS plug-in for Abaqus, where a mapped mesh is used and the models comparable to the models of compmech.conecyl.

The frame and nodes input types are described in Abaqus Scripting Reference Manual and they can be obtained through:

>>> frame = session.odbs['odb_name.odb'].steps['step_name'].frames[0]
>>> nodes = mdb.models['model_name'].parts['part_name'].nodes
Parameters:
frameOdbFrame

The frame from where the field output will be taken from.

fieldOutputKeystr

The field output key to be used. It must be available in frame.fieldOutputs.keys(). This function was tested with 'UT' and 'U' only.

vecstr

The displacement vector to be plotted: 'u', 'v' or 'w'.

nodesMeshNodeArray

The part nodes.

numel_cirint

The number of elements around the circumference.

elem_typestr, optional

The element type. The elements 'S4R', 'S4R5' where tested.

ignorelist, optional

A list with the node ids to be ignored. It must contain any nodes outside the mapped mesh included in parts['part_name'].nodes.

axAxesSubplot, optional

When ax is given, the contour plot will be created inside it.

figsizetuple, optional

The figure size given by (width, height).

savebool, optional

Flag telling whether the contour should be saved to an image file.

aspectstr, optional

String that will be passed to the AxesSubplot.set_aspect() method.

cleanbool, optional

Clean axes ticks, grids, spines etc.

plot_typeint, optional

See plot().

outpathstr, optional

Output path where the data from Abaqus and the plots are saved (see notes).

filenamestr, optional

The file name for the generated image file.

npznamestr, optional

The file name for the generated npz file.

pynamestr, optional

The file name for the generated Python file.

num_levelsint, optional

Number of contour levels (higher values make the contour smoother).

Returns:
outtuple

Where out[0] and out[1] contain the circumferential and meridional grids of coordinates and out[2] the corresponding field output.

Notes

The data is saved using np.savez() into outpath as abaqus_output.npz with an accompanying script for plotting abaqus_output_plot.py, very handy when Matplotlib is not importable from Abaqus.

save()[source]#

Save the ConeCyl object using pickle

Notes

The pickled file will have the name stored in ConeCyl.name followed by a '.ConeCyl' extension.

static(NLgeom=False, silent=False)[source]#

Static analysis for cones and cylinders

The analysis can be linear or geometrically non-linear. See Analysis for further details about the parameters controlling the non-linear analysis.

Parameters:
NLgeombool

Flag to indicate whether a linear or a non-linear analysis is to be performed.

silentbool, optional

A boolean to tell whether the msg messages should be printed.

Returns:
cslist

A list containing the Ritz constants for each load increment of the static analysis. The list will have only one entry in case of a linear analysis.

Notes

The returned cs is stored in the cs parameter of the ConeCyl object. The actual increments used in the non-linear analysis are stored in the increments parameter.

strain(c, xs=None, ts=None, gridx=300, gridt=300, inc=1.0)[source]#

Calculates the strain field

Parameters:
cnp.ndarray

The Ritz constants vector to be used for the strain field calculation.

xsnp.ndarray, optional

The \(x\) coordinates where to calculate the strains.

tsnp.ndarray, optional

The \(\theta\) coordinates where to calculate the strains, must have the same shape as xs.

gridxint, optional

When xs and ts are not supplied, gridx and gridt are used.

gridtint, optional

When xs and ts are not supplied, gridx and gridt are used.

incfloat, optional

Load increment, necessary to calculate the full set of Ritz constants using calc_full_c().

stress(c, xs=None, ts=None, gridx=300, gridt=300, inc=1.0)[source]#

Calculates the stress field

Parameters:
cnp.ndarray

The Ritz constants vector to be used for the strain field calculation.

xsnp.ndarray, optional

The \(x\) coordinates where to calculate the strains.

tsnp.ndarray, optional

The \(\theta\) coordinates where to calculate the strains, must have the same shape as xs.

gridxint, optional

When xs and ts are not supplied, gridx and gridt are used.

gridtint, optional

When xs and ts are not supplied, gridx and gridt are used.

incfloat, optional

Load increment, necessary to calculate the full set of Ritz constants using calc_full_c().

uvw(c, xs=None, ts=None, gridx=300, gridt=300, inc=1.0)[source]#

Calculates the displacement field

For a given full set of Ritz constants c, the displacement field is calculated and stored in the parameters u, v, w, phix, phit of the ConeCyl object.

Parameters:
cfloat

The full set of Ritz constants

xsnp.ndarray

The \(x\) positions where to calculate the displacement field. Default is None and method _default_field is used.

tsnp.ndarray

The theta positions where to calculate the displacement field. Default is None and method _default_field is used.

gridxint

Number of points along the \(x\) axis where to calculate the displacement field.

gridtint

Number of points along the \(theta\) where to calculate the displacement field.

incfloat, optional

Load increment, necessary to calculate the full set of Ritz constants using calc_full_c().

Returns:
outtuple

A tuple of np.ndarrays containing (u, v, w, phix, phit).

Notes

The returned values u`, v, w, phix, phit are stored as parameters with the same name in the ConeCyl object.

Non-linear analysis (compmech.conecyl.non_linear)#

Cone / Cylinder Database (compmech.conecyl.conecylDB)#

This database is a group of three dictionaries defines as following:

  • laminaprops: the material properties for an orthotropic lamina

  • allowables: the allowables for an orthotropic lamina

  • ccs: the geometric, stacking sequence and material data for a conical / cylindrical structure

One could import these dictionaries doing:

from compmech.conecyl.conecylDB import laminaprops, allowables, ccs

Models’ Database (compmech.conecyl.modelDB)#

Used to configure the main parameters for each implemented model.

compmech.conecyl.modelDB.get_linear_matrices(cc, combined_load_case=None)[source]#

Obtain the right functions to calculate hte linear matrices for a given model.

The model parameter of the ConeCyl object is used to search for the functions fG0, fG0_cyl, fkG0, fkG0_cyl, and the matrix k0edges is calculated, when applicable.

Parameters:
cccompmech.conecyl.ConeCyl

The ConeCyl object.

combined_load_caseint, optional

As explained in the ConeCyl.lb() method, the integer indicating which combined load case should be used. Default is None.

Returns:
outtuple

A tuple containing (fk0, fk0_cyl, fkG0, fkG0_cyl, k0edges).

Imperfections (compmech.conecyl.imperfections)#

Convenient routines to transform discrete measured data into continuous functions applicable in semi-analytical analyses.

The implemented geometric imperfection of this module is the one representing the normal displacement of the mid-surface and will be called Measured Geometric Imperfection (MGI).

The non-linear analysis using a MGI require the calculation of the initial imperfection field, called \(w_0\), and the corresponding partial derivatives \({w_0}_{,x}\) and \({w_0}_{,\theta}\).

The function calc_c0 described below is implemented to find the best-fit for a given imperfection data. Three different approximation functions can be selected using the funcnum parameter.

The imperfection file should be in the form:

theta1  height1  imp1
theta2  height2  imp2
...
thetan  heightn  impn

where height is measured from the bottom to the top of the cylinder or cone, parallel to the axial axis, and theta is the circumferential coordinate \(\theta\) measured as shown in the semi-analytical model.

When implementing a non-linear analysis algorithm, see for example fsdt_donnell_bc1_nonlinear.pyx , the functions to calculate the partial derivatives of the geometric imperfection are accessible using Cython

from compmech.conecyl.imperfection.mgi cimport cfw0x, cfw0t

The cfw0x and cfw0t function headers are:

cdef void cfw0x(double x, double t, double *c0, double L,
                int m, int n, double *w0xref, int funcnum) nogil

cdef void cfw0t(double x, double t, double *c0, double L,
                int m, int n, double *w0tref, int funcnum) nogil

where c0 is the array containing the coefficients of the approximation function, L is the meridional length of the cone or cylinder (as shown here), x and t the \(x\) and \(\theta\) coordinates, m and n the number of terms of the approximation series as described above, and finally, w0xref and w0tref are pointers to double variables.

compmech.conecyl.imperfections.imperfections.calc_c0(path, m0=40, n0=40, funcnum=2, sample_size=None, maxmem=8, save=True, offset_w0=None)[source]#

Find the coefficients \(c_0\) that best fit the \(w_0\) function.

The measured data will be fit using one of the following functions, selected using the funcnum parameter:

funcnum=1

\[w_0 = \sum_{i=1}^{m_0}{ \sum_{j=0}^{n_0}{ c_{ij}^a sin{b_x} sin{b_\theta} +c_{ij}^b sin{b_x} cos{b_\theta}}}\]

funcnum=2 (default)

\[w_0 = \sum_{i=0}^{m_0}{ \sum_{j=0}^{n_0}{ c_{ij}^a cos{b_x} sin{b_\theta} +c_{ij}^b cos{b_x} cos{b_\theta}}}\]

funcnum=3

\[w_0 = \sum_{i=0}^{m_0}{ \sum_{j=0}^{n_0}{ c_{ij}^a sin{b_x} sin{b_\theta} +c_{ij}^b sin{b_x} cos{b_\theta} +c_{ij}^c cos{b_x} sin{b_\theta} +c_{ij}^d cos{b_x} cos{b_\theta}}}\]

where:

\[ \begin{align}\begin{aligned}b_x = i \pi \frac x L_{points}\\b_\theta = j \theta\end{aligned}\end{align} \]

where \(L_{points}\) represents the difference between the maximum and the height values in the imperfection file divided by the cosine of the semi-vertex angle:

\[L_{points} = \frac{H_{max} - H_{min}}{cos(\alpha)} = \frac{H_{points}}{cos(\alpha)}\]

In this form \({}^x/_{L_{points}}\) will vary from \(0.\) (at the top) to \(1.\) (at the bottom).

Note

Note that if the measured sample does not cover all the height, it will be stretched.

The approximation can be written in matrix form as:

\[w_0 = [g] \{c\}\]

where \([g]\) carries the base functions and \({c}\) the respective amplitudes. The solution consists on finding the best \(\{c\}\) that minimizes the least-square error between the measured imperfection pattern and the \(w_0\) function.

Parameters:
pathstr or numpy.ndarray

The path of the file containing the data. Can be a full path using r"C:\Temp\inputfile.txt", for example. The input file must have 3 columns: \(\theta\), \(height\), \(imp\); expressed in Cartesian coordinates.

This input can also be a numpy.ndarray object, with \(\theta\), \(height\), \(imp\) in each corresponding column.

m0int

Number of terms along the meridian (\(x\)).

n0int

Number of terms along the circumference (\(\theta\)).

funcnumint, optional

As explained above, selects the base functions used for the approximation.

sample_sizeint or None, optional

Specifies how many points of the imperfection file should be used. If None all points will be used in the computations.

maxmemint, optional

Maximum RAM memory in GB allowed to compute the base functions. The scipy.interpolate.lstsq will go beyond this limit.

savebool, optional

If True saves the calculated coefficients in the compmech/conecyl/imperfections/c0 folder.

Returns:
outnumpy.ndarray

A 1-D array with the best-fit coefficients.

compmech.conecyl.imperfections.imperfections.fw0(m0, n0, c0, xs_norm, ts, funcnum=2)[source]#

Calculates the imperfection field \(w_0\) for a given input.

Parameters:
m0int

The number of terms along the meridian.

n0int

The number of terms along the circumference.

c0numpy.ndarray

The coefficients of the imperfection pattern.

xs_normnumpy.ndarray

The meridian coordinate (\(x\)) normalized to be between 0. and 1..

tsnumpy.ndarray

The angles in radians representing the circumferential coordinate (\(\theta\)).

funcnumint, optional

The function used for the approximation (see the calc_c0 function)

Notes

The inputs xs_norm and ts must be of the same size.

If funcnum==1 or funcnum==2 then size=2, if funcnum==3 then size=4 and the inputs must satisfy c0.shape[0] == size*m0*n0.