A block-coupled Finite Volume methodology for problems of large strain and large displacement
AA block-coupled Finite Volume methodology forproblems of large strain and large displacement
L.R. Azevedo a, ∗ , P. Cardiff d , F.J. Galindo-Rosales c , M. Schafer b a Graduate School of Computational Engineering, Technische Universitat Darmstadt,Dolivostrasse 15, 64293 Darmstadt, Germany b Chair of Numerical Methods in Mechanical Engineering, Technische UniversitatDarmstadt, Dolivostrasse 15, 64293 Darmstadt, Germany c Centro de Estudos de Fen´omenos de Transporte (CEFT), Dept. Engenharia Qu´ımica,Faculdade de Engenharia da Universidade do Porto, 4200-465 Porto, Portugal d University College Dublin, School of Mechanical and Materials Engineering, Belfield,Ireland
Abstract
A nonlinear block-coupled Finite Volume methodology is developed for largedisplacement and large strain regime. The new methodology uses the samenormal and tangential face derivative discretisations found in the original fullycoupled cell-centred Finite Volume solution methodology for linear elasticity,meaning that existing block-coupled implementations may easily be extendedto include finite strains. Details are given of the novel approach, including useof the Newton-Raphson procedure on a residual functional defined using thelinear momentum equation. A number of 2-D benchmark cases have shownthat, compared with a segregated procedure, the new approach exhibits errorswith many orders of magnitude smaller and a much higher convergence rate.
Keywords:
Cell-centred Finite Volume method, Finite Area method, Finiteelasticity, Block-coupled, Solid mechanics OpenFOAM
1. Introduction
The Finite Volume Method (FVM) has been been successfully used forcomputational solid mechanics (CSM) since late 1980s. For a detailed his-torical review, see e.g. [1]. At present, the typically employed formulationis known as Segregated (SEG). This methodology closely resembles the pro-cedures commonly used in fluid dynamics where memory-efficient segregatedsolution algorithms are used in conjunction with iterative linear solvers. Inpractice, the linear momentum vector equation is temporarily decoupled intothree scalar component equations that are independently solved, where outer ∗ Corresponding author
Email address: [email protected] (L.R. Azevedo)
Preprint submitted to Computer Meth. Appl. Mech. Engineering September 15, 2020 a r X i v : . [ c s . C E ] S e p ixed-Point/Picard iterations provide the required coupling [2]. It is a flexiblemethod of discretisation, in the sense that it does not constraint the constitu-tive equations. But, its major drawback is that it can present poor convergencewhenever there is a strong coupling between displacement components [2]. Toovercome such inadequacy, it was recently proposed by Cardiff et al. [2] ablock-coupled solution methodology, where inter-component coupling is implic-itly included as coefficients in a block matrix; hereafter named BC, it has shownitself to be much faster than SEG for those strongly coupling test cases (by afactor of 2.5-6 times [2]). Furthermore, the BC solver resulted in less executiontime and memory requirements than a finite element software for the set of casestested (in fact it was almost 6 times faster and used 8 times less memory). Nev-ertheless, the current BC formulation is tied to only one constitutive equationand linear elasticity. The current article presents the first attempt to generalizesuch methodology in order to add support for large strain and large displace-ment. The article is constructed as follows: Section 2 outlines the mathematicalmodel, derived from the governing momentum equation and neo-Hookean con-stitutive relation. The novel nonlinear FV discretisation is presented in Section3. Subsequently, in Section 4 it is presented the application of the new approachto five representative benchmark test cases, where accuracy of the method iscompared with that of Segregated approach. Finally, the main findings of thecurrent investigation, and suggestions for future works, are given in Section 5.
2. Mathematical model
Neglecting inertia and body forces for clarity, the conservation of linear mo-mentum for an arbitrary body of volume Ω bounded by surface ∂ Ω with outwardfacing unit normal N is given in strong integral form as: (cid:90) Ω ∇ · P dV = (cid:73) ∂ Ω P · N dS = (1)The first Piola-Kirchhoff P is given by P = F T · Σ , (2)where F is the deformation gradient tensor and Σ is the second Piola-Kichhoffstress tensor. This work adopts the the compressible and isotropic neo-Hookeanhyperelastic model, i.e. Σ = µ ( I − C − ) + λ (ln J ) C − , (3)where C = F T · F is the right Cauchy-Green deformation tensor and { µ, λ } arethe Lam constants. The elasticity tensor C = 2 ∂ Σ /∂ C of this model has theright-minor symmetry (see Appendix A).The discretisation of the new methodology requires a new tensor, say T d ( d = 1 , , transformed lasticity tensor M = F · C · (3) F T (the operator · (3) is a contraction at thethird index) and the face normal N , and is defined as T d = M aJdL N J e a ⊗ e L = F aI C IJKL F dK N J e a ⊗ e L = F aI C IJKL f dK N J e a ⊗ e L ( f d = f dK e K = F dK e K ) ≡ ( F · C ) : (2 , ( N ⊗ f d ) . . (4)Considering the neo-Hookean model, T d is: T d = λ ( A · N ) ⊗ ( C − · f d ) + ( µ − λ ln J ) (cid:104) ( A · f d ) ⊗ b + ( b · f d ) A (cid:105) , (5)where A = F · C − and b = N · C − = C − · N .
3. Numerical method
The mathematical model presented in the preceding section is now discre-tised using a semi-implicit coupled manner and a cell-centred-based FV ap-proach, providing a discrete approximation of the previously presented exactintegral. The discretisation procedure is separated into two distinct parts: dis-cretisation of the solution domain and discretisation of the governing equations.If the temporal effects were considered, time would also be discretised into afinite number of time increments, where the mathematical model is solved in atime-marching manner.
The starting point for a FV discretisation is to decompose the solution spatialdomain B , which is usually approximated by arbitrary and finite number n C ofcontiguous convex polyhedral cells (also known as finite volume) Ω C ’s boundedby faces that do not overlap. But this work adopts a specific polyhedral: therectangular cuboid. The reason for choosing rectangular cuboids is to avoidnon-conformal (skewed and/or non-conjunctional) mesh [3] and the complexitiesthat arise from it. This way, investigation efforts focus only on the core (i.e.minimal structure to be fully usable) of the NLBC methodology. Non-essentialextensions can be added to NLBC after an extensive investigation of the core.The approximation mentioned above is written as B ≈ B d = n C (cid:91) C =1 Ω C , (6)i.e. the continuous body B is approximated by the computational domain B d which is the union of n C cells. The Figure 1 shows (for two-dimensional case)the configuration of one cell Ω C ∈ B d . Before proceeding, note the geometric3 igure 1: Discretisation of a body B into cells Ω C ’s. Every cell Ω C has a boundary ∂ Ω C . Notethat because of the rectangular cuboid restriction, the boundary domain ∂B is approximatedin a castellated staircase manner. parameters shown in the Figures (2a-b) which are needed in the FV discreti-sation process of the governing equations. The Figure (2a) shows a cell withits neighbours F ’s and their face centroids f ’s. The other image (Fig. 2b)exemplifies a typical cuboid cell Ω C , having volume V C and the centroid, orcomputational node, located at the point C. To solve the governing equation (Eq. (1)), it is first rewritten as R ( ∇ U ) = (cid:73) ∂ Ω (cid:101) P ( ∇ U ) · N dS = (7)where R can be called residual function , since R ( ∇ U ) is the so-called, inFinite Element Analysis terminology, residual or out-of-balance force [5]. Thesolution of this equation is sought using the Newton-Raphson iterative processwhereby, given a solution estimate ∇ U n − at iteration n −
1, a new value ∇ U n = ∇ U n − + ∇ ◦ δ u · F ◦ is obtained by establishing the linear approximation: R ( ∇ U n ) ≈ R ( ∇ U n − ) + ∂ R ( ∇ U n − ) ∂ ∇ U : ( ∇ ◦ δ u · F ◦ ) = (8)where ∇ ◦ δ u is the incremental displacement gradient and F ◦ = F n − (seeAppendix C). Using a simplified notation, the equation to be solved is ∂ R ( ∇ U ◦ ) ∂ ∇ U : ( ∇ ◦ δ u · F ◦ ) = − R ( ∇ U ◦ ) (9)4 igure 2: a) A cell Ω C , with centroid C of a 2D discretised domain, and its neighbours F i ; b)A cell Ω C with its geometric parameters used in the finite volume discretisation. In style of[4]. which corresponds to (cid:73) ∂ Ω (cid:34) ∂ (cid:101) P ◦ ∂ ∇ U : ( ∇ ◦ δ u · F ◦ ) (cid:35) · N dS (cid:124) (cid:123)(cid:122) (cid:125) surface force increment = − (cid:73) ∂ Ω (cid:101) P ◦ · N dS (cid:124) (cid:123)(cid:122) (cid:125) old surface force (10)where the simplified notation (cid:101) P ◦ = (cid:101) P ( ∇ U ◦ ) has been used. The integral of the surface force increment term is approximated as: (cid:73) ∂ Ω C (cid:34) ∂ (cid:101) P ◦ ∂ ∇ U : ( ∇ ◦ δ u · F ◦ ) (cid:35) · N dS = (cid:88) Γ f ∈ ∂ Ω C (cid:90) Γ f (cid:34) ∂ (cid:101) P ◦ ∂ ∇ U : ( ∇ ◦ δ u · F ◦ ) (cid:35) · N dS ( ∂ Ω C is a polyhedral) ≈ (cid:88) f S f (cid:34)(cid:40) ∂ (cid:101) P ◦ ∂ ∇ U : ( ∇ ◦ δ u · F ◦ ) (cid:41) · N (cid:35) f (mid-point rule integration) . (11)Substituting for A = ∇ ◦ δ u · F ◦ into equation (see Appendix A for derivation) ∂ (cid:101) P ∂ ∇ U : A = A · (cid:101) Σ + M : A , ∀ A ∈ V , (12)5ields (cid:34)(cid:40) ∂ (cid:101) P ◦ ∂ ∇ U : ( ∇ ◦ δ u · F ◦ ) (cid:41) · N (cid:35) f = (cid:20)(cid:18) ∇ ◦ δ u · F ◦ · (cid:101) Σ ◦ (cid:19) · N (cid:21) f + (cid:20)(cid:26) M ◦ : (cid:18) ∇ ◦ δ u · F ◦ (cid:19)(cid:27) · N (cid:21) f , (13)where (cid:101) Σ ◦ = (cid:101) Σ ( ∇ U ◦ ) and M ◦ = M ( ∇ U ◦ ). Letting NN = N ⊗ N , the firstterm on the right-hand side of Equation (13) is (cid:20) ∇ ◦ δ u · F ◦ · (cid:101) Σ ◦ · N (cid:21) f = (cid:20) ( ∇ ◦ δu ) ab F ◦ bc Σ ◦ cd N d e a (cid:21) f = (cid:20) ( ∇ ◦ δu ) ab v ◦ b e a (cid:21) f ( v ◦ = v ◦ b e b = F ◦ bc Σ ◦ cd N d e b )= (cid:20) ∇ ◦ δ u · v ◦ (cid:21) f = (cid:20) ∇ ◦ δ u · v ◦ n + ∇ ◦ δ u · v ◦ t (cid:21) f ( v ◦ n = ( v ◦ · N ) N , v ◦ t = ( I − NN ) · v ◦ )= (cid:20) ( v ◦ · N ) ∇ ◦ δ u · N + ∇ ◦ δ u · v ◦ t (cid:21) f . (14)Note the projection of v ◦ onto the face normal direction and onto the face plane.This step creates the opportunity to apply the same discretisation procedures,employed by the BC method, to calculate the normal and tangential derivativeterms (i.e. ∇ ◦ δ u · N and ∇ ◦ δ u · v ◦ t , respectively). The second term on the6ight-hand side of Equation (13) is computed as: (cid:20)(cid:26) M ◦ : (cid:18) ∇ ◦ δ u · F ◦ (cid:19)(cid:27) · N (cid:21) f = (cid:20) M ◦ abcd ( ∇ ◦ δu ) ce F ◦ ed N b e a (cid:21) f = (cid:20) M ◦ acd ( ∇ ◦ δu ) ce F ◦ ed e a (cid:21) f ( M ◦ acd = M ◦ abcd N b )= (cid:20) (cid:88) d M ◦ acd ( ∇ ◦ δu ) ce g ◦ ed e a (cid:21) f ( g ◦ d = g ◦ ed e e = F ◦ ed e e )= (cid:20) (cid:88) d T ◦ d · ∇ ◦ δ u · g ◦ d (cid:21) f ( T ◦ d = M ◦ acd e a ⊗ e c )= (cid:20) (cid:88) d T ◦ d · ∇ ◦ δ u · (( g ◦ d · N ) N + ( I − NN ) · g ◦ d )) (cid:21) f (project g ◦ d )= (cid:20) (cid:88) d ( g ◦ d · N ) T ◦ d · ( ∇ ◦ δ u · N ) (cid:21) f + (cid:20) (cid:88) d T ◦ d · ( ∇ ◦ δ u · h ◦ d ) (cid:21) f ( h ◦ d = ( I − NN ) · g ◦ d ) . (15)Using (14) and (15), the term (13) is given as: (cid:20)(cid:26) ∂ (cid:101) P ◦ ∂ ∇ U : ( ∇ ◦ δ u · F ◦ ) (cid:27) · N (cid:21) f = (cid:20) ( v ◦ · N ) ∇ ◦ δ u · N + ∇ ◦ δ u · v ◦ t (cid:21) f + (cid:20) (cid:88) d ( g ◦ d · N ) T ◦ d · ( ∇ ◦ δ u · N ) (cid:21) f + (cid:20) (cid:88) d T ◦ d · ( ∇ ◦ δ u · h ◦ d ) (cid:21) f = (cid:20)(cid:110) ( v ◦ · N ) I + (cid:88) d ( g ◦ d · N ) T ◦ d (cid:111) · ( ∇ ◦ δ u · N ) (cid:21) f (cid:124) (cid:123)(cid:122) (cid:125) (1) + (cid:20) ∇ ◦ δ u · v ◦ t + (cid:88) d T ◦ d · ( ∇ ◦ δ u · h ◦ d ) (cid:21) f (cid:124) (cid:123)(cid:122) (cid:125) (2) . (16)The underlined terms (1) and (2) from the equation before are approximatedusing the same approach employed by the BC method [2], i.e. the normal7erivative term (1) is discretised using the central differencing method as (cid:20) (cid:110) ( v ◦ · N ) I + (cid:88) d ( g ◦ d · N ) T ◦ d (cid:111)(cid:124) (cid:123)(cid:122) (cid:125) H ◦ n · (cid:32) δ u C − δ u F | d CF | (cid:33) (cid:21) f , (17)where the vector connecting the centroids of the cells sharing the common face d CF = X F − X C , and the tangential face derivative term (2) above is discretisedusing the face-Gauss/Finite Area method as (cid:20) S (cid:88) e L e ( M e · v ◦ t ) δ u e + (cid:88) d T ◦ d · (cid:110) | S | (cid:88) e L e ( M e · h ◦ d ) δ u e ) (cid:111)(cid:21) f = (cid:34) S (cid:88) e L e (cid:20) ( M e · v ◦ t ) I + (cid:88) d ( M e · h ◦ d ) T ◦ d (cid:21)(cid:124) (cid:123)(cid:122) (cid:125) H ◦ t · δ u e (cid:35) f . (18) The discretisation of this term uses the mid-point integration approximationas (cid:73) ∂ Ω C (cid:101) P ◦ · N dS = (cid:88) Γ f ∈ ∂ Ω C (cid:90) Γ f (cid:101) P ◦ · N dS ( ∂ Ω C is a polyhedral) ≈ (cid:88) f (cid:20) (cid:101) P ◦ · S (cid:21) f (mid-point rule integration) , (19)where (cid:101) P ◦ is the last known value of the first Piola-Kirchhof stress tensor. The boundary conditions are handled in the same way as in the BC method,except by the fact that, instead of U f , U ◦ f + δ u f is used. Thus, the originalEquation I · U = U b in BC, for example, becomes I · (cid:0) U ◦ f + δ u f (cid:1) = U b = ⇒ I · δ u f = U b − U ◦ f , (20)and Equation (18) in BC becomes T b = (cid:20) (cid:101) P ( ∇ U ◦ ) · N + (cid:26) ∂ (cid:101) P ( ∇ U ◦ ) ∂ ∇ U : (cid:18) ∇ ◦ δ u · F ◦ (cid:19)(cid:27) · N (cid:21) f (21)to be discretised using the same processes applied to the surface force incre-ment and to the old surface force terms described before. The symmetry planeboundary condition is discretised analogously to BC’s approach.8 .3. Solution procedure Assembling Equation (10) using (17), (18) and (19) along with the neo-Hookean material model equations produces a linear algebraic equation withthe same structure seen in the linear Block-Coupled method. However, insteadof solving for the total displacement U , it is solved for incremental displacement δ u . That is, for each control volume C , the final discretised form of the mo-mentum equation can be arranged in the form of N i linear algebraic equations: A C · δ u C + (cid:88) F A F · δ u F = R C (22)where the summation is over the control volume faces. The boundary discreti-sation creates an additional N b linear equations with the same structure as Eq.(22), one for each boundary face centre. These two sets of linear equations arethen assembled forming a linear system of equations:[ A ][ δ u ] = [ R ] (23)where [ A ] is a sparse N × N matrix with the tensorial coefficients A C on thediagonal and the tensorial coefficients A F form the matrix off-diagonal. Thetotal number of computational points being N = N i + N b .Just as in BC, the tangential derivative terms contribute solely to the off-diagonal coefficients, thus [ A ] is not diagonally dominant, in contrast to thesegregated methodology. Therefore, the standard preconditioned ConjugateGradient (CG) methods may not guarantee convergence. As an alternative,the system of linear equations can be solved using, e.g. Bi-Conjugate GradientStabilised (BiCGStab), Generalised Minimal Residual (GMRes) or even directmethods [2].
4. Method verification
In this section, the accuracy and robustness of the novel nonlinear block-coupled methodology is examined for five separate representative test cases andcomparing the numerical prediction to the available analytical solutions. Themethods SEG, BC and NLBC were implemented as a Matlab toolbox callednFVM to generate the results presented in this section. Note that, for all testcases examined here, a solution is considered converged when the residual fallsbelow 10 − . It can easily be shown that when the NLBC method is restricted to thelinearised elasticity framework, it reduces to BC formulation. Thus, the lattercan be seen as a special case of the former. The results from the next test case,that of a slender 2-D cantilever undergoing bending, show this fact by means ofnumerical simulation. This case was used by Cardiff et al. [2] in their seminalwork on the BC method. 9 L Figure 3: Geometry and boundary conditions for the slender cantilever beam in bending testcase. Figure 4: Deformed profile (scaled by factor of 10) for mesh 100x5 cells.
The geometry of the test case, shown in Figure (3), consists of a rectanglebeam 2 x 0.1 m with a Youngs modulus E of 200 GPa and a Poissons ratio ν of 0.3. Three uniform quadrilateral meshes were considered: 60x3, 100x5and 300x15 cells. The mesh with 100x5 cells is shown in Figure (4). Thebeam is fixed at the left end, by imposing the boundary displacement condition U = [0 0] T m, and is subjected to a uniform distributed traction at the otherend, by imposing the boundary traction condition T = [0 1] T MPa. The topand bottom boundaries are traction-free, i.e. T = . Plane strain conditionsare assumed.This problem has analytical solution and the deflection on the right-end ofthe beam is given as [6]:∆ = P L (cid:16) E − ν (cid:17) I = 14 . × − m (24)where P = 0 . × N is the applied load, L = 2 m is the length of the beam, and I = bh = . m m is the second moment of area of the beam about its bendingaxis. A metric defined as the difference between the predicted displacementand the analytical solution shows that both results from BC and NLBC matchconsistently (Fig. 5), reflecting the analytical proof of equivalence between theformulations inside the boundaries of the linearised elasticity framework.The NLBC method converged with only one correction step. Finally, justfor comparison, the SEG method needs more than 23000 correction steps formesh 60x3 cells. 10 B ea m de f e c t i on e rr o r ( % ) Mesh spacing (normalized by beam length)BCNLBC Figure 5: Error in cantilever end-deflection for different mesh refinements. Clearly the ap-proaches match consistently.Figure 6: The Coarsest (3 × ×
64 cells).
Using finite elasticity framework, which allows simulation of accurate largedisplacement-large strain models, it is presented here the comparison of NLBCwith the SEG solution procedure. All test cases use the unit square domain andits five uniform discretisation levels. In particular, five Cartesian meshes wereconsidered: 3 ×
3, 8 ×
8, 16 ×
16, 32 ×
32 and 64 ×
64 cells. The coarsest and finestmeshes are shown in Figure (6).The following metrics were defined to quantify the difference between thepredicted displacement and the analytical solution: e abs Mean error = n cells n cells (cid:80) i =1 r i Max error = max { r , r , ..., r n cells } Min error = min { r , r , ..., r n cells } , (25)11here n cells is the total number of cells composing the mesh, the sum is overall cells and considering a cell C a , r a = | U C a calculated − U C a analytic | . Every test casewas split into two versions: one for displacement-only (Dirichlet) boundary con-ditions and another for traction-only boundary (Neumann) conditions (exceptfor one boundary, which is set to zero-displacement in order to avoid rigid-bodymotions). This split scheme isolates patterns which arise due to different bound-ary condition discretisations employed by NLBC and errors from each one canbe investigated individually.All test cases were created using the accepted standard of verification testing,the Method of Manufactured Solutions (MMS), which allows validation againstanalytical solution [7]. A MMS test prescribes the deformation map ϕ , or anyother map that allows one to recover it.The density ρ was set to 216 kg/m ; the Youngs Moduli E and Poissonsratio ν were set to 0.02 GPa and 0.3, respectively. Two homogeneous uniaxial strain MMS were simulated. The deformationgradient F is the mapping prescribed for these cases and it is given as: F = φ ( t ) 0 00 1 00 0 1 , where φ ( t ) = 1 + (Λ − t and 0 ≤ t ≤ . (26)Note that F is homogeneous, i.e. does not depend on a material point X . Thedeformation map is defined as: x = ϕ ( X ) ≡ F · X and it is used to set thedisplacement boundary condition by imposing U = x − X (27)at the boundary face centroids. A traction boundary counterpart can be set bynoting that a traction T acting on the face with unit normal N is T = (cid:101) P ( ∇ U ) · N = (cid:101) P ( F − I ) · N . (28) Compression for displacement boundary
A variation of the homogeneous uniaxial strain test case described in [7] ispresented in this section. However, instead of traction, displacement boundarycondition was adopted. Two compression levels were investigated by assigningdifferent values for the compression factor Λ, in particular, Λ = 0 .
65 and = 0 . e abs < − , regardless the method, for Λ = 0 .
65 (see Fig. 9). Theconvergence in all scenarios was achieved with only one correction step, i.e. n corr = 1. When Λ is decreased to 0 .
1, the SEG method produces e abs < − .The errors for NLBC also increase when Λ get smaller, but they are still rela-tively small ( e abs < − ) and only one correction is needed, considering anymesh. 12 igure 7: The final deformed domain at com-pression level Λ = 0 .
65. Figure 8: The final deformed domain at com-pression level Λ = 0 . -18 -16 -14 -12 -10 -8 e ab s Mesh refnement (number o f FVs) Mean
Error (SEG)Max.
Error (SEG)Min.
Error (SEG)Mean
Error (NLBC)Max
Error (NLBC)Min
Error (NLBC)
Figure 9: Errors from compression for displacement boundary test case using Neo-Hookeanmaterial. The missing data corresponds to when the difference between the solutions is belowmachine precision. igure 10: The final deformed domain for tensile strain case and for displacement boundarycondition. Compression for traction boundary
Just changing from Dirichlet to Neumann makes the convergence a challengefor both methods, in particular, they are not able to simulate big compression.The summarized results gathered from simulations are: • The SEG method converges only when using the 3 × ≥ .
8, but with relatively high errors ( e abs > − ). • The NLBC method also converges only for Λ ≥ . ×
16. For thesescenarios, e abs < − . Tension for displacement boundary
The cases above were repeated, but with Λ > e abs < − (see Fig. 11). Note that NLBCproduces significantly smaller errors. Both methods converge with only onecorrection step. Tension for traction boundary
Once again, when traction is introduced, SEG does have convergence prob-lems. In fact, it does not converge for Λ much greater than one. And even whenΛ is close to one, e.g. 1.2, the errors are relatively high (either with nFVMor S4F). Regarding NLBC’s results, they show good agreement with analyticalsolution. The method converges for all meshes, for any Λ ∈ (1 ,
2] and the errorsare relatively small ( e abs < − ), but much higher than the corresponding testwhich uses displacement boundaries (see Fig. 12 and compare with Fig. 11).14 -18 -16 -14 -12 -10 -8 e ab s Mesh refnement (number o f FVs)Mean
Error (SEG)Max.
Error (SEG)Min.
Error (SEG)Mean
Error (NLBC)Max
Error (NLBC)
Figure 11: Errors from tension for displacement boundary test case. The missing data corre-sponds to when the difference between the solutions is below machine precision. -10 -9 -8 -7 -6 e ab s n c o rr Mesh refnement (number o f FVs) Mean
Error (NLBC)Max
Error (NLBC)Min
Error (NLBC)n corr (NLBC)
Figure 12: Error in tension for traction boundary test case as mesh is refined. The n corr asmesh is refined is also shown. Results are only for NLBC, since SEG method could not handlethis case. igure 13: Deformed profile for mesh 16 ×
16 in shear test case.
This test case consists of a simple shear [5]. The deformation gradient forthis manufactured solution is similar to that of the uniaxial test case and isgiven by (being the shear factor ω = 0 .
45 chosen arbitrarily): F = φ ( t ) 00 1 00 0 1 , where φ ( t ) = ωt and 0 ≤ t ≤ . (29)The Figure (13) shows the deformed profile for mesh 16 ×
16. The boundarycondition is imposed in the same manner as it was done in uniaxial test cases.
Shear for displacement boundary
The results from simulations were qualitatively similar to that of the uniaxialcompression, or tension, for displacement boundary (compare Fig. 9 and 11with Fig. 14). The nFVM’s SEG and NLBC needed only one correction forall meshes. The S4F’s SEG was also tested and the output shows that as meshgets refined, it needs more corrections to achieve convergence ( n corr = 23 , ×
3, 8 ×
8, 16 ×
16 respectively). Besides, it did not convergefor meshes finer than 16 × Shear for traction boundary
Once more, when displacement boundaries are replaced by traction bound-aries, the SEG method has convergence problems (both in nFVM and in S4F).In fact, convergence is achieved, however with relatively high errors (see Fig.15). The SEG approach needs more than 100 correction steps to converge andfor the finer the mesh, more correction steps are necessary for convergence (seeFig. 16).The results from the NLBC method were in good agreement with analyticalsolutions and only one correction was needed in order to achieve convergence(see Fig. 16) using any mesh. 16 -18 -16 -14 -12 -10 -8 -6 e ab s n c o rr Mesh refnement (number o f FVs) Mean
Error (SEG)Max.
Error (SEG)Min.
Error (SEG)Mean
Error (NLBC)Max
Error (NLBC)Mean
Error (SEG-S4F)Max.
Error (SEG-S4F)Min.
Error (SEG-S4F)n corr (SEG-S4F)
Figure 14: Error in shear for displacement boundary test case as mesh is refined. The numberof correction n corr as mesh is refined for S4F’s SEG is also shown. The SEG and NLBCimplementations in nFVM needed only one correction to achieve convergence. The final deformation domain (Fig. 17) for SEG is clearly “warped” (andrefining the mesh does not reduce this spurious artifact).
5. Discussion & conclusions
It has been presented a novel nonlinear block-coupled FV methodology whichgeneralises the work of Cardiff et al. [2]. The developed methodology has beeninvestigated by means of numerical simulations, i.e. one test case for infinites-imal elasticity and four for finite elasticity. The accuracy of the methodologyhas been shown through detailed comparison with analytical solutions and nu-merical benchmarks. For all the test cases analysed, NLBC has shown to bean efficient and accurate alternative to SEG method for the analysis of 2-Dproblems in finite elasticity.The novel methodology does not assume small strain or small displacementduring the discretisation process and only requires the presence of a right-minorsymmetric elasticity tensor. Thus, it defines a class of officially supported solidmodels. In fact, it can be demonstrated that a large set of important solid mod-els (those that are hyperelastic, frame-indifferent, homogeneous and isotropic)belong to this class. As a matter of fact, frame-indifference should be requiredindependently of the FVM methodology, since it is required in finite elasticity,otherwise different observers could collect different results [5]. This symmetric-related elasticity restriction should be subjected to investigation in order toestablish if it can be removed or at least weakened.As regards mesh support, NLBC assumes that finite volumes are rectangu-lar cuboids. However, by judging how other FV methodologies handle convex17 -16 -14 -12 -10 -8 -6 -4 -2 e ab s Mesh refnement (number o f FVs) Mean
Error (SEG)Max.
Error (SEG)Min.
Error (SEG)Mean
Error (NLBC)Max
Error (NLBC)Min
Error (NLBC)Mean
Error (SEG-S4F)Max.
Error (SEG-S4F)Min.
Error (SEG-S4F)
Figure 15: Error in shear for traction boundary test case as mesh is refined. n c o rr Mesh refnement (number o f FVs) SEG (nFVM)SEG (S4F)NLBC
Figure 16: Number of corrections in shear for traction boundary test case as mesh is refined.The n corr for NLBC is 1. igure 17: Deformed profiles for SEG using S4F (lower-left corner) and nFVM (lower-rightcorner). On top the result from NLBC. Traction boundary was used. polyhedral, the modification to add support to it should be relatively straight-forward.In conclusion, it has been presented the first attempt to generalise the BCsolution methodology to finite elasticity, for which a Newton-Raphson methodwas employed similar as in finite element analysis.In conclusion, it has been presented the first attempt to generalise the BCsolution methodology to finite elasticity, which closely resembles the Finite Ele-ment Methodology in the sense that all displacement components are solved atthe same time in a big linear system generated by applying the Newton-Raphsonprocedure on an out-of-balance force function.
6. Acknowledgments
This work was supported by the program Cincia sem Fronteiras (Grant233309/2014-4, CNPq, Brazil) and the “Excellence Initiative” of the GermanFederal and State Governments and the Graduate School of ComputationalEngineering at Technische Universitt Darmstadt.
Appendix A. Elastic body
An elastic body can be defined through the following elastic body axiom [8]: a continuum body with reference configuration B is elastic if ∃ (cid:98) σ : V × B →V such that σ m ( X , t ) = (cid:98) σ ( F ( X , t ) , X ) , ∀ X ∈ B, t ≥ (cid:98) σ ( F , X ) T = (cid:98) σ ( F , X ) , ∀ X ∈ B, F ∈ V , det F > . (A.1)Since this work considers only homogeneous bodies, the stress response func-tion (cid:98) σ is considered independent of X , thus σ m ( X , t ) = (cid:98) σ ( F ( X , t )). Because19f this axiom, there are two functions (cid:98) P : V → V and (cid:98) Σ : V → V such that P ( X , t ) = (cid:98) P ( F ( X , t )) and Σ ( X , t ) = (cid:98) Σ ( F ( X , t )) , (A.2)in particular, they must satisfy the relations (cid:98) P ( F ) = (det F ) (cid:98) σ ( F ) · F − T and (cid:98) Σ ( F ) = F − · (cid:98) P ( F ) . (A.3)The axiom of material frame-indifference implies that: ∃ Σ : V → V such that (cid:98) P ( F ) = F · Σ ( C ) and (cid:98) Σ ( F ) = Σ ( C ) (A.4)where C = F T · F is the right Cauchy-Green deformation tensor.Let A ∈ V , then ∂ C ∂ F : A = A T · F + F T · A , (A.5)thus (using the chain rule and ∂ F /∂ ∇ U = I , i.e. the fourth-order identity tensor ) ∂ C ∂ ∇ U : A = ∂ C ∂ F : (cid:18) ∂ F ∂ ∇ U : A (cid:19) = ∂ C ∂ F : (cid:18) I : A (cid:19) = ∂ C ∂ F : A = A T · F + F T · A ≡ · sym( F T · A ) , (A.6)where sym( · ) denotes the symmetric component of a tensor . The field ∂ C /∂ ∇ U is used next.The Green-Lagrange strain tensor E = ( C − I ) and the new function ( Σ ( E ) = Σ ( C ( E )) can be used to find the elasticity tensor C = ∂ ( Σ /∂ E in terms of ∂ Σ /∂ C as ∂ ( Σ ∂ E : A = ∂ Σ ∂ C : (cid:18) ∂ C ∂ E : A (cid:19) (using again the chain rule)= ∂ Σ ∂ C : (cid:18) I : A (cid:19) = 2 ∂ Σ ∂ C : A . (A.7)The arbitrariness of A implies that C = ∂ ( Σ ∂ E = 2 ∂ Σ ∂ C . (A.8) The definition is I ≡ δ ac δ bd e a ⊗ e b ⊗ e c ⊗ e d which implies that I : A = A . The last equation shows that sym( B ) ≡ ( B + B T ). (cid:101) P ( ∇ U ) = (cid:98) P ( F ( ∇ U )) = (cid:98) P ( I + ∇ U ) and (cid:101) Σ ( ∇ U ) = Σ ( C ( ∇ U )), the derivative of the stressresponse function (cid:101) P is given as ∂ (cid:101) P ∂ ∇ U : A = ∂ ( F · (cid:101) Σ ) ∂ ∇ U : A (using Equation (A.4))= (cid:18) ∂ F ∂ ∇ U : A (cid:19) · (cid:101) Σ + F · (cid:18) ∂ (cid:101) Σ ∂ ∇ U : A (cid:19) = A · (cid:101) Σ + F · (cid:18) ∂ (cid:101) Σ ∂ ∇ U : A (cid:19) (using ∂ F /∂ ∇ U = I )= A · (cid:101) Σ + F · (cid:20) ∂ Σ ∂ C : (cid:18) ∂ C ∂ ∇ U : A (cid:19)(cid:21) (using the chain rule)= A · (cid:101) Σ + F · (cid:20) ∂ Σ ∂ C : (cid:18) · sym( F T · A ) (cid:19)(cid:21) (using Eq. (A.6))= A · (cid:101) Σ + F · (cid:20) C : (cid:18) sym( F T · A ) (cid:19)(cid:21) (using Eq. (A.8)) . (A.9)The last equation needs to be extended by taking the C ’s right-minor symmetryinto consideration (2.51) as ∂ (cid:101) P ∂ ∇ U : A = A · (cid:101) Σ + F · (cid:20) C : (cid:18) sym( F T · A ) (cid:19)(cid:21) = A · (cid:101) Σ + F · (cid:20) C : (cid:18) F T · A (cid:19)(cid:21) (using C ’s symmetric property)= A · (cid:101) Σ + F · (cid:20) C αβγδ F aγ A aδ e α ⊗ e β (cid:21) = A · (cid:101) Σ + (cid:18) F · C · (3) F T (cid:19) : A = A · (cid:101) Σ + M : A , ∀ A ∈ V , (A.10)where M is the transformed elasticity tensor defined in the paragraph precedingEquation (2.53). Note that the right-minor symmetry restriction, which couldnot be overcome, creates a class of supported materials. Appendix B. Neo-Hookean model
This is a compressible isotropic hyperelastic material model and its strain-energy function is defined as [5, 9] W ( C ) = (cid:99) W ( I C ) = µ I ( C ) − − µ ln J + λ J ) , (B.1)21here I ( C ) = det C = J . The Lam (material) coefficients λ and µ relatingto the Youngs modulus E and Poisson’s ratio, ν , are given respectively as: µ = E ν ) ; νE (1 + ν )(1 − ν ) for plane stress; and νE (1 + ν )(1 − ν ) for plane strainand 3-D. The second Piola-Kirchhoff stress tensor is obtained from Equation(B.1) as Σ = Σ ( C ) = 2 ∂W ( C ) ∂ C = µ ( I − C − ) + λ (ln J ) C − . (B.2)The elasticity tensor can be obtained by differentiation of Equation (B.2) withrespect to the components of E to give, after some algebra using ∂I ( C ) /∂ C = J C − , C as C = ∂ ( Σ ∂ E = 2 ∂ Σ ∂ C = λ C − ⊗ C − + 2( µ − λ ln J ) J , (B.3)where the fourth-order tensor J is defined as J = − ∂ C − ∂ C ⇐⇒ J IJKL = 12 (cid:104) ( C − ) IK ( C − ) JL + ( C − ) IL ( C − ) JK (cid:105) . (B.4)It is straightforward to show that J , C − ⊗ C − (using ( C − ) KL = ( C − ) LK )and therefore the elasticity tensor C above has right-minor symmetry, i.e. C IJKL = C IJLK . (B.5)The full expression for T d is obtained by substituting Equation (B.3) intoEquation (4) resulting in T d = λ ( F · C − · N ) ⊗ ( f d · C − )+ ( µ − λ ln J ) (cid:104) ( F · C − · f d ) ⊗ ( N · C − ) + ( N · C − · f d )( F · C − ) (cid:105) = λ ( A · N ) ⊗ ( f d · C − ) + ( µ − λ ln J ) (cid:104) ( A · f d ) ⊗ b + ( b · f d ) A (cid:105) = λ ( A · N ) ⊗ ( C − · f d ) + ( µ − λ ln J ) (cid:104) ( A · f d ) ⊗ b + ( b · f d ) A (cid:105) ( C − is symmetric) , (B.6)where A = F · C − and b = N · C − = C − · N . Appendix C. Mathematical framework for incremental description
To describe the incremental approach, let the following maps be defined: ϕ : X ∈ B → B (cid:48) (cid:51) x , φ : X ∈ B → B ◦ (cid:51) y and χ : y ∈ B ◦ → B (cid:48) (cid:51) x , (C.1)22 igure C.18: A deformation is illustrated by considering the reference configuration B , theold configuration B ◦ and the current deformed configuration B (cid:48) . The black dots representone and the same material particle. where B ◦ can be thought as an intermediate (also labeled as old) body statebetween the reference body state B and the current body state B (cid:48) (see Fig.C.18). Then, by using the composition χ ◦ φ , it is derived the relation betweenthe deformation gradients associated with the mappings as x = ϕ ( X ) = χ ( φ ( X )) = ⇒ ∂ ϕ ∂ X (cid:124)(cid:123)(cid:122)(cid:125) F = ∂ χ ∂ y (cid:124)(cid:123)(cid:122)(cid:125) δ F · ∂ φ ∂ X (cid:124)(cid:123)(cid:122)(cid:125) F ◦ , (C.2)thus obtaining the relation between the deformation gradients as F = δ F · F ◦ . (C.3)The symbol δ F is the well known [10] incremental (or relative ) deforma-tion gradient , and its relation with the so-called incremental (or relative ) displacement gradient ∇ ◦ δ u is found using the incremental displacement field δ u : B ◦ → B (cid:48) (Fig. C.18) as: δ u ( y ) = χ ( y ) − y = ⇒ ∂ δ u ∂ y (cid:124) (cid:123)(cid:122) (cid:125) ∇ ◦ δ u = ∂ χ ∂ y − ∂ y ∂ y = δ F − I , (C.4)therefore I + ∇ ◦ δ u = δ F . (C.5)The intermediate (or old) displacement field U ◦ ( X ) = φ ( X ) − X gives riseto the intermediate (or old) displacement gradient ∇ U ◦ = ∇ φ − I = ⇒ ∇ U ◦ = F ◦ − I . (C.6)The gradient increment ∇ ◦ δ u · F ◦ is finally found using (C.3), (C.5) and (C.6)23s F = δ F · F ◦ I + ∇ U = ( I + ∇ ◦ δ u ) · F ◦ = ⇒∇ U = F ◦ + ∇ ◦ δ u · F ◦ − I = ⇒∇ U = ∇ U ◦ + ∇ ◦ δ u · F ◦ . (C.7) References [1] P. Cardiff, I. Demirdi, Thirty years of the finite volume method for solidmechanics (2018). arXiv:1810.02105 .[2] P. Cardiff, ˇZ. Tukovi, H. Jasak, A. Ivankovi, A block-coupled finite volumemethodology for linear elasticity and unstructured meshes, Computers &Structures 175 (15) (2016) 100–122. doi:10.1016/j.compstruc.2016.07.004 .[3] F. Moukalled, L. Mangani, M. Darwish, The Finite Volume Method inComputational Fluid Dynamics: An Advanced Introduction with Open-FOAM and Matlab, 1st Edition, Springer, 2015.[4] I. de Oliveira, Using foam-extend to assess the influence of fluid-structureinteraction on the rupture of intracranial aneurysms, Ph.D. thesis (082017).[5] J. Bonet, R. D. Wood, Nonlinear Continuum Mechanics for Finite ElementAnalysis, 2nd Edition, Cambridge University Press, 2008.[6] S. Timoshenko, J. N. Goodier, Theory of elasticity, 3rd Edition, McGraw-Hill, 1970.[7] K. Kamojjala, R. Brannon, A. Sadeghirad, J. Guilkey, Verification testsin solid mechanics, Engineering with Computers 31 (2013) 193–213. doi:10.1007/s00366-013-0342-xdoi:10.1007/s00366-013-0342-x