Computational Stability Analysis of Diffusion‐Deformation Processes at Finite Strains

In the present work, we investigate the onset of structural instabilities occuring in soft hydrophilic elastomers undergoing geometrically constrained swelling, which is characterized by large elastic deformations coupled with transient fluid diffusion. We model this phenomenon using a minimization‐based variational formulation having the deformation map and the fluid flux as the independent variables [1,3]. The space‐time‐discrete form of the variational problem is implemented using a conforming Raviart–Thomas‐type finite element formulation, which yields a symmetric and positive definite global tangent matrix for a stable deformation state. The onset of instability is captured as the point where the tangent matrix loses its positive definiteness [4,5]. To validate the proposed theoretical framework, we consider the surface wrinkling phenomenon observed in constrained hydrogel bilayers subjected to fluid diffusion. Variations in the critical buckling load and the corresponding mode number with respect to the geometry and material parameters of the bilayer are studied. These studies could, for instance, serve as potential design guidelines for applications of such materials in microsensors or in the fabrication of soft microgears for actuators, where the buckling modes are exploited to produce a desired mechanical output [8].


Minimization-based Formulation of Diffusion-driven Finite Elasticity
Diffusion-driven finite elasticity is modeled as a two-field problem consisting of the deformation ϕ and fluid volume flux H as global primary fields, with F := ∇ϕ. Additionally, we have the swelling-fraction field s representing the local relative change in volume of the continuum due to species diffusion ϕ : Given a constitutive response in terms of the free energyψ (F , s) and dissipation potentialφ (H; F , s) functions, a rate-type potential functional [1,3] is formulated for the coupled problem in terms of the deformation rate and fluid-volume flux fields For a pure Dirichlet problem, the contribution P ext (φ, H) vanishes. Applying an implicit Euler time integration scheme to the governing balance law for fluid concentrationṡ = −Div [H] and to (2) over the time step τ := t n+1 − t n , the primary fields ϕ and H are obtained by a two-field minimization of the resulting incremental potential

Finite Element Implementation and Structural Stability Analysis
In order to obtain a conforming global finite element design, the variational structure (3) requires the ansatz functions for the volume flux H to be in the space H (Div, B). To fulfill this, we make use of the Q 1 RT 0 finite element for the spatial discretization. Here, the ansatz for H constitutes scalar flux degrees of freedom at edge centers interpolated by vectorial Raviart-Thomas shape functions. This ensures the continuity of the normal component of H across element boundaries and thus fulfills the criterion for H (Div, B) conformity [7]. A chemo-mechanical state ω 1 := {ϕ 1 , H 1 } is considered to be stable if it satisfies the criterion Π τ (ω 2 ) − Π τ (ω 1 ) > 0 for a state ω 2 := ω 1 + δω in the neighbourhood of ω 1 [5]. A Taylor-series expansion of Π τ (ω 2 ) about ω 1 leads to the inequality In a space-discrete setting, (4) translates to the inequality δd T K (d 1 ) δd > 0 demanding the positive definiteness of the global tangent matrix K for a stable state of deformation d 1 [4]. The system becomes structurally unstable and bifurcates to an alternate deformed state when K loses its positive definiteness.
As constitutive inputs, we adopt a free-energy functionψ =ψ mech (F ) +ψ chem (s) +ψ coup (det F , s) and a convex dissipation potentialφ (H; C n , s n ) = (2M s n ) −1 C n : H ⊗ H , refer [1]. The mechanical part of the free energy is modeled using a standard neo-Hookean function and for the chemical part, a Flory-Rehner type function based on statistical molecular thermodynamics is adopted [2]. The coupling term is a penalty enforcement of the incompressibility constraint det [F ] = 1 + s .
To validate the proposed theoretical framework, we investigate the onset of sinusoidal wrinkles on the surface of a representative hydrogel bilayer, having length L and constituting a stiff film of thickness w bound to a compliant substrate of thickness H. Fluid is injected into the bilayer uniformly and at a constant flow rate from the top film surface. Mechanical boundary conditions are imposed such that the bilayer, upon absorbing the fluid, is free to swell only along the direction normal to the film surface. We define the critical fluid volume v c as the total volume of fluid accumulated inside the bilayer specimen for which the longitudinal compression due to the mechanical constraints is sufficient enough to trigger the formation of sinusoidal wrinkles on the film surface. This together with the number of wrinkles N c is investigated for a wide range of shear moduli ratios γ f /γ s and film apect ratios w/L. Following [1], a set of realistic material parameters are adopted. Plots a) and b) correspond to a fixed modulus ratio γ f /γ s = 8 whereas a film aspect ratio w/L = 0.02 is used for plots c) and d).
The plot in Fig. 1a is the curve enveloping a family of isomodes and bears resemblence to the dependence of the buckling load of a plate under uniaxial compression on its aspect ratio [6]. There is a continuous shift between the isomodes along this curve owing to the existence of an alternate buckled configuration at a lower critical load. Increasing the relative film stiffness causes surface wrinkling to be activated at a lower critical volume as seen in Fig. 1c. The wrinkle count N c is observed to be a decreasing function with respect to both w/L and γ f /γ s . This is to be expected as the film tends to develop more resistance to being bent into wrinkles when its thickness and shear modulus are increased. hh Fig. 2: Buckling modes of the hydrogel bilayer for a few selected film aspect ratios. The trend depicted in Fig. 1b is clearly evident.