Bardell’s hierarchical functions#
Every domain of panels, a Shell, a panel of a
MultiDomain assembly or the flange of a BladeStiff2D,
approximates its displacement fields with the Ritz method, as series of
products of one-dimensional functions along \(x\) and \(y\). The functions are
the hierarchical Legendre-based polynomials in Rodrigues form applied by
Bardell [bardell1991] in the hierarchical finite element method, and used
as approximation functions of semi-analytical models of plates, shells and
stiffened panels by Castro and Donadon [castro2017Multidomain]. They are
coded in C++ up to 30 terms and accessed from Python with
panels.bardell, where their expressions are given.
The functions#
The functions \(P_i(\xi)\), \(i = 0, 1, \dots, 29\), are defined over the natural coordinate \(-1 \le \xi \le 1\), which maps the length of the domain, \(\xi = 2x/a - 1\) along \(x\) and \(\eta = 2y/b - 1\) along \(y\). The first four terms are cubic Hermite polynomials and the terms \(i \ge 4\) are polynomials of degree \(i\) that vanish together with their first derivative at both ends. Only the first four terms therefore carry a value or a slope at the ends, each one at one end only:
term \(i\) |
\(P_i(-1)\) |
\(P_i'(-1)\) |
\(P_i(+1)\) |
\(P_i'(+1)\) |
flag |
|---|---|---|---|---|---|
0 |
1 |
0 |
0 |
0 |
|
1 |
0 |
\(1/2\) |
0 |
0 |
|
2 |
0 |
0 |
1 |
0 |
|
3 |
0 |
0 |
0 |
\(1/2\) |
|
\(\ge 4\) |
0 |
0 |
0 |
0 |
where the prime is the derivative with respect to \(\xi\), and
\(\partial/\partial x = (2/a) \, \partial/\partial \xi\). The set is
hierarchical: increasing the number of terms, m along \(x\) and n
along \(y\), adds functions without changing the previous ones. Evaluated in
floating point, the terms of high degree carry a round-off that grows with
the degree, below \(10^{-4}\) of the largest value of each term at the ends,
as asserted by tests/test_bardell.py:
def test_end_values_and_flags():
flags = (1., 1., 1., 1.) # x1t, x1r, x2t, x2r
xis = np.linspace(-1, 1, 2001)
fmax = np.abs([calc_vec_f(xi, *flags) for xi in xis]).max(axis=0)
fpmax = np.abs([calc_vec_fp(xi, *flags) for xi in xis]).max(axis=0)
# only the first four terms have a value or a slope at the ends
expected = {
# term: (f(-1), f'(-1), f(+1), f'(+1))
0: (1., 0., 0., 0.),
1: (0., 0.5, 0., 0.),
2: (0., 0., 1., 0.),
3: (0., 0., 0., 0.5),
}
for i in range(NUM):
ends = (calc_f(i, -1., *flags), calc_fp(i, -1., *flags),
calc_f(i, +1., *flags), calc_fp(i, +1., *flags))
# round-off of the evaluation in floating point, which grows with
# the degree of the polynomial
tol = 1e-4*np.array([fmax[i], fpmax[i], fmax[i], fpmax[i]])
assert np.all(np.abs(np.subtract(ends, expected.get(i, (0., 0., 0.,
0.)))) <= tol), i
# a flag 0 removes its term: here a clamped end at xi = -1
x1t, x1r = 0., 0.
for xi in np.linspace(-1, 1, 7):
f = calc_vec_f(xi, x1t, x1r, 1., 1.)
assert f[0] == 0. and f[1] == 0.
assert f[2] == calc_f(2, xi, 1., 1., 1., 1.)
assert np.allclose(calc_vec_f(-1., x1t, x1r, 1., 1.), 0., atol=1e-4*fmax)
assert np.allclose(calc_vec_fp(-1., x1t, x1r, 1., 1.), 0.,
atol=1e-4*fpmax)
# any other non-zero flag only scales its term
assert np.isclose(calc_fp(1, -1., 1., 2., 1., 1.), 1.)
Boundary flags#
The four flags multiply the first four terms: x1t and x1r the terms
of the value (translation) and of the slope (rotation) at \(\xi = -1\),
x2t and x2r those at \(\xi = +1\). A flag 0. removes its term,
which makes the value or the slope of every function of the series zero at
that end, an essential boundary condition imposed exactly on the whole
edge; a flag 1. keeps the term, and the corresponding natural boundary
condition follows from the energy. Any other non-zero flag only scales the
term.
Each field of a domain has its own flags along \(x\) and along \(y\),
attributes of Shell named after the edge, the field and the kind
of term:
x1u,x1ur,x2u,x2ur: the flagsx1t,x1r,x2t,x2rof the functions of \(u\) along \(x\), i.e. at the edges \(x = 0\) and \(x = a\);y1u,y1ur,y2u,y2ur: the same along \(y\), at the edges \(y = 0\) and \(y = b\);likewise for \(v\) and \(w\), and for the rotations \(\phi_x\) and \(\phi_y\) of the models based on shear deformation theories, e.g.
x1phix,y2phiyr.
For example, x1w = 0 imposes \(w = 0\) and x1wr = 0 imposes \(w_{,x} =
0\) on the edge \(x = 0\). The combinations of the flags that give the usual
supports are listed in Boundary conditions. A removed term keeps its
place in the vector of the Ritz constants, see Shell.get_size(), with
null rows and columns in the matrices, which the solvers of structsolve
remove before solving.
Integrals and structural matrices#
An entry of a structural matrix is the integral over the domain of the product of two terms of the series, \((i, j)\) and \((k, l)\), or of their derivatives, multiplied by the constitutive coefficients. When these coefficients are constant over the domain and the integration domain is the whole shell, the double integral separates in the two directions:
with \(p, q, r, s \le 2\), where \(f\) and \(g\) carry the flags of the fields of
each term along \(x\) and \(y\). The six one-dimensional integrals,
calc_integral_ff(), calc_integral_ffp(),
calc_integral_ffpp(), calc_integral_fpfp(),
calc_integral_fpfpp() and calc_integral_fppfpp(), are
tabulated in closed form, generated with SymPy in rational arithmetic and
exact. They give the closed-form (analytical) matrices of the models, e.g.
the linear constitutive, geometric, mass and aerodynamic matrices of a
Shell, Shell.calc_kC(), Shell.calc_kG(),
Shell.calc_kM() and Shell.calc_kA(), and the penalty
connections along the edges, see panels.multidomain.connections.
When the integrand does not separate, i.e. for the non-linear matrices and
internal forces, for a stiffness given per integration point by
ABDnxny, for the follower pressure loads, or for a partial integration
domain, see Shell.is_partial_domain(), the matrices are integrated
numerically with nx by ny Gauss-Legendre points,
Shell.integration_points(), evaluating the functions and their first
and second derivatives at each point, see calc_vec_f(),
calc_vec_fp() and calc_vec_fpp(). The same functions give the
field outputs, Shell.uvw(), Shell.strain() and
Shell.stress(), through the field modules
panels.models.clpt_field and panels.models.fsdt_tsdt_field.
The values of the functions at the integration points can be tabulated once
per kernel, see Approximation functions tabulated at the integration points.
The closed-form integrals agree with a Gauss-Legendre quadrature of 30 points, exact for the products of two terms up to the round-off of the functions, and Eq. (1) gives the physical integrals:
def test_integrals_against_quadrature():
# products of two terms are polynomials of degree up to 58, integrated
# exactly by 30 Gauss-Legendre points
xis, ws = leggauss(30)
fx = (1., 0., 1., 1.) # x1t, x1r, x2t, x2r of the term i
fy = (0., 1., 1., 0.) # x1t, x1r, x2t, x2r of the term j
F = np.array([calc_vec_f(xi, *fx) for xi in xis])
F_xi = np.array([calc_vec_fp(xi, *fx) for xi in xis])
F_xixi = np.array([calc_vec_fpp(xi, *fx) for xi in xis])
G = np.array([calc_vec_f(xi, *fy) for xi in xis])
G_xi = np.array([calc_vec_fp(xi, *fy) for xi in xis])
G_xixi = np.array([calc_vec_fpp(xi, *fy) for xi in xis])
for func, A, B in ((calc_integral_ff, F, G),
(calc_integral_ffp, F, G_xi),
(calc_integral_fpfp, F_xi, G_xi),
(calc_integral_fppfpp, F_xixi, G_xixi)):
quad = (A*ws[:, None]).T @ B
closed = np.array([[func(i, j, *fx, *fy) for j in range(NUM)]
for i in range(NUM)])
assert np.allclose(closed, quad, rtol=1e-6,
atol=1e-6*np.abs(quad).max()), func.__name__
# physical integral over 0 <= x <= a of the product of the derivatives
# with respect to x of two terms
a = 0.7
i, j = 5, 7
F_x = F_xi*2/a
ref = np.sum(ws*a/2*F_x[:, i]*F_x[:, j])
assert np.isclose((2/a)**2*a/2*calc_integral_fpfp(i, j, *fx, *fx), ref)