# The Hamiltonian System Method for the Stress Analysis in Axisymmetric Problems of Viscoelastic Solids

Authors: W. X. Zhang; Y. Bai; F. Yuan
Journal: Journal of Applied Mathematics
 (2012)

---

## Abstract

With the use of the Laplace integral transformation and state space formalism, the classical axial symmetric quasistatic problem of viscoelastic solids is discussed. By employing the method of separation of variables, the governing equations under Hamiltonian system are established, and hence, general solutions including the zero eigensolutions and nonzero eigensolutions are obtained analytically. Due to the completeness property of the general solutions, their linear combinations can describe various boundary conditions. Simply by applying the adjoint relationships of the symplectic orthogonality, the eigensolution expansion method for boundary condition problems is given. In the numerical examples, stress distributions of a circular cylinder under the end and lateral boundary conditions are obtained. The results exhibit that stress concentrations appear due to the displacement constraints, and that the effects are seriously confined near the constraints, decreasing rapidly with the distance from the boundary.


---

## Body

## 1.  Introduction

In the modern engineering designs, the knowledge in material behavior is necessary to obtain predictive numerical simulations. This is particularly true for viscoelastic materials [1–3]. However, analytical solutions for viscoelasticity are difficult to be found due to the time-dependent property of its stress-strain relations, and hence numerical approaches are often used. Among the computational techniques, one of the most popular approaches is the finite element method, which is an appreciable and widespread tool for researches of solids and structures. The main feature of finite elements is the modularity in their usage and the easy framework for integration. Dureisseix and Bavestrello focused on solving coupled problems, like those arising from multiphysics models, where different meshes are used for different physics [4]. Bottoni et al. developed a finite element model for the time-dependent analysis of pultruded thin-walled beams subject to long-term loadings [5]. The other important tool is the boundary element method. Zhang et al. discussed crack problems in linear viscoelastic materials by generalizing the Heaviside function to represent the displacement discontinuity across the crack surface [6]. Based on differential constitutive relations, Mesquita and Coda provided the important algebraic equations and presented a method for the treatment of two dimensional coupling problems between the finite element method and the boundary element method by discussing Kelvin and Boltzmann models [7]. Compared with the finite element method, this method can reduce the dimension of the problem and provide an attractive idea for the viscoelastic research.Due to the use of single variable (displacement or stress), the above mentioned methods inevitably produce high order partial differential equations. Zhong developed the Hamiltonian system method for deriving exact analytical solutions to some basic problems in elastic mechanics [8, 9]. This method is developed on the basis of the mathematical theory on Hamiltonian geometry, by which the method of separation of variables can be applied by introducing dual variables. Xu et al. and Zhang and Xu discussed Saint-Venant problems and Saint-Venant principle of elasticity and viscoelasticity, the decay property of non-zero eigensolutions, and local effect near the boundary [10, 11]. Zhang and Deng presented the Lie group integration method for the constrained generalized Hamiltonian system [12]. Yao and Li derived the plane magnetoelectroelastic solids problem by means of the generalized variable principle [13]. Since the difficulty for solving high order differential equations in the traditional methods, such as the semi-inverse method, is overcome, the Hamiltonian system method gained much attention in recent years, and has been applied successfully in applied mechanics [14].The Hamiltonian system method is further introduced in this paper to analyse quasistatic axisymmetric viscoelastic problems. Since viscoelastic constitutive equations can be transformed into a set of corresponding elastic ones using the Laplace integral transform, the potential energy of viscoelastic body belongs to energy conservative system in the Laplace domain, and the Hamiltonian system method can be applied. Based on the adjoint symplectic relationships of the general solutions in the time domain, the eigensolution expansion method is introduced to satisfy the boundary conditions. Using this method, various boundary conditions, such as displacement conditions, stress conditions, and mixed conditions of displacement and stress, can be conveniently described by the combination of the eigensolutions. Accordingly, the assumption of stress or displacement, usually used in the semi-inverse method, is not needed. In numerical example, stress distributions of the boundary condition problems of a circular cylinder are discussed. The results show that stress concentrations occur in the region near the boundary when the cylinder is subjected to displacement constraint. It is verified that the local effect plays an important role in the analysis of stresses. However, in the region far from the boundary, the effect vanished and usually can be neglected.


## 2. The Governing Equations in the Hamiltonian System

An isotropic viscoelastic circular cylinder, whose radius and length areR and l, respectively, is considered. As Figure 1 shows, the origin of the cylindrical polar coordinates (r,θ,z) is located at the center of the bottom end, with the z-axis pointing the top end along the axial direction. For axisymmetric problems, the non-zero components are stresses σr, σθ, σz, and τrz and displacements u and w along r and z coordinates, respectively. The constitutive relations of viscoelastic media can be expressed in an integral form:
(2.1)σm(w,t)=3∫0tK(t-τ)dεm(w,τ)dτdτ,sij(w,t)=2∫0tG(t-τ)deij(w,τ)dτdτ,
where w is a position vector; K(t) and G(t) are the relaxations of bulk modulus and shear modulus; σm and εm are the volumetric stress and strain, respectively, and the deviatoric components of the stress and strain tensors sij=σij-σmδij/3,eij=εij-εmδij/3. The viscoelastic constitutive equations can be transformed into a set of corresponding elastic governing equations using the Laplace transform. As a result, the transformed constitutive equation (2.1) in the Laplace domain is written as
(2.2)σ-m(w,s)=3Kε-m(w,s),s-ij(w,s)=2G-(s)e-ij(w,s).
Take the Burgers model, for example, which consists of a Kelvin model and a Maxwell model connected in series. As Figure 2 shows, the parameters G1 and G2 denote the spring constants, and η1 and η2 denote the viscosity coefficients of the dashpots. In general, the Poisson’s ratio ν is a function of time, but the time dependence is much weaker than other modului. Here we assume that the Poisson’s ratio is a constant, independent of time. Thus, the relaxation of shear modulus and the Young’s modulus can be written as
(2.3)G-(s)=a1s+a2s22(1+b1s+b2s2),E-(s)=2(1+ν)G-(s),
where a1=η1, a2=η1η2/G2, b1=η1/G1+η1/G2+η2/G1, and b2=η1η2/(G1G2), s is the Laplace transform parameter, and a bar over a modulus denotes its inverse Laplace transform. In order to derive the final governing equations of the Hamiltonian system, the Lagrange function of strain energy density in the Laplace domain is introduced as
(2.4)L-=G-νr(1-2ν)(∂u-∂r+u-r+w-˙)2+G-r[(∂u-∂r)2+u-2r2+w-˙2+12(∂w-∂r+u-˙)2],
in which an overdot on a variable denotes its partial derivative with respect to coordinate z. Write the displacement variables in vector form
(2.5)Q-=[u-,w-]T.
Then the dual vector can be obtained as
(2.6)P-=[∂L-∂u-˙,∂L-∂w-˙]T=[P1P2]=[rτ-rzrσ-z].
The Hamiltonian function can be introduced as
(2.7)H(Q-,P-)=P-TQ--L-(Q-,P-).
Applying the variational method
(2.8)δ∫0R∫0l[H(Q-,P-)-P-TQ-]drdz=0,
we get the governing equations under the Hamiltonian system
(2.9)ψ-˙=Hψ-,
where ψ-=[u-,w-,rτ-rz,rσ-z]T, and the Hamiltonian operator matrix
(2.10)H=[0-α11G-r0-ν1-να200(1-2ν)2G-r(1-ν)2G-1-να300-ν1-να400-α10],
in which the partial differential operators α1=r∂/∂r, α2=r∂/∂r+1, α3=-r2∂2/∂r2-r∂/∂r+1, and α4=r∂/∂r-1. The other two non-zero stress components can be obtained accordingly as
(2.11)[rσ-rrσ-θ]=11-ν[νrσ-z+2G-r∂u-∂r+2G-νu-νrσ-z+2G-νr∂u-∂r+2G-u-].Figure 1
The coordinate system of the circular cylinder.Figure 2
The Burgers viscoelastic model.


## 3. General Solutions

To derive the general solutions of the governing equations (2.9), we suppose the lateral boundary of the cylinder is stress free. The condition for this case is
(3.1)τrz=σr=0(r=R).
Using the method of separation of variables, the solution is written as  ψ-(r,z)=ω-(z)X-(r). Considering (2.9), we get  ω-(z)=eμz and the eigenequation
(3.2)HX-(r)=μX-(r),
where μ is an eigenvalue, and X- is its corresponding eigenvector: X-=[q-p-]=[q-1,q-2,p-1,p-2]T. The displacement components for the eigenequation (3.2) can be expressed as
(3.3)[q-1q-2]=[4(1-ν)ξr-ddr(ξ0+rξr)-μ(ξ0+rξr)],
in which the functions ξr and ξ0 satisfy
(3.4)[rd2dr+ddr+rκ2rd2dr+d2dr+rκ2-1r][ξ0ξr]=0.
The solutions are
(3.5)ξ0=β1J0,ξr=β2J1,
where J0 and J1 are Bessel functions: J0=J0(μr), J1=J1(μr), and β1 and β2 are integral constants. Substituting (3.3) into (3.2), we get
(3.6)q-1=h1μJ1β1+h1[4(1-ν)J1-μrJ0]β2,q-2=h1μJ0β1+h1μrJ1β2,p-1=h2μrJ1β2,p-2=h3μ2rJ0β1+h3μr(μrJ1+2νJ0-4ν2J0)β2,
in which h1=1/(4E--4E-ν),h2=1/(2+2ν),andh3=1/[4(1-ν2)(1-2ν)]. Using the lateral boundary condition (3.1), we get
(3.7)[βij]β=0(i,j=1,2),
in which β=[β1β2], and the components are β11=(2ν-1)J1(μR)/R+μJ0(μR), β21=0, β12=4(2ν-1)(1-ν)J1(μR)/(μR)-ν[J0(μR)-μRJ1(μR)]+4(1-ν)2J0(μR), and β22=J1(μR). To ensure the existence of non-zero eigensolutions, the integral constants β1 and β2 cannot be zeros, simultaneously. Thus, the equation about the eigenvalues can be constructed as
(3.8)|βij|=0.For the case ofμ=0, by solving (3.2) directly, we can get the fundamental zero eigensolution:
(3.9)X-0=[q-1=0,q-2=1,p-1=0,p-2=0]T.
Besides solution (3.9), there is a corresponding solution of Jordan form. The Jordan form solution should satisfy
(3.10)HX--0=X-0.
The solution is
(3.11)X--0=[q-1=-νr,q-2=0,p-1=0,p-2=E-r]T.
It should be mentioned that (3.10) and (3.11) are not solutions of the governing equations (2.9). However, we can get the general solutions from these zero eigensolutions. The corresponding solutions are
(3.12)ψ-0=X-0=[0,1,0,0]T,ψ--0=X--0+zX-0=[-νr,z,0,E-r]T,
which are solutions of the rigid translation and the simple extension along z direction, respectively. Noticing that both the zero and non-zero eigensolutions mentioned above are expressed in concise analytical forms, we can discuss the problem in the time domain directly simply by using the inverse Laplace transform.


## 4. Boundary Conditions

It is well known that boundary conditions can be displacement conditions or stress conditions, and they also can be the mixed conditions of displacements and stresses. In the Hamiltonian system, the boundary conditions just correspond to the fundamental variables (displacements and stresses). Therefore, it is very convenient to discuss boundary condition problems.


### 4.1. Lateral Boundary Conditions

The above discussion of general solutions is in condition that the lateral boundary conditions (r=R) are homogeneous. For nonhomogeneous case, these eigensolutions can not be used directly. However, simply by employing the variable substitution method, we can transform the nonhomogeneous lateral boundary conditions into the homogeneous ones. As an example, we discuss the stress lateral boundary conditions:
(4.1)σ-r=σ-0,τ-rz=τ-0(r=R).
Equation (4.1) can be described as
(4.2)E-1-ν2∂u-∂r+E-ν1-ν2u-r+ν1-νP-2r=σ-0(r=R),P-1=Rτ-0(r=R).
To homogenize boundary conditions (4.2), we introduce new variables
(4.3)ψ-*=[u-*,w*,P-1*,P-2*]T=ψ--φ-,
where
(4.4)φ-=[(1-ν)σ-0rE-0r2τ-0R0].
It can be verified that ψ-* satisfy homogeneous condition
(4.5)E-1-ν2∂u-*∂r+E-ν1-ν2u-*r+ν1-νP-2*r=0(r=R),P-2*=0(r=R).
However, the governing equation (2.9) is changed into nonhomogeneous one
(4.6)ψ-˙*=Hψ-*+g-*,
where
(4.7)g-*=Hφ--φ-˙.
The components are g-1*=2(1+ν)rτ-0/(E-R)-(1-ν)σ-˙0r/E-, g-2*=-2v-σ-0/E-, g-3*=(ν-1)σ-˙0r/E-, and g-4*=(ν-1)σ-0/E-. Thus, the problem of nonhomogeneous lateral boundary conditions is transformed into finding a particular solution of (4.6). To derive a particular solution, we define the symplectic integral product based on the property of the Hamiltonian operator matrix H:
(4.8)〈Xi,Xj〉=∫0RXiTJXjdr,
where Xi and Xj are arbitrary eigensolutions of the time domain, and J is a unit rotational matrix [11]. For the convenience of discussion, we redescribe the eigensolutions as
(4.9)Xm=X|μ=μm,X-m=ηmX|μ=-μm,
in which 1/ηm=〈X|μ=μm,X|μ=-μm〉. Thus, the eigensolutions satisfy the adjoint symplectic relationships
(4.10)〈Xi,X-j〉=-〈X-i,Xj〉=1(j=0,1,2,…),〈Xi,Xj〉=0(i,j=0,±1,±2,…;i≠-j).
Since the eigensolutions constructed a complete solution space, the nonhomogeneous term of (4.6) can be developed as
(4.11)g*=∑j=0∞[bj(z)Xj(r)+b-j(z)X-j(r)],
where bj(z)=〈ψ,X-j〉 and b-j(z)=-〈ψ,Xj〉. Suppose that ψp is a particular solution
(4.12)ψp=∑j=0∞[cj(z)Xj(r)+c-j(z)X-j(r)].
Therefore,
(4.13)c˙j(z)=μjcj(z)+bj(z),c˙-j(z)=-μjc-j(z)+b-j(z).
The solutions are
(4.14)cj(z)=∫0zbj(t)eμj(z-t)dtc-j(z)=∫0zb-j(t)e-μj(z-t)dt.


### 4.2. End Conditions

To explain the eigensolution expansion method for the satisfaction of end conditions, we suppose the end conditions:(4.15)u=u~,w=w~,(z=0)τrz=τ~rz,σz=σ~z,(z=l).
Equation (4.15) can be expressed as
(4.16)Pz=l(r)=Pl,Qz=0(r)=Q0.
The complete solution of the problem should be the linear combination of the eigensolutions and a particular solution
(4.17)ψ=∑n=0∞(aneμnzXj+bne-μnzX-n)+ψp,
where ψp is the particular solution of the problem given in (4.12), and an and bn are coefficients to be determined. Using the solution (4.17), the boundary conditions (4.16) can be rewritten as
(4.18)Pl(r)∑nanpn(r)eμnl+∑nbnp-n(r)e-μnl+Pp(z=l),Q0(r)∑nanqn(r)+∑nbnq-n(r)+Qp(z=0).
or
(4.19)ψ|z=0[Q0(r)∑nanpn+∑nbnp-n+Pp]ψ|z=l[∑nanqneμnl+∑nbnq-ne-μnl+QpPl(r)].
Here Qp and Pp are components of the particular solution ψp. Based on the adjoint symplectic relationships (4.10), an infinite set of equations about the coefficients an and bn can be obtained
(4.20)∫0Rq-j·(Pl-Pp)dr=∫0R∑n=0∞[aneμnlq-j·pn+bn-μnlq-j·p-n]dr(z=l),∫0R(Q0-Qp)·pjdr=∫0R∑n=0∞[anqn·pj+bnq-n·pj]dr(z=0),
where j=0,1,2,…. In the numerical calculations, we usually take the first N terms in (4.20). Thus, there are 2N undetermined coefficients and 2N algebra equations. So the combination of the eigensolution expansion is solely determined.


## 4.1. Lateral Boundary Conditions

The above discussion of general solutions is in condition that the lateral boundary conditions (r=R) are homogeneous. For nonhomogeneous case, these eigensolutions can not be used directly. However, simply by employing the variable substitution method, we can transform the nonhomogeneous lateral boundary conditions into the homogeneous ones. As an example, we discuss the stress lateral boundary conditions:
(4.1)σ-r=σ-0,τ-rz=τ-0(r=R).
Equation (4.1) can be described as
(4.2)E-1-ν2∂u-∂r+E-ν1-ν2u-r+ν1-νP-2r=σ-0(r=R),P-1=Rτ-0(r=R).
To homogenize boundary conditions (4.2), we introduce new variables
(4.3)ψ-*=[u-*,w*,P-1*,P-2*]T=ψ--φ-,
where
(4.4)φ-=[(1-ν)σ-0rE-0r2τ-0R0].
It can be verified that ψ-* satisfy homogeneous condition
(4.5)E-1-ν2∂u-*∂r+E-ν1-ν2u-*r+ν1-νP-2*r=0(r=R),P-2*=0(r=R).
However, the governing equation (2.9) is changed into nonhomogeneous one
(4.6)ψ-˙*=Hψ-*+g-*,
where
(4.7)g-*=Hφ--φ-˙.
The components are g-1*=2(1+ν)rτ-0/(E-R)-(1-ν)σ-˙0r/E-, g-2*=-2v-σ-0/E-, g-3*=(ν-1)σ-˙0r/E-, and g-4*=(ν-1)σ-0/E-. Thus, the problem of nonhomogeneous lateral boundary conditions is transformed into finding a particular solution of (4.6). To derive a particular solution, we define the symplectic integral product based on the property of the Hamiltonian operator matrix H:
(4.8)〈Xi,Xj〉=∫0RXiTJXjdr,
where Xi and Xj are arbitrary eigensolutions of the time domain, and J is a unit rotational matrix [11]. For the convenience of discussion, we redescribe the eigensolutions as
(4.9)Xm=X|μ=μm,X-m=ηmX|μ=-μm,
in which 1/ηm=〈X|μ=μm,X|μ=-μm〉. Thus, the eigensolutions satisfy the adjoint symplectic relationships
(4.10)〈Xi,X-j〉=-〈X-i,Xj〉=1(j=0,1,2,…),〈Xi,Xj〉=0(i,j=0,±1,±2,…;i≠-j).
Since the eigensolutions constructed a complete solution space, the nonhomogeneous term of (4.6) can be developed as
(4.11)g*=∑j=0∞[bj(z)Xj(r)+b-j(z)X-j(r)],
where bj(z)=〈ψ,X-j〉 and b-j(z)=-〈ψ,Xj〉. Suppose that ψp is a particular solution
(4.12)ψp=∑j=0∞[cj(z)Xj(r)+c-j(z)X-j(r)].
Therefore,
(4.13)c˙j(z)=μjcj(z)+bj(z),c˙-j(z)=-μjc-j(z)+b-j(z).
The solutions are
(4.14)cj(z)=∫0zbj(t)eμj(z-t)dtc-j(z)=∫0zb-j(t)e-μj(z-t)dt.


## 4.2. End Conditions

To explain the eigensolution expansion method for the satisfaction of end conditions, we suppose the end conditions:(4.15)u=u~,w=w~,(z=0)τrz=τ~rz,σz=σ~z,(z=l).
Equation (4.15) can be expressed as
(4.16)Pz=l(r)=Pl,Qz=0(r)=Q0.
The complete solution of the problem should be the linear combination of the eigensolutions and a particular solution
(4.17)ψ=∑n=0∞(aneμnzXj+bne-μnzX-n)+ψp,
where ψp is the particular solution of the problem given in (4.12), and an and bn are coefficients to be determined. Using the solution (4.17), the boundary conditions (4.16) can be rewritten as
(4.18)Pl(r)∑nanpn(r)eμnl+∑nbnp-n(r)e-μnl+Pp(z=l),Q0(r)∑nanqn(r)+∑nbnq-n(r)+Qp(z=0).
or
(4.19)ψ|z=0[Q0(r)∑nanpn+∑nbnp-n+Pp]ψ|z=l[∑nanqneμnl+∑nbnq-ne-μnl+QpPl(r)].
Here Qp and Pp are components of the particular solution ψp. Based on the adjoint symplectic relationships (4.10), an infinite set of equations about the coefficients an and bn can be obtained
(4.20)∫0Rq-j·(Pl-Pp)dr=∫0R∑n=0∞[aneμnlq-j·pn+bn-μnlq-j·p-n]dr(z=l),∫0R(Q0-Qp)·pjdr=∫0R∑n=0∞[anqn·pj+bnq-n·pj]dr(z=0),
where j=0,1,2,…. In the numerical calculations, we usually take the first N terms in (4.20). Thus, there are 2N undetermined coefficients and 2N algebra equations. So the combination of the eigensolution expansion is solely determined.


## 5. Numerical Examples

In order to verify the proposed approach, numerical investigations are carried out in this section. In the numerical calculations, the geometrical data and the parameter are selected asl/R=3, G1=G2=G, and η1=η2=η, and the Poisson’s ratio ν is taken to be 0.25.We, firstly, discuss the end condition problem. The computational model is supposed to be as follows: the top (z=l) end of the cylinder is subjected to stress condition, while the bottom end (z=0) is clamped. The boundary conditions for this case can be expressed as follows.End conditions:(5.1)τrz=0,σz=σ0(z=l),u=0,w=0(z=0).Lateral conditions:(5.2)σr=0,τrz=0(r=R).
Based on (4.20), algebra equations about the coefficients in the eigensolution expansion are established, and hence the distributions of stresses are obtained numerically. According to the results, the effects of σr and τrz are much weaker than the other stresses. Thus, only σz and σθ are considered here. As Figures 3 and 4 show, stress concentrations appear clearly near the clamped end of the domain, especially in the region around the original point. It is obvious that the effects decrease rapidly, and the stresses tend to be constant with the distance from the bottom end. This character indicates that zeroeigensolutions can approximately describe solutions near the top end. However, the solution is not accurate enough near the fixed end where local effects appear, and nonzero-eigenvalue eigensolutions are required to describe the stress concentrations.Figure 3
Distribution of the stressσz/(σ0) in the end condition example.Figure 4
Distribution of the stressσθ/(σ0) in the end condition example.Now, consider secondly a lateral condition problem. The boundaries are given as follows.Lateral conditions:(5.3)σr=-σ0,τrz=0(r=R).End conditions:(5.4)τrz=0,σz=0(z=l),u=0,w=0(z=0).
Based on the technique of lateral boundary conditions discussed in Section 4.1, we obtained the numerical result. Among the stresses, σr and σθ have more important effect during the compression for this case. Figures 5 and 6 show exhibit that stress concentrations appear near the bottom end similarly, and the local effects decay rapidly with coordinate z. The results of local effects can be well explained by the famous Saint-Venant principle. According to Saint-Venant principle, local effects of stresses and displacements must appear near the region, where the external displacement or force boundary conditions are given.Figure 5
Distribution of the stressσr/(σ0) in the lateral condition example.Figure 6
Distribution of the stressσθ/(σ0) in the lateral condition example.


## 6.  Conclusion

Because of the existence of the energy nonconservation, the Hamiltonian system method can not be applied directly for viscoelasticity. However, the potential energy of viscoelastic body has conservative form in the Laplace domain, and the problem can be transformed into problem of conservative system in which the formulation is applicable. The Hamiltonian system is a direct method by which the order of differential governing equations can be reduced. With the direct method, all the general solutions of the governing solutions, including zero eigensolutions and non-zero eigensolutions, are obtained analytically. In fact, the solution of the problem should be composed of zero eigensolutions and non-zero eigensolutions. In the Hamiltonian system, the boundary conditions just correspond to the fundamental variables, and therefore it is very convenient to discuss nonhomogeneous problems and boundary condition problems.



---

*Source: 945238-2012-10-02.xml*