Application of non-classical operational calculus to indicate hazards in numerical solutions of engineering problems

The differential elastic deflection line of a beam axis is widespread in engineering practice. Non-uniform settlement of a tank yields additional stresses in its structure, possibly leading to failure, that is deterioration or damage, triggering a significant threat to petroleum storage tanks. In order to capture these phenomena, the tank bottom plate is modelled by a simply supported beam resting either on stiff bedding [7] or on elastic bedding [18]. Other models are applied to assess internal forces in the storage tank structure [7], [17], [18]. The US Standard [3] makes use of working out design criteria of these models for large storage tanks.

Deflection of a rectangular plate, assumed a preliminary small cylindrical curvature [16], is expressed by a differential equation: (1)d2ydx2+k2y=−k2f(x),{{{d^2}y} \over {d{x^2}}} + {k^2}y = - {k^2}f(x), where f(x) is a preliminary deflection function and k reflects physical (mechanical) properties of the plate.

Equations of a similar sort are applied to model the intercomponent interaction of a vehicle/railway track/track bed/subsoil system. In railway engineering domain, the beam models on elastic foundation are applied in the range stated above, see [4], [5], [13], [14], and in railway track geometry issues [11]. The differential equation (1) is also applied in quantum mechanics to deal with the particle vibration in potential well analysis. The relevance of this equation is available.

A wide application range makes it reasonable to propose a new solution method of the equations, keeping in mind critical solution aspects, especially in numerical approach.

The modern design engineering practice requires a fluent use of CAD (Computer-aided design) methods in computational mechanics. Contemporary constitutive theories bring about non-linearity in the system of differential equations defining the initial boundary value problem. The equivalent non-linear equation systems require complex numerical procedures.

A widespread use of CAD methods linked with finite element method (FEM) or finite difference method (FDM) provides a vehicle to develop commercial software. The users, however, are frequently banned from numerical procedures; thus their interaction with the software is limited to operate the so-called black box, that is, they enter the input data and receive the final solution. The problem appears when the results of the computing process are not the ones expected by the designer.

It is a generally non-trivial task to relevantly interpret the computational, discretisation-based results. Whilst attention is paid to the finite element dimensions, marking the FE (finite element) mesh sensibility, the computational results may strongly differ from the observations, for example, in the steel plate cyclic load test, involving shear degradation [2].

Solution ambiguity may appear in the simplest possible cases, for example, beam deflection or other engineering issues incorporating similar mathematical models.

Most boundary problems make it extremely difficult to fulfil the conditions of a well-stated boundary problem, that is, to fulfil Hadamard's conditions. These conditions concern solution existence, uniqueness and stability, the second one is the focus of this article.

Non-classical operational calculus is basically aimed at creating uniform mathematical methods to solve specified forms of differential, difference or integral equations. It naturally generalises classified mathematical issues of a Heaviside idea or Mikusinski operators.

In classification of Mathematical Reviews^® Mathematical Subject Classification (MSC2010), they are described under No. 44 Integral Transforms, Operational Calculus, and then in detail under No. 44A40 – Calculus of Mikusinski and other operational calculi.

The non-classical operational calculus allows to uniformly describe and analyse a variety of engineering (mechanical) problems concerning continuum or discrete fields, including methods to solve differential or difference equations and their systems. Three linear operations, S, T_q and s_q, and two linear spaces, L¹ and L⁰, are defined, assuming L¹ ⊂ L⁰, to read as (2)S:L1→L0, Tq:L0→L1, sq:L1→L0.\matrix{{S:{L^1} \to {L^0},} \hfill & {{T_q}:{L^0} \to {L^1},} \hfill & {{s_q}:{L^1} \to {L^0}.} \hfill \cr}

The following conditions hold: (3)STq=id, S(L1)=L0, TqS=id−sq,\matrix{{S{T_q} = id,} \hfill & {S({L^1}) = {L^0},} \hfill & {{T_q}S = id - {s_q},} \hfill \cr} where q ∈ Q is the index set. The operation characteristics S, T_q and s_q are included in [8, 9].

Note that many operations exist to fulfil the conditions stated before; their examples are shown in [8]. In the article, it is sufficient to mention two representations only.

In a continuous field case, the following conditions hold: (4)S=ddx, Tx0=∫x0x, Sx0=|x=x0, x0∈<a,b>\matrix{{S = {d \over {dx}},} \hfill & {{T_{{x_0}}} = \int_{{x_0}}^x,} \hfill & {{S_{{x_0}}} = {|_{x = {x_0}}},} \hfill & {{x_0} \in < a,b >} \hfill \cr} and (5)L0=C0(<a,b>,R), L1=C1(<a,b>,R), q=x0∈<a,b>=Q.\matrix{{{L^0} = {C^0}(< a,b >,R),} \hfill & {{L^1} = {C^1}(< a,b >,R),} \hfill & {q = {x_0} \in < a,b > = Q}.}

In a discrete field, that is, the space C(N) of a series x={x_k}=(x₀, x₁, x₂, ...), it can be considered as: (6)S{xk}={xk+1−xk}=Δ{xk}S\left\{{{x_k}} \right\} = \left\{{{x_{k + 1}} - {x_k}} \right\} = \Delta \left\{{{x_k}} \right\}(7)Tk0{xk}={−xk0−1−xk0−2−…−xkdlak0>k0dlak=k0xk0+xk0+1+…+xk−1dlak0<k{T_{{k_0}}}\left\{{{x_k}} \right\} = \left\{{\matrix{{- {x_{{k_0} - 1}} - {x_{{k_0} - 2}} - \ldots - {x_k}} & {dla} & {{k_0} > k} \cr 0 & {dla} & {k = {k_0}} \cr {{x_{{k_0}}} + {x_{{k_0} + 1}} + \ldots + {x_{k - 1}}} & {dla} & {{k_0} < k} \cr}} \right.(8)Sk0{xk}={xk0}.{S_{{k_0}}}\left\{{{x_k}} \right\} = \left\{{{x_{{k_0}}}} \right\}.

Whilst there exist an unambiguously defined solution to the equation, (9)Sx=Rx, x∈L1,\matrix{{Sx = Rx,} \hfill & {x \in {L^1},} \hfill \cr} with a condition (10)sqx=c∈KerS,{s_q}x = c \in KerS, where R is the endomorphism of spaces L¹ and L⁰ commutative with operations S and s_q, the solution may be called a generalised exponential function, which is defined as (11)x=eRtqc.x = {e^{R{t_q}}}c.

Whilst the model features an operation S=ddxS = {d \over {dx}} , a simple exponential function may be eventually obtained.

In the case of endomorphism R=a∂∂xR = a{\partial \over {\partial x}} , a ≠ 0 and a two-variable function f={f(t,x)} with operations S=∂∂tS = {\partial \over {\partial t}} , Tq=∫0t{T_q} = \int_0^t , s_q = |_t=0, where L⁰ and L¹ are the function spaces of relevant classes [9], the character of a generalised exponential function is (12)ea∂∂xtqc={c(x+at)},c={c(x)}∈Ker∂∂t.\matrix{{{e^{a{\partial \over {\partial x}}{t_q}}}c = \left\{{c(x + at)} \right\},} \hfill & {c = \left\{{c(x)} \right\}} \hfill \cr} \in Ker{\partial \over {\partial t}}.

It is easy to check that the function (13)eiRtqc−e−iRtqc2i,{{{e^{iR{t_q}}}c - {e^{- iR{t_q}}}c} \over {2i}}, fulfil the conditions stated by the equation: (14)S2x=−R2x,{S^2}x = - {R^2}x, whilst s_qx=0 and s_qSx=c. The function (13) is also called a generalised trigonometric function, defined as sinRt_qc.

Given a discrete field with the operations, Δ,Tk0,sk0, k0=0\matrix{{\Delta,{T_{{k_0}}},{s_{{k_0}}},} \hfill & {{k_0} = 0} \hfill \cr} (see eq. 7, 8, 9) and endomorphism R, which means multiplication by a number a and continuous series {x₀}, the exponential function is a series defined as (15)eatq{x0}={(1+a)kx0}.{e^{a{t_q}}}\left\{{{x_0}} \right\} = \left\{{{{(1 + a)}^k}{x_0}} \right\}.

Figure 1 shows the diagram of series (15) for different values of number a.

Figure 1

The shape and characteristics of a generalised exponential function e^Rt_qc are predictable. Given a complex number domain, a well-known exponential function is periodic in contrast to the real domain.

Equation (13) leads to the following: (16)sinatq{x0}{(1+ai)k−(1−ai)k2ix0},\sin a{t_q}\left\{{{x_0}} \right\}\left\{{{{{{(1 + ai)}^k} - {{(1 - ai)}^k}} \over {2i}}{x_0}} \right\}, which is next transformed into the series form: (17)sinatq{x0}={(1+a2)ksin(arctga)x0}==x0{(1+a2)ksin(k arctga)},\matrix{{\sin a{t_q}\left\{{{x_0}} \right\} = \left\{{{{\left({\sqrt {1 + {a^2}}} \right)}^k}\sin ({\rm{arctg}}a){x_0}} \right\} =} \cr {= {x_0}\left\{{{{\left({\sqrt {1 + {a^2}}} \right)}^k}\sin (k\,{\rm{arctg}}a)} \right\},} \cr} which is presented in Figure 2 in a constant sequence case {x₀}={1}.

Figure 2

The diagrams in Figures 1 and 2 mark sequential solutions to some equations, possibly to represent stable and unstable problems on both modelling and solution phases, including the use of approximate methods [1].

The equation for the deflection of a beam of length l and variable cross section is subjected to axial compressive force P and distributed lateral load q(x). The differential equations of the beam are (18)M″(x)+PB(x)M(x)=q(x),M^{''}(x) + {P \over {B(x)}}M(x) = q(x),(19)B(x)y″=M(x),B(x)y^{''} = M(x), where M(x) means bending moment and B(x)=EJ(x) is a variable beam bending stiffness. This coordination of axis direction is accepted in soil mechanics; however, considering classical mechanics, the minus sign cannot to be omitted, and then B(x) y″= −M(x). The case of the operational calculus model with an operation S=ddtS = {d \over {dt}} makes it possible to redefine the problem: (20)S2M+PB(x)M=q(x){S^2}M + {P \over {B(x)}}M = q(x)(21)B(x)S2y=M(x).B(x){S^2}y = M(x).

Equations (20) and (21) correspond to equations (18) and (19), respectively. Whilst variable stiffness is assumed, equation (18) is analysed by operational methods described in [11] by applying a chain-type connection of two generalised inertial parts of the first order.

Owing to a specific supporting variant, the equations are complemented by relevant boundary conditions. In the simply supported case, the following conditions hold: (22)M(0)=0 M(l)=0.\matrix{{M(0) = 0} & {M(l) = 0.} \cr}

It is a considerable task to analytically find a deflection function of the beam, thus approximative methods are open, however, excluding a ‘blind user’ routine.

A related case concerns neutral axis definition in the cross-sectional plane of a compressed bar, for example, a concrete post covered with steel or composite coating (reinforced by carbon or glass fibres). The neutral axis passes through the cross section of the beam (post); it is the locus of zero normal stresses and axial strains, that is, the undeflected region. Whilst no axial forces act in a non-curved bar made of homogeneous material, the neutral axis at each cross section intercepts the centroid. The presence of cross-sectional axial forces makes the neutral axis shift. The location of neutral axis is linked with structural state, possibly detected in every time instant based on the deformation measurement.

The application of numerical method to assess deflections of bars brings about a number of disadvantages. Whilst a differential equation is represented by its approximation, for example, a difference ratio, it is probable to get quantitatively and qualitatively incorrect results even after entering a simple input differential equation form.

Equation (18) or its equivalent form (20) may be transformed, assuming (23)PB(x)=1, q(x)=0.\matrix{{{P \over {B(x)}} = 1,} \hfill & {q(x) = 0.} \hfill \cr}

Here the assumptions are still general. Whilst in operational calculus, the operation S can be replaced with formulas α=ddt\alpha = {d \over {dt}} and α2=PB(x){\alpha ^2} = {P \over {B(x)}} of operational origin; the coefficient at M(x) is unity. A simple equation can be written as (24)M″(x)+M(x)=0.M^{''}(x) + M(x) = 0.

The derivative ddx{d \over {dx}} is represented by its finite approximation, the difference ratio, by assuming a small h parameter: (25)Δhy(x)=y(x+h)−y(x)h;{\Delta _h}y(x) = {{y(x + h) - y(x)} \over h}; subsequently, the second derivative d2dx2{{{d^2}} \over {d{x^2}}} is approximated by a doubled Δ_h operation, that is, Δh2\Delta _h^2 .

Stating x=kh and M(kh)=M_k, the differential equation (24) turns into difference equation (recurrence equation); thus a continuous–discrete transformation is conducted here: (26)1h2(Mk+2−2Mk+1+Mk)+Mk=0.{1 \over {{h^2}}}\left({{M_{k + 2}} - 2{M_{k + 1}} + {M_k}} \right) + {M_k} = 0.

The above equation may be expressed by means of operation S stated by (6): (27)S2Mk+h2Mk=0,{S^2}{M_k} + {h^2}{M_k} = 0, including the conditions (28)M0=0, Mk0=0,\matrix{{{M_0} = 0,} \hfill & {{M_{{k_0}}} = 0,} \hfill \cr} where k₀ is the assumed natural number, fulfilling condition l=k₀h, complied with the second condition (22).

Let l=k₀h≠nπ and n=1,2,...,n₀, where n₀ is assumed. Then the equation (24) considered in continuous field with conditions (22) has one solution, and the solution is equal to zero.

The solution of the recurrence equation (27) with condition M₀=0 according to (17) can be written as (29)Mk=α(1+h2)ksin(k arctgh),α∈ℝ.\matrix{{{M_k} = \alpha {{{\sqrt {\left(1 + {h^2}\right)^k}}}}\sin (k\,{\rm{arctg}}h),} \hfill & {\alpha \in {\mathbb R}} \hfill \cr}.

Introducing the second condition M_k0=0 leads to (30)Mk0=α(1+h2)k0sin(k0 arctgh),α∈ℝ.\matrix{{{M_{k_0}} = \alpha {{{\sqrt {\left(1 + {h^2}\right)^{k_0}}}}}\sin ({k_0}\,{\rm{arctg}}h),} \hfill & {\alpha \in {\mathbb R}} \hfill \cr}.

If it happens, (31)k0 arctgh=nπ,n=1,2,…,\matrix{{{k_0}\,{\rm{arctg}}h = n\pi,} \hfill & {n = 1,2, \ldots,} \hfill \cr}equation (27) with conditions (28) shows an infinity of solutions.

If condition (31) holds, the difference equation approximates the differential equation with conditions (28) and shows an infinite solution range, keeping in mind that the differential equation shows exactly a single zero solution.

The presented numerical methodology brings an essential quantitative and qualitative problem reclassification, that is, the solution passes from a differential unique form to an ambiguous sort. The variation between the exact and the approximate solutions (the latter presents values at nodes) is larger whilst taking a large factor α.

The work [8] shows the operational calculus model, detecting the form S²M+M=0 without a solution.

Both Kerr's model analysis and Winkler–Pasternak foundation concept regarded in [15] showed another practical application of (27); it used a shear layer (soil) with geosynthetic membranes. The presented model highly resembles the real geosynthetic-soil systems, possible road or railway flexible pavements [10].

Kerr introduced a model, even a large spectrum of models made of Winkler, Pasternak and Héteny layers. The equation of Kerr's foundation model is similar to the form presented in [10] and Figure 3 whilst not regarding a geosynthetic layer, that is, T_p=0 and G_t=G_b=G; it is written as follows: (32)q(1+kskt)−Gkt∇2q=ksw−G∇2w,q\left({1 + {{{k_s}} \over {{k_t}}}} \right) - {G \over {{k_t}}}{\nabla ^2}q = {k_s}w - G{\nabla ^2}w, where q=q(x,y) is a model surface load, w=w(x,y) is a total surface deflection, k_t, k_s are upper and lower elastic foundation constants, and G is the shear modulus.

Figure 3

In a one-dimensional case, the model equation (32) modifies the Laplace operator ∇², which becomes a common second derivative d2dx2{{d^2} \over {d{x^2}}} . Whilst no load occurs, the differential equation (24) appears. Whilst a complete equation (32) is taken, its numerical solution is possible.

Sequences (15) are applied to determine responses of discrete two-dimensional systems [9]. Whilst discretising a continuous system of distributed parameters, the sequences can be applied to approximate the system response. The systems cover the performance of civil engineering plate and shell structures, wave propagation in soils and so on [12].

The presented discrete time signals are relevant to approximately determine the response of some class of differential problems adjusted to mechanics, for example, axially compressed columns. The signals may also be used to approximate harmonic and inharmonic signals and modulated signals.

The application of properly generalised exponential and trigonometric function makes it possible to work numerically stable or unstable procedures. A continuous–discrete engineering problem transformation causes a hazard of its quantitative and qualitative re-classification. In the course of continuous–discrete re-modelling of an engineering system, the initial unique solution turns into ambiguous solution or no solution at all. The numerical procedure user should be aware of such situations and take them into account in the design process.

Even a highly selected constitutive model is not enough for the considered boundary problem to relevantly reflect the real analysed case. The static or dynamic load feature possibly complements the differential equation by inertial effects, effects of elastic and shock wave propagation, loading rate and so on. Numerical modelling is also highly influenced by the specific material instability or imperfections, that is, density and temperature variations, microcracks showcasing ductile fracture patterns out of yield zones.

All the enlisted issues consequently increase the algorithm sophistication and computational effort [6]. On the other hand, they trigger more complex mathematical forms, making the errors in numerical procedures more difficult to detect.

The design quality progress is chiefly an effect of progress in numerical domain. However, the engineers should stay critical to numerical results of the ‘black box’ software, excluding the algorithm tracing.