Flat-Joint Model

The flat-joint model can be installed at both ball-ball and ball-facet contacts, and is referred to in commands and FISH by the name flatjoint.

Introduction

A flat-joint contact and its corresponding flat-jointed material are shown in Figure 1. The flat-joint contact model provides the macroscopic behavior of a finite-size, linear elastic, and either bonded or frictional interface that may sustain partial damage. A flat-jointed material mimics the microstructure of angular, interlocked grains that is similar to marble. The model formulation is given in this document. The initial two- and three-dimensional flat-joint models are described in [Potyondy2012a], [Potyondy2012b], and [Potyondy2013]. The present description defines both the 2D and 3D flat-joint models. The creation and laboratory testing of a flat-jointed material is described in the section “Example Materials 2: Flat-Jointed Material Example” in [Potyondy2017], which is provided in the material-modeling support package. Test problems that examine the behaviors of a single flat-jointed contact and an interlocked grain are provided in [Potyondy2016], which is included in the documentation of the material-modeling support package.

The formulation begins with a description of a flat-joint contact and its corresponding flat-jointed material, and is followed by a description of the flat-jointed interface. The activity-deletion criterion, force-displacement law, energy partitions, properties, methods, and callback events of the flat-joint contact model are then presented, followed by the stiffnesses required to determine a stable timestep. Expressions for element normal force and bending moment that are used in the force-displacement law are provided in the final subsection.

Behavior Summary

The behavior summary consists of a description of the flat-jointed material, followed by a description of the flat-jointed interface.

Flat-Jointed Material

A flat-joint contact and its corresponding flat-jointed material are shown in Figure 1. The flat-joint contact model provides the macroscopic behavior of a finite-size, linear elastic, and either bonded or frictional interface that may sustain partial damage. A flat-jointed material is defined in [Potyondy2017] as a granular assembly in which the flat-joint contact model exists at all grain-grain contacts with a gap less than or equal to the installation gap at the end of the material-finalization phase; all other grain-grain contacts as well as new grain-grain contacts that may form during subsequent motion are assigned the linear contact model. A flat-jointed material mimics a microstructure of angular, interlocked grains that is similar to marble.

../../../../../_images/cmflatjoint_fig1.png

Figure 1: Flat-joint contact (left) and flat-jointed material (right).

The flat-joint contact model provides the macroscopic behavior of a finite-size, linear elastic and either bonded or frictional interface that may sustain partial damage (see Figure 2). The interface is discretized into elements. Each element is either bonded or unbonded, and the breakage of each bonded element contributes partial damage to the interface. The behavior of a bonded element is linear elastic until the strength limit is exceeded and the bond breaks, making the element unbonded; the behavior of an unbonded element is linear elastic and frictional, with slip accommodated by imposing a Coulomb limit on the shear force. Each element carries a force and moment that obey the force-displacement law described below, while the force-displacement response of the flat-joint interface is an emergent behavior that includes evolving from a fully bonded state to a fully unbonded and frictional state.

../../../../../_images/cmflatjoint_fig8.png

Figure 2: Behavior and rheological components of the flat-joint model.

A flat-joint contact simulates the behavior of an interface between two notional surfaces, each of which is connected rigidly to a piece of a body. A flat-jointed material consists of bodies (balls, clumps, or walls) joined by flat-joint contacts such that the effective surface of each body is defined by the notional surfaces of its pieces, which interact at each flat-joint contact with the notional surface of the contacting piece. The notional surfaces are called faces, which are {lines in 2D; disks in 3D}.

The following description of a flat-jointed material applies to the case in which bodies are balls; however, the flat-joint model can be installed at both ball-ball and ball-facet contacts. We refer to the balls of a flat-jointed material as faced grains, each of which is depicted as a {circular in 2D; spherical in 3D} core and a number of skirted faces. The faced grains are created when the flat-joint model is installed at the ball-ball contacts of a packed ball assembly (see Figure 3). An interface exists between each set of adjoining faces and is discretized into elements, with each element being either bonded or unbonded. The breakage of each bonded element contributes partial damage to the interface, and each breakage event is denoted as a crack (see Figure 4)[1]. If the relative displacement at a flat-joint contact becomes larger than the flat-joint diameter, then the adjoining faces may be removed (because the contact may be deleted), making the associated balls locally {circular in 2D; spherical in 3D}. If these balls come back into contact, the behavior will be that of an interface between {circular in 2D; spherical in 3D} surfaces (if the linear contact model is assigned to the new contact).

../../../../../_images/cmflatjoint_fig2.png

Figure 3: Creation of faced grains showing packed ball assembly (left) and initial faced-grain assembly (right).

../../../../../_images/cmflatjoint_fig3.png

Figure 4: Partially damaged flat-jointed material showing faced grains with cracks colored red/blue for tensile/shear failure.

Flat-Jointed Interface

The interface of the flat-joint model remains centered on, and rotates with, the contact plane such that the interface coordinate system coincides with the contact-plane coordinate system. The interface gap (g,g>0isopen) is expressed below in terms of the surface gap and the relative bend-rotation vector.[2] The relative motion of the notional surfaces varies over the interface, and is expressed by the symbols δ and θ, which are related to the relative motion of the piece surfaces at the contact location:

(1)Δδδ=Δδnˆnc+ΔδδsΔθθΔθθwithΔδδs=Δδδs+(Δθθt×r)=˙δδsΔt+(˙θθtΔt×r)[Δδss=Δδδsˆsc,Δδst=Δδδsˆtc]

where the relative displacement and rotation increments (Δδδ and Δθθ) are from Equation (12) in the “Contact Resolution” section. It is only the relative shear displacement (Δδδs) of the 3D model that varies over the interface (when the relative twist rotation is nonzero). Thus, the circumflex will be used only for this quantity. In the remainder of this section, we describe the kinematic variables and interface discretization, first for the 2D model and then for the 3D model.

The interface of the 2D flat-joint model is a rectangle of width 2R and unit-thickness depth. The location of a point x that lies on the interface is expressed by its relative position r=xxc. The interface discretization (see Figure 5) is controlled by the number of equal-length elements in the radial direction (Nr). The area and centroid location of each element are denoted by A(e) and x(e), respectively. These quantities are given by

(2)A(e)=2RtNr(t=1)x(e)=xc+r(e),e=1,2,,Nrwithr(e)=ρ(e)ˆtcρ(e)=R(2e+1+NrNr).

The quantities A(e) and ρ(e) do not change during a simulation. They are set and stored during the first cycle after flat-joint installation.

../../../../../_images/cmflatjoint_fig4.png

Figure 5: Interface discretization of the 2D flat-joint model showing element-numbering convention (left) and parameterization of a typical element (right).

The interface gap of the 2D flat-joint model (see Figure 6) is given by

(3)g(r)=g(γ)=gs+(γR)θbk(smallangleapprox.θbktanθbk)withθbk=θθbˆk,γ[0,2R]

where gs is the surface gap and θθb is the relative bend-rotation vector of Equation (12) in the “Contact Resolution” section. When expressed in the γ system, the interface gap varies only with γ, and is related to gs and θbk. The interface gap is not affected by the relative shear motion (which is a small-strain approximation).

../../../../../_images/cmflatjoint_fig5.png

Figure 6: Interface gap of the 2D flat-joint model corresponding with motion at top of image.

For the 2D flat-joint model, the following components of the interface relative motion are tracked:

(4)gs:=gs+Δδnδst:=δst+Δδδsˆtcθbk:=θbk+Δθθbˆk

The interface of the 3D flat-joint model is a disk of diameter 2R. The location of a point x that lies on the interface is expressed by its relative position r=xxc. The interface discretization (see Figure 7) is controlled by the number of elements in the radial (Nr) and circumferential (Nα) directions. The area and centroid location of each element are denoted by A(e) and x(e), respectively. These quantities are given by

(5)A(e)=dA=12(α2α1)(r22r21),r1<r2,α1<α2,e=1,2,,NrNαx(e)=xc+r(e)withr(e)=r(e)sˆsc+r(e)tˆtc=s(e)ˆsc+t(e)ˆtcs(e)=sdAdA=23(sinα2sinα1α2α1)(r32r31r22r21)t(e)=tdAdA=23(cosα1cosα2α2α1)(r32r31r22r21)

where the superscript (e) notation has been dropped for the r and α terms, which are given by

(6)r(e)1=IΔr,e=1,2,,NrNαr(e)2=(I+1)Δrα(e)1=JΔαα(e)2=(J+1)ΔαwithΔr=R/RNrNr,Δα=2π/2πNαNαI=(e1)/(e1)NαNα,J=e1INα

where floor(x)=x is the largest integer not greater than x. The quantities A(e), s(e), and t(e) do not change during a simulation. They are set and stored during the first cycle after flat-joint installation.

../../../../../_images/cmflatjoint_fig6.png

Figure 7: Interface discretization of the 3D flat-joint model showing element-numbering convention (left) and parameterization of a typical element (right).

There are two coordinate systems associated with the interface of the 3D flat-joint model (see Figure 8): the nst system, which coincides with the contact plane coordinate system, and the nξη system, which is dependent on the contact plane coordinate system as follows. The relative bend-rotation vector (θθb) orients the nξη coordinate system such that ξ is aligned with θθb (at an angle βs from the positive s axis) and ˆnc=^ξξ×^ηη. The relative position of a point xx that lies on the interface can be expressed in either coordinate system by the relations:

(7)r=xxc  =r(s,t)=rsˆsc+rtˆtc  =r(ξ,η)=rξˆξ+rη^ηηwithrs=rˆsc=s,rt=rˆtc=trξ=r^ξξ=ξ,rη=r^ηη=η.

The mapping between these two coordinate systems of a vector S that lies on the interface is given by the relations:

(8){SξSη}=[tgl]{SsSt},{SsSt}=[tgl]T{SξSη}[tgl]=[cosβssinβssinβscosβs]ξs=^ξξˆsc=cosβs,ξt=^ξξˆtc=sinβsηs=^ηηˆsc=sinβs,ηt=^ηηˆtc=cosβs.

The interface gap of the 3D flat-joint model (see the next figure) is given by

(9)g(r)=g(ξ,η)=gs+ηθbξ(smallangleapprox.θbξtanθbξ)withθbξ=θθb^ξξ=θθb,{ξ,η}[R,R]

where gs is the surface gap and θθb is the relative bend-rotation vector of Equation (12) in the “Contact Resolution” section. When expressed in the nξη system, the interface gap varies only with η and is related to gs and θbξ. The interface gap is not affected by the relative shear displacement or the relative twist rotation (which is a small-strain approximation). The interface gap is expressed in terms of s, t, and βs (by substituting the mapping of Equation (8) into (9)):

(10)g(r)=g(s,t,βs)=gs+[tcosβsssinβs]θbξ=gs+[tξssξt]θbξwithθbξ=θθb^ξξ=θθb,{s,t}[R,R]βs[0,2π),ξs=^ξξˆsc=cosβs,ξt=^ξξˆtc=sinβs.
../../../../../_images/cmflatjoint_fig7.png

Figure 8: Interface gap of the 3D flat-joint model corresponding with motion at top of image.

For the 3D flat-joint model, the following components of the interface relative motion are tracked:

(11)gs :=gs+Δδnδss:=δss+Δδδsˆsc,  δst:=δst+Δδδsˆtcθbs:=θbs+Δθθbˆsc,  θbt:=θbt+Δθθbˆtcξs={0, θθb1×1012θbsθθb,  otherwise,  ξt={0, θθb1×1012θbtθθb,  otherwise

Activity-Deletion Criteria

A force may arise at a flat-joint element if it is either bonded or has a negative gap. The contact remains active until the distance between the centers of the notional surfaces is greater than the flat-joint diameter (g2s+(δss)2+(δst)2>2R). Subsequently, the contact is deemed inactive and may be deleted. An inactive flat-joint contact may become active if the pieces overlap. In this case the notional surfaces are initialized.

Force Displacement Law

The force-displacement law for the flat-joint model updates the contact force and moment (Fc=˜F and Mc=˜M) that act at the contact location in an equal and opposite sense on the two pieces (see Figure 1 in the “Contacts Model Framework” topic). This is done by discretizing the interface into elements and providing a force-displacement law for each element so that the force-displacement response of the flat-joint interface is an emergent behavior. The force-displacement law for a flat-joint element is described in this section.

Each element carries a force and moment (F(e) and M(e)) that act at the element centroid in an equal and opposite sense on the notional surfaces. The element forces and moments produce a statically equivalent force and moment at the center of the interface (which coincides with the contact location—see Figure 9 and Figure 10) given by

(12)˜F=eF(e),˜M=e{(r(e)×F(e))+M(e)}

where r(e)=x(e)xc is the relative position and x(e) is the centroid location of element (e). For the 2D model, Equation (12) becomes

(13)˜F=eF(e),˜M=e{(ρ(e)F(e)nˆk)+M(e)}(2Dmodel)

by substituting r(e) from Equation (2) into Equation (12) and expressing the element force in terms of its normal and shear components (defined below).

../../../../../_images/cmflatjoint_fig9.png

Figure 9: Force and moment acting on a typical 2D element (right) and its centroid relative position vector used to determine their contribution to force and moment at the center of the interface (left).

../../../../../_images/cmflatjoint_fig10.png

Figure 10: Force and moment acting on a typical 3D element (right) and its centroid relative position vector used to determine their contribution to force and moment at the center of the interface (left).

The force-displacement law for a flat-joint element updates the element force and moment (F(e) and M(e)), and may modify the element bond state (B(e)). The element force is resolved into a normal and shear force, and the element moment is resolved into a twisting and bending moment:

(14)F(e)=F(e)nˆnc+F(e)sM(e)=M(e)tˆnc+M(e)b(2Dmodel:M(e)t0)

where F(e)n>0 is tension and M(e)t>0 is shown in Figure 8 in the “Contact Resolution” section. The shear force and bending moment lie on the interface and are expressed in the interface coordinate systems:

(15)F(e)s=F(e)ssˆsc+F(e)stˆtc(2Dmodel:F(e)ss0)M(e)b=M(e)b{ˆk,  2D^ξξ,  3D.

For the 3D model, a simplifying assumption[3] gives M(e)t0. F(e)n and M(e)b are updated by integrating the normal stress acting over the element, and F(e)s is updated incrementally based on the effective portion of the relative shear-displacement increment at the element centroid.

The element normal and shear stresses are given by

(16)σ(e)=F(e)nA(e),τ(e)=F(e)sA(e)

where σ(e)>0 is tension, and the stresses act at the element centroid. The interface normal stress (σ,σ>0istension) is given by

(17)σ(r)={0,unbondedandg(r)0kng(r),otherwise

where kn is the normal stiffness and g is the interface gap of {Equation (3) in 2D; Equation (10) in 3D}. The interface normal stress is proportional to the gap, and the unbonded region with a positive gap carries no load. A tensile load is carried only in a bonded region with a positive gap, and a compressive load is carried wherever the gap is negative (regardless of its bonding state). The gap and normal stress are continuous over the interface, and the gap may change sign within an element; however, the bonding is not continuous over the interface because each element is either bonded or unbonded (see below).

../../../../../_images/cmflatjoint_fig11.png

Figure 11: Variation of gap, normal stress, and bond state over a four-element, 2D flat-joint model corresponding with motion at top of image.

When the flat-joint model is installed at a contact, the following quantities are initialized:

(18)e:{F(e)n,F(e)ss,F(e)st,M(e)b}=E(e)μ=0.

When the first cycle occurs after the flat-joint model is installed at a contact, the following properties are fixed:

(19){Nr,Nα,go},R=λ{  min(R(1),R(2)),  ballball  R(1),  ballfacet

and initialized:

(20){gs=go,θbs=θbt=0}and{δss=δst=0}.

The first condition establishes the initial location and orientation of the notional surfaces relative to the pieces. The initial surface gap (go, go>0 is open) is the distance between the finite-size notional surfaces measured along the dotted line in Figure 2. The surface gap (gs=go+Δδn, gs>0 is open) is the cumulative relative normal displacement of the piece surfaces.

The force-displacement law for a flat-joint element consists of the following steps (see Figure 12 and Figure 13).

  1. Update F(e)n:

    (21)F(e)n=eσdA

    where σ is the interface normal stress of Equation (17), and the integration is performed over element (e). Analytical expressions for this integral (given in the section “Expressions for Element Normal Force and Bending Moment” below) are used to update F(e)n, which in turn, is used to update σ(e)=F(e)n/F(e)nA(e)A(e). If the element is bonded and the tensile-strength limit is exceeded (σ(e)>σc), then break the bond in tension (B(e)=1,{F(e)n,F(e)ss,F(e)st}=0), trigger the bond-break callback event, and skip the next step. If the element is unbonded and σ(e)0, then set the element shear force to zero ({F(e)ss,F(e)st}=0) and skip the next step.

  2. Update F(e)s as follows. Compute a trial element shear force:

    (22)F(e)s=F(e)sksA(e)Δδδ(e)s[F(e)ss=F(e)ssksA(e)Δδ(e)ss,F(e)st=F(e)stksA(e)Δδ(e)st]with  Δδδ(e)s=αΔδδ(e)s[Δδ(e)ss=αΔδ(e)ss,Δδ(e)st=αΔδ(e)st]α={g(e)1g(e)1g(e)0,unbondedandg(e)0>0andg(e)1<01,otherwise

    where Δδδ(e)s is the effective portion of the relative shear-displacement increment at the element centroid. If the element is bonded, then the entire increment is effective; if the element is not bonded, then only the portion of the increment that occurs while g(e)<0 is effective; Δδδ(e)s is Δδδs of Equation (1) evaluated at the element centroid and g(e) at the beginning and end of the timestep is denoted by g(e)0 and g(e)1, respectively. The trial element shear stress:

    (23)τ(e)=F(e)s/F(e)sA(e)A(e).

    The procedure now differs based on the bond state as follows.

    • Unbonded. The shear strength τ(e)c=crμσ(e). If |τ(e)|τ(e)c, then F(e)s=F(e)s; otherwise, enforce the shear-strength limit (F(e)s=τ(e)cA(e)(F(e)s/F(e)sF(e)sF(e)s)), which implies that slip has occurred and the slip energy is updated via Equation (28). If the slip state has changed, then trigger the slip-change callback event.
    • Bonded. The shear strength τ(e)c=cσ(e)tanϕ. If |τ(e)|τ(e)c, then F(e)s=F(e)s; otherwise, the shear-strength limit has been exceeded. Break the bond in shear (by setting B(e)=2, reevaluating F(e)n as in step 1 and triggering the bond-break callback event). If Mr=0 (shear drop to zero), then set {F(e)ss,F(e)st}=0. If Mr=1 (shear drop to residual), then set F(e)s=(crμσ(e))A(e)(F(e)s/F(e)sF(e)sF(e)s).
  3. Update M(e)b:

    (24)M(e)b=erσdA

    where σ is the interface normal stress of Equation (17), r is the moment arm with respect to the element centroid, and the integration is performed over element (e). Analytical expressions for this integral (given in the section “Expressions for Element Normal Force and Bending Moment” below) are used to update M(e)b.

../../../../../_images/cmflatjoint_fig12.png

Figure 12: Force-displacement law for an unbonded flat-joint element: (a) normal stress versus element gap; (b) shear stress versus relative shear displacement; and (c) slip envelope.

../../../../../_images/cmflatjoint_fig13.png

Figure 13: Force-displacement law for a bonded flat-joint element: (a) normal stress versus element gap; (b) shear stress versus relative shear displacement; and (c) failure envelope.

Energy Partitions

The flat-joint model provides two energy partitions:

  • strain energy, Ek, stored in the springs; and
  • slip energy, Eμ, defined as the total energy dissipated by frictional slip.

The strain energy in the flat joint is obtained by summing the strain energy in each element:

(25)Ek=eE(e)k.

The strain energy in each element is updated after the force-displacement law:

(26)E(e)k=12((F(e)n)2knA(e)+F(e)s2ksA(e)+(M(e)b)2knI(e)+(M(e)t)2ksJ(e))withI(e)=eη2dA={23tˉr3,2D(t=1)14πr4e,3DJ(e)=e(ξ2+η2)dA={0,2D(t=1)12πr4e,3Dre=A(e)/A(e)ππ

where I(e) is the moment of inertia of the element cross-section (about the line passing through x(e) and in the direction of θθb), J(e) is the polar moment of inertia of the element cross-section (about the line through x(e) and in the direction of ˆnc), ˉr is the element half-length, and re is the element effective radius. The integrals in Equation (26) are {evaluated exactly in 2D; and approximated in 3D} (by assuming that the element is a disk with an area equal to that of the element).

The slip energy in the flat joint is obtained by summing the slip energy in each element:

(27)Eμ=eE(e)μ.

The slip energy in each element is updated during the force-displacement law whenever the shear-strength limit has been exceeded via:

(28)E(e)μ:=E(e)μ+τ(e)cA(e)Δδδ(e)s

where the quantities are defined in the “Force Displacement Law” section.

Table 1: Flat-Joint Model Energy Partitions
Keyword Symbol Description Range Accumulated
energy‑strain Ek strain energy [0.0,+) NO
energy‑slip Eμ total energy dissipated by slip [0.0,+) YES

Properties

The properties defined by the flat-joint contact model are listed in the table below as a concise reference; see the “Contact Properties” section for a description of the information in the table columns.

Table 2: Flat-Joint Model Properties
Keyword Symbol Description Type Range Default Modifiable Inheritable
flatjoint   Model name          
fj_nr Nr Number of elements in radial direc. (total number of elements in 2D) INT [1, ∞) 2 NO [1] NO
fj_nal (3D) Nα Number of elements in circumferential direc. INT [3, ∞) 4 NO [1] NO
fj_rmul λ Radius multiplier FLT (0.0, ∞) 1.0 NO [1] NO
fj_gap0 g0 Initial surface gap FLT [0.0, ∞) 0.0 NO[1] ({go>0}{e:B(e)=0}) NO
fj_kn kn Normal stiffness [stress/disp.] FLT [0.0, ∞) 0.0 YES NO
fj_ks ks Shear stiffness [stress/disp.] FLT [0.0, ∞) 0.0 YES NO
fj_fric μ Friction coefficient FLT [0.0, ∞) 0.0 YES NO
fj_ten σc Tensile strength [stress] FLT [0.0, ∞) 0.0 YES NO
fj_coh c Cohesion [stress] FLT [0.0, ∞) 0.0 YES NO
fj_cohres cr Residual cohesion [stress] FLT [0.0, ∞) 0.0 YES NO
fj_fa ϕ Friction angle [degrees] FLT [0.0, 90.0) 0.0 YES NO
fj_resmode Mr Shear-drop residual mode {0,sheardroptozero1,sheardroptoresidual} INT {0,1} 0 YES NO
fj_elem e Element ( e ) — accesses element-based data (element numbering in Figure 5 and Figure 16) INT [1, NrNα] 1 YES NO
fj_emod E Effective modulus FLT [0.0, ∞) 0.0 NO NO
fj_kratio κ Normal-to-shear stiffness ratio, κkn/ks FLT [0.0, ∞) 0.0 NO [2] NO
fj_slip s(e) Slip state of element (e) {true,slippingfalse,notslipping} BOOL {true,false} false NO NO
fj_state B(e) Bond state of element (e) {0,unbonded1,unbonded&brokeintension2,unbonded&brokeinshear3,bonded} INT {0,1,2,3} 0 NO NO
fj_radius R Flat-joint radius FLT (0, ∞) NA NO (set via λ) NO
fj_gap gs Surface gap FLT R go NO NO
fj_relbr θθb Relative bend-rotation (θbs,θbt) (2D model: θbk=θbs,θbt0) VEC2 NA 0 NO NO
fj_cen x(e) Centroid location of element (e) VEC NA NA NO NO
fj_area A(e) Area of element (e) FLT (0, ∞) NA NO NO
fj_egap g(e) Gap at centroid of element (e) FLT R NA NO NO
fj_shear τ(e)c Shear strength [stress] at centroid of element (e) FLT [0.0, ∞) 0.0 NO (set via c and σc) NO
fj_sigma σ(e) Normal stress at centroid of element (e) FLT R NA NO NO
fj_tau τ(e) Shear stress at centroid of element (e) FLT [0.0, ∞) NA NO NO
fj_force ˜F Interface force VEC NA 0 NO NO
fj_moment ˜M Interface moment VEC NA 0 NO NO
fj_track   Tracking state for use when plotting contacts with the specific numeric value BOOL {true,false} false YES NO
fj_mtype   Initial microstructural type {1,bonded2,gapped3,slit4,undefined} INT {1,2,3,4} NA NO NO
user_area A Constant area [length*length] FLT (0.0,+) 0.0 YES NO

[1] Specify before cycling; cannot modify thereafter.

[2] If either the normal or shear stiffness is zero, then κ is zero.

Methods

Table 3: Flat-Joint Model Methods
Method Argument Symbol Type Range Default Description
area Set user_area to the area
deformability Set deformability
  emod E FLT [0.0,+ ∞) NA Effective modulus
  kratio κkn/ks FLT [0.0,+ ∞) NA Normal-to-shear stiffness ratio
initialize Initialize the forces and moment.
  force   VEC     Total local force
  Moment   VEC     Total local moment
bond Bond element ( e ) or all elements
  element         Element number ( e ) [1]
  gap G VEC2 R2 (,0] Gap interval
unbond Unbond element ( e ) or all elements.
  element         Element number ( e ) [1]
  gap G VEC2 R2 (,0] Gap interval

[1] Element numbering is shown in figures above.

Area

Set the user_area property via the current contact area. This operation means that the contact area stays constant and is fixed independent of changes to the piece sizes/geometries. In order for the stiffnesses to be recomputed accounting for this area, one should subsequently call the deformabilty method.

Deformability

The deformability provided by the flat-joint interface can be specified with the deformability method which sets:

(29)kn:=E/ELL,ks:=kn/knκκwithL={R(1)+R(2),ballballR(1),ballfacet.

The first term in this expression is obtained by equating the normal stiffness to the axial stiffness of the volume of material shown in Figure 6 of the “Linear Model” description.

The deformability of a homogeneous, isotropic, and well-connected granular assembly experiencing small-strain deformation can be fit by an isotropic material model, which is described by the elastic constants of Young’s modulus (E) and Poisson’s ratio (ν). E and ν are emergent properties that can be related to the effective modulus (E) and the normal-to-shear stiffness ratio (κkn/knksks) at the contact as follows: E is related to E with E increasing as E increases, and ν is related to κ with ν increasing up to a limiting positive value as κ increases. These relationships are obtained by specifying E and κ as the arguments of the deformability method.

Initialize

Initialize the local force and moment presuming all elements are bonded. This can be useful when inserting a flat-joint model in a preexisting contact in a system that has been cycled to equilibrium. The gaps of each element are set to provide the desired force. A separate moment is stored and added to each bonded element to produce the desired moment.

Bond

Bond element (e) if the contact gap between the pieces is within the bonding-gap interval. If no gap is specified, then the element is bonded if the pieces overlap. A single value can be specified with the gap keyword corresponding to the maximum gap. If the element number (e) is not specified, then all elements are bonded. Only elements with a zero gap over their entire surface can be bonded. If an element becomes bonded, then the element bond state becomes bonded (B(e)=3). The element force and moment are unaffected and will be updated during the next cycle. One can ensure the existence of contacts between all pieces with a contact gap less than a specified bonding gap (gb) by specifying gb with the proximity in the contact cmat default command of the Contact Model Assignment Table (CMAT).

Unbond

Unbond element (e) if the contact gap between the pieces is within the bonding-gap interval. If no gap is specified, then the element is unbonded if the pieces overlap. A single value can be specified with the gap keyword corresponding to the maximum gap. If an element number (e) is not specified, then unbond all elements. If the element becomes unbonded, then the element bond state becomes unbonded (B(e)=0). The element force and moment are unaffected and will be updated during the next cycle.

Callback Events

Table 4: Flat-Joint Model Callback Events
Event Array Slot Value Type Range Description
contact_activated Contact has become active
  1 C_PNT N/A Contact address
slip_change Slip state of element ( e ) has changed.
  1 C_PNT NA Contact address
  2     Element number ( e ) [1]
  3 INT {0,1} Slip change mode {0: slip has initiated1: slip has ended
bond_break Bond of element ( e ) has broken.
  1 C_PNT NA Contact address
  2     Element number ( e ) [1]
  3 INT {1,2} Break mode {1: broke in tension2: broke in shear
  4 FLT [0.0,+ ∞) Strain energy at onset of failure
broken All bonds are broken.
  1 C_PNT NA Contact address

[1] Element numbering is shown in figures above.

Stiffnesses for Timestep Estimation Scheme

Expressions for Element Normal Force and Bending Moment

Analytical expressions for the element normal force and bending moment (F(e)n and M(e)b) are presented here. These expressions are obtained by integrating the interface normal stress over the element. The element normal force and bending moment satisfy:

(30)F(e)n=eσdAM(e)b=erσdA

where σ is the interface normal stress of Equation (17), r is the moment arm with respect to the element centroid, and the integration is performed over element (e). These integrals are {evaluated exactly in 2D; and approximated in 3D}.

For the 2D model, the integrals are evaluated exactly as follows. The gap across the element (see Figure 14) is given by

(31)g(ζ)=g(e)a+(g(e)bg(e)a2ˉr(e))ζ,ζ[0,2ˉr]

where g(e)a and g(e)b are the values of the gap at the element ends (obtained by evaluating Equation (3) at these locations) and ˉr(e) is the element half-length.

../../../../../_images/cmflatjoint_fig14.png

Figure 14: Variation of interface gap over a typical 2D element.

The element is mapped into one of the following three cases (in order): (1) the gap changes sign within the element (g(e)ag(e)b<0); (2) the gap remains positive or zero (g(e)a0andg(e)b0); or (3) the gap remains negative (g(e)a0 and g(e)b0).

The integrals in Equation (30) are expressed in the element system by

(32)F(e)n=2ˉr0σtdζ(t=1)M(e)b=2ˉr0rσtdζ(t=1),r=ζˉr(e).

where the g(r) term of the interface normal stress of Equation (17) is expressed via Equation (31). These integrals are evaluated analytically for each case to yield (with the superscript (e) notation dropped for the ga, gb, and ˉr terms):

(33)case1:F(e)n=knˉrt{ga+gb,bondedg2agagb,unbondedandga<0g2bgbga,unbondedandga>0
(34)case1:M(e)b=knˉr2t{gagb3,bondedg2a(ga3gb)3(gbga)2,unbondedandga<0(gbga)3g2a(ga3gb)3(gbga)2,unbondedandga>0
(35)case2:F(e)n=knˉrt{ga+gb,bonded[3mm]0,unbondedM(e)b=knˉr2t{gagb3,bonded0,unbonded
(36)case3:F(e)n=knˉrt(ga+gb)M(e)b=knˉr2t(gagb3).

For the 3D model, it is assumed that the normal stress is constant over the element and equal to its value at the element centroid[4] such that

(37)F(e)n=eσdA={0,unbondedandg(e)0kng(e)A(e),otherwiseM(e)b=erσdA=0

where g(e) is the gap at the element centroid obtained by evaluating Equation (10) at the element centroid.

Usage and Verification Examples

The “Rock Testing” example application demonstrates compressive and tensile strength tests applied to simple BPMs using the linear parallel bond and flat joint contact models.

Endnotes

[1]A crack is a disk for the 3D model and a linear segment of unit-thickness depth for the 2D model. A crack has failure mode (tensile or shear) and geometric information (size, position, normal direction and gap). The size is two times the element radius (which is the radius giving the same area as the element), or the element length for the 2D model. The position, normal direction and gap are updated to correspond with material motion subsequent to bond breakage. Cracks are displayed by the crack-monitoring package (see section “Crack Monitoring” in [Potyondy2017]).
[2]The kinematic formulation approximates the motion of two notional surfaces, each of which is connected rigidly to a piece of a body. A limitation of the current formulation is that when two bodies joined by a flat-joint contact are moved tangentially to one another, the interface gap increases (because gs increases). It may be possible to remove this limitation by introducing an alternative interface coordinate system that does not rotate under the above condition.
[3]It is assumed that the shear stress arising from relative twist rotation is constant over the element and equal to its value at the element centroid. An inaccuracy introduced by this assumption is that the twisting moment is zero with respect to the element centroid. The interface shear stress arising from relative twist rotation actually varies linearly over the element, which makes the twisting moment nonzero with respect to the element centroid. As a result, the twisting moment at the center of the interface converges to the correct value from below as the discretization is refined. A similar assumption and associated inaccuracy exists for the bending moment of the 3D flat-joint model.
[4]An inaccuracy introduced by this assumption is that the bending moment is zero with respect to the element centroid. The interface normal stress actually varies linearly over the element, which makes the bending moment nonzero with respect to the element centroid. As a result, the bending moment at the center of the interface converges to the correct value from below as the discretization is refined. A similar assumption and associated inaccuracy exists for the twisting moment of the 3D flat-joint model.

Model Summary

An alphabetical list of the flat-joint contact model methods is given here. An alphabetical list of the flat-joint contact model properties is given here.

References

[Potyondy2017](1, 2, 3) Potyondy, D. O. “Material-Modeling Support in PFC [fistPkg25],” Itasca Consulting Group, Inc., Technical Memorandum ICG7766-L (March 16, 2017), Minneapolis, Minnesota.
[Potyondy2016]Potyondy, D. (2016) “Flat-Joint Contact Model [version 1],” Itasca Consulting Group, Inc., Minneapolis, MN, Technical Memorandum 5-8106:16TM47, October 12, 2016.
[Potyondy2013]Potyondy, D. (2013) “PFC3D Flat-Joint Contact Model (version 1),” Itasca Consulting Group, Inc., Minneapolis, MN, Technical Memorandum ICG7234-L, June 25, 2013.
[Potyondy2012a]Potyondy, D. (2012a) “PFC2D Flat-Joint Contact Model,” Itasca Consulting Group, Inc., Minneapolis, MN, Technical Memorandum ICG7138-L, July 26, 2012.
[Potyondy2012b]Potyondy, D.O. (2012b) “A Flat-Jointed Bonded-Particle Material for Hard Rock,” paper ARMA 12-501 in Proceedings of 46th U.S. Rock Mechanics/Geomechanics Symposium, Chicago, USA, 24–27 June 2012.