Analysis
The file src/Analysis.jl and files contained in the src/analysis/ directory are referred to as the analysis module. All finite element analysis is carried out here, including
- determination of the membrane geometry and dynamics
- calculation of the residual vector and tangent matrix
- implementation of Newton–Raphson iteration, and time stepping
- calculation of the pull force
Overview
Though several different mesh motions are implemented, we provide details for only the most involved one: the $\texttt{ALE-vb}$ mesh motion. At any time step, we begin with the vector of known degrees of freedom $[\mathbf{u} (t)]$, and seek to calculate the unknown degrees of freedom $[\mathbf{u} (t + \Delta t)]$. In practice, we work primarily with the degrees of freedom corresponding to the various unknowns—in this case, the velocities $[\mathbf{v} (t)]$, mesh velocities $[\mathbf{v}^{\text{m}} (t)]$, surface tensions $[\bm{\lambda} (t)]$, and mesh pressures $[\mathbf{p}^{\text{m}} (t)]$.
The direct Galerkin expression for the $\texttt{ALE-vb}$ mesh motion is provided in §2.7.3 of our manuscript (Sahu, 2024). Upon discretization, the direct Galerkin expression is satisfied by solving the set of equations
\[[\mathbf{r}^\lambda] \, = \, [\mathbf{0}] ~, \qquad [\mathbf{r}^{\text{v}}] \, = \, [\mathbf{0}] ~, \qquad [\mathbf{r}^{\text{m}}] \, = \, [\mathbf{0}] ~, \qquad \text{and} \qquad [\mathbf{r}^{\text{p}}] \, = \, [\mathbf{0}] ~,\]
where each vector corresponds to the portion of the residual vector $[\mathbf{r}]$ associated with a particular unknown. The unknown degree of freedom vector $[\mathbf{u} (t + \Delta t)]$ then satisfies
\[[\mathbf{r} ( [\mathbf{u} (t + \Delta t)] ) ] \, = \, [\mathbf{0}] ~.\]
The above equation is solved with the Newton–Raphson method, in which our initial approximation for $[\mathbf{u} (t + \Delta t)]$, denoted $[\mathbf{u} (t + \Delta t)]_0$, is chosen to be $[\mathbf{u} (t)]$. Successive approximations are determined according to
\[[\mathbf{u} (t + \Delta t)]_{j+1} \, = \, [\mathbf{u} (t + \Delta t)]_j \, - \, [\mathbf{K}]_j^{-1} [\mathbf{r} ( [\mathbf{u} (t + \Delta t)]_j ) ] ~,\]
where the global tangent diffusion matrix at the $j^{\text{th}}$ iteration is given by
\[[\mathbf{K}]_j \, := \, \dfrac{\partial [\mathbf{r}]}{\partial [\mathbf{u}]} \bigg\rvert_{[\mathbf{u} (t + \Delta t)]_j} ~.\]
Note that the diffusion matrix is calculated numerically by extending $[\mathbf{r}]$ and $[\mathbf{u}]$ into the complex plane [Lyness and Moler (1967), Lyness (1968)]. Details of each portion of our calculation can be found below.
GeoDynStress.jl
To calculate the residual vector, integrals over the parametric domain are evaluated as the sum of integrals over each finite element. The latter are numerically calculated via Gauss point integration, for which (see GaussPoint.jl and GpBasisFn.jl)
\[\int_\Omega f(\zeta^\alpha) ~ \text{d}\Omega \ = \ \sum_{e=1}^{\texttt{nel}} \int_{\Omega^e} f(\zeta^\alpha) ~ \text{d}\Omega \ = \ \sum_{e=1}^{\texttt{nel}} \sum_{k=1}^{\texttt{ngp}} w_k f(\zeta^{\alpha, e}_k) ~.\]
To simplify such calculations, the struct GeoDynStress contains all quantities at a single Gauss point that are used in the determination of the residual vector. Many of the quantities are self-explanatory; here we describe the few that are not.
To begin, we use Voigt notation to represent 2×2 matrices as 3×1 vectors, where we take advantage of symmetries in the stresses, couple stresses, and arbitrary variations. Thus,
𝞂is a3×1vector containing $(σ^{1 1}, σ^{2 2}, 2 σ^{1 2})$
𝞂mis a3×1vector containing $(σ^{1 1}_{\text{m}}, σ^{2 2}_{\text{m}}, 2 σ^{1 2}_{\text{m}})$
𝗠is a3×1vector containing $(M^{1 1}, M^{2 2}, M^{1 2} + M^{2 1})$
In addition, we introduce the matrices 𝗕a and 𝗕b containing products of the basis functions with surface geometry. More precisely, we have
\[[\mathbf{B}^a] \, = \, \bigg[ ~ [\mathbf{B}^a_1] ~ [\mathbf{B}^a_2] ~ \ldots ~ [\mathbf{B}^a_{\text{nen}}] ~ \bigg] \qquad \text{and} \qquad [\mathbf{B}^b] \, = \, \bigg[ ~ [\mathbf{B}^b_1] ~ [\mathbf{B}^b_2] ~ \ldots ~ [\mathbf{B}^b_{\text{nen}}] ~ \bigg] ~,\]
where
\[[\mathbf{B}^a_j] \, := \, \begin{bmatrix} [\bm{a}_1]^{\text{T}} N_{j, 1} \\[4pt] [\bm{a}_2]^{\text{T}} N_{j, 2} \\[4pt] ( [\bm{a}_1]^{\text{T}} N_{j, 2} \, + \, [\bm{a}_2]^{\text{T}} N_{j, 1} ) / 2 \\[2pt] \end{bmatrix} \qquad \text{and} \qquad [\mathbf{B}^b_j] \, := \, \begin{bmatrix} [\bm{n}]^{\text{T}} N_{j; 1 1} \\[4pt] [\bm{n}]^{\text{T}} N_{j; 2 2} \\[4pt] [\bm{n}]^{\text{T}} \, ( N_{j; 1 2} \, + \, N_{j; 2 1} ) / 2 \\[2pt] \end{bmatrix}\]
As we will see, determining these quantities greatly simplifies our later finite element calculations.
MembraneAleFem.GeoDynStress — TypeGeoDynStress(xms, cps, dofs, N, ∂Nα, ∂∂Nαβ, p::Params)Container for all geometric, dynamic, and stress-related quantities at a single point $(ζ^1, ζ^2)$.
Accessible fields:
| code | symbol |
|---|---|
𝘅 | $\bm{x}$ |
𝗮_α | $\bm{a}_α$ |
∂∂𝘅αβ | $\bm{x}_{, α β}$ |
a_αβ | $a_{α β}$ |
a△αβ | $a^{α β}$ |
𝗮△α | $\bm{a}^α$ |
JΩ | $J_{\Omega}$ |
Γ△μ_αβ | $\Gamma^\mu_{\alpha \beta}$ |
𝗻 | $\bm{n}$ |
b_αβ | $b_{α β}$ |
H | $H$ |
K | $K$ |
𝘃 | $\bm{v}$ |
∂𝘃α | $\bm{v}_{, α}$ |
∂∂𝘃αβ | $\bm{v}_{, α β}$ |
λ | $λ$ |
pm | $p^{\text{m}}$ |
𝘃m | $\bm{v}^{\text{m}}$ |
∂𝘃mα | $\bm{v}^{\text{m}}_{, α}$ |
∂∂𝘃mαβ | $\bm{v}^{\text{m}}_{, α β}$ |
𝞂 | $\langle \bm{\sigma} \rangle$ |
𝞂m | $\langle \bm{\sigma}^{\text{m}} \rangle$ |
𝗠 | $\langle \bm{M} \rangle$ |
𝗕a | $[\mathbf{B}^a]$ |
𝗕b | $[\mathbf{B}^b]$ |
MembraneAleFem.get_dof_cps — Functionget_dof_cps(cps, dofs, dof_list)Get control point values for dof_list, which may not be contained in dofs.
Allows one to evaluate various physical quantities (e.g. $\bm{v}$, $λ$) even if some of the associated control points are not degrees of freedom—and are thus not solved for.
FiniteElement.jl
MembraneAleFem.time_step! — Functiontime_step!(mesh, xms, cps, Δt, p::Params; args...)Given a known state of the membrane, determine the state a time Δt later.
The Newton–Raphson method is employed to advance forward in time, which uses the result of calc_r_K.
MembraneAleFem.calc_r_K — Functioncalc_r_K(mesh, xms, cps, time, Δt, p::Params; args...)Calculate the residual vector r and diffusion matrix K.
The global tangent diffusion matrix K is calculated via numerical differentiation, in which the residual vector and its argument are extended into the complex plane Lyness and Moler (1967), Lyness (1968).
MembraneAleFem.calc_elem_residual — Functioncalc_elem_residual(mesh, el_id, xms_el, cps_el, p::Params)Calculate the area element residual vector by looping over Gauss points.
MembraneAleFem.calc_elem_dof_residuals — Functioncalc_elem_dof_residuals(mesh, el_id, xms_el, cps_el, p::Params)Calculate the area element residual vectors for each degree of freedom by looping over Gauss points.
The discretized form of the residual vectors can be read directly from this function.
MembraneAleFem.calc_bdry_element_residual — Functioncalc_bdry_element_residual(mesh, bdry, ntype, nval, ..., p::Params; args...)Calculate the boundary element residual vector by looping over Gauss points.
MembraneAleFem.update_xms! — Functionupdate_xms!(motion, xms, cps, Δt, dofs)Update membrane positions according to the mesh velocity, using forward Euler.
MembraneAleFem.calc_τ_ν — Functioncalc_τ_ν(bdry::Boundary, 𝗮_α::Matrix{ComplexF64}, 𝗻::Vector{ComplexF64})Determine the in-plane unit vectors τ and ν on the boundary.
PullForce.jl
Specific calculations related to determining the force during tether pulling, as detailed in Appendix C of Sahu (2024).
MembraneAleFem.get_pull_el_id — Functionget_pull_el_id(numel::Int64)Return the id of the element that will be pulled, given the number of elements.
MembraneAleFem.get_adj_maps — Functionget_adj_maps(num1el::Int64, numel::Int64, IX::Matrix{Int64})For each element adjacent to the pulled one, return local ids of pulled nodes.
These mappings are required to determine the residual vector of pulled nodes, even though the nodes are treated as Dirichlet boundary conditions and thus excluded from the list of degrees of freedom.
MembraneAleFem.calc_pull_force — Functioncalc_pull_force(mesh, xms, cps, adj_el_ids, adj_node_map, p::Params)Return the calculated pull force 𝓕, as a 3×1 vector.
Analysis.jl
MembraneAleFem.run_analysis! — Functionrun_analysis!(mesh, xms, cps, p::Params; args...)Execute the finite element analysis and step forward in time through args[:Δts].
At every time step, carry out Newton–Raphson iteration to advance to the next time step. The global residual vector and tangent diffusion matrix are generated at every iteration, and mesh positions are updated according to the mesh velocities.