Asymptotically exact theory of fiber-reinforced composite beams
AAsymptotically exact theory of fiber-reinforcedcomposite beams
K. C. Le a,b , T. M. Tran a,b a Materials Mechanics Research Group, Ton Duc Thang University, Ho Chi Minh City,Vietnam b Faculty of Civil Engineering, Ton Duc Thang University, Ho Chi Minh City, Vietnam
Abstract
An asymptotic analysis of the energy functional of a fiber-reinforced com-posite beam with a periodic microstructure in its cross section is performed.From this analysis the asymptotically exact energy as well as the 1-D beamtheory of first order is derived. The effective stiffnesses of the beam are cal-culated in the most general case from the numerical solution of the cell andhomogenized cross-sectional problems.
Keywords: fiber-reinforced, composite, beam, periodic microstructure,variational-asymptotic method.
1. Introduction
Fibre reinforced composite (FRC) beams are widely used in civil, me-chanical, and aerospace engineering due to their low weight, high strength,and good damping properties (10; 9; 19). The exact treatment of FRCbeams within the 3D elasticity is only possible in a few exceptional casesdue to their complicated microstructure (see, e.g., (24; 26) and the refer-ences therein). For this reason, different approaches have been developeddepending on the type of beams. If FRC beams are thick, no exact one-dimensional theory can be constructed, so only the numerical methods orapproximate semi-analytical methods applied to three-dimensional elastic-ity theory make sense. However, if the FRC beams are thin, the reductionfrom the three-dimensional to the one-dimensional theory is possible. This Email: [email protected]
Preprint submitted to Composite Structures March 30, 2020 a r X i v : . [ phy s i c s . c l a ss - ph ] M a r imension reduction can be made asymptotically exact in the limit whenthe thickness-to-length ratio of the beam goes to zero. The rigorous deriva-tion of the asymptotically exact one-dimensional beam theory based on thevariational-asymptotic method (VAM) developed by Berdichevsky (3) wasfirst performed in (4). His asymptotic analysis shows that the static (or dy-namic) three-dimensional problem of the beam can be split into two problems:(i) the two-dimensional cross-sectional problem and (ii) the one-dimensionalvariational problem whose energy functional should be found from the so-lution of the cross-sectional problem. The latter has been solved both foranisotropic homogeneous beams and for inhomogeneous beams with the con-stant Poisson’s ratio (4). In addition to these findings, Berdichevsky (4) hasshown that the average energy as well as the extension and bending stiffnessesof FRC beams with piecewise constant Poisson’s ratio must be larger thanthose of FRC beams with constant Poisson’s ratio (see also (23)). However,to our knowledge, the question of how the corrections in energy and stiffnessesdepend on the difference in Poisson’s ratio of fibers and matrix of FRC beamsremains still an issue. It should be noted that VAM has been further devel-oped in connection with the numerical analysis of cross-sectional problemsfor geometrically nonlinear composite beams by Hodges, Yu and colleaguesin a series of papers (11; 28; 29; 30). Note also that VAM has been used,among others, to derive the 2D theory of homogeneous piezoelectric shells(13), the 2D theory of purely elastic anisotropic and inhomogeneous shells(5), the 2D theory of purely elastic sandwich plates and shells (6; 7), thetheory of smart beams (25), the theory of low and high frequency vibrationsof laminate composite shells (17; 18), and more recently, the theory of smartsandwich and functionally graded shells (15; 16).For FRC beams that have the periodic microstructure in the cross section,an additional small parameter appears in the cross-sectional problems: Theratio between the length of the periodic cell and the characteristic size ofthe cross section. In this case the finite element code VABS developed inthe above mentioned papers (11; 28; 29; 30) cannot be applied to the cross-sectional problem because it requires a large computational capacity. Thepresence of this small parameter allows however an additional asymptoticanalysis of the cross-sectional problems to simplify them. By solving the cellproblems imposed with the periodic boundary conditions according to thehomogenization technique (8; 22), one finds the effective coefficients in thehomogenized cross-sectional problems, which can then be solved analyticallyor numerically (cf. also (20)). The aim of this paper is to derive and solve2he cell and homogenized cross-sectional problems for unidirectional FRCbeams whose cross section has the periodic microstructure. For simplicity,we will assume that both matrix and fibers are elastically isotropic but havedifferent Poisson’s ratio. The solution of the cell problems found with thefinite element method is used to calculate the asymptotically exact energyand the extension and bending stiffnesses in the 1-D theory of FRC beams.Thus, we determine the dependence of the latter quantities on the shape andvolume fraction of the fibers and on the difference in Poisson’s ratio of fibersand matrix that solves the above mentioned issue.The paper is organized as follows. After this short introduction the vari-ational formulation of the problem is given in Section 2. Sections 3 and 4are devoted to the multi-scaled asymptotic analysis of the energy functionalof FRC beams leading to the cross-sectional and cell problems. In Section5 the cell problems are solved by the finite element method. Section 6 pro-vides the solutions of the homogenized cross-sectional problems. Section 7present one-dimensional theory of FRC beams. Finally, Section 8 concludesthe paper.
2. Variational formulation for FRC beams x x Figure 1: A cross section of a FRC beam
Consider a FRC beam that occupies the domain B = A × (0 , L ) of the3-D euclidean space in its undeformed state. Let x ≡ x be the coordinatealong the beam axis. The cross section of the beam in the ( x , x )-plane, A ,consist of two separated 2-D sub-domains A m and A f such that the matrixoccupies the domain A m × (0 , L ), while the uni-directional fibers occupy the3omain A f × (0 , L ). We choose the origin of the ( x , x )-coordinates so thatit matches the centroid of A . We assume that the fibers are periodicallysituated in the matrix and the bond between the fibers and the matrix isperfect (see Fig. 1 representing the cross section of the beam where A f isthe set of gray circles). We first consider the beam in equilibrium, whosedeformation is completely determined by the displacement field u ( x ). Let,for simplicity, the edge x = 0 of the beam be clamped, while at the otheredge x = L of the beam the traction t ( x , x ) be specified. Gibbs’ variationalprinciple (5) states that the true displacement ˇ u ( x ) minimizes the energyfunctional I [ u ( x )] = (cid:90) B W ( x , ε ) d x − (cid:90) A t ( x , x ) · u ( x , x , L ) d x (1)among all admissible displacements u ( x ) satisfying the boundary conditions u ( x , x ,
0) = 0 . Here, d x = d x d x d x is the volume element, d x = d x d x the area element,and the dot denotes the scalar product. Function W ( x , ε ), called storedenergy density, reads W ( x , ε ) = 12 λ ( x )(tr ε ) + µ ( x ) ε : ε , where ε is the strain tensor ε = 12 ( ∇ u + ( ∇ u ) T ) . The problem is to replace the three-dimensional energy functional by anapproximate one-dimensional energy functional for a thin FRC beam, whosefunctions depend only on the longitudinal co-ordinate x . The possibilityof reduction of the three- to the one-dimensional functional is related to thesmallness of two parameters: (i) the ratio between the thickness h of the crosssection and the length L of the beam, and (ii) the ratio between the size c of the periodic cell and h . By using the variational-asymptotic method, theone-dimensional energy functional will be constructed below in which termsof the order h/L and (cid:15) = c/h are neglected as compared with unity (thefirst-order or “classical” approximation). Formally, this corresponds to thelimits h → (cid:15) →
0. 4n order to perform the variational-asymptotic analysis, it is convenient touse the index notation, with Greek indices running from 1 to 2 and denotingthe component of vectors and tensors in the ( x , x )-plane. Summation overrepeated Greek indices is understood. Index 3 of the coordinate x , thedisplacement u , and the traction t is dropped for short. To fix the domainof the transverse co-ordinates in the passage to the limit h →
0, we introducethe dimensionless co-ordinates y α = x α h , y α ∈ ¯ A , and transform the energy functional of the beam to I = (cid:90) L (cid:90) ¯ A h W ( y α , ε ) d y d x − (cid:90) ¯ A h t ( y α ) · u ( y α , L ) d y . (2)The transverse coordinates y α play the role of the “fast” variables as opposedto the slow variable x . Regarding u as function of y α and x , we separate thefast variables from the slow one. Now h enters the action functional explicitlythrough the components of the strain tensor ε ε αβ = 1 h u ( α ; β ) , ε α = 1 h u ; α + u α,x , ε = u ,x . (3)Here and below, the semicolon preceding Greek indices denotes the deriva-tives with respect to y α , while the parentheses surrounding a pair of indicesthe symmetrization operation.Besides y α there are also much faster variables z α = x α c = y α (cid:15) associated with the periodicity of the elastic moduli of this composite materi-als in the transverse directions leading to the fast oscillation of the stress andstrain fields in these directions. In order to characterize these faster changesof the stress and strain fields in the transverse directions we will assume thatthe displacement field u is a composite function of y and x such that u = u ( y , y /(cid:15), x ) = u ( y , z , x ) , where y , z are two-dimensional vectors with the components y α , z α , respec-tively, and where u ( y , z , x ) is a doubly periodic function in z with the period1. This is a typical multi-scale Ansatz to composite structures. The asymp-totic analysis must therefore be performed twice, first in the limit h → (cid:15) → . Dimension reduction Before starting the asymptotic analysis of the energy functional in thelimit h → W ( z , ε ),the derivatives w ( α ; β ) /h in ε αβ and w ; α /h in ε α are the main ones in theasymptotic sense. Therefore it is convenient to single out the components ε αβ and ε α in the stored energy density. We represent the latter as the sumof three quadratic forms W (cid:107) , W ∠ , and W ⊥ according to W (cid:107) = min ε αβ ,ε α W, W ∠ = min ε αβ ( W − W (cid:107) ) , W ⊥ = W − W (cid:107) − W ∠ . The “longitudinal” energy density W (cid:107) depends only on ε and coincideswith W when the stresses σ αβ and σ α vanish; the “shear” energy density W ∠ depends only on ε α ; the remaining part W ⊥ is called the “transverse”energy density. From the definitions of W (cid:107) , W ∠ and W ⊥ one easily finds outthat W (cid:107) = 12 E ( ε ) , E = µ (3 λ + 2 µ ) λ + µ ,W ∠ = 2 µε α ε α ,W ⊥ = 12 λ ( ε αα + 2 νε ) + µ ( ε αβ + νδ αβ ε )( ε αβ + νδ αβ ε ) , where E is the Young modulus and ν = λ/ λ + µ ) the Poisson ratio. Notethat E , ν as well as Lame’s constants λ , µ are doubly periodic functions ofthe fast variables z . Note also the following identities W ∠ = σ α ε α ,W ⊥ = 12 ( σ αβ ε αβ + σ ε − E ( ε ) ) . We could start the variational-asymptotic analysis in the limit h → N according to its general scheme (14).As a result, it would turn out that, at the first step, the function u does notdepend on the transverse co-ordinates y (and z ): u = v ( x ); at the secondstep the function u (cid:63) becomes a linear function of y and involves one moredegree of freedom ϕ ( x ) representing the twist angle; and at the next step u (cid:63)(cid:63) is completely determined through u and ϕ . Thus, the set N according to6he variational-asymptotic scheme consists of functions v ( x ) and ϕ ( x ). Wewill pass over these long, but otherwise standard, deliberations and make achange of unknown functions according to u α ( y , z , x ) = v α ( x ) − he αβ ϕ ( x ) y β + hw α ( y , z , x ) ,u ( y , z , x ) = v ( x ) − hv α,x ( x ) y α + hw ( y , z , x ) , (4)where e αβ are the two-dimensional permutation symbols ( e = e = 0 , e = − e = 1). By redefining v ( x ) and ϕ ( x ) if necessary we can impose onfunctions w α and w the following constraints (cid:104) w α (cid:105) = 0 , e αβ (cid:104) w α ; β (cid:105) = 0 , (cid:104) w (cid:105) = 0 , (cid:104) . (cid:105) ≡ (cid:90) ¯ A . d y . (5)According to these constraints v α ( x ) and v ( x ) describe the mean displace-ments of the beam, while ϕ ( x ) corresponds to the mean rotation of its crosssection. Equations (4) and (5) set up a one-to-one correspondence between u α , u and the set of functions v α , v, ϕ, w α , w and determine the change in theunknown functions { u α , u } → { v α , v, ϕ, w α , w } .Based on the Saint-Venant principle for elastic beams, we may assumethat the domain occupied by the beam consists of the inner domain and twoboundary layers near the edges of the beam with width of the order h wherethe stress and strain states are essentially three-dimensional. Then functional(2) can be decomposed into the sum of two functionals, an inner one for whichan iteration process will be applied, and a boundary layer functional. Whensearching for w α and w the boundary layer functional can be neglected inthe first-order approximation. Therefore, the dimension reduction problemreduces to finding the minimizer ˇ w α and ˇ w of the inner functional that, inthe limit h →
0, can be identified with the functional (2) without the lastterm.We now fix v α , v, ϕ and seek w α , w . Substituting (4) into the action func-tional (2) with the last term being removed, we will keep in it the asymp-totically principal terms containing w α , w and neglect all other terms. Theestimations based on Eqs. (3) and (4) lead to the asymptotic formulas ε αβ = w ( α ; β ) , ε α = w ; α − he αβ Ω y β , ε = γ + h Ω α y α , (6)where γ , Ω α , and Ω are the measures of elongation, bending, and twist definedby γ = v ,x , Ω α = − v α,xx , Ω = ϕ ,x . w α , w with respect to x do not enter the energy functional. As x becomes the formal parameter, wemay drop the integral over x and reduce the determination of w α , w to theuncoupled minimization problems for every fixed x of the functionals I ⊥ [ w α ] = h (cid:104) λ ( z )[ w α ; α + 2 ν ( z )( γ + h Ω σ y σ )] + µ ( z )[ w ( α ; β ) + ν ( z ) δ αβ ( γ + h Ω γ y γ )][ w ( α ; β ) + ν ( z ) δ αβ ( γ + h Ω δ y δ )] (cid:105) , (7) I ∠ [ w ] = h (cid:104) µ ( z )( w ; α − h Ω e αβ y β )( w ; α − h Ω e αγ y γ ) (cid:105) . (8)The minima are searched among all admissible functions w α , w satisfyingthe constraints (5). Note that the decoupling of problems (7) and (8) holdstrue in the most general case of anisotropy (5). This can be seen from theasymptotically main terms containing the unknown functions w α and w inthe transverse and shear strain energy densities: w ; α does not enter W ⊥ while w α ; β does not enter W ∠ . Functionals (7) and (8) represent the transverse andshear strain energies, integrated over the cross section of the beam. Theyare positive definite and convex, so the existence of their minimizers ˇ w α , ˇ w is guaranteed. We shall see in the next Sections that the minimum of (7) isequal to zero if ν is equal for both matrix and fiber, while that of (8) is equalto 1 / C Ω , with C the torsional rigidity.
4. Homogenization
Consider now the other limit (cid:15) →
0. In this limit y plays the role ofthe “slow” variable, while z = y /(cid:15) becomes the fast variable. We start withthe cross-sectional problem of minimizing functional (8) among w ( y , z , x )satisfying the constraint (5) , where µ ( z ), expressed in terms of the fastvariable z = y /(cid:15) , is a doubly periodic function with period 1. Since x isfixed in this cross-sectional problem, we shall drop this formal variable of w in this Section. Following the homogenization technique based also on thevariational-asymptotic method (5), we look for the minimizer in the form w ( y , z ) = ψ ( y ) + (cid:15)χ ( z ) , (9)where χ ( z ) is a doubly periodic function with period 1. Note that χ maydepend on the slow variable y , but this dependency is suppressed for short.8n addition to the constraint (5) we may impose the following constraint on χ ( z ) (cid:104)(cid:104) χ (cid:105)(cid:105) ≡ (cid:90) C χ ( z ) d z = 0 , (10)where C = (0 , × (0 ,
1) is the unit periodic cell. In this case ψ ( y ) can beinterpreted as the average value of w ( y , z ) over the cell. Although the secondterm in (9) is small compared to ψ ( y ) and goes to zero in the limit (cid:15) → w ( y , z ) has the sameorder as the gradient of ψ ( y ) as seen from the asymptotic formula w ; α = ψ ; α + χ | α , α = 1 , . (11)Here, the vertical bar followed by an index α denotes the partial derivativewith respect to the corresponding fast variable z α .Inserting (11) into the energy functional (8), we get I ∠ [ ψ ( y ) , χ ( y /(cid:15) )] = h (cid:104) µ ( y /(cid:15) )( χ | α + ξ α )( χ | α + ξ α ) (cid:105) , with ξ α = ψ ; α − h Ω e αβ y β being the function of the “slow” variable y . Wereplace the double integral (cid:104) . (cid:105) by the sum of double integrals over the cells.Then I ∠ [ ψ ( y ) , χ ( y /(cid:15) )] = h N (cid:88) n =1 (cid:90) C n µ ( y /(cid:15) )( χ | α + ξ α )( χ | α + ξ α ) d y , (12)where N is the total number of cells and C n is the cell (( i − (cid:15), i(cid:15) ) × (( j − (cid:15), j(cid:15) ). We minimize this functional in two steps: (i) Fix ψ ( y ) and minimizethe functional among doubly periodic χ ( y /(cid:15) ), (ii) Minimize the obtainedfunctional among ψ ( y ). Because χ ( y /(cid:15) ) is doubly periodic with respect to y , we can minimize each integral in the sum independently. Besides, as ξ α ( y ) change slowly in one cell, we may regard them in each cell integral asconstants equal to their value in the middle of the cell. It is convenient tochange the variable y in the cell integrals to z = y /(cid:15) − ( i − , j − (cid:104)(cid:104) µ ( z )( χ | α + ξ α )( χ | α + ξ α ) (cid:105)(cid:105) (13)among doubly periodic functions χ ( z ) satisfying the constraint (10), where,as before, (cid:104)(cid:104) . (cid:105)(cid:105) is the double integral in z over the unit cell. We denote the9btained minimum by ¯ W ∠ ( ξ α ) and call it the average shear energy density.Then, replacing the sum in (12) by the integral in the limit (cid:15) →
0, we arriveat the following homogenized cross-sectional problem: Minimize the averagefunctional ¯ I ∠ [ ψ ( y )] = h (cid:104) ¯ W ∠ ( ψ ; α − h Ω e αβ y β ) (cid:105) (14)among ψ ( y ) satisfying the condition (cid:104) ψ (cid:105) = 0 . The standard calculus of variations shows that the minimizer of (13)satisfies the equation [ µ ( z )( χ | α + ξ α )] | α = 0 (15)in the unit cell, where ξ can be regarded as the constant vector in the cell(since it does not depend on the fast variable z ). This equation is subjectedto the periodic boundary condition and the constraint (10). Besides, as theshear modulus suffers a jump at the boundary ∂ i between the fiber and thematrix, the continuity of the traction[[ µ ( z )( χ | α + ξ α )]] ν α = 0 (16)should be fulfilled at this internal boundary ∂ i , where ν is the unit normalvector outward to the fiber and [[ . ]] denotes the jump.After finding the minimizer ˇ χ ( z ) as a solution to the boundary-valueproblem (15)-(16), we substitute it back into functional (13) to calculate theaverage shear energy density¯ W ∠ ( ξ α ) = (cid:104)(cid:104) µ ( z )( ˇ χ | α + ξ α )( ˇ χ | α + ξ α ) (cid:105)(cid:105) = (cid:104)(cid:104) µ ( z )( ˇ χ | α + ξ α ) ˇ χ | α (cid:105)(cid:105) + (cid:104)(cid:104) µ ( z )( ˇ χ | α + ξ α ) ξ α (cid:105)(cid:105) . (17)As the minimizer ˇ χ satisfies (15)-(16), the first integral must vanish. Usingthe constancy of ξ α , we reduce the second integral to¯ W ∠ ( ξ α ) = 12 ξ α (cid:104)(cid:104) µ ( z )( ˇ χ | α + ξ α ) (cid:105)(cid:105) . The integrand is the shear stress σ α , so the integral gives the average shearstress ¯ σ α . On the other side, due to the constraint (10), (cid:104)(cid:104) ˇ χ | α + ξ α (cid:105)(cid:105) = ξ α , ξ α is the average engineering shear strain ¯ ε α . We define the effectiveelastic shear moduli µ ∗ αβ by the linear equation¯ σ α = (cid:104)(cid:104) µ ( z )( ˇ χ | α + ξ α ) (cid:105)(cid:105) = µ ∗ αβ ξ β . (18)Note that the homogenized material must not necessarily be isotropic evenif the components of the composite are. So, in general µ ∗ αβ is a tensor ofsecond order. With this we get the average shear energy density in terms of ξ α = ψ ; α − h Ω e αβ y β ¯ W ∠ ( ξ α ) = 12 µ ∗ αβ ξ α ξ β . (19)We want to show that the tensor µ ∗ αβ in (18) is symmetric. Indeed, since¯ W ∠ ( ξ α ) in (19) is a quadratic form, we can replace µ ∗ αβ there by the symmetrictensor µ ∗ ( αβ ) = ( µ ∗ αβ + µ ∗ βα ). If we substitute (19) into (17) and differentiateit with respect to ξ α , we get µ ∗ ( αβ ) ξ β = (cid:104)(cid:104) µ ( z )( ˇ χ | α + ξ α ) (cid:105)(cid:105) = µ ∗ αβ ξ β , which proves the symmetry of µ ∗ αβ .We turn now to the cross-sectional problem (7). Before doing the asymp-totic analysis for it in the limit (cid:15) → ν isan equal constant for both fiber and matrix (4). Indeed, let us choose theminimizer in the formˇ w α = − νγδ αβ y β − νa αβγ ( y β y γ − (cid:104) y β y γ (cid:105) / | ¯ A| ) , (20)where | ¯ A| is the area of ¯ A and a αβγ = h ( δ αβ Ω γ + δ αγ Ω β − δ βγ Ω α ) . Note that the tensor a αβγ is symmetric with respect to the last two indices,and ˇ w α ; β = − νγδ αβ − a αβγ y γ . Due to our choice of the origin, (cid:104) y α (cid:105) = 0. Therefore the chosen field (20)satisfies the constraints (5) , . It is easy to check that the transverse energyevaluated at ˇ w α vanishes identically. Since functional (7) is non-negativedefinite, its minimum is obviously zero in this case.11he asymptotic analysis of problem (7) for changeable ν is quite similarto that of problem (8). In this case the minimizer is sought in the form w α ( y , z ) = ψ α ( y ) + (cid:15)χ α ( z ) , where χ α satisfy the constraints (cid:104)(cid:104) χ α (cid:105)(cid:105) = 0 . (21)To determine χ α we need to solve the following cell problem: Minimize thefunctional (cid:104)(cid:104) λ ( z )( χ α | α + ¯ ε αα + 2 ν ( z ) ξ ) + µ ( z )( χ ( α | β ) + ¯ ε αβ + ν ( z ) δ αβ ξ )( χ ( α | β ) + ¯ ε αβ + ν ( z ) δ αβ ξ ) (cid:105)(cid:105) (22)among doubly periodic functions χ α ( z ) satisfying the constraint (21), where¯ ε αβ = ψ ( α | β ) , ξ = γ + h Ω α y α . Note that this problem is quite similar to the problem of determining the ef-fective thermal expansion of composite material with periodic microstructure,where ξ plays the role of the temperature increase (27). Let the minimumof functional (22) be ¯ W ⊥ (¯ ε αβ , ξ ). Then the determination of ψ α reduces tominimizing the following average functional¯ I ⊥ [ ψ α ( y )] = h (cid:104) ¯ W ⊥ ( ψ ( α ; β ) , γ + h Ω α y α ) (cid:105) (23)among ψ α ( y ) satisfying the constraints (cid:104) ψ α (cid:105) = 0 , e αβ (cid:104) ψ α ; β (cid:105) = 0 . (24)The minimizer of functional (22) satisfy the equilibrium equations[ λ ( z )( χ β | β + ¯ ε ββ + ξ )] | α + [2 µ ( z )( χ ( α | β ) + ¯ ε αβ )] | β = 0in the unit cell, where ¯ ε αβ and ξ can be regarded as the constant tensorand scalar in the cell (since they do not depend on the fast variable z ).These equations are subjected to the periodic boundary conditions and theconstraints (21). Besides, as the elastic moduli suffers jumps at the boundary ∂ i between the fiber and the matrix, the continuity of the traction[[ λ ( z )( χ β | β + ¯ ε ββ + ξ ) ν α + 2 µ ( z )( χ ( α | β ) + ¯ ε αβ ) ν β ]] = 012hould be fulfilled there. Since the functional (22) is quadratic, its minimummust be a quadratic form of ¯ ε αβ and ξ ¯ W ⊥ (¯ ε αβ , ξ ) = 12 C ∗ αβγδ ¯ ε αβ ¯ ε γδ + D ∗ αβ ¯ ε αβ ξ + 12 F ∗ ξ . (25)Using the same arguments as in the previous case, we can prove the followingsymmetry properties of the effective moduli C ∗ αβγδ = C ∗ βαγδ = C ∗ αβδγ = C ∗ γδαβ ,D ∗ αβ = D ∗ βα . We also introduce the effective tensor of Poisson’s ratios by ν ∗ αβ = C ∗ ( − αβγδ D γν , where C ∗ ( − αβγδ is the tensor of elastic compliances defined by C ∗ ( − αβγδ C γδζη = δ α ( ζ δ βη ) . It is easy to see that ν αβ is symmetric.
5. Numerical solution of the cell problems
It is convenient to rewrite functional (22) by changing the sign of χ α ( y )12 (cid:104)(cid:104) C αβγδ ( z )( − χ ( α | β ) + ¯ ε αβ + α αβ ( z ) ξ )( − χ ( γ | δ ) + ¯ ε γδ + α γδ ( z ) ξ ) (cid:105)(cid:105) , (26)where C αβγδ ( z ) = λ ( z ) δ αβ δ γδ + µ ( z )( δ αγ δ βγ + δ αδ δ βγ ) , α αβ ( z ) = ν ( z ) δ αβ . The minimizer of functional (26) satisfies the variational equation (cid:104)(cid:104) C αβγδ ( z ) χ ( γ | δ ) δχ α | β (cid:105)(cid:105) = (cid:104)(cid:104) C αβγδ ( z )(¯ ε γδ + α γδ ( z ) ξ ) δχ α | β (cid:105)(cid:105) (27)for all doubly periodic functions δχ α . Eq. (27) will be solved by the finiteelement method (12; 1). For this purpose it is convenient to change fromtensor notation to matrix notation (12), in which Eq. (27) becomes (cid:104)(cid:104) ε T ( ψ ) C ( z ) ε ( χ ) (cid:105)(cid:105) = (cid:104)(cid:104) ε T ( ψ ) C ( z )(¯ ε + α ( z ) ξ ) (cid:105)(cid:105) , (28)13here ψ = δ χ , ε ( χ ) is the strain vector given by ε ( χ ) = ( χ | , χ | , χ | + χ | ) T , while C ( z ) = λ ( z ) + µ ( z ) , α ( z ) = ν ( z ) . (29)The discretization of this equation based on the bilinear isoparametric ele-ments is standard. The periodic boundary conditions are imposed by iden-tifying the nodes on opposite sides of the unit cell. This is implemented byusing the matrix edofMat for a full, regular grid to index into a periodicversion of the grid (1). The global stiffness matrix is K = N A e =1 ( k e ) , k e = (cid:90) A e B Te C e B e d z , where A is the assembly operator taken over the total number N of finiteelements. The matrix B e is the element strain-displacement matrix, A e isthe domain of element e , and C e is the constitutive matrix for the element.The indicator matrix is introduced that specifies whether the element is in A m ( n e = 1) or in A f ( n e = 2). The piecewise constant C ( z ) from (29) takesvalue C e in the element in accordance with this indicator matrix.The discretization of the right hand side of (28) yields the loads f ( i ) f ( i ) = N A e =1 ( f ( i ) e ) , f ( i ) e = (cid:90) A e B Te C e ¯ ε ( i ) d z , i = 1 , , ε (1) = (1 , , T , ¯ ε (2) = (0 , , T , ¯ ε (3) = (0 , , T , (30)and f (4) = N A e =1 ( f (4) e ) , f (4) e = (cid:90) A e B Te C e α e d z , which correspond to ξ = 1, with α e being the constitutive vector from (29)that takes values in the element in accordance with the indicator matrix.14he displacement fields are computed by solving the corresponding linearequations with four load-cases K χ ( i ) = f ( i ) , i = 1 , , , . When the displacements are obtained, the elements of the effective matrix C ∗ are found as: C ∗ ij = N (cid:88) e =1 (cid:90) A e ( χ i ) e − χ ( i ) e ) T k e ( χ j ) e − χ ( j ) e ) d z , i, j = 1 , , , where χ i ) e are the three displacement fields corresponding to the unit strainsfrom (30), and χ ( i ) e contains three columns corresponding to the three dis-placement fields resulting from globally enforcing the unit strains (30). Theindices in parentheses refer to the column number. The components of theeffective moduli D ∗ i ( D ∗ αβ ) due to ξ are computed according to D ∗ i = N (cid:88) e =1 (cid:90) A e ( α e − B e χ (4) e ) T C e (¯ ε i − B e χ ( i ) e ) d z , i = 1 , , . Having D ∗ i we can also compute the effective “Poisson” ratios ν ∗ i ( ν ∗ αβ ) asfollows ν ∗ i = ( C ∗ ij ) − D ∗ j . (31)Finally, the effective coefficient F ∗ is given by F ∗ = N (cid:88) e =1 (cid:90) A e ( α e − B e χ (4) e ) T C e ( α e − B e χ (4) e ) d z . The Matlab-code homogenizecs.m to solve the cell problem and compute theeffective moduli which is a modification of the code homogenize.m writtenby Andreasen (1) is presented in the Appendix.For the FRC bar it will be shown in the next Section that the minimumof functional (23) is ( F ∗ − ν ∗ αβ C ∗ αβγδ ν ∗ γδ ) h (cid:104) ξ (cid:105) . Therefore it makes senseto investigate the quantity H ∗ = F ∗ − ν ∗ αβ C ∗ αβγδ ν ∗ γδ giving the correction tothe stiffnesses on extension and bending. We take the microstructure of thecomposite in such a way that the cross sections of the fibers are ellipses of15 H * (GPa) r a=1a=2a=4a=8 Figure 2: The plot of H ∗ versus the largest half-axis of fibers of elliptical cross sections,where E = 1 GPa, E = 2 GPa, ν = 0 . ν = 0 .
4: (i) a = 1 (blue), (ii) a = 2 (orange),(iii) a = 4 (red), (iv) a = 8 (black). half-axes r and r/a placed in the middle of the unit quadratic periodic cells.Then the region occupied by one fiber in the unit cell is given by the equation( z − / r + ( z − / ( r/a ) ≤ . Note that for the circle with a = 1 the effective Poisson ratio tensor is ν ∗ αβ = ν ∗ δ αβ , however this property is no longer valid for ellipses with a > E = 1 GPa, E = 2 GPa, ν = 0 . ν = 0 .
4. The plot of H ∗ as function of the largest half-axis r for four different aspect ratios of the ellipses a = 1 , , , H ∗ varies in the range (0 , .
03) GPa whichis not so large compared with the mean value of the Young moduli. Since thevolume fraction of the fibers is s = πr /a , we can also plot H ∗ as functionof this volume fraction.Let us now fix the largest half-axis of elliptical cross section of the fibersto be r = 0 .
4, choose E = 1 GPa, E = 2 GPa, ν = 0 .
3, vary the Poissonratio of the fiber ν in its admissible range between − .
5, and plot H ∗ as function of ν . Looking at this plot shown in Fig. 3 we see that the largerthe difference between the Poisson ratios, the larger is the correction H ∗ . For ν = 0 . H ∗ vanishes as expected. For ν near − H * (GPa) a=1a=2a=4a=8 ν Figure 3: The plot of H ∗ versus ν of fibers of elliptical cross sections, where r = 0 . E = 1 GPa, E = 2 GPa, ν = 0 .
3: (i) a = 1 (blue), (ii) a = 2 (orange), (iii) a = 4 (red),(iv) a = 8 (black). modulus which can no longer be neglected.The solution of the cell problem (13) by the finite element method is quitesimilar. Its solution can be obtained from the previous problem if we put λ ( z ) = 0. Below we present the numerical solution for the fibers with circularcross sections periodically embedded in the matrix. The unit cell is a squareof length 1 with the circle of radius r in the middle representing the crosssection of the fiber. In this case the effective tensor is µ ∗ αβ = µ ∗ δ αβ . Theplot of µ ∗ as function of s = πr is presented in Fig. 4. For the numericalsimulation we choose µ = 1 GPa, µ = 2 GPa.
6. Solution of the homogenized cross-sectional problems
Let us first consider the homogenized cross-sectional problem (23) withthe average transverse energy being given by (25). It is convenient to presentthe latter formula in the form¯ W ⊥ (¯ ε αβ , ξ ) = 12 C ∗ αβγδ (¯ ε αβ + ν ∗ αβ ξ )(¯ ε γδ + ν ∗ γδ ξ ) + 12 H ∗ ξ , (32)where ξ = γ + h Ω α y α . Note that the second term in (32) does not dependon ψ α ( y ), so we need to minimize functional (23) that contains just thefirst term. We want to show that the minimum of the functional (23) thatcontains just the first term among ψ α ( y ) satisfying the constraints (24) is17 μ * (GPa) s Figure 4: The plot of µ ∗ versus the volume fraction s = πr of fibers of circular crosssection, where µ = 1 GPa, µ = 2 GPa. equal to zero. Indeed, following the same line of reasoning as in Section 4 wechoose the minimizer in the formˇ ψ α = − ν ∗ αβ γy β − a αβγ ( y β y γ − (cid:104) y β y γ (cid:105) / | ¯ A| ) , where a αβγ = h ( ν ∗ αβ Ω γ + ν ∗ αγ Ω β − ν ∗ βγ Ω α ) . It is easy to check that ˇ ψ α satisfy the constraints (24) and that the first termevaluated at these functions vanishes identically, somin ψ α ( y ) ∈ (24) ¯ I ⊥ [ ψ α ( y )] = 12 H ∗ h (cid:104) ξ (cid:105) . Adding the average longitudinal energy density h (cid:104) W (cid:107) (cid:105) = ¯ E (cid:104) ξ (cid:105) to thisaverage transverse energy density and integrating ξ = ( γ + h Ω α y α ) overthe cross section, we find the final energy density of extension and bendingof FRC beam in the formΦ( γ, Ω α ) = 12 ( ¯ E + H ∗ )( |A| γ + I αβ Ω α Ω β ) , where ¯ E = (cid:104)(cid:104) E (cid:105)(cid:105) = s E + s E , I αβ = h (cid:104) y α y β (cid:105) . Thus, H ∗ is the correction to the stiffnesses on extension and bending ofFRC beam. To see how this correction changes the stiffnesses on extension18nd bending of FRC beam we plot both ¯ E and ¯ E + H ∗ as function of thevolume fraction of fibers with circular cross sections placed in the middleof the quadratic periodic cell (see Fig. 5). The chosen elastic moduli are: E = 1 GPa, E = 2 GPa, ν = 0 . ν = 0 .
4. The correction is about 3percent of ¯ E in this case. Note that this correction becomes one third of ¯ E if ν is near − without correctionwith correction E+H * (GPa) s _ Figure 5: The plot of ¯ E and ¯ E + H ∗ versus the volume fraction s = πr of fibers ofcircular cross section, where E = 1 GPa, E = 2 GPa, ν = 0 . ν = 0 . We now turn to the minimization problem (14). First of all it is easyto see that the constraint (cid:104) ψ (cid:105) = 0 does not affect the minimum value of(14), because the latter is invariant with respect to the change of unknownfunction ψ → ψ + c , with c being a constant. By such the change one canalways achieve the fulfillment of the constraint (cid:104) ψ (cid:105) = 0. The minimizer ˇ ψ of¯ I ∠ is proportional to h Ω. Therefore the torsional rigidity C can be calculatedby C = h inf ¯ ψ (cid:104) µ ∗ δ αβ ( ¯ ψ | α − e ασ y σ )( ¯ ψ | β − e βκ y κ ) (cid:105) , (33)where ¯ ψ = ψ/h Ω. The solution to this problem is well-known for the ellipticaland rectangular cross sections (14). For the FRC beam with the ellipticalcross section described by the equation c αβ y α y β ≤ , where c αβ are the components of a positive definite symmetric second-rank19ensor, the torsional rigidity reads C = 4 µ ∗ c c µλ c λτ I µτ , (34)where I µτ = h (cid:104) y µ y τ (cid:105) are the moments of inertia of the cross section.In the co-ordinate system associated with the principal axes of the ellipse c = h b , c = h b , c = 0 ,I = 14 πb b , I = 14 πb b , I = 0 , where b , b are the half-lengths of the major and minor axes. Substitutingthese formulas into (34) we get finally C = πµ ∗ b b b + b = 4 µ ∗ ( I − ) αα , where ( I − ) αβ is a tensor inverse to I αβ .It is easy to see that (34) is also the solution to the problem (33) for ahollow elliptical cross section λ ≤ c αβ y α y β ≤ , λ < . Similar calculations give C = πµ ∗ (1 − λ ) b b b + b . For the rectangular cross section of width a and height 1, where a < C = µ ∗ ch a , where c = 13 − aπ ∞ (cid:88) n =0 tanh( π (2 n + 1) / a )(2 n + 1) . For other cross sections the problem (33) can be solved by the finite elementmethod. 20 . 1-D beam theory
Summarizing the results obtained in Sections 3-6, we can now reduce the3-D problem of equilibrium of the FRC beam to the following 1-D variationalproblem: Minimize the energy functional J [ v , ϕ ] = (cid:90) L Φ( γ, Ω α , Ω) d x − f · v ( L ) + m α v α,x ( L ) − mϕ ( L )among functions v and ϕ satisfying the kinematic boundary conditions v (0) = 0 , v α,x (0) = 0 , ϕ (0) = 0 . In this functional the stored energy density Φ( γ, Ω α , Ω) is given byΦ( γ, Ω α , Ω) = 12 ( ¯ E + H ∗ )( |A| γ + I αβ Ω α Ω β ) + 12 C Ω , while the resultant forces and moments acting at x = L are equal to f = (cid:90) A t ( x , x ) d x ,m α = (cid:90) A x α t ( x , x ) d x , m = (cid:90) A e αβ x α t β ( x , x ) d x . The standard calculus of variations shows that v and ϕ satisfy the equi-librium equations T ,x = 0 , M α,xx = 0 , M ,x = 0 , (35)where T = ∂ Φ ∂γ = ( ¯ E + H ∗ ) |A| γ, M α = ∂ Φ ∂ Ω α = ( ¯ E + H ∗ ) I αβ Ω β , M = ∂ Φ ∂ Ω = C Ω . Besides, the following boundary conditions must be fulfilled at x = L : T ( L ) = f, M α ( L ) = m α , M α,x = f α , M = m. (36)Using the technique of Gamma-convergence (8; 22), one can prove that thesolution of (35)-(36) converges to the minimizer of (1) in the energetic normas h/L → (cid:15) = c/h →
0. 21t is easy to extend this one-dimensional theory to the dynamics of FRCbeam, where functions v and ϕ depend on x and t . One need just to includeinto the one-dimensional functional the kinetic energy density which, afterthe dimension reduction and homogenization in accordance with (4), takesthe form Θ = 12 ¯ ρ ( |A| ˙ v · ˙ v + I αα ˙ ϕ ) , (37)with ¯ ρ = (cid:104)(cid:104) ρ (cid:105)(cid:105) = s ρ + s ρ being the mass density averaged over the unitcell. Hamilton’s variational principle for 1-D beam theory states that, amongall admissible functions v ( x, t ) and ϕ ( x, t ) satisfying the initial and end con-ditions as well as the kinematic boundary conditions, the true displacementˇ v and rotation ˇ ϕ are the extremal of the action functional J [ v , ϕ ] = (cid:90) t t (cid:90) L (Θ − Φ) d x d t + (cid:90) t t ( f · v ( L ) − m α v α,x ( L ) + mϕ ( L )) d t . The Euler equations become¯ ρ |A| ¨ v = T ,x , ¯ ρ |A| ¨ v α = M α,xx , ¯ ρI αα ¨ ϕ = M ,x . These equations are subjected to the boundary conditions (37) and the initialconditions. Similar to the static case, it can be proved that the solution ofthis 1-D dynamic theory converges to the solution of the 3-D theory in thelimits h/L → (cid:15) = c/h →
8. Conclusion
It is shown in this paper that the rigorous first order approximate 1-Dtheory of thin FRC beams can be derived from the exact 3-D elasticity theoryby the variational-asymptotic method. The developed finite element code canbe used to solve cross-sectional problems with arbitrary elastic moduli of thefibers and matrix as well as arbitrary distributions and shapes of fibers. Theextension of this multi-scaled asymptotic analysis to curved and naturallytwisted FRC beams is straightforward (14). As seen from (14), the extension,bending, and torsion modes of beams are coupled in that case. The VAMcombined with the finite element cross-sectional analysis can also be appliedto the nonlinear FRC beams with periodic microstructure in the spirit of(28; 29; 31; 20). Finally, let us mention the FRC beams with randomlydistributed fibers under torsion analyzed in (2). The analysis of extensionand bending of such beams requires the solution of the plane strain problemwhich is quite challenging. 22 ppendix: Matlab-code function [CH ,DH ,FH] = homogenizecs (lx , ly , lambda , mu , phi , x) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % lx = Unit cell length in x- direction . % ly = Unit cell length in y- direction . % lambda = Lame ’s first parameter for both materials . Two entries . % mu = Lame ’s second parameter for both materials . Two entries . % phi = Angle between horizontal and vertical cell wall . Degrees % x = Material indicator matrix . Size used to determine nelx / nely %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% INITIALIZE % Compute contraction ratios nu = [ lambda (1) /(2*( lambda (1) + mu (1) )) ... lambda (2) /(2*( lambda (2) + mu (2) ))]; % Deduce discretization [nely , nelx ] = size (x); % Stiffness matrix consists of two parts , one belonging to lambda and % one belonging to mu. Same goes for load vector dx = lx/ nelx ; dy = ly/ nely ; nel = nelx * nely ; [ keLambda , keMu , feLambda , feMu , feXiLambda , feXiMu ] = elementMatVec (dx /2 ,dy /2 , phi ); % Node numbers and element degrees of freedom for full ( not periodic ) mesh nodenrs = reshape (1:(1+ nelx ) *(1+ nely ) ,1+ nely ,1+ nelx ); edofVec = reshape (2* nodenrs (1: end -1 ,1: end -1) +1 , nel ,1) ; edofMat = repmat ( edofVec ,1 ,8) + repmat ([0 1 2* nely +[2 3 0 1] -2 -1],nel ,1) ; %% IMPOSE PERIODIC BOUNDARY CONDITIONS % Use original edofMat to index into list with the periodic dofs nn = ( nelx +1) *( nely +1) ; % Total number of nodes nnP = ( nelx )*( nely ); % Total number of unique nodes nnPArray = reshape (1: nnP , nely , nelx ); % Extend with a mirror of the top border nnPArray ( end +1 ,:) = nnPArray (1 ,:) ; % Extend with a mirror of the left border nnPArray (: , end +1) = nnPArray (: ,1) ; % Make a vector into which we can index using edofMat : dofVector = zeros (2* nn , 1); dofVector (1:2: end ) = 2* nnPArray (:) -1; dofVector (2:2: end ) = 2* nnPArray (:) ; edofMat = dofVector ( edofMat ); ndof = 2* nnP ; % Number of dofs %% ASSEMBLE STIFFNESS MATRIX % Indexing vectors iK = kron ( edofMat , ones (8 ,1) ) ’; jK = kron ( edofMat , ones (1 ,8) ) ’; % Material properties in the different elements lambda = lambda (1) *(x ==1) + lambda (2) *(x ==2) ; mu = mu (1) *(x ==1) + mu (2) *(x ==2) ; nu = nu (1) *(x ==1) + nu (2) *(x ==2) ; % The corresponding stiffness matrix entries sK = keLambda (:) * lambda (:) .’ + keMu (:) *mu (:) . ’; K = sparse (iK (:) , jK (:) , sK (:) , ndof , ndof ); %% LOAD VECTORS AND SOLUTION % Assembly three load cases corresponding to the three strain cases sF = feLambda (:) * lambda (:) . ’+ feMu (:) *mu (:) . ’; lambda_nu = lambda .* nu; mu_nu = mu .* nu; sF =[ sF; feXiLambda (:) * lambda_nu (:) .’ + feXiMu (:) * mu_nu (:) . ’]; iF = repmat ( edofMat ’ ,4 ,1); jF = [ ones (8 , nel ); 2* ones (8 , nel ); 3* ones (8 , nel ); 4* ones (8 , nel )]; F = sparse (iF (:) , jF (:) , sF (:) , ndof , 4); % Solve ( remember to constrain one node ) chi (3: ndof ,:) = K (3: ndof ,3: ndof )\F (3: ndof ,:) ; chi_4 = chi (: ,4) ; %% HOMOGENIZATION % The displacement vectors corresponding to the unit strain cases chi0 = zeros (nel , 8, 3); % The element displacements for the three unit strains chi0_e = zeros (8 , 4); ke = keMu + keLambda ; % Here the exact ratio does not matter , because fe = feMu + feLambda ; % it is reflected in the load vector fe = [fe feXiLambda + feXiMu ]; chi0_e ([3 5: end ] ,:) = ke ([3 5: end ] ,[3 5: end ])\fe ([3 5: end ] ,:); chi0_4_e =nu (:) * chi0_e (: ,4) ’; % epsilon0_11 = (1 , 0, 0) chi0 (: ,: ,1) = kron ( chi0_e (: ,1) ’, ones (nel ,1) ); % epsilon0_22 = (0 , 1, 0) chi0 (: ,: ,2) = kron ( chi0_e (: ,2) ’, ones (nel ,1) ); % epsilon0_12 = (0 , 0, 1) chi0 (: ,: ,3) = kron ( chi0_e (: ,3) ’, ones (nel ,1) ); CH = zeros (3) ; DH = [0; 0; 0]; cellVolume = lx*ly; sumLambda = (( chi0_4_e (: ,:) -chi_4 ( edofMat )) *... keLambda ) .*( chi0_4_e (: ,:) -chi_4 ( edofMat )); sumMu = (( chi0_4_e (: ,:) -chi_4 ( edofMat ))* keMu ) .*... ( chi0_4_e (: ,:) -chi_4 ( edofMat )); sumLambda = reshape ( sum ( sumLambda ,2) , nely , nelx ); sumMu = reshape ( sum ( sumMu ,2) , nely , nelx ); FH = 1/ cellVolume * sum ( sum ( lambda .* sumLambda +mu .* sumMu )); for i = 1:3 sumLambda = (( chi0_4_e (: ,:) -chi_4 ( edofMat )) *... keLambda ) .*( chi0 (: ,: ,i)-chi ( edofMat +(i -1) * ndof )); sumMu = (( chi0_4_e (: ,:) -chi_4 ( edofMat ))* keMu ) .*... ( chi0 (: ,: ,i)-chi ( edofMat +(i -1) * ndof )); sumLambda = reshape ( sum ( sumLambda ,2) , nely , nelx ); sumMu = reshape ( sum ( sumMu ,2) , nely , nelx ); DH(i) = 1/ cellVolume * sum ( sum ( lambda .* sumLambda +mu .* sumMu )); for j = 1:3 sumLambda = (( chi0 (: ,: ,i) - chi ( edofMat +(i -1) * ndof ))* keLambda ) .*... ( chi0 (: ,: ,j) - chi ( edofMat +(j -1) * ndof )); sumMu = (( chi0 (: ,: ,i) - chi ( edofMat +(i -1) * ndof ))* keMu ) .*... ( chi0 (: ,: ,j) - chi ( edofMat +(j -1) * ndof )); sumLambda = reshape ( sum ( sumLambda ,2) , nely , nelx ); sumMu = reshape ( sum ( sumMu ,2) , nely , nelx ); % Homogenized elasticity tensor CH(i,j) = 1/ cellVolume * sum ( sum ( lambda .* sumLambda + mu .* sumMu )); end end disp (’--- Homogenized elasticity CH ---’); disp (CH) disp (’--- Homogenized elasticity DH ---’); disp (DH) disp (’--- Homogenized elasticity FH ---’); disp (FH) %% COMPUTE ELEMENT STIFFNESS MATRIX AND FORCE VECTOR ( NUMERICALLY ) function [ keLambda , keMu , feLambda , feMu , feXiLambda , feXiMu ] =elementMatVec (a, b, phi ) % Constitutive matrix contributions CMu = diag ([2 2 1]) ; CLambda = zeros (3) ; CLambda (1:2 ,1:2) = 1; % Two Gauss points in both directions xx =[ -1/ sqrt (3) , 1/ sqrt (3) ]; yy = xx; ww =[1 ,1]; % Initialize keLambda = zeros (8 ,8) ; keMu = zeros (8 ,8) ; feLambda = zeros (8 ,3) ; feMu = zeros (8 ,3) ; feXiLambda = zeros (8 ,1) ; feXiMu =zeros (8 ,1) ;
L = zeros (3 ,4) ; L (1 ,1) = 1; L (2 ,4) = 1; L (3 ,2:3) = 1; for ii =1: length (xx) for jj =1: length (yy) % Integration point x = xx(ii); y = yy(jj); % Differentiated shape functions dNx = 1/4*[ -(1 - y) (1 -y) (1+ y) -(1+y)]; dNy = 1/4*[ -(1 - x) -(1+x) (1+ x) (1 -x)]; % Jacobian
J = [ dNx ; dNy ]*[ -a a a +2* b/ tan ( phi *pi /180) 2*b/ tan ( phi *pi /180) -a; ... -b -b b b] ’; detJ = J (1 ,1) *J (2 ,2) - J (1 ,2) *J (2 ,1) ; invJ = 1/ detJ *[J (2 ,2) -J (1 ,2) ; -J (2 ,1) J (1 ,1) ]; % Weight factor at this point weight = ww(ii)*ww(jj)* detJ ; % Strain - displacement matrix
G = [ invJ zeros (2) ; zeros (2) invJ ]; dN = zeros (4 ,8) ; dN (1 ,1:2:8) = dNx ; dN (2 ,1:2:8) = dNy ; dN (3 ,2:2:8) = dNx ; dN (4 ,2:2:8) = dNy ;
B = L*G*dN; % Element matrices keLambda = keLambda + weight *(B’ * CLambda * B); keMu = keMu + weight *(B’ * CMu * B); % Element loads feLambda = feLambda + weight *(B’ * CLambda * diag ([1 1 1]) ); feMu = feMu + weight *(B’ * CMu * diag ([1 1 1]) ); feXiLambda = feXiLambda + weight *(B’ * CLambda * [1; 1; 0]) ; feXiMu = feXiMu + weight *(B’ * CMu * [1; 1; 0]) ; end end