Fractal mechanics

Physica D 220 (2006) 54–68 www.elsevier.com/locate/physd

Fractal mechanics b ´ Marcelo Epstein a,∗ , J˛edrzej Sniatycki a Department of Mechanical and Manufacturing Engineering, University of Calgary, Calgary, Alberta T2N 1N4, Canada b Department of Mathematics and Statistics, University of Calgary, Calgary, Alberta T2N 1N4, Canada

Received 30 July 2005; received in revised form 19 June 2006; accepted 26 June 2006 Available online 1 August 2006 Communicated by S. Kai

Abstract A mechanical theory of fractals and of non-smooth objects in general is developed on the basis of the theory of differential spaces of Sikorski. Once the (generally infinite dimensional) configuration space is identified, an extended form of the principle of virtual work is used to define the concept of generalized force and stress. For the case of self-similar fractals, an appropriate integration based on the Hausdorff measure is introduced and applied to the numerical formulation of stiffness matrices of some common fractals, which can be used in a finite element implementation. c 2006 Elsevier B.V. All rights reserved.

Keywords: Solid mechanics; Differential geometry; Differential spaces; Fractals; Finite elements

1. Introduction Considered not just as geometrical entities but also as collections of real-life material particles, can fractals possess elasticity or display any other signs of material response? And, if so, how can one go about formulating a consistent mechanical theory of fractals? These questions, and others of this kind, can and should be raised and answered, since the (often informal) use of fractals has become a powerful technique for modelling a large variety of natural phenomena: from snowflakes to crack formation in rocks, from plant growth to coastal erosion, from the mechanics of porous media to the innervation of muscle tissue. Our aim in this paper is twofold. On the one hand, we hope to achieve a formulation of Mechanics that includes both Classical Particle Mechanics and the Mechanics of Continuous Media as particular cases of a single theory. On the other hand, and more importantly, we expect to develop a mechanical theory of deformable entities which are neither discrete sets of points nor differentiable manifolds. Questions such as “What is the stiffness of a Koch beam?”, and “How does a Sierpi´nski carpet deform under stress?” make sense within this new framework. Foundational work in this field has been carried out, to the best of our knowledge, either by means of fractional calculus [3,4] or through deeper techniques of functional analysis [2,8]. Our proposal herein is to use the geometrical theory of differential spaces [10–12] in order to characterize the material body and its configurations. Among the potential advantages of this approach we may mention the following: (1) a differential space is a geometric object whose mathematical treatment closely resembles that of a differentiable manifold, a terrain where the techniques of traditional Continuum Mechanics are well established; (2) the notion of differential space is not limited to fractals, but is applicable also to any subset of Rn ; (3) the groundwork laid down in works such as [3,4,2,8], when reinterpreted within the framework of differential spaces, is likely to acquire a renewed geometrical meaning. The present work is by no means definitive. On the contrary, it should be regarded as the outline of a plan for future work. On the other hand, even at this initial level of treatment, we obtain specific results of great practical utility, such as exact stiffness matrices ∗ Corresponding author. Tel.: +1 403 220 5791; fax: +1 403 282 8406.

´ E-mail addresses: [email protected] (M. Epstein), [email protected] (J. Sniatycki). c 2006 Elsevier B.V. All rights reserved. 0167-2789/$ - see front matter doi:10.1016/j.physd.2006.06.008

´ M. Epstein, J. Sniatycki / Physica D 220 (2006) 54–68

55

to be used in the Finite Element analysis of fractals. After a brief introduction to fractals and differential spaces (Section 2) and to a mechanical theory based on the Principle of Virtual Work (Section 3), specific results are developed for the mechanics of fractals (Section 4), including Finite Element evaluations for the Sierpi´nski gasket and carpet. A whole section (Section 5) is devoted to the bending of a Koch beam, a topic that raises interesting and challenging questions as to the appropriate differential structure to assign to a given material body as part of its very definition. 2. Material bodies 2.1. General ideas Our point of departure is a material body B, which we conceive as an arbitrary subset of R3 . In Classical Mechanics of systems of point particles, the body B is a discrete subset {x1 , . . . , xn } of R3 . In Continuum Mechanics, B is a smooth submanifold of R3 (possibly with boundary). These two examples can be seen as extreme cases, and their detailed study by the respective disciplines has been amply justified by the wide range of their scientific and technological applications. It is clear, however, that these very applications are in need of a mechanical theory of rather more general bodies. The perfect example of this need is provided by the increasing popularity of the use of fractals in modelling natural phenomena. A precise definition of the term fractal does not exist (cf. [6]), nor is it necessary for our purposes. We will treat the geometry of fractals mainly as a rich source of examples of non-conventional material bodies, all of which can be treated by means of the proposed mechanical theory. 2.2. Fractals Some of the most common and useful examples of fractals are: the Koch curve, the Sierpi´nski gasket and carpet, and the Menger sponge. All of these fractals enjoy the property of self-similarity. Roughly speaking, a subset of Rn is said to be self-similar if it is a union of a number of smaller similar copies of itself. We will limit ourselves, therefore, to briefly review some of the essential notions of the theory of fractals as it applies to self-similar sets. A previous acquaintance with at least one of the sets mentioned above is useful but not necessary. Given an arbitrary subset B of Rn , a cover of B is a family U of sets Uα ⊂ Rn such that: [ B⊂ Uα . (2.1) α

If the family U is countable (or finite) the cover is said to be countable (or finite). It is customary in that case to indicate the members of the family with a Latin subscript, namely, Ui , where i ranges over the natural numbers. The diameter of a subset of Rn is the least upper bound (i.e., the supremum) of the distance between pairs of points in the subset. By convention, the diameter of the empty set is zero. A δ-cover of B ⊂ Rn is a (countable) cover such that, for every i, diam(Ui ) ≤ δ, where δ is a positive real number and where diam(.) is the diameter function on subsets of Rn . The members of a δ-cover are indicated by Uiδ . The s-dimensional Hausdorff measure Hs (B) of B is defined as: Hs (B) = lim inf δ→0

∞ X

[diam(Uiδ )]s ,

(2.2)

i=1

where the inf extends over all δ-covers. It can be shown that this definition indeed provides an outer measure for all non-negative values of s and for any subset of Rn . The Hausdorff dimension dim H B of B is defined as: dim B = inf{s ≥ 0 : Hs (B) = 0}. H

(2.3)

The Hausdorff dimension of a subset of Rn cannot exceed n. If dim H B > 0, then it can be shown that for all values of s strictly smaller (respectively, larger) than the dimension, the s-dimensional Hausdorff measure of B is infinite (respectively, zero). In this sense, the Hausdorff dimension represents a critical value of discontinuity of the s-dimensional Hausdorff measure, thus making Eq. (2.3) meaningful. Of particular interest is the value of the s-dimensional Hausdorff measure for s equal to the Hausdorff dimension of the set. We shall call this the Hausdorff measure of the set. The Hausdorff measure of a set may turn out to have any value, including zero or infinite. Remark 2.1. The Hausdorff dimension reduces to the ordinary one when the set under consideration is a manifold. Moreover, for integer n, the n-dimensional Hausdorff measure is a constant multiple of the n-dimensional Lebesgue measure (i.e., the ordinary notion of n-dimensional volume). Remark 2.2. Although the notion of s-dimensional Hausdorff measure is not difficult to grasp, its evaluation is in general very difficult. Conceptually (and, perhaps, numerically), what needs to be done is the following: For some δ > 0 one constructs a δ-cover. The (finite or infinite) sum of the s-powers of the diameters of this cover is then evaluated and recorded. The same is to

´ M. Epstein, J. Sniatycki / Physica D 220 (2006) 54–68

56

be done for all possible δ-covers, for the same value of δ. In this way, an infinite set of real numbers is obtained, whose largest lower bound (i.e., its infimum) is calculated. By construction, this infimum is never negative, but it may certainly be infinite. One is to repeat this process for smaller and smaller values of δ so as to evaluate the limit of these infima as δ tends to zero. The result (which may be infinite) is the desired s-Hausdorff measure. Notice that as δ decreases, the available covers become scarcer, so that the infimum necessarily increases and the limit is approached from below. The determination of the Hausdorff measure of fractals is an active area of research. For the unit-sided Sierpi´nski triangle, for example, it is only known [13,7] that the Hausdorff measure lies between 0.5 and 0.81. A similar lack of knowledge applies to the Koch curve [14]. An important feature of the s-dimensional Hausdorff measure is the scaling property, namely, the way it changes under similarity transformations of Rn . Recall that a transformation S : Rn → Rn is called a similarity if there exists a real scale factor λ > 0 such that, for all x, y ∈ Rn , the following equation is satisfied: |S(y) − S(x)| = λ|y − x|,

(2.4) Rn .

Rn

where |.| denotes the Euclidean distance in It is not difficult to prove that for all subsets B ⊂ and for all values of s, under a similarity transformation S with scale factor λ the s-dimensional Hausdorff measure transforms according to the formula: Hs (S(B)) = λs Hs (B).

(2.5)

Another important property of the s-dimensional Hausdorff measure is that it is preserved under Euclidean isometries (translations, rotations, reflections). In other words, congruent sets have the same s-dimensional Hausdorff measure. The behaviour of the Hausdorff measure under general affine transformations, on the other hand, cannot be captured under the umbrella of a simple formula. The properties just described of the s-dimensional Hausdorff measure can be used to obtain a straightforward evaluation of the Hausdorff dimension of self-similar fractals. Indeed, let B be the union of m copies of λ-scaled mutually congruent copies of itself (with λ < 1). If these copies are disjoint Borel sets, we clearly have: Hs (B) = mλs H(B),

(2.6)

by virtue of the scaling property (2.5). Assume now that the Hausdorff measure of B (namely, the value of Hs (B) for s = dim H (B)) is finite and positive. In that case, it follows form (2.6) that: dim(B) = − log(m)/ log(λ). H

(2.7)

There are two problems with this simple “derivation”. The first one is the assumption that the Hausdorff measure of the fractal is finite and positive. The second difficulty arises from the fact that, in many examples of interest, the self-similar fractal is not quite the union of strictly disjoint sets. The first difficulty is eliminated by the fact that one can actually prove that the Hausdorff measure of a self-similar set is necessarily finite and positive. As far as the second problem is concerned, the additivity property of a measure usually extends beyond the case of disjoint Borel sets. Intuitively, the sets under consideration can intersect, but not too much. More technically, the additivity will hold if the sets satisfy a so-called open set condition (see, e.g., [6]). This condition is satisfied by the Koch curve and by the Sierpi´nski and Menger fractals mentioned above. Example 2.3 (Koch Curve). Starting from a (horizontal) line segment (of, say, unit length), the middle third is replaced by two pieces of equal length, to form a zigzag made of 4 pieces of length 1/3. Each of these pieces is subjected to the same procedure, and so on ad infinitum. To apply our formula (2.7), we note that, by construction, m = 4 and λ = 1/3, with the result that dim H B = log(4)/ log(3). In the following we shall need an observation that the resulting Koch curve B is the topological closure of the set of all vertices of all zigzags. Moreover, for each vertex p ∈ B, there are two sequences of vertices, pn and pn0 , in B, which converge √ √to p and are contained in two distinct lines through p. If p is an interior √ point of √B, we can choose lines with slopes 0, 3 or − 3. For the left (right) end point, we can choose lines with slopes 0 and 3/2 (− 3/2). Example 2.4 (Sierpi´nski Carpet). The Sierpi´nski carpet is generated from an initial square (with unit side, say) from which the central square (of side 1/3) is removed. The resulting figure can be seen as the union of 8 equal squares (of side 1/3), any two of which intersect at most along a side. The fractal is generated by repeating this process for each of these 8 squares, and so on ad infinitum. We have now that m = 8 and λ = 1/3 yielding: dim H (B) = log(8)/ log(3). In the following we shall need an observation that the resulting Sierpi´nski carpet B is the topological closure of the union of all sides of all squares appearing in the process. Moreover, for each point p ∈ B, contained in one of the sides, there are two sequences of vertices xn and yn in B, which converge to p with all xn contained in the horizontal line through p and all yn contained in the vertical line through p. Example 2.5 (Sierpi´nski Gasket). The Sierpi´nski gasket (or triangle) is generated by removing an inverted central (1/2)-sided equilateral triangle from an initial unit-sided equilateral triangle. Its Hausdorff dimension is log(3)/ log(2). As before, the Sierpi´nski gasket B is the closure of all sides of all the triangles appearing in the process, Moreover, for each point p ∈ B contained in one of

´ M. Epstein, J. Sniatycki / Physica D 220 (2006) 54–68

57

the sides, there√are two sequences of vertices xn and yn in B, which converge to p and are contained in two distinct lines through p √ with slopes 0, 3, or − 3. Fractals are geometrical, not physical, objects. But so are differentiable manifolds. When using such geometric concepts to model physical entities, care should be exercised in terms of the physical interpretation. Thus, for instance, a “real” piece of matter does not strictly have a volume or an area. We talk, nevertheless, of forces per unit volume, pressure, mass density, etc., and we have become so used to these abuses that we no longer question them. Analogously, if we wish to model a body by means of a fractal, we should be willing to accept a new conceptual framework. Since the dimension of a fractal is not an integer, its Hausdorff measure will be expressed in terms of units of length raised to a fractional power. We should be prepared, therefore, to accept that the material properties of such an entity be expressed and understood in similar terms. Thus, for instance, if the fractal is made of a linearly elastic material, we should not be surprised if its modulus of elasticity is expressed in units of force divided by a fractional power of units of length. The physical context of each application will be the guide to interpret this (only apparently) curious fact accordingly. 2.3. Differential spaces Fractals are not the only possible, or even interesting, generalization of the concept of material body. Manifolds with corners and unions of manifolds of different dimensions (including zero) are two examples that immediately come to mind. In fact, there is no reason to exclude a priori any subset of R3 as a potential model for a material entity on which to formulate a mechanical theory.1 In order to formulate a mechanical theory that is valid for bodies that are arbitrary subsets of R3 , we will only require that the notion of smooth function over such sets be available and well defined. This is the domain of the theory of differential spaces [10–12], which we now briefly review. A differential space consists of a topological space B and a collection, denoted as C ∞ (B), of real-valued continuous functions. That is, out of the class C 0 (B) of all the continuous functions f : B → R, a particular subclass C ∞ (B) is singled out to be considered as deserving the name of “smooth”. This subclass, however, cannot be arbitrarily chosen, but must satisfy the following three conditions: (i) it must be closed under composition with smooth functions in Rn ; (ii) it must be consistent with localized restrictions; (iii) the family of sets { f −1 ((a, b))| f ∈ C ∞ (B), a, b ∈ R} must generate the topology of B. More specifically, let f : B → R be a smooth function on B and let φ : R → R be a smooth function (in the ordinary sense) on R. Then the composition φ ◦ f must be a smooth function on B. In other words: f ∈ C ∞ (B), φ ∈ C ∞ (R) ⇒ φ ◦ f ∈ C ∞ (B). Similarly, if for any n and any smooth functions f 1 , . . . , f n we form the function F : B → the composition of F with any smooth function Φ : Rn → R belongs to C ∞ (B):

(2.8) Rn

given by x 7→ ( f 1 (x), . . . , f n (x)),

f i ∈ C ∞ (B) (i = 1, . . . , n), Φ ∈ C ∞ (Rn ) ⇒ Φ ◦ ( f 1 , . . . , f n ) ∈ C ∞ (B).

(2.9)

As far as condition (ii) is concerned, it expresses the fact that if a continuous function g : B → R has the property that for each x ∈ B there exists a neighbourhood U and a smooth function f ∈ C ∞ (B) such that g|U = f |U , then g must be in C ∞ (B). In fact, given any initial choice of functions generating the topology of B, these two properties (being closed under composition and being consistent with restrictions) allow us to complete a class of smooth functions on B by adding all other functions that fulfil these conditions. In this way, the set B is said to have acquired a differential structure and becomes a differential space. This structure enables us to use calculus to study the geometry of B even though B itself need not be smooth in the usual sense. Let B1 and B2 be differential spaces and let κ be a map from B1 to B2 . The map κ is said to be smooth if for each smooth function f ∈ C ∞ (B2 ), the pull-back (to B1 ) of f by κ is smooth. In other words, κ is smooth if the following implication holds: f ∈ C ∞ (B2 ) ⇒ κ ∗ f = f ◦ κ ∈ C ∞ (B1 ).

(2.10)

The map κ : B1 −→ B2 is a diffeomorphism if it is smooth with a smooth inverse κ −1 : B2 −→ B1 . Example 2.6. It is not difficult to verify that a differentiable manifold together with the class of smooth functions is a differential space. More striking, however, is the fact that any subset B of Rn can be endowed with a natural differential structure by adopting the topology induced by the metric topology of Rn and defining the smooth functions C ∞ (B) to consist of all functions on B that locally extend to smooth functions on Rn . In other words, f ∈ C ∞ (B) if, for each x ∈ B ⊂ Rn , there exists a neighbourhood U of x in B and a smooth function F on Rn such that f restricted to U coincides with the restriction of F to U . Loosely speaking, we say that in this case the smooth functions in C ∞ (B) are defined as restrictions to B of smooth functions on Rn . We will call this the induced differential structure on B. 1 Even this restriction to subsets of R3 may be too stringent, as the various theories of media with internal structure, whereby the material body is identified with some fibre bundle over a standard body-manifold, clearly demonstrates.

´ M. Epstein, J. Sniatycki / Physica D 220 (2006) 54–68

58

Using functions Φ : R2 → R such as (x, y) 7→ x + y and (x, y) 7→ x y, and invoking property (2.9), one can see that C ∞ (B) is an (associative and commutative) algebra (with pointwise sum and product of functions). A derivation at p ∈ B is a linear map v : C ∞ (B) → R satisfying the product (or Leibnitz) rule: v( f g) = v( f )g( p) + f ( p)v(g),

∀ f, g ∈ C ∞ (B).

(2.11)

The collection of all derivations at p forms a vector space (under the obvious operations) which, by analogy with the case of differentiable manifolds, is called the tangent space to B at p, denoted by T p B. Note that, while in the case of differentiable manifolds the tangent spaces at all points have the same dimension, in the case of differential spaces the dimension of T p B may depend on p. The set of derivations is never empty, since the zero operator is linear and satisfies the Leibnitz rule. The collection of all tangent vectors at all points of B is denoted T B and called the tangent bundle space of B (in the literature the term tangent pseudobundle is also used). Example 2.7. Consider the case in which B is an “angle”, namely, the union of two (closed) non-collinear line segments with one common end point (the corner). Seen as a subset of R3 and with the induced differential structure, this is a well-defined differential space. A smooth function f on B is the restriction to B of some smooth function on R3 . On the other hand, given f ∈ C ∞ (B) there exists an infinite number of extensions to smooth functions on R3 , all of which have, point by point, the same directional derivative in the direction of the leg to which the point belongs. At the corner, all these extensions will automatically share 2 directional derivatives. These directional derivatives constitute a basis for all possible derivations at the point. The dimension of T p (B) is, therefore, 1 at all points p (including the ends of the angle) except at the corner, where the dimension is 2. Proposition 2.8. Let B be the Koch curve, the Sierpi´nski carpet or the Sierpi´nski gasket with the differential structure induced by its embedding in R2 . Then, for each p ∈ B, dim T p B = 2. Proof. By construction, B is a closed subset of R2 . Hence, C ∞ (B) consists of restrictions to B of smooth functions on R2 . Consider first the case when B is the Koch curve, and let p be a vertex in B. For every function f ∈ C ∞ (R2 ), its gradient at p is uniquely determined by the values of f at the points of the sequences pn and pn0 introduced in Example 2.3. Since these sequences are contained in different lines through p, it follows that there are two linearly independent derivations of C ∞ (B) at p. Hence, dim T p B = 2. Suppose now that p ∈ B is a limit point of the set of vertices √ in B. Let qn be a sequence of vertices in B convergent to p. For every f ∈ C ∞ (R2 ), and each line through p with slope 0 or ± 3, the partial derivative of f at p in the direction of this line is the limit, as n → ∞, of partial derivatives of f at qn in the directions of the lines through qn with the same slope. Hence, the space of derivations of C ∞ (B) at p is spanned by two linearly independent derivations. This implies that dim T p B = 2. The argument for the Sierpi´nski carpet and gasket is similar.  A map κ from B to R3 is a smooth embedding if it is an embedding2 and if, for every φ ∈ C ∞ (R3 ), the composition φ ◦ κ is in Note that, under a smooth embedding κ, B is diffeomorphic to its image κ(B). A smooth embedding κ of B into R3 gives rise to a map T κ from T B to T R3 . For every x ∈ B and v ∈ Tx B, T κ(v) is a vector tangent to R3 at κ(x) such that, for every function φ ∈ C ∞ (R3 ),

C ∞ (B).

T κ(v)(φ) = v(φ ◦ κ).

(2.12)

Remark 2.9. For the case of interest, in which B is a compact subset of R3 endowed with the induced differential structure, it is sometimes intuitively useful to think of a smooth embedding κ : B → R3 as simply the restriction to B of a diffeomorphism Φ of R3 . Similarly, the tangent map T κ at p ∈ B can be thought of as the restriction of T Φ to the subspace T p B. Remark 2.10. For p ∈ B the dimension of the tangent space T p B is invariant under diffeomorphisms. In other words, if κ : B1 −→ B2 is a diffeomorphisms between differential spaces, then necessarily: dim(T p B1 ) = dim(Tκ( p) B2 )

∀ p ∈ B1 .

(2.13)

This result follows directly form the definitions. 3. Mechanics The most fundamental kinematic notion is that of configuration space, namely, the collection of all possible configurations of a mechanical system. In Classical Mechanics of a system of particles, a configuration is a 1–1 map from B = {x1 , . . . , xn } to R3 . 2 Recall that a map κ : T −→ T between the topological spaces T and T is an embedding (or, more precisely, a topological embedding) if it is a 1 2 1 2 homeomorphism of T1 onto its image κ(T1 ) ⊂ T2 in the subspace topology.

´ M. Epstein, J. Sniatycki / Physica D 220 (2006) 54–68

59

Under the assumption of smoothness (or absence) of constraints among particles, the configuration space of Classical Mechanics is a smooth m-dimensional manifold Q, where m is the number of degrees of freedom. In Continuum Mechanics, on the other hand, a configuration is a smooth embedding of B into R3 . If B is a compact manifold with boundary (say, a closed ball in R3 ), the configuration space Q is an infinite dimensional manifold of maps. If B is not compact, then there are difficulties with defining a manifold structure on Q. Nevertheless, even in this case it is convenient to use the notation and terminology of the theory of manifolds of maps. As we have shown above, the notion of smooth embedding of B into R3 makes sense for arbitrary subsets of R3 . Accordingly, we define a configuration of any given body B ⊂ R3 as a smooth embedding κ : B → R3 . The space Q of all smooth embeddings of B into R3 (satisfying, if necessary, any given kinematic constraints) is the configuration space of our system. We do not know if Q is a manifold of maps. Nevertheless, following the tradition of continuum mechanics for non-compact bodies, we use the manifold notation and terminology. Given a configuration κ, the derived map T κ : T B → T R3 is called the deformation gradient. Virtual displacements are tangent vectors to Q. A virtual displacement of a configuration κ ∈ Q is a vector δκ ∈ Tκ Q, which can be visualized as vector field in R3 defined on the range of κ. Remark 3.1. Intuitively, we may think of a tangent vector δκ to a configuration κ ∈ Q as an equivalence class of “curves” in Q through κ. Somewhat more precisely, the construction of a tangent vector is the following: We start from a smooth one-parameter (β, say) family of configurations of B, such that the given κ corresponds to β = 0. As β varies, each point of B describes an ordinary smooth curve in R3 with a well-defined tangent vector at β = 0. It is clear that the tangent vector δκ corresponding to this one-parameter family can be identified with the collection of all these “ordinary” tangent vectors, that is, with a vector field on κ(B). The collection of all tangent vectors δκ at the configuration κ constitutes the (finite or infinite dimensional) tangent space Tκ Q. In Newtonian Mechanics, the concept of force is a primary notion. Using the underlying metric structure of Euclidean space, forces are defined as vectors at the same level as displacements, velocities and accelerations. The point of view of Analytical Mechanics, on the other hand, considers forces (or “generalized forces”) to be linear operators on virtual displacements. In other words, forces are elements of the cotangent bundle of the configuration space. The evaluation of a co-vector (force) on a vector (virtual displacement) is called virtual work. No metric structure is, therefore, invoked to define forces. The Newtonian (vectorial) conditions of equilibrium of a system are replaced with a single scalar identity stating the vanishing of the virtual work for all virtual displacements. This approach, which assigns a primary status to the kinematics of a system, is not in standard use in Continuum Mechanics, which, with few exceptions (e.g. [5,9]), tends to follow the Newtonian point of view. Nevertheless, its usefulness should be evident whenever we are confronted with a non-standard mechanical system for which the kinematical setting is well defined. As we have shown, this is the case of the mechanics of differential spaces and, in particular, of fractals. Following the lead of Analytical Mechanics, we define a generalized force f at a configuration κ ∈ Q as a bounded linear operator on Tκ Q. The evaluation h f, δκi of a generalized force f at a configuration κ on a virtual displacement δκ at the same configuration is called the virtual work of f on δκ. Remark 3.2. This notion of generalized force is wide enough to encompass all the classical notions of external force, internal stress, non-local interactions, etc. A distinction between internal and external forces can be postulated as part of a given mechanical theory (for instance by demanding that the internal forces perform no virtual work on a class of virtual displacements). Such distinction, however, is not essential for the development of the general theory. To completely define the response of a mechanical system, the generalized forces acting on it must be specified by means of a so-called constitutive equation. In principle, the generalized forces may depend on the whole past history of the configurations of the system. Excluding all such memory effects, we will consider the case in which the generalized forces depend only on the present configuration of the system, namely, f = f (κ). From a differential geometric point of view, we observe that this constitutive equation is tantamount to choosing a particular cross-section of the cotangent bundle T ∗ Q. Remark 3.3. Under conditions of smoothness, a constitutive equation is simply a differential 1-form on Q. In other words: f : Q → T ∗ Q,

(3.14)

with the condition that π ◦ f = idQ , where π : T ∗ Q → Q is the cotangent-bundle projection and idQ : Q → Q is the identity map. A system is called conservative if this 1-form is exact, namely, if there exists a scalar function W : Q → R such that f = dW . A configuration of a system is said to be an equilibrium configuration if the equation: h f (κ), δκi = 0 is satisfied identically for all virtual displacements.

(3.15)

´ M. Epstein, J. Sniatycki / Physica D 220 (2006) 54–68

60

Remark 3.4. We are confining out attention to Statics alone. In this case, the identity (3.15) can be considered as an equation (or system of equations) to be solved for the unknown equilibrium configuration κ. The nature of this equation (algebraic, ODE, PDE, integro-differential, etc.) depends on the particular theory at hand. Example 3.5. In the framework of Continuum Mechanics, the standard theory of elasticity corresponds to the case in which the 1-form can be represented as an integral of the form: Z Z i I i ti δκ i dA, (3.16) h f, δκi = [ f i δκ − f i δκ,I ]dV + ∂B

B

where the usual coordinate notation is used and where dV and dA denote, respectively, the ordinary Cartesian volume and area element. The quantities f i , f iI and ti represent, respectively, the body force, the (Piola) stress and the surface traction. A material is of the first grade (or simple) if the stress at a point p is a function only of the local value of the deformation gradient at that point, namely: j

f iI = f iI (κ,J ( p), p).

(3.17)

Dependence on higher gradients of the deformation, particularly the second gradient, is also of practical interest. 4. Fractal mechanics 4.1. Introduction In formulating a theory for the mechanical response of fractals, one is tempted to proceed by a tedious passing to the limit, whereby at each step of the generation of the fractal the body is understood as an ordinary body (a collection of rods or of plates, etc.). We propose to show how this procedure may be replaced by a more rigorous one formulated directly on the limit entity, that is, on the fractal itself. We will confine our attention to self-similar fractals. As already pointed out, the self-similar fractals that we have considered have a constant (integer) tangent dimension. To further simplify the description, we will assume that this tangent dimension is equal to the dimension of the ambient space (R2 or R3 ). For example, for the Sierpi´nski carpet, we will limit the analysis to in-plane deformations. As a consequence of this assumption, and in accordance with the definition of tangent map, Eq. (2.12), the deformation gradient at each point of the fractal looks exactly as an ordinary deformation gradient at a point of a continuum. What this means is that a constitutive equation such as (3.16) becomes available and perhaps meaningful, provided we interpret the integration appropriately, namely, with respect to the Hausdorff measure. 4.2. Integration over a fractal domain We have already noted that the Hausdorff measure is notoriously difficult to calculate and, in fact, even for some of the most frequently used fractals the exact value of the Hausdorff measure is not known. It is only by virtue of the property of self-similarity that we can succeed in performing approximate Hausdorff integration over a fractal domain, provided that the function to be integrated is reasonably well behaved in the underlying Euclidean domain (see [1]). For example, we may assume that the function satisfies a Lipschitz condition. As we will presently see, integration over a self-similar fractal domain will always provide a result which is proportional to the (usually unknown) Hausdorff measure of the fractal. We are interested, therefore, in obtaining as good an approximation as possible of the constant of proportionality. We start by recalling the scaling property (2.5) of the Hausdorff measure. Another useful result is the following [6]: Let B ⊂ Rn and let f : B → Rn be a Lipschitz mapping, namely, there exists a constant c such that: | f (x) − f (y)| ≤ c|x − y|

∀x, y ∈ B,

(4.18)

where |.| denotes the Euclidean distance. Then Hs ( f (B)) ≤ cs Hs (B).

(4.19)

We already pointed out that the Hausdorff measure is invariant under Euclidean translations and rotations. Unfortunately, it does not seem possible to capture the behaviour of the Hausdorff measure under general affine transformations. Consider a fractal (such as the Koch curve or the Sierpi´nski carpet) obtained from a so-called initiator (zero-th iteration) and a generator (first iteration). Let the ordinary (Euclidean) measure of the initiator be a. Let the generator consist of m congruent scaled-down versions of the generator, each of ordinary measure λa (λ < 1 being the scaling factor). Example 4.1. For the Koch curve the initiator is a segment of some length a, and the generator consists of 4 segments, each of length a/3. For the Sierpi´nski gasket the initiator is an equilateral triangle, while the generator consists of 3 equal equilateral triangles scaled down from the initiator by a factor of 1/2.

´ M. Epstein, J. Sniatycki / Physica D 220 (2006) 54–68

61

Because of the self-similarity of the generating process, we can subdivide the fractal B into m J mutually congruent pieces Bi (i = 1, . . . , m J ), where J is an arbitrary positive integer. Each of these pieces has a Hausdorff measure equal to the Hausdorff measure of the fractal multiplied by the factor λs J . Let H0 denote the Hausdorff measure of the fractal. Recall that this measure is finite and positive. Let φ : B → R be a “well-behaved” (for example, Lipschitz) function, and denote by φi the value of φ at some point within the piece Bi . Then, by definition, the following expression gives the value of the Hausdorff integral of φ over the fractal domain B: Z B

φdH(B) = lim

J →∞

mJ X

φi λs J H0

(4.20)

i=1

where s is the Hausdorff dimension. We omit the proof of the fact that, when this limit exists, it is independent of the choice of the points within the individual pieces. As a practical demonstration of the use of this definition to calculate Hausdorff integrals over a fractal set, we list a simple computer code written in BASIC to effect this task on a Sierpi´nski carpet of sides of length a aligned with Cartesian coordinates x, y whose origin is at the centre of the carpet. The method is based on the so-called ternary representation technique. The original domain is subdivided chessboard-like into squares of side a/3 J and the algorithm determines which squares actually belong to the J -th generation step by representing its (integer) coordinates in base 3 and choosing those squares for which it never happens that they have a digit 1 in the same position. The function φ is then evaluated at the centres of the surviving squares and Eq. (4.20) is applied. Function integrate(a, j) Sum = 0 For k = 1 To 3ˆj For l = 1 To 3ˆj kx = k − 1 lx = l − 1 For m = 1 To j res k = kx − Int(kx/3) ∗ 3 kx = Int(kx/3) res l = lx − Int(lx/3) ∗ 3 lx = Int(lx/3) If (res k= 1 And res l= 1) Then GoTo 100 Next m x = −a/2 + (k − 0.5) ∗ a/3ˆj y = −a/2 + (l − 0.5) ∗ a/3ˆj Sum = Sum + phi(x, y) ∗ (1/3)ˆ ( j ∗ Log(8)/Log(3)) 100 Next l Next k integrate = Sum End Function Notice that, for clarity, we have not written the code in a most efficient way (for example, we could have replaced (1/3)ˆ (j ∗ Log(8)/Log(3)) with (1/8)ˆj, and we could have avoided repeated calculations of residuals within the loops). Be that as it may, using this simple code with a = 1 and φ(x, y) = x 2 + y 2 , we obtained (as an example) a value of 0.1875H0 . Convergence to 3 significant digits was achieved already at the third generating step. For later use, Table 1 provides a few “moments” of the Sierpi´nski carpet over the unit square, where the moment of order i, j is defined as: Z Ii j = x i y j dH(B). (4.21) B

4.3. Reduced configuration spaces and the Rayleigh–Ritz method So far we have considered a configuration space Q consisting of all the possible smooth embeddings of the body B into the physical space. There are several reasons that justify the consideration of a mechanical theory in other configuration spaces as well. One of these reasons could be that the space of smooth embeddings may turn out to be too restrictive for certain applications. Among the many instances where this might be the case, we can mention that of layered bodies containing sharp surfaces of discontinuity between adjacent layers. For such bodies, C 0 -maps may be required. Another example is provided by the propagation of wave fronts.

´ M. Epstein, J. Sniatycki / Physica D 220 (2006) 54–68

62 Table 1 Moments i

j

Ii j /H0

0 0 2 2 0 4

0 2 0 2 4 0

1 0.09375 0.09375 0.008008 0.014648 0.014648

A completely different set of reasons may tip the balance to the side of restricting, rather than enlarging, the configuration space. One of these reasons is the presence of external or internal constrains, such as incompressibility or inextensibility. Another reason is not from physical but from numerical origin, and this is the motivation behind the treatment of this section. One of the many approximation techniques for the solution of problems in mechanics goes back to an idea of Lord Rayleigh, perfected later by his disciple Ritz and further enriched by Galerkin, culminating in the modern theory of finite elements. According to the main principle underlying this technique, the infinite dimensional configuration space of a body is curtailed to a finite dimensional subspace, controlled by a finite number of generalized coordinates. Within this reduced subspace, the principle of virtual work can be enforced to obtain the “best” approximation to the actual solution. When properly formulated, the method can be applied systematically to increasingly larger subspaces and shown to lead to convergence (in some norm) to the solution of the original problem. We will exemplify the application of the Rayleigh–Ritz method to the approximate solution of the problem of the uniaxial extension of a Sierpi´nski carpet of side a. Let the Cartesian coordinate system X, Y have its origin at the centre of the carpet, such that the corners have coordinates whose absolute values are a/2. The two opposite edges X = ±a/2 are subjected to a uniformly distributed tension q, measured in units of force per unit length. Identifying the reference space and its coordinate system with their physical counterpart, we introduce the displacement vector with components u = x − X and v = y − Y , where x and y are the spatial coordinates of the point with referential coordinates X and Y , respectively. If we assume these displacements and their gradients to be very small, the so-called infinitesimal theory can be used. The strain measures are defined as: X =

∂u , ∂X

Y =

∂v , ∂Y

γX Y =

∂u ∂v + , ∂Y ∂X

(4.22)

where γ denotes the shear strain, while the ’s denote the longitudinal strains. As an example of a restriction of the configuration space of this body, we adopt for each of the displacement fields a truncated double power-series expression, namely: u(X, Y ) =

M,N X

Ci j X i Y j ,

v(X, Y ) =

i, j=0

M,N X

Di j X i Y j ,

(4.23)

i, j=0

where Ci j and Di j are constants to be determined. By the symmetry of our problem, for the coefficients Ci j we need only consider odd i’s and even j’s, whereas the reverse is true for the coefficients Di j . It is convenient to introduce the non-dimensional coordinates ξ = X/a and η = Y /a and to rewrite (4.23) as: u(X, Y ) =

M,N X

ci j ξ i η j ,

v(X, Y ) =

i, j=0

M,N X

di j ξ i η j ,

(4.24)

i, j=0

with an obvious notation. We obtain: 1X X = ci j iξ i−1 η j a i, j Y =

1X di j i jξ i η j−1 a i, j

γX Y =

 1 X ci j jξ i η j−1 + di j iξ i−1 η j . a i, j

(4.25)

If we assume the material to abide by Hooke’s law of isotropic linear elasticity in a plane-stress regime, we can write the stress components as: σX =

E ( X + νY ) 1 − ν2

´ M. Epstein, J. Sniatycki / Physica D 220 (2006) 54–68

E (Y + ν X ) 1 − ν2 E γX Y , = 2(1 + ν)

63

σY = τX Y

(4.26)

where E and ν are Young’s modulus and Poisson’s ratio, respectively. According to our philosophy, the expression for the virtual work of the stresses is borrowed from the ordinary theory of continuous media as: Z IVW = h (4.27) (σ X δ X + σY δY + τ X Y δγ X Y ) dH(B), B

where h is the known thickness of the carpet. Note that, in order for this expression to be consistent in terms of physical units, Young’s modulus has to be interpreted in units of force divided by units of length to the power s, the Hausdorff dimension of the fractal (recall that in the particular case of the Sierpi´nski carpet, s = log(8)/ log(3)). The variations of the strain fields are to be taken from Eq. (4.26) by varying the coefficients c and d. As far as the virtual work of the applied external forces is concerned, it can be calculated as: Z a/2 qδu(a/2, Y )dY, (4.28) EV W = 2 −a/2

where the displacement variations can be obtained by varying the coefficients in (4.24). Equating (identically for all independent variations of the coefficients) the last two expressions, we obtain a system of linear equations for the unknown coefficients. We illustrate such a calculation for the case in which the series for the displacements are truncated at the total power of 3, namely: u = c10 + c12 ξ η2 + c30 ξ 3 v = d01 + d21 ξ 2 η + d03 η3 . Carrying out the calculations just explained, we obtain the following system of equations:      c10  I00 I02 3I20 ν I00 ν I20 3ν I02  1                 c12  I04 + 2(1 − ν)I22 3I22 ν I02 (2 − ν)I22 3ν I04   1/12     I02                   3I20  3I22 9I40 3ν I20 3ν I40 9ν I22  c30 1/4  qa 3−s (1 − ν 2 )  = .   d01   ν I00  Eh ν I02 3ν I20 I00 I20 3I02  0                     d21  (2 − ν)I22 3ν I40 I20 I40 + 2(1 − ν)I22 3I22   0       ν I20             d03 3ν I02 3ν I04 9ν I22 3I02 3I22 9I04 0 Using the values in Table 1 for the various moments, the solution of this system for ν = 0.3 is obtained as:     c10   1.47241       c12      −1.61966           c30 −0.788152 qa 3−s (1 − ν 2 ) H0 , = d01  −0.505637   Eh              d21     1.31088   d03 0.188697

(4.29)

(4.30)

(4.31)

where H0 is still the Hausdorff measure of the unit-sided carpet. Figs. 1–3 display these results. For graphical purposes, we have arbitrarily assumed that the multiplier on the right-hand side of Eq. (4.31) is equal to 0.2. It is not our purpose to propound this kind of approximation, since the convergence of this type of scheme may be very slow and erratic, and the addition of degrees of freedom (in this case, by adding higher powers of the coordinates) somewhat awkward to implement. The Finite Element Method addresses precisely these issues by substituting for higher-order approximations a domain subdivision into smaller and smaller pieces, within each of which a coarse approximation of the type just exhibited is assumed. In other words, the configuration space is severely curtailed within each element, but expanded in the sense that the global fields are no longer required to be smooth at the inter-element boundaries. Remark 4.2. When used in conventional Continuum Mechanics, the Rayleigh–Ritz method just illustrated leads usually to a violation of the boundary conditions of stress. In the case of a fractal, it is not even clear what such “boundary” conditions are. In fact, in a certain sense, the Sierpi´nski carpet is nothing but an infinite union of boundaries of squares. One of the merits of the weak formulation of Mechanics is precisely the possibility of obtaining approximate solutions in this way, without a need to specify what exactly a Neumann-type boundary condition is.

64

´ M. Epstein, J. Sniatycki / Physica D 220 (2006) 54–68

Fig. 1. u-displacement of the uniaxially loaded Sierpi´nski carpet.

Fig. 2. v-displacement of the uniaxially loaded Sierpi´nski carpet.

Fig. 3. Deformed shape of the uniaxially loaded Sierpi´nski carpet.

4.4. Sierpi´nski finite elements The simplest possible finite element for plane-stress analysis in Continuum Mechanics is the constant-strain triangle. Its name derives form the fact that the interpolation functions (or shape functions) for the displacement fields within the element are assumed to be linear, so that the deformation gradient is constant. The element is thus endowed with 6 degrees of freedom (the displacements of the vertices, or nodes). We propose to generalize this element to Fractal Mechanics by using the Sierpi´nski gasket (also known as the Sierpi´nski triangle) as the underlying body. We need only develop the stiffness matrix for a typical gasket of, say, side b (and thickness t) with its base lying on the X -axis, symmetrically about the origin of coordinates. Rotated versions of this stiffness matrix

´ M. Epstein, J. Sniatycki / Physica D 220 (2006) 54–68

65

can be obtained by standard means. The calculation for arbitrary triangles can be conducted in an analogous fashion. Unfortunately, though, since we do not have a formula for the transformation of the Hausdorff measure under general affine transformations, we have no exact means to relate √ the stiffness matrices of non-similar triangles. For the equilateral triangle described above, we number the nodes (b/2, 0), (0, 3b/2), (−b/2, 0) as 1, 2, 3, respectively, and order the degrees of freedom correspondingly as: u 1 , v1 , u 2 , v2 , u 3 , v3 . Defining the following linear shape functions: √ 1 3 1 X− Y+ , b 3b 2 √ 2 3 Y, N2 (X, Y ) = 3b N1 (X, Y ) =

(4.32) (4.33)

and √ 3 1 1 N3 (X, Y ) = − X − Y+ , b 3b 2

(4.34)

the displacement fields are interpolated within the element as follows: u(X, Y ) = u 1 N1 (X, Y ) + u 2 N2 (X, Y ) + u 3 N3 (X, Y ),

(4.35)

v(X, Y ) = v1 N1 (X, Y ) + v2 N2 (X, Y ) + v3 N3 (X, Y ).

(4.36)

and

Assuming the same material behaviour we used for the Sierpi´nski carpet example, the stiffness matrix with respect to the above listed degrees of freedom is obtained by (Hausdorff) integration as: 7−ν   √ 6  − 3(1 + ν)  6   −1 + ν    √3 [K ] =   2 3ν   3   −5 − ν   √ 6  3(1 − 3ν) 6 

√

3(1 + ν) 6 5 − 3ν √ 6 3(1 − ν) 3 2 − 3 √ 3(3ν − 1) 6 −1 + 3ν 6

−

−1 + ν √ 3 3(1 − ν) 3 2(1 − ν) 3

√ 2 3ν 3 2 − 3 0 4 3 √

0 −1 + ν √ 3 3(−1 + ν) 3

−

2 3ν 3 2 − 3

−5 − ν √ 6 3(3ν − 1) 6 −1 + ν 3 √ 2 3ν − 3 7−ν √ 6 3(1 + ν) 6

√  3(1 − 3ν)  6  −1 + 3ν     √ 6 3(−1 + ν)    t EH(b)  3 ,   (1 − ν 2 )b2 2  −  √ 3  3(1 + ν)    6  5 − 3ν 

(4.37)

6

where H(b) is the Hausdorff measure of the Sierpi´nski gasket of side b. Recalling that the Hausdorff dimension of the gasket is log(4)/ log(3), we can write H(b) = blog(4)/ log(3) H(1). In the case of the triangular element just discussed, due to the constancy of the strain, the numerical coefficients in the stiffness matrix are identical to their continuum counterparts. It is clear, on the other hand, that this will not be the case for the bilinear Sierpi´nski carpet element, since the various integrations must be weighted by the Hausdorff measure (just as in the case of the calculation of moments in the previous subsection). Indeed, let a Sierpi´nski carpet of side 2b, with centre at the origin, be endowed with the following bilinear interpolation functions: Ni (ξ, η) =

1 (1 + ξi ξ )(1 + ηi η) 4

(i = 1, . . . , 4)

(4.38)

in terms of the non-dimensional coordinates ξ = X/b and η = Y /b. The nodes (1, 2, 3, 4) occupy the positions (ξi , ηi ), i = 1, . . . , 4, given, respectively, by: (1, 1), (−1, 1), (−1, −1) and (1, −1). Following the same procedure as for the triangular element, the element stiffness matrix is given by:   [M] [N ] t EH(2b) [K ] = [N ] [M] (1 − ν 2 )b2

(4.39)

´ M. Epstein, J. Sniatycki / Physica D 220 (2006) 54–68

66

where:  A (1 − ν)B +  16 32  (1 − ν)A ν A  +   32 16 [M] =  C (1 − ν)D  − +  16 32  ν A (1 − ν)A − 16 32

(1 − ν)A ν A + 32 16 (1 − ν)C B + 32 16 (1 − ν)A ν A − 32 16 (1 − ν)C D − + 32 16

C (1 − ν)D + 16 32 (1 − ν) A ν A − 32 16 C (1 − ν)B + 16 32 (1 − ν)A ν A − − 32 16

ν A (1 − ν)A − 16 32 D (1 − ν)F − − 16 32 ν A (1 − ν)A − 16 32 B (1 − ν)F − + 16 32

F (1 − ν)B − 16 32 ν A (1 − ν)A − 16 32 F (1 − ν)D − − 16 32 ν A (1 − ν)A + 16 32

−

ν A (1 − ν)A  −  16 32 (1 − ν)C D  − +  32 16  (1 − ν)A ν A   − −  32 16  (1 − ν)C B  + 32 16

(4.40)

ν A (1 − ν)A  + 16 32  B (1 − ν)F   − +   16 32 . ν A (1 − ν)A   +   16 32 D (1 − ν)F  − − 16 32

(4.41)

and (1 − ν)D F −  16 32  ν A (1 − ν)A  − −  16 32 [N ] =  F (1 − ν)B  −   16 32  ν A (1 − ν)A + − 16 32 

−

−

−

The constants A, B, C, D, F in these expressions are related to the moments Ii j as follows: A = I00 ,

B = I00 + I20 ,

C = I00 + I02 ,

D = I00 − I20 ,

F = I00 − I02 ,

(4.42)

where the moments have to be taken from Table 1. The symbol H(2b) stands now for the Hausdorff measure of the Sierpi´nski carpet with side equal to 2b. Notice that the modulus of elasticity E carries now different units of “area” than before. The numerical entries in the stiffness matrix are different from their counterparts in classical Continuum Mechanics in that the moments Ii j to be used are already affected by a Hausdorff integration over the fractal domain (rather than a Riemann integration over the square domain). One may therefore say, perhaps too simplistically, that the Hausdorff measure provides a means for calculating homogenized elastic constants in a fractal. Unlike the case of Continuum Mechanics, in Fractal Mechanics the subdivision into (fractal) finite elements can be conceived in two parts, the first of which carries a physical, rather than just numerical, meaning. Indeed, the body under consideration may itself consist of a tiled collection of fractals of a definite size. This first subdivision has an intrinsic physical meaning dictated by the physical situation at hand (for instance, the maximum pore size of a porous medium). The subdivision into finite elements, therefore, cannot be coarser than this initial tiling, but it can certainly be finer, always respecting the natural self-similarity of the original fractals. 5. The Koch beam The differential structure induced on a subset of Rn by restrictions of smooth functions is not the only possible structure of interest in mechanics. In fact, this induced differential structure, which we have been using successfully to model the mechanical behaviour of the Sierpi´nski carpet and triangle, may prove to be inadequate for the mechanical modelling of the Koch beam, namely, of a Koch curve endowed with bending resistance. The reason for this inadequacy can be easily understood when one considers that the bending stiffness of a beam must be somehow connected with its change of curvature, a concept that, in the case of a line-generated fractal, will not necessarily emerge automatically from the smooth embeddings already considered. Let γ denote a Koch curve generated upon the interval [0, 1] × {0} ⊂ R2 , and let H(γ ) denote its Hausdorff measure. From the process of generation of the Koch curve it is clear that every point x ∈ γ divides γ − {x} into two disjoint subsets. Let, therefore, a function f : γ → [0, H(γ )] be defined such that, for every x ∈ γ , f (x) is the Hausdorff measure of the subset of γ between (0, 0) and x. The arguments that follow are based on the following conjecture: Conjecture 5.1. f is continuous in the topology of γ induced by its inclusion in R2 and has a continuous inverse. The above conjecture implies that γ is a one-dimensional topological manifold with boundary. We have no proof of this conjecture, although its various aspects are stated in the literature as facts. We now endow γ with a differential structure, which we will call the Hausdorff differential structure, by declaring as smooth functions the compositions φ ◦ f of f with arbitrary smooth functions φ on I and completing the collection of smooth functions as required by the definition of a differential space. Remark 5.2. The function f can be interpreted as an “arc length” defined in terms of the Hausdorff measure. With this interpretation, the Hausdorff differential structure is given by smooth functions of this arc length.

´ M. Epstein, J. Sniatycki / Physica D 220 (2006) 54–68

67

One of the attractive features of this differential space is that on it we can establish a version of the fundamental theorem of Calculus, in terms of Hausdorff integration. Indeed, let g be a smooth function on γ . Then g ◦ f −1 is a smooth (and, therefore, integrable) function on I , whose Riemann definite integral coincides with our previous definition of Hausdorff integration over a fractal. On the other hand, we can now consider the Hausdorff integral with variable upper limit and, clearly, the fundamental theorem of Calculus applies. The same pull-back technique can be used to define derivatives of smooth functions with respect to the Hausdorff measure. In other words, the whole apparatus of Calculus becomes available in a natural fashion. In particular, consider a configuration of the Koch beam, namely, a smooth embedding κ : γ → R2 . One way to specify such an embedding is by prescribing two smooth functions, u and v say, representing the horizontal and vertical displacements of the beam.3 Considering the running parameter r of I as measuring Hausdorff length in γ , we have at our disposal all the tools necessary to define changes of slope and of curvature at any point of the Koch beam and to mimic the bending energy expression from the classical theory. Under the usual assumption that the beam is practically inextensible, so that the internal virtual work is due mainly to the bending moment, and assuming the latter to be constitutively related to the local change in curvature of the beam axis, we write the following expression for the virtual work in a Koch beam with one end clamped (at the origin, say) while the other is free and subjected to a couple M1 : Z dθ dθ K δ dH − M1 δθ1 , VW = (5.43) dr dr γ where K is a bending stiffness constant and θ represents the change of slope. Notice that the bending stiffness constant must be expressed in appropriate units of force divided by units of length raised to the power (3 − s), where s is the Hausdorff dimension of the Koch curve. Clearly, the identical vanishing of the virtual work will be achieved when: K

dθ = M1 , dr

(5.44)

so that the internal bending moment is constant and equal to the externally applied couple and θ is a linear function of r , as expected. In this particular case, it would not have been difficult to obtain this result (or similar results for arbitrary end loadings) on the basis of an ad hoc heuristic argument in which the classical beam problem is solved at a typical step in the generation of the curve and then a passage to the limit is invoked. Nevertheless, if such a technique were used, a rule for modifying the bending stiffness as the generation of the curve progresses would have to be specified. For the limit to be finite, the stiffness would have to be modified in step with the Hausdorff dimension of the curve. A detailed analysis reveals that at each successive step of the generation, the bending stiffness must necessarily increase by the factor 32−s for the limiting rotation of the free end to remain finite and non-zero. The physical justification of this law should be found within the particular context in which the problem arises (for example, it could be due to the exposure of more and more lateral area, which may require extra energy to create). When working in the more rigorous context of the Hausdorff differential structure, on the other hand, this supposed law of stiffness evolution is not needed, since the bending stiffness already incorporates the Hausdorff measure in its very definition. To complete our treatment of the Koch beam, we will explicitly calculate the stiffness matrix of a beam element with respect to its end degrees of freedom. Let the end points A and B of the undeformed Koch curve be located, respectively, at (0, 0) and (b, 0) in R2 . Let the left end A be fully fixed, as before, while the right end B is subjected to, say, a unit horizontal force that produces an internal bending moment M(r ). According to the principle of virtual work we must have that: Z dθ (5.45) Mδ dH = 1δu B , dr γ identically for all variations. Notice that the relationship between δθ and δu B is a very subtle one. Specifically: Z δu B = δθ sin(α)dH, γ

(5.46)

where α represents the local “slope” of the Koch curve. Although this slope is not necessarily smooth, it is certainly integrable. Indeed, the (Hausdorff) integral between two points of the curve is simply the difference between the corresponding y-coordinates of the points. To obtain our stiffness coefficients, however, we may bypass this relation by using the so-called “unit-force” method, valid for infinitesimal deformations. Indeed, we may consider as a virtual deformation the real deformation due to, successively, a unit horizontal force, a unit vertical force and a unit couple applied at B. The bending moments due to these actions are easily evaluated (from standard equilibrium considerations), respectively, as: y, b − x and 1. By the linearity of the material response, the corresponding rates of rotation are obtained by simply dividing these values by the stiffness constant K . It follows that the

3 In fact, a more detailed analysis of this case should be based on C k (rather than C ∞ ) differential structures.

´ M. Epstein, J. Sniatycki / Physica D 220 (2006) 54–68

68

horizontal displacements at point B due to these actions are given, according to (5.45), respectively by: Z 2 y d11 = dH, γ K Z y(b − x) dH, d12 = K γ

(5.47) (5.48)

and Z d13 =

γ

y dH, K

(5.49)

These definite Hausdorff integrals can be evaluated by implementing numerically the definition of integration already used in the Sierpi´nski carpet, namely, by subdividing the fractal into smaller self-similar pieces. The result (to 4 accurate significant digits) is: b2+s b1+s b2+s H(1) d12 = 0.048113 H(1) d13 = 0.096225 H(1), (5.50) K K K where H(1) is the Hausdorff measure of the unit-based Koch curve and s its Hausdorff dimension. It is not difficult, following this procedure, to complete the 3 × 3-symmetric flexibility matrix with the entries: d11 = 0.016667

b1+s bs b2+s H(1) d23 = 0.500000 H(1) d33 = 1.000000 H(1). (5.51) K K K The inverse of the matrix {di j } thus constructed occupies the lower right corner of the desired 6 × 6-stiffness matrix, which can now be completed by elementary means. d22 = 0.316667

6. Concluding remarks We close by noting that there are many open problems in Fractal Mechanics, some of them stemming from the somewhat unchartered mathematical territory covered. In particular, within the Finite Element context, recall that we do not know how the Hausdorff measure changes under affine transformations. Since most domains cannot be closely tiled with, say, equilateral triangles, an obvious problem arises. We hope, nevertheless, to have provided a setting for the formulation and numerical solution of problems in the mechanics of fractals and non-smooth bodies in general. Acknowledgment This work was partially supported by the Natural Sciences and Engineering Research Council of Canada. References [1] L. Ambrosio, N. Fusco, D. Pallara, Functions of Bounded Variation and Free Discontinuity Problems, Oxford Science Publishers, Clarendon Press, Oxford, 2000. [2] R. Capitanelli, M.R. Lancia, Nonlinear energy forms and Lipschitz spaces on the Koch curve, Journal of Convex Analysis 9 (2002) 245–257. [3] A. Carpinteri, B. Chiaia, P. Cornetti, Static–kinematic duality and the principle of virtual work in the mechanics of fractal media, Computer Methods in Applied Mechanics and Engineering 191 (2001) 3–19. [4] A. Carpinteri, P. Cornetti, A fractional calculus approach to the description of stress and strain localization in fractal media, Chaos, Solitons and Fractals 13 (2002) 85–94. [5] M. Epstein, R. Segev, Differentiable manifolds and the principle of virtual work in continuum mechanics, Journal of Mathematical Physics 21 (1980) 1243–1245. [6] K. Falconer, Fractal Geometry, John Wiley & Sons, 2003. [7] B.G. Jia, Z.L. Zhou, Z.W. Zhu, A lower bound for the Hausdorff measure of the Sierpi´nski gasket, Nonlinearity 15 (2002) 393–404. [8] U. Mosco, Energy functionals on certain fractal structures, Journal of Convex Analysis 9 (2002) 581–600. [9] G. Rodnay, R. Segev, Cauchy’s flux theorem in light of geometric integration theory, Journal of Elasticity 71 (2003) 183–203. [10] R. Sikorski, Abstract covariant derivative, Colloquium Mathematicum 18 (1967) 251–272. [11] R. Sikorski, Differential modules, Colloquium Mathematicum 24 (1971) 45–79. [12] R. Sikorski, Wst¸ep do geometrii Ró˙zniczkowej, PWN (Polish Scientific Publishers), Warszawa, 1972. [13] Z.L. Zhou, L. Feng, A new estimate of the Hausdorff measure of the Sierpi´nski gasket, Nonlinearity 13 (2000) 479–491. [14] Z.W. Zhu, Z.L. Zhou, B.G. Jia, On the lower bound of the Hausdorff measure of the Koch curve, Acta Mathematica Sinica, English Series 19 (2003) 715–728.