Stress‐controlled weakly periodic boundary conditions: Axial stress under varying orientations

The accuracy of multiscale modeling approaches for the analysis of heterogeneous materials hinges on the representativeness of the micromodel. One of the issues that affects this representativeness is the application of appropriate boundary conditions. Periodic boundary conditions are the most common choice. However, when localization takes place, periodic boundary conditions tend to overconstrain the microscopic problem. Weakly periodic boundary conditions have been proposed to overcome this effect. In this study, the effectiveness of weakly periodic boundary conditions in restoring transverse isotropy of representative volume elements (RVE) for a fiber‐reinforced composite with elastoplastic matrix is investigated. The formulation of weakly periodic boundary conditions is extended to allow for force‐controlled simulations where a uniaxial stress can be applied. A series of simulations is performed where the orientation of applied stress is gradually varied and the influence of this orientation on the averaged response is examined. An original method is presented to test the correlation between the ultimate principal stress and average localization angle of shear bands within an RVE. It is concluded that weakly periodic boundary conditions alleviate anisotropy in the RVE response but do not remove it.

which linear displacement boundary conditions, periodic boundary conditions, and uniform traction (Neumann) boundary conditions are the most straightforward options. Of these, displacement boundary conditions give an upper bound for the stiffness, while the uniform traction boundary conditions provide a lower bound and periodic boundary conditions are generally most accurate.
Strong periodic boundary conditions (PBC) are applied by defining a linear constraint for displacements on the edge of the RVE such that the difference in displacement at opposite edges exactly matches the applied strain. A periodic fluctuation is then permitted inside the RVE on top of the applied strain. PBC tend to become less accurate when localized deformations arise in the microscale problem. The fact that the orientation of localization bands is enforced to respect the periodicity constrains the problem and results in an overly stiff response.
Recently developed weakly periodic boundary conditions (WPBC) serve as a transition between strong periodicity and Neumann boundary conditions. 8,10 Periodicity is enforced in a weak sense with Lagrange multipliers rather than directly with linear constraints. By coarsening the discretization associated with the Lagrange multipliers, the periodic constraint is released, thereby allowing some nonperiodic localization and improving the stress and stiffness estimates. In the limit, when only a single Langrange multiplier is used per edge, WPBC become equivalent to Neumann boundary conditions. An additional advantage of WPBC is that they allow for nonmatching discretization at opposite edges of the RVE. This also gives rise to the option to define boundary conditions where periodic boundary conditions are aligned with a preferential direction. 10 If the orientation of the localization band is known, it is possible to apply a shift to the coordinates of one of the opposite edges in the definition of the boundary conditions such that localization at the specified orientation is supported.
The objective of this work is to investigate how effective WPBC are in restoring the representativeness of micromodels with localized deformation. RVEs are generated for fiber-reinforced composites where the constitutive model for the matrix includes perfect plasticity. The heterogeneous model with perfectly plastic limit behavior is prone to localization. The robustness of the plasticity model with respect to softening models makes it a suitable candidate for a systematic study with multiple RVEs and multiple load cases. The representativeness of the micromodels is examined in terms of transverse isotropy of the response. In absence of directional bias in the generated fiber distributions, the RVE response should be equal in all directions perpendicular to the fiber direction. Whether this is the case is tested by applying uniaxial stress to an RVE under different directions and comparing the averaged stress-strain response. A batch of different RVEs with random fiber distributions is used to eliminate effects that may be specific for a given fiber arrangement.
In order to be able to apply uniaxial stress conditions, this work introduces a novel modification to the WPBC framework which allows using WPBC in "stress-controlled" procedures. Uniaxial stress is applied to the RVE instead of strain, because the macroscale strain that would need to be applied to obtain the required stress is not known a priori. The stress-driven procedure presented herein provides a direct way to investigate the effectiveness of WPBC in alleviating anisotropy of micromodels with localized deformation without having to iterate on the macroscale strain until the sought-after stress is obtained.
The influence of load direction on the RVE response is statistically assessed. An original postprocessing procedure is introduced to detect the orientation of a localization band inside an RVE and the results of this procedure are used to explain the observed trends in averaged stress-strain response.
This article is organized as follows. In Section 2, the formulation of weakly periodic boundary conditions is reviewed after which it is extended in Section 3 to allow for force-controlled analysis. Section 4 presents the methodology of the study into the transverse isotropy of the RVE response. Results of this study are presented in Section 5.

Boundary value problem
This article is restricted to single scale analysis, focusing on what would be the subscale model in a computational homogenization approach. The subscale BVP is described by pointwise equilibrium on the domain Ω □ , with appropriate subscale boundary conditions on Ω □ ≡ Γ □ . The weak form of the equilibrium equation is given by: Find u and t that solves where is the first Piola-Kirchhoff stress, ∇ is the spatial gradient with respect to x in the reference configuration, f is the vector of body forces, u is the displacement field, and t ≜ ⋅n is the vector of boundary tractions. Finding a unique solution to Equation (1) requires removing translational and rotational rigid body motions (RBM). While u and t are a priori unknowns, the displacement u can be decomposed into: 8 where u is the macroscale displacement at x-the center of the RVE, typically coincident with a macroscale integration point, (u ⊗ ) ⋅ (x − x) is the first-order macroscale deformation component, and u s is the subscale fluctuation. Rigid body translations do not affect the stress resultant, and thus they are removed. The resulting weak form is as follows: Find u ′ ∈ U □ and t ′ ∈ T □ that solves where u ′ represents the non-RBM component of the displacement vector found in H 1 (Ω □ ), or the space of functions with square integrable gradients on Ω □ . The vector x 0 is an arbitrarily chosen point in Ω □ . In this work, x 0 is placed at the lower left corner of the RVE for ease of implementation. Likewise, t ′ represents the self-equilibrating traction forces on the boundaries of the RVE and belongs to L 2 (Γ □ ), which is the space of square integrable functions on Γ □ . Here, d denotes the dimension of the problem.

Strong periodic boundary conditions
In order to formulate strong periodic boundary conditions, Γ □ is split into an image boundary (Γ + □ ) and a mirror boundary Figure 1). The mapping per ∶ Γ + □ → Γ − □ mirrors any x + ∈ Γ + □ onto x − ∈ Γ − □ . That is, x − = per (x + ). Strong periodicity postulates a simultaneous variation of u ′ and t ′ , which allows the subscale fluctuation u s to exist on Γ □ : 2 where ⟦u ′ ⟧ □ ≜ u ′ (x) − u ′ ( per (x)). Equation (4) removes rigid body rotations and implies ⟦u s ⟧ □ = 0 on Γ + □ . In other words, the macroscale deformation gradient is prescribed on Γ □ in an average sense due to the strong periodicitiy requirement of u s on Γ □ .

Weakly periodic boundary conditions
Weakly periodic boundary conditions present a variational (weak) form of the constraint of periodicity of displacements presented in (4): 8 where T + □ is the trace of functions in T □ on the image boundary Γ + □ . Thence, all trial and test functions reside on Γ + □ . Using criterion (5), it follows that Using (7), Equation (3) is reformulated for weak microperiodicity as: Find u ′ ∈ U □ and t ′ ∈ T + □ for a given macroscale strain (u ⊗ ) that solves Note that Equations (6) and (8) must be solved simultaneously. The following notation changes are adopted for the sake of brevity: u ′ → u, t ′ → t, and T + □ → T □ . The system of Equations (6) and (8) can be interpreted as a minimization of total potential energy Π(u,t): where (u ⊗ ) is the elastic energy in the material, for which it holds that the derivative with respect to deformation gradient is the stress: / (u⊗∇) = .
Aligned weakly periodic boundary conditions are obtained by introducing a shift in the mapping x − = per (x + ).

Galerkin approximation
The Galerkin method entails the construction of finite-dimensional spaces of weighted residual approximations of the weak form: Both the weight and trial functions can be found within these spaces. Therefore, the BVP with weak microperiodicity is reduced to the Galerkin approximation: The reduction to finite-dimensional spaces enables the numerical computation of u h and t h as linear combinations of the basis vectors in U h □ and T h □ , respectively. In this article, the Bubnov-Galerkin approximation is used, which stipulates that u h and t h also originate from U h □ and T h □ .

STRESS-CONTROLLED WEAKLY PERIODIC BOUNDARY CONDITIONS
Unlike strong periodicity, weak periodicity is currently limited to strain-controlled procedures. Here, the WPBC formulation is expanded such that stress can be applied. Consider an RVE with strong periodic boundary conditions subject to prescribed strain. The macroscale stress tensor is given by: where f cor are the reaction forces at the corner nodes (x cor ) whose imposed displacement is given by (u ⊗ ) ⋅ (x cor − x). Instead of prescribing (u ⊗ ), it is possible to apply f cor to represent . In this case, relation (12) is employed in the formulation of the governing system of equations as shown in: 2 whereû + cor is the displacement of a corner x + cor on Γ + □ andû − cor = u( per (x + cor )) on Γ − □ . This concept is adapted for WPBC: apply f cor and solve iteratively forû ± cor . Note that it is not possible to apply the stress as constant tractions on the boundaries, because this would prevent prescribing periodicity of displacements. Instead, corner forces f cor are applied andû ± cor is solved for iteratively. Antiperiodicity of tractions (5) is enforced a priori, thus (11) may be evaluated over Γ + □ . The tying relations have to be in equilibrium, thus t + = −t + cor .
The vector ⟦x − x⟧ □ = ⟦x + cor − x⟧ □ remains constant over each image boundary of the RVE, thus: where f + cor are the corner forces that lie on Γ + □ . The sum of the corner forces on Γ + □ equals f − cor at x − cor , which equals x 0 herein. The potential energy (9) is updated by incorporating (12) and adding the potential energy term contributed by the corner forces in (14): Minimizing the potential energy yields the weak form of the governing equation. The Galerkin approximation follows: The Galerkin approximation is solved by discretizing the domain and boundary (Ω □ and Γ □ ) into finite-sized elements with known test functions ( u h and t h ). This involves the creation of a displacement mesh ( h ) and a traction mesh ( h ) with linearly independent shape functions n i (x) and h j (x) for each node i in  h and j in  h . The vectors u h and t h are now represented by linear combinations of the shape functions times their corresponding discrete nodal values, or where N and H are the matrices of shape functions associated with  h and  h , andû andt are vectors of nodal displacements and nodal tractions, respectively. The Bubnov-Galerkin approximation requires that the weight and trial functions originate from the same finite-dimensional space, hence The vectorû includes the boundary displacementsû + andû − on Γ + □ and Γ − □ . There are elements on the boundary of  h with shape functions contained in N + and N − that map said boundary displacements onto u h+ and u h-, thus The second-order tensors in (16) are rewritten as vectors using Voigt notation, such that ( u h ) ≡ u⊗∇ and (u h ) ≡ . The deformation matrix B and tangent stiffness operator D mapû onto anḋontȯ, respectively: The aforementioned equations are used to linearize (16) so it can be solved iteratively using Newton-Raphson methods. The following system of equations ensues: where f andû int represent the vectors of externally applied corner forces and displacements of internal nodes, respectively; and Note that when the displacement in the corner nodes is prescribed, elimination of Δû ± cor from the system of equations leads to the system as given in Section 2.
Weakly periodic boundary conditions are only effective if  h -the mesh associated with the Lagrange multipliers-is more coarsely discretized than  h on Γ □ . In this study,  h is constructed for piecewise linear and continuous tractions following the procedure in Reference 8. The standard traction mesh is created by projecting all nodes on Γ − □ and Γ + □ onto Γ + □ . Then,  h is coarsened by removing nodes that are closer than a given tolerance (dx) to other nodes until the desired dimension of T h □ has been achieved ( Figure 2). The level of coarsening is adjusted using coarsening factors (cf_). Here, cf_ is defined as the ratio of the smallest element dimension in the standard mesh (dx0) to dx.

METHODOLOGY
A methodology is presented to compare the performance of WPBC to PBC: Uniaxial stress is applied to a batch of RVEs at 46 angles between 0 • and 90 • to generate curves of ultimate principal stress ( 1 ) vs orientation angle ( ). Fibers are  Figure 3A shows the uniaxial stress ( x'x' ) being applied to an RVE at an angle from the global x-axis. The uniaxial stress x'x' is transformed into global coordinates and then into the corner forces shown in Figure 3A via where Δx and Δy are the RVE dimensions along the x and y axes, respectively. The BVP is solved and the resulting stresses and strains are transformed back into stresses and strains along the local x ′ -axis. The ultimate principal stress is denoted 1 . Figure 4 shows an example stress-strain curve with 1 highlighted. At each orientation angle , the value of 1 is recorded to generate a 1 -curve as in Figure 5. Each RVE yields a different 1 -curve due to the randomness in the fiber distribution.
This randomness is overcome by averaging the curves from the batch of RVEs, whose random fiber distributions have constant density. This results in a curve of average ultimate principal stress ( 1 ) vs . If the model maintains transverse isotropy, then 1 should be independent of .
In addition, a method for estimating the orientation of shear bands is presented in order to study their influence on 1 . An area-weighted average of the localization angles of shear bands ( ) is estimated using the orientation angles of plastic zones in an image of the RVE subject to 1 , as shown in Figure 6.

F I G U R E 8 Ultimate principal
stress vs orientation angle, strong periodic boundary conditions. Individual lines are from individual RVEs, the thick red line is the mean from all RVEs, and the error bars indicate the standard deviation MATLAB's Image Processing Toolbox is used to extract the blue channel from the RGB image of an RVE. An appropriately chosen threshold yields a black and white image whose blobs accurately portray the location of shear bands. The localization angle is defined as the smallest angle between a shear band and the x-or y-axis, such that lies between 0 • (vertical or horizontal bands) and 45 • (diagonal bands). The blobs' orientation angles ( or ∈ (−90 • ,90 • )) are transformed into localization angles ( loc ∈ (0 • ,45 • )) via The angle is estimated using the area-weighted average of loc of the blobs: where A i and loc,i are the area and localization angle of each blob, respectively. Figure 7 shows an example RVE subject to axial stress at an angle of = 24 • .

Strong periodic boundary conditions
The methodology was first executed on a batch of 48 fiber-matrix composite RVEs, with elastoplastic matrix and strong periodic boundary conditions. Uniaxial stress was applied to the RVE batch at 46 different angles between 0 • and 90 • with respect to the orientation of the RVE boundaries. Figure 8 shows the 1 -curves for each RVE. The curve of average ultimate principal stress ( 1 ) vs (shown in red) shows a dependency between 1 and .
For scrutiny, null hypothesis significance testing was used to verify whether the observed difference in 1 at = 44 • and = 0 • is truly significant. Indeed, the two-tailed Welch's t-test rejects the hypothesis that the sample means of 1 at = 44 • and = 0 • stem from populations with equal means with a p-value of 7.9⋅10 −19 .

Weakly periodic boundary conditions
The performance of WPBC is assessed against PBC. A subset of 38 RVEs * was subject to uniaxial stress applied at 46 different angles between 0 • and 90 • to generate a 1 -curve like the red curve on Figure 8. This process was repeated with 10 different coarsening factors (cf_), whose 1 -curves are shown in Figure 9. It is observed that WPBC fail to remove the dependency between 1 and . Even at the lowest cf_, the two-tailed Welch's t-test rejects the hypothesis that the sample mean of 1 at = 44 • and = 0 • stem from populations with equal means with a p-value of 3.7⋅10 −10 . However, low coarsening factors do reduce the difference between max( 1 ) and min( 1 ). It is also worth noting that increasing cf_ elevates the values of 1 throughout the range of , but more so toward = 0 • , by imposing more exacting periodicity conditions.
Most notably, RVEs tend to develop single vertical or horizontal shear bands ( ≈0 • ) when subject to axial loads oriented at ≈45 • (where 1 is lowest). Figure 10 shows the deformed shape and average localization angles of the two RVEs * 12 RVEs were removed because they diverged before reaching the ultimate principal stress.

F I G U R E 13
Averaged principal stress vs orientation angle with aligned WBPC and strong PBC with lowest 1 in (A) and the two RVEs with highest 1 in (B). This seems to indicate that vertical and horizontal shear bands that respect periodicity lead to lower 1 . That is, values of closer to 0 • lead to lower 1 .
The correlation between and 1 was tested using various coarsening factors. Figure 11 shows the average localization angle of each RVE plotted against its associated ultimate principal stresses 1 at the lowest coarsening factor, cf_ = 0.016. The strength of the correlation between 1 and is measured by Pearson's coefficient ( ) and its corresponding p-value, where p < = .05 disproves H 0 : = 0. For instance, Figure 11 shows = 0.395 and p = .014 at = 44 • , indicating a highly significant positive correlation between and 1 .
This process was repeated for all 10 coarsening factors, all of which displayed very similar results to those observed in Figure 11. This shows that load oriented at = 45 • is bound to yield lower 1 by favoring periodicity-abiding vertical and horizontal shear bands. It is gathered that: F I G U R E 14 Principal stress vs orientation angle with aligned WBPC and strong PBC from single RVE • The positive correlation between and 1 at = 44 • is highly significant, meaning that fiber distributions that accommodate horizontal or vertical shear bands tend to yield lower values of 1 .
• The correlation between and 1 at = 0 • is less significant, meaning that diagonal shear bands do not necessarily come with lower 1 values.
• Overall, the lower values of 1 around = 45 • can be attributed to the formation of localization bands under shear which respect periodicity. Figure 12 presents both 1 (in red) and (in blue) against for a single RVE. Three coarsening factors are selected to show their difference. The shear bands due to the lowest coarsening factor are displayed atop for some orientation angles. It seems that the location of shear bands is dictated by both and the fiber arrangement, as evidenced by their sudden transitions. For example, = 0 • to = 8 • displays shear bands at the corners whose wedges slide horizontally outward resulting in a quasi-periodic deformation. Likewise, = 16 • to = 32 • shows diagonal shear bands and = 40 • to = 48 • shows a single vertical shear band whose ensuing deformation respects periodicity.
The lowest coarsening factor leads to lower stresses by allowing the formation of shear bands, which do not respect periodicity, particularly toward = 0 • , where axial stress led to the formation of quasi-periodic corner shear bands.
Increasing cf_ leads to more stringent periodicity requirements. As such, the quasi-periodic deformation near = 0 • and the diagonal shear bands between = 16 • and = 32 • start to favor plastification of the entire matrix to satisfy periodicity. This yields higher 1 . Au contraire, the shear band near = 44 • already respects periodicity. Thus, increasing cf_ does not raise 1 considerably.
All in all, low values of 1 are obtained if both the fiber distribution and boundary conditions allow the formation of localization bands which are favorably aligned with respect to the load orientation angle (eg, = 45 • favors horizontal or vertical shear bands, which respect periodicity). Summarizing: • Shear bands form only if both the fibers and boundary conditions are favorably aligned with respect to the load orientation angle.
• Very low cf_ allows nonperiodic shear bands with low 1 throughout . However, single shear bands respecting periodicity tend to yield lower 1 .
• Increasing cf_ raises 1 toward = 0 • due to more stringent periodicity, which favors plastification of the entire matrix to satisfy periodicity.
• Increasing cf_ does not raise 1 as much near = 45 • because that load angle favors shear bands which abide periodicity of displacements.
• Increasing cf_ limits the range of over which localization with a single shear band is permissible. Few RVEs' fibers are aligned such that falls within this narrower range.

Aligned weakly periodic boundary conditions
Finally, simulations are performed where the WPBC are aligned to the preferential direction for localization, that is, at an angle of +45 • . A uniform traction mesh with 20 nodes per edge is used. Averaged results from all 50 RVEs are shown along with earlier results with strong and weakly periodic boundary conditions (cf_ = 0.25) in Figure 13. It is observed that the aligned WPBC are not improving the performance. Strong fluctuations appear in the curve, which, as the curve represents the mean from 50 different RVEs, indicates an effect that is repeated for a significant number of fiber distributions.
To elucidate the cause for this effect, a single RVE geometry is taken for which the fluctuations are particularly pronounced. Ultimate stress as a function of the applied load direction is plotted in Figure 14. It is observed that for some angles, the ultimate stress is indeed closer to that for = 0 • with aligned WPBC. Deformations are shown for = 58 • and = 66 • . The low ultimate stress at = 58 • indeed conforms with the desired nonperiodic localization band. For the high ultimate stress, no localization takes place. In Figure 15 the RVE geometry is visualized as surrounded with shifted copies of itself, which is a way to illustrate the microstructure that is implied with the aligned boundary conditions. The reason for the high ultimate stress at = 66 • is found in the fact that a geometric artifact is formed along the top/bottom boundary: neighboring fibers on one side of the boundary are interconnected by the positioning of the fibers on the other side of the boundary. Effectively, the microstructure is now changing as a function of the load orientation and, for some angles, an artificial reinforcement appears when fibers cut by the boundary are coincidentally grouped together. The significant fluctuation in the averaged curve shows that such artificial reinforcement is more likely to appear at some angles than at others, which is related to the fiber diameter and the RVE size. It is likely that this effect could be removed by using nonperiodic microstructures, then the averaged curve may become more flat, but individual RVEs will still suffer from angle dependence.
For plasticity simulations as performed here, there are in fact two preferred orientations for a localization band, that is, at +45 • and −45 • . Only one of the two can be supported with the aligned periodic boundary conditions, there that is a localization band at an angle of +45 • . For ∈ [0,90] this means that the shift is always performed along the top/bottom boundary. As a consequence, a horizontal localization band remains supported for all applied angles .
For use in full multiscale analysis, an additional issue with aligned WPBC is that the alignment angle should be known a priori. In this study, uniaxial loading under a predefined angle is used, for which case the desired orientation of localization is known a priori, but in general this is not the case.

CONCLUSIONS
A novel modification to weakly periodic boundary conditions-which allows for force-controlled simulations-was tested on a batch of 38 fiber-reinforced composite RVEs with elastoplastic matrix to verify whether they alleviate the dependency between 1 and observed when applying strong periodic boundary conditions. A methodology for the automated analysis of the localization angles provided further insight into the behavior of the RVEs.
Overall, WPBC alleviate transverse anisotropy in the RVE response but do not remove it. Low coarsening factors did not remove the dependency between 1 and , but they did reduce the difference between max( 1 ) and min( 1 ). It is also observed the low cf_ led to lower 1 throughout the range of The statistical study shows ample evidence of a high positive correlation between 1 and at = 44 • for all coarsening factors, indicating that values of near 0 • (vertical or horizontal shear bands) lead to lower 1 . The evidence of a weaker negative correlation at = 0 • was not as substantial. This shows that load oriented at = 45 • is bound to yield lower 1 by favoring periodicity-abiding vertical and horizontal shear bands.
Furthermore, increasing cf_ vastly reduces the range of over which localization with a single shear band is permissible. Thus, increasing cf_ raises 1 more toward = 0 • , where stricter periodicity favors plastification of the entire matrix, as opposed to near = 45 • , where the shear bands already tend to abide periodicity.
Finally, aligned weakly periodic boundary conditions have been explored as a solution to the direction dependence. It is found that this sometimes leads to a much lower ultimate stress, when a favorable localization band is supported. However, results are not consistently improved. For some alignment angles, geometric artifacts are formed along the boundary that increase the ultimate stress.