Sliding frame-solid interaction using BEM/FEM coupling

^IFederal University of Sergipe, Department of Civil Engineering, Av. Marechal Rondon, s/n, 49100]000 Sao Cristovão-SE, Brazil. University of São Paulo, School of Engineering of São Carlos, Department of Structural Engineering, Av. Trabalhador SãoCarlense, 400, 13566-590 São Carlos-SP, Brazil. [e-mail]

^IIUniversity of São Paulo, School of Engineering of São Carlos, Department of Structural Engineering, Av. Trabalhador SãoCarlense, 400, 13566-590 São Carlos-SP, Brazil. [e-mail]

^IIIUniversity of São Paulo, School of Engineering of São Carlos, Department of Structural Engineering, Av. Trabalhador SãoCarlense, 400, 13566-590 São Carlos-SP, Brazil. [e-mail]

ABSTRACT

In this work it is presented a coupling between the Boundary Element Method (BEM) and the Finite Element Method (FEM) for two-dimensional elastostatic analysis of frame-solid interaction. The BEM is used to model the matrix while the reinforcement is modeled by the FEM. Regarding the coupling formulation a third degree polynomial is adopted to describe the displacement and rotations of the reinforcement, while a linear polynomial is used to describe the contact force among the domain (matrix) and the reinforcement. Perfect bounding contact forces are improved by means of redundant equations and Least squares method. Slip-bounding with two and three paths written as function of relative displacement are used to calculate the transmitted contact forces. Examples are used to demonstrate that the proposed slip-bounding procedure regularizes the contact force behavior.

Keywords: Boundary Element Method, Finite Element Method, BEM/FEM coupling, adherence models.

1 INTRODUCTION

The combination of the finite element method (FEM) and boundary element method (BEM) to solve structural analysis is attractive because it allows for an optimal exploitation of the respective advantage of the methods (Zienkiewicz et al., 1977). The main strength of the BEM for boundary-value problems governed by linear, homogenous, and elliptic differential equations with constant coefficients is the reduction of the dimensionality of the problem by one unit for linear constitutive relations (Brebbia, 1978; Brebbia, 1980). Particularly, BEM is useful to model special situations such as very large or unbounded domains, geometrical singularities (e. g. cracks) or to obtain very accurate results in regions of complicated shape (Aliabadi, 1997; Bonnet, 1999; Frangi et al., 2002). Thus, coupling the BEM and the FEM allows exploiting their complementary advantages. By the other hand, the FEM is appropriate to solve a lot of problems, including e. g. those with heterogeneous or non-linear constitutive properties, or finite deformations.

The idea of combining these two methods goes back to Zienkiewicz et al. (1977). One branch of BEM/FEM coupling is the iterative coupling in which the individual sub-domains are treated independently by either method. The procedure starts with an initial guess of the interface unknowns that will be improved by solving each sub-domain and returned to interface. This procedure repeats until an error tolerance is achieved. Although this iterative coupling is very attractive to software design, its convergence commonly depends on relaxation parameters which are rather empirical (Estorff and Hagen, 2005). For this reason, a direct coupling approach is adopted in this work.

Standard BEM formulations to deal with solids stiffened by bars or fibers are derived by combining the BEM and FEM algebraic equations. The domain (continuum or matrix) is analyzed by BEM, while finite elements are used to represent inclusions (bars for instance). The coupling is always done by enforcing displacement compatibility and traction equilibrium at interface nodes. Practical applications for non-slipping or slipping coupling using BEM/FEM coupling are presented, for instance, in Beer and Watson (1996), Coda and Venturini (1995), Coda (2001), Leite et al. (2002), Botta and Venturini (2005, 2003), Leonel (2009), Rocha (2010).

As contribution, the present work is able to capture bending effects in the analysis of frame-solid interaction, which includes bending stiffness in reinforcements immersed in 2D continuum. For this purpose, in this study it is adopted a Rissner-Mindlin type frame finite element with third order of approximation for both displacement and rotations resulting in four nodes and 12 degrees of freedom. However, the contact force is modeled by linear approximation resulting in only 4 independent values which leads to a not square force matrix. These approximations (displacements and tractions) were used because it is the lowest order that satisfies the differential equation governing the problem. In this approach, the boundary element force lines are build in a compatible way, that is there are 4 source points generating a third order approximation for displacements, but the contact force approximation is linear, also resulting in a not square force matrix.

The least Squares technique is used to eliminate the dependent equations due to the above mentioned difference in approximation order for displacement and tractions. Moreover some authors, Botta and Venturini (2005, 2003) and Leonel (2009), claim that this procedure reduces contact force oscillations.

This paper is organized as follows. It is presented in section 2 the FEM formulation to model the frame structure, which is shown the kinematics adopted. In section 3 is shown the BEM formulation adapted to domain modeling. In section 4 it is presented the proposed coupling formulation between BEM and FEM considering both perfect bonding and debonding cases. This section is divided in subsection 4.1 in which there are presented the basic equations to perfect bonding and an example to verify this formulation. In sub-section 4.3 it is presented the coupling formulation considering the slip between reinforcement and domain. Debonding models, basic equations and the non-linear formulation to solve slipping are presented in sub-sections 4.3.1, 4.3.2 and 4.3.3, respectively. The sub-section 4.4 presents two examples. The first simulates a pullout test and the second solve a soil-structure interaction case. Finally, in section 5 the final remarks and conclusions are given.

2 THE FRAME ELMENT MODELING - FEM FORMULATION

As mentioned before, the FEM is used to model frame elements and structures. Here, elements which have three degrees of freedom per node and cubic approximation for displacement and rotation are employed. This way, the elements have four nodes and these nodes present two translations (vertical and horizontal) and one rotation. Moreover, the distributed applied forces will follow linear approximation.

2.1 Kinematics

For any point on the structure, the horizontal and vertical components of displacements are given by, respectively:

where and are the reference system in the center of the layer, as shown in figure 1 and and are the displacements of point .

From equation (1), one can apply the differential operator to obtain the linear strain components:

Applying the constitutive law for the isotropic materials, the stress components at the point "p" are obtained:

where E and G are the longitudinal and shear elastic moduli, respectively.

To write the equilibrium equation it is used the Principle of Minimum Total Potential Energy. Using equations (2) and (3) one writes the Total Potential Energy equation as,

where and are the internal and potential energy of external forces, and are the components of the distributed loading (contact tractions) applied to the structure, and are the length and cross sectional area of the frame element, respectively. For approximate unknowns (, e ) cubic independent approaches were used, as shown:

with , e being the nodal values (unknown) . Since and are shape functions:

Minimizing the energy functional, equation (4), one finds the algebraic equilibrium system given by:

where: is stiffness matrix, is equivalence force matrix, unknown vector with displacements and rotations, is the vector containing the nodal values of the distributed load, is the concentrated force. Labels and are node numbers for the displacements (and rotations), four per element, and node numbers for forces (two per element), respectively

3 THE DOMAIN MODELING - BEM FORMULATION

Let us consider the domain and its boundary . For an elastic body defined by the domain , the equilibrium equation, written in terms of displacements, is given by:

where represents the displacement vector, is the shear modulus and is the Poisson's ratio.

For a domain with boundary , the integral representation of displacements is derived by applying reciprocity theorem (or Green's second identity).

Where the symbols is used to indicate fundamental solution (equation 10) and represent boundary traction values.

where with being the distance between the source and field points, , are components of the vector , is the boundary normal unit vector and is the second order identity tensor (or Kronecker delta).

For numerical solution, and are approximate by polynomial functions over boundary elements, in this work linear polynomials are used for boundary elements (internal force lines are further discussed), and the integral equation (9) is converted into an equivalent algebraic system as follows:

where matrix is obtained from the left terms in equation (9) and matrix from the terms on the right side. is a vector which contains the nodal values of displacements for all boundary nodes and is the nodal traction vector.

After substitution of the prescribed boundary conditions, the algebraic equations may be written as

The vector contains all the unknown boundary displacements and tractions, is a coefficient matrix which is usually non-symmetric and densely populated, and is a matrix which contains the coefficients corresponding to the prescribed boundary conditions .

One can differentiate equation (9) to derive the integral representation of strains and then apply the Hooke's law to obtain the stress integral equation, written for internal points, as follows,

Where and are well known third-order tensor for the stress equation obtained by applying the Hooke's law on the fundamental solution at the source point (Brebbia and Domingues, 1992).

4 BEM/FEM COUPLING FORMULATION

In this section, additional terms inserted in the classical formulations of BEM and FEM as well as the coupling between BEM and FEM are shown.

4.1 Basic equations - perfect bonding

In this subsection the perfect bounding situation is described. It implies a direct compatibility between displacements and contact force equilibrium (or continuity). Thus:

Where and are vectors containing nodal displacements for frame element and domain, respectively; and are nodal distributed force vector applied on the frame finite elements and on the force line boundary element (in the 2D domain), respectively. Once vector contains three components (two translations and one rotation) and vector contains two components (two translations), the correlation between these two vectors is done by matrix.

The components of the vectors , , and are shown in figure 2.

For 3D solids the unknown forces, , would be distributed over internal surfaces. For 2D solids this forces appears distributed along internal lines, which, roughly, work as internal "boundaries" dedicated only for tractions. Thus, the integral equation (9) is modified to include the additional term:

where is the internal force acting along the interface, , between the two materials and represents the fiber effect applied in the domain. Similarly, the integral equation (13) including this additional effect is written as:

Selecting a proper number of collocation points (source points) at the boundary and at the frame element, two sets of algebraic equations are, respectively, written from equation (16), as follows:

where index denotes boundary and index denotes frame, respectively. Additionally, the first index denotes source points and the second index denotes field points, which, are on boundary or load line.

The superimposed frame element (reinforcement) and load line representations can be seen in figure 3 in which elements (finite and line) are employed. The displacement nodes for each finite element are represented by squares and crosses represent the collocation points of equation (19) where displacements are calculated. It is observed that for fiber ends BEM equations are not written in the same position of the finite element nodes to avoid singularities; however displacements are extrapolated to the nodal positions in order to make them compatible. Moreover, this modeling allows the reinforcement (frame) ends to reach the body boundary without interfering in basic BEM equations.

Therefore, the displacement compatibility is rewritten as:

Where is the matrix that relates the nodal positions and direction of the vector with the vector. Furthermore, the matrix has dimension columns by rows, for which is the amount of BEM interface points.

According to figure 3, the amount of finite element nodes is equal to the amount of boundary interface nodes, i. e. , therefore the determination of the coupling parameters and the boundary values can be summarized by the following system of equations:

The values of and values of are known on and respectively, hence there are only unknowns in the system of equations (21). As usual, to introduce these boundary conditions into (21) one has to rearrange the system by moving columns of , with , from one side to the other, respectively. Once all unknowns are passed to the left-hand side and applying conditions of equations (14) and (20) one can write the new system of equations as:

Or in matrix form,

As described above, at the fiber element level, the number of displacement values is larger than the number of bonding (or contact) force nodal values. This occurs because in equation (19) the number of algebraic relations is much larger than the number of force values in . To reduce the number of equations to the same as the number of unknowns one can apply the Least Square Method (LSM). In this work the LSM is applied over equation (21) or (19), as follows:

Or,

Where the matrix defined for each fiber and is the transpose matrix of , i.e. .

Therefore, the system of equations (23) turns into a square system as:

4.2 Numerical example - perfect bonding

In this section, one numerical example is analyzed to check the performance and accuracy of the proposed BEM/FEM coupling for two-dimensional reinforced solids.

The reinforced simple supported beam subjected to homogeneous transversal surface load (), shown in figure 4, is analyzed. This plane structure is five meter length (), one meter height () and the position of the reinforcement is 25 centimeter () from the lower part of the matrix and is four meters long (). The following properties for domain () and reinforcement () are considered: Elasticity modulus and , Poisson's ration and , inertia moment and cross-sectional area .

The beam is discretized by 120 linear boundary elements with the same length. Three discretizations are employed to model the reinforcement, with 25, 50 and 100 finite elements. The same number of force line BEM elements is employed to model the interface. The results are compared with the commercial software ANSYS employing 100 BEAM3 elements to model the reinforcement and 2800 PLANE42 2D solid elements to model the domain. Figures 5, 6 and 7compare results (axial, transverse displacements and rotations) achieved using the BEM/FEM coupling discretizations and ANSYS.

The above results make evident that even the poorest BEM/FEM coupling discretization has good accuracy when compared with the reference result.

Regarding tractions at the interface, the BEM/FEM results can be seen in Figures 8 and 9.

According to Figures 8 and 9, for perfect bonded elastic reinforcement, it's observed that for the three BEM/FEM discretizations there is a perturbation in the traction values at the ends of the fiber. The behavior of axial forces, as shown in Figure 8, had been observed by Botta and Venturini (2005), however, the behavior shown in Figure 9 has not been reported before. Moreover, the use of LSM partially smoothes the behavior of contact force when compared with results that do not apply this strategy for the coupling, Leite et al. (2003). At this point an important result of this study should be advanced; when debondig is allowed (always occurs in a small vicinity of infinity contact stress) the above perturbation disappears.

When using nonlinear constitutive relations for matrix, there are evidences of smooth solutions for this kind of problems which were reported by Coda (2001). Moreover, one may note, in Figure 9, that the extension of the contact force perturbation reduces as discretization increases.

4.3 Coupling formulation - debonding

In this section, the previous coupling is extended to accomplish debonding. The adopted models to represent debonding as well as the nonlinear formulation to simulate the slip between domain and frame elements are presented. The feasibility of this formulation is shown through numerical examples.

4.3.1 Debonding models

Fibers embedded in the domain play an important role to improve solid stiffness and loading capacity if enough internal forces along the interface can be sustained. Sliding along the interface may be allowed when a certain amount of strength is preserved. The ideal situation for which perfect bonding is assumed, as shown in previous section, is impossible in practice; at least in the vicinity of fiber ends, as the interface forces approach the infinity (Radtke et al, 2011). Therefore, a certain amount of sliding occurs according to the bonding carrying capacity.

To model the slip that may occur in the fiber-domain interface, a debonding criterion should be considered. In this work, two models were implemented together with the proposed BEM-FEM coupling.

The curves shown in Figure 10 and 11 represent the debonding criterion that relates the bonding force with the relative displacement at interface (slip) . The following parameters define the two models: model 1 depends only on the maximum bonding force and model 2 depends on the maximum force , residual bonding force and slip characteristic values and .

From the Figure 10 and 11, the following relationships are written for the adopted models, respectively:

Model 1

Model 2

4.3.2 Basic equations

The slip consideration introduces a new variable in equations (20). This new variable represents the relative displacement, , between the domain and frame elements. The compatibility equations are now expressed by:

Where vector contains the nodal relative displacement values, relates the nodal positions with . Furthermore, matrix has dimension columns by rows. The other terms in equations (31) and (32) have already been explained, but now tractions on interface depend on relative displacements . In this work, and have been approximated by cubic polynomial and by linear polynomial.

From the introduction of relative displacements, the equilibrium equation (19) of the boundary element method may be rewritten as:

Therefore, the BEM coupling equation can be rewritten as,:

Applying boundary conditions, the equation (34) results:

If is known, the matrix form of (35) is:

By the other hand, if is known, one writes:

As can be seen, in equations (36) and (37), there are more equations than unknowns, since . Thus, to reduce the number of equations to be equal to the number of unknowns the Least Square Method (LSM) is applied over internal point equations (35). This way:

Where is the transpose matrix of .

Therefore, equations (36) and (37) become square, as follows:

Or:

4.3.3 Non-linear formulation

As one can see, equations (39) and (40) are non linear regarding slip . To solve them, one has to take into account the non-linear relationship described by the debonding model presented in item 4.3.1, in which the relation between the debonding force and the slip is established. The equilibrium equation (39) is then rewritten in terms of the variable increments, as follows:

Isolating and in first and second equation, respectively, of (41), one has:

which can be replaced in equation (41), resulting:

With

Where is the identity matrix.

Equation (44) represents a non-linear system of equations given in terms of the slip increment . It can be solved by applying the iterative Newton-Raphson scheme. Then, from the iteration the next try, , for the time increment is given by:

Linearizing equation (44) and using the first term of the Taylor's expansion, results:

The derivative that appears in equation (47) is directly obtained from equation (44) using the debonding model relationships given by equations (27)-(30). Then, one has:

The matrix , in equation (48), is the consistent tangent operator of the proposed algorithm. The derivatives on the right hand side of equation (48) depend on the updated slip value, computed appropriately according to the adopted model defined in equations (27)-(30). These derivatives are locally defined by:

To model 1

To model 2

Reaching the convergence in equation (44) for the time increment after iterations, one has to compute the slip variable to start the next increment, as follows:

After finding , other variables are directly obtained. The internal displacements, and boundary tractions and displacements are computed from equation (41). The debonding forces are computed from the constitutive relation .

4.4 Numerical example - debonding

In this item, two numerical examples are analyzed to examine the performance and accuracy of the proposed BEM/FEM combination using two models to consider the debonding between a straight bar and a two-dimensional solid.

4.4.1 Example 1

In this example the capability of the formulation to model the bonding shear contact force distribution along the bar-matrix interface during a classical pulling test is analyzed. In figure 12, a bar is partially embedded into a 2D domain and a small part to the bar is not immersed to allow applying the pulling force. The adopted geometric dimensions are , and , see figure 12. Null displacements are prescribed along the left vertical side of the two-dimensional domain. Whereas at the opposite side the load is applied by prescribing the displacement, , at the bar extremity; the 2D domain right end is free to move. As is applied, its conjugate force is calculated.

The adopted domain elastic properties are: Young's modulus, and Poisson's ratio . The bar properties are: Young's modulus, , inertia moment, and cross sectional area . For this example it is considered the debonding model 2, with the following parameters:, , and .

A boundary mesh with 120 linear elements is adopted to approximate the matrix, while 100 uniform cubic finite elements were adopted to model the single bar. Finer meshes have been tested to confirm that the discretization adopted was enough fine to give accurate results.

In the Figure 13 it is presented the traction curves along the interface (at bar) to nine different imposed displacements, which shows the evolution of the bar pullout. As one can see, the debonded region keeps the constant value, , at all increments. When the pulling out is completed the final load is exactly the expected.

Figure 14 shows the displacement curves, in meter, at the bar-domain interface. It is possible to verify a decrease on the slope of displacements for the points located near the load application, as the loading progresses.

Figure 15 shows the relative displacements between bar and domain. According to the definition given in equation (32), the bar displacements are obtained by subtracting the results of figure 14 from figure 15.

The evolution of decoupling is illustrated in figure 16 wherein the domain displacements are unequal to the bar displacements. At the first displacement increment it is verified that almost all nodes are perfectly coupled, except the ends nodes. As the free end displacement is increased, other nodes begin to decouple until the ninth increment situation in which the bar displacements continue increasing and domain displacement decreasing.

One important aspect shown by this example, more evident in figure 13, is that the limitation of adherence force by the debonding process regularizes the shear contact force and other variables.

4.4.2 Example 2

The structure analyzed in this example is shown in Figure 17, it is important to note that this example cannot be solved without considering bending stiffness and transverse contact forces as did in this work. This structure is a 2D deep foundation, i.e., a pile embedded in an infinite soil. In order to simulate the presence of nearby structures and a rigid supporting rock mass, some soil displacement restrictions are imposed. The pile is considered inclined at an angle of . A soil region, with 5 meters length () and 20 meters depth (), is considered. The pile has 4 meters in length () and a prescribed load of , and at its top. Loads increase from zero to the reference values in 50 equal increments.

Figure 17

A mesh of 300 linear elements is used to model the soil and 100 finite elements (line elements) are used to model the pile. The soil elastic properties are: and . The pile properties are: , , and . Both adherence models are considered. The model 2 parameters are: , , and . For model 1 it is adopted: .

Figures 19-21 show the displacement in the axial and transverse directions, as well as the rotation along the pile for each reference node.

As can be seen in figures ^19-21, the pile behavior for model 1 presents severe changes. These changes occur when the maximum load capacity is reached (increment 35); then slip occurs without further gain of resistance. It is important to mention that over the load capacity the solution is unstable and the system loss objectivity.

For model 2 as the interface forces do not reach the residual part of the model the maximum load capacity is not achieved and no abrupt change occurs in the pile behavior. However, if the total load capacity is reached (all contact points reach the residual part of the model) the collapse occurs.

Figure 22 (a) and (b) show the displacements in the direction, in meters, for the soil internal points considering the models 1 and 2, respectively. Figures 23 (a) and (b) show the stress values, , for the soil internal points for models 1 and 2.

5 CONCLUSIONS

In this work a BEM/FEM coupling among frame bars and 2D continuum is successfully developed and implemented. The domain is modeled by BEM and reinforcement is modeled by FEM. The combination of the two methods is made by writing displacement compatibility and interface equilibrium. The most important feature of the formulation is the consideration of sliding at frame/continuum interface. This procedure is able to model, for example, the progressive failure of pile-soil interaction until reaching the collapse load, or the progressive failure of fiber reinforced bodies considering the influence of shear and normal contact forces.

Regarding the solution behavior an important conclusion should be stated. As boundary elements are able to model high stress concentrations, the values of contact forces at the beginning or ending of any perfect bounding region present a strong perturbation. The use of redundant algebraic equations and the least squares method are tested here and result in a small improvement of this phenomenon. The increasing of discretization reduces the extension of perturbation but increases the near singular stresses.

The complete solution for this problem results when using a more realistic model that allows the natural stress relaxation at singularities. The developed non linear behavior of contact forces, considering the sliding or decoupling between reinforcement and continuum, completely regularizes the contact force behavior, leading to reliable solutions for low or high load situations. Further developments are the consideration of non-linear behavior for both continuum and frame media.

Acknowledgements

Authors would like to acknowledge CNPq (National Counsel of Technological and Scientific Development) and FAPESP (São Paulo Research Foundation) for the financial support.

• Aliabadi, M. H., (1997). Boundary element formulations in fracture mechanics. Apply Mechanics Review 50:83 - 96.
• Beer, G., Watson, J., (1992). Introduction to finite and boundary element methods for engineers. Wiley (New York).
• Bonnet, M. (1999).Boundary integral equation method for solids and fluids. Wiley .
• Botta, A. S. (2003).Método dos Elementos de Contorno para análise de corpos danificados com ênfase no fenômeno da localização de deformações, Ph.D. Thesis (in Portuguese), University of Sao Paulo, Brazil
• Botta, A. S., Venturini, W. S., (2005). Reinforced 2d domain analysis using BEM and regularized BEM/FEM combination. Computer Modeling in Engineering and Sciences 8(1):15 - 28.
• Brebbia, C. A.(1978). The boundary element method for engineers. Pentech Press (London).
• Brebbia, C. A., Domingues, J.(1992). Boundary elements: an introductory course. Comp. Mech. Publ. (Southampton).
• Brebbia, C. A., Walker, S. (1980). Boundary element techniques in engineering. Newnes-Butterworths (London).
• Coda, H. B., (2001). Dynamic and static non-linear analysis of reinforced media: a BEM/FEM coupling approach. Computers and Structures 79:2751 - 2765.
• Coda, H. B., Venturini, W. S., (1999). On the coupling of 3D BEM and FEM frame model applied to elastodynamic analysis. International Journal of Solids and Structures 36(31/32):4789 - 4804.
• Estorff, O., Hagen, C., Iterative coupling of FEM and BEM in 3D transient elastodynamics, Eng. Anal. Boundary Elements 29:775 - 787
• Frangi, A., Novati, G., Springhetti, R., Rovizzi, M. (2002). 3D fracture analysis by symmetric Galerking BEM. Comput Mechanic 28:220 - 232
• Leite, L. G. S., Coda, H. B., Venturini, W. S., (2003). Two-dimensional solids reinforced by thin bars using the boundary element method. Engineering Analysis with Boundary Elements 27:193 - 201
• Leonel, E. D. (2009).Nonlinear boundary element models to analyze fracture problems and reliability / optimization models applied to structures submitted to fatigue, Ph.D. Thesis (in Portuguese), University of Sao Paulo, Brazil
• Pozrikidis, C., (2006). A note on the relation between the boundary- and finite-element method with application to Laplace's equations in two dimensions. Engineering Analysis with Boundary Elements 30:143 - 147.
• Radtke, F. K. F., Simone, A., Sluys, L. J., (2011). A partition of unity finite element method for simulating nonlinear debonding and matrix failure in thin fiber composites. International Journal for Numerical Methods in Engineering 86:453 - 476.
• Rocha, FC., (2010). Analysis of domains reinforced considering models of adherence (in Portuguese). Caderno de Engenharia de Estruturas, 12: 93 - 105.
• Zienkiewicz, O., Kelly, D., Bettess, P., (1977). The coupling of the finite element method and boundary solution procedures. International Journal for Numerical Methods in Engineering 11:355 - 375.

Sliding frame-solid interaction using BEM/FEM coupling

Fabio Carlos da Rocha^I; Wilson Sergio Venturini^II; Humberto Breves Coda^III

Publication Dates

• Publication in this collection
  06 May 2014
• Date of issue
  Dec 2014