{'fraction_non_alphanumeric': 0.03840604026845638, 'fraction_numerical': 0.0173993288590604, 'mean_word_length': 4.38060846799675, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A finite membrane element formulation for surfactants

Surfactants play an important role in various physiological and biomechanical applications. An example is the respiratory system, where pulmonary surfactants facilitate the breathing and reduce the possibility of airway blocking by lowering the surface tension when the lung volume decreases during exhalation. This function is due to the dynamic surface tension of pulmonary surfactants, which depends on the concentration of surfactants spread on the liquid layer lining the interior surface of the airways and alveoli. Here, a finite membrane element formulation for liquids is introduced that allows for the dynamics of concentration-dependent surface tension, as is the particular case for pulmonary surfactants. A straightforward approach is suggested to model the contact line between liquid drops/menisci and planar solid substrates, which allows the presented framework to be easily used for drop shape analysis. It is further shown how line tension can be taken into account. Following an isogeometric approach, NURBS-based finite elements are used for the discretization of the membrane surface. The capabilities of the presented computational model is demonstrated by different numerical examples - such as the simulation of liquid films, constrained and unconstrained sessile drops, pendant drops and liquid bridges - and the results are compared with experimental data.


Introduction
Surface tension, which is due to the attraction of molecules at a fluid interface by the bulk of the fluid, leads to a minimization of the fluid surface area. However, if molecules of a surface active agent, usually called surfactant, are added to the fluid interface, the surface tension decreases. Surfactants, which are amphiphilic compounds, have a hydrophobic tail, which is oil-soluble and allow them to stay on the surface, and a hydrophilic head, which is water-soluble and is able to reduce the surface tension by disrupting hydrogen bonds. This function is particularly crucial for the lung biomechanics as pulmonary surfactants facilitate the breathing and reduce the possibility of airways being blocked by lowering the surface tension as the lung volume decreases (Nkadi et al., 2009). By adding lung surfactants to pure water with a surface tension of about 70 mN/m, the surface tension can decrease to values lower than 2 mN/m. Pulmonary surfactants are composed of roughly 80% phospholipids (PL), 5-10% neutral lipids (NL), mainly cholesterol, and 8-10% surfactants proteins (SP) (Goerke, 1998). Around half of the PL content is dipalmitoylphosphatidylcholine (DPPC), which is the strongest surfactant compound in the mixture of the pulmonary surfactants and has the largest contribution to the reduction of surface tension. Although a solution of pure DPPC can reduce the surface tension to zero during monolayer compression, it does not respread properly during expansion and adsorbs very slowly to the interface; therefore, the DPPC molecules are lost during cyclic compression and expansion of the liquid-air interface. In fact, other components of the pulmonary surfactants also contribute to the dynamic characteristics of lung surfactants. Surfactant proteins, particularly SP-B and SB-C, among other functionalities, improve the adsorption rate of DPPC molecules during the expansion of the interface (Goerke, 1998 andVeldhuizen andHaagsman, 2000) and prevent the collapse of the surfactants mixture film (Ingenito et al., 1999). There are many references on the link between deficiencies in the components of the pulmonary surfactants and respiratory distress syndrome (RDS) in adults and infants (e.g Gregory et al., 1991;Lewis et al., 1993;Hallman et al., 2001 andMa andMa, 2012). The pulmonary surfactants have crucial contributions to the physiology and biomecanics of the respiratory system: They accelerate the oxygen diffusion through the air-water interface (Olmeda et al., 2010); they enhance lung immunity and its resistance against pulmonary infections (Han and Mallampalli, 2015) and foremost they play the key role in lung stability by increasing the stability of alveolar surfaces .
Computational models of the respiratory system can provide new insights to these mechanisms and offer a path to developing better treatment of lung diseases such as RDS. Therefore, it is necessary to develop a computational framework to study the dynamics of pulmonary surfactants. With this motivation in mind, the objective of this paper is to present a new finite membrane element formulation that allows for concentration-dependent surface tension like in pulmonary surfactants. Although the constitutive models presented here are developed for pulmonary surfactants, other fluids with dynamic surface tension can be modeled similarly.
Compared to earlier works, this paper presents a dynamic finite membrane element formulation for liquids that: 1) allows for different concentration-dependent constitutive laws for dynamic surface tension of the liquid-gas interface; 2) can be coupled with structural finite membrane, shell and solid elements; 3) can model thin films, bubbles, liquid droplets and menisci; 4) is able to treat different external forces (like gravity), contact constraints and various boundary conditions like contact lines with specified contact angles; 5) can include the effect of line tension if needed.
Furthermore, the presented formulation is combined with the computational biological membrane model of Roohbakhshan et al. (2016) in order to model surfactant-lined alveolar tissue in Roohbakhshan (2018).
The rest of this paper is organized as follows: In Sec. 2, the existing empirical, theoretical and computational approaches towards describing the dynamic surface tension are briefly reviewed. Sec. 3 summarizes the liquid membrane theory developed by Sauer et al. (2014) and Sauer (2014), which includes the treatment of contact lines for quasi static droplets and menisci on the rigid planar substrate, and extends it by dynamic surface tension and contact line tension . In Sec. 4, the two existing dynamic surface tension models, proposed by Otis et al. (1994) and Saad et al. (2010), are adapted to general membranes. These models allow for dynamic changes of surface tension due to the variation of surfactants concentration. Sec. 5 is devoted to the finite element solution, where the time discretization and integration of the dynamic model is discussed. Sec. 6 presents numerical examples, which show the performance and robustness of the formulation. Sec. 7 concludes the paper.

Dynamic surface tension
In this section, first a concise history of the empirical methods for dynamic measurement of the surface tension is presented; then, the existing theories and computational models are briefly reviewed.

Experimental methods
In general, the existing techniques for dynamic measurement of surface tension can be grouped into three main categories: 1) Force method, 2) drop shape and pressure method and 3) flow method (Franses et al., 1996). In the force method, a solid is inserted in the liquid-gas interface, which creates a meniscus confined by a boundary with a specific contact angle on the solid surface. By balancing the forces at the contact line, the surface tension of the liquid-gas interface can be found. The typical examples for the force method are the Wilhelmy plate method and the du Noüy ring method (Butt et al., 2006). In the shape and pressure method, the well known Young-Laplace equation, ∆p = 2 H γ (1) which relates the pressure jump across the interface ∆p to the interface surface tension γ and the interface mean curvature H, is inversely solved based on the captured shape of an interface. The contact angle is included as a boundary condition for Eq. (1). If gravity is not neglected, ∆p and accordingly H vary across the interface. The Bond number, which defines the ratio of gravitational to interfacial forces, is usually used to evaluate the importance of gravity. Here ∆ρ is the density difference of the media on the two sides of the interface, g is the gravitational constant and L is a characteristic length, which can be for example the radius of a drop. In order to ignore gravity, Bo 1. In general, the drop shape method can be further grouped into (Saad and Neumann, 2016) a) Volume-radius limited, e.g. pendant drop (PD) and constrained sessile drop (CSD); b) Volume-angle limited, e.g. sessile drop (SD) and captive bubble (CB); c) Volume-radius-radius limited, e.g. 2-edges-constrained liquid bridge 2 (CLB), and d) Volume-radius-angle limited, e.g. 1-edge-constrained liquid bridge (LB).
In the flow method, the interface is controlled in a way that the shape analysis in not required to find the interface profile and numerical methods are not needed for solving the Young-Laplace equation (1). The well know examples for the flow method are the maximum bubble pressure method (MBPM), growing drop method and oscillating jet method. For instance, in the MBPM and growing drop method, the interface is restricted to be sections of a sphere with radius R, so that H = 2/R and Eq. (1) can be solved analytically as γ = ∆p R/2. To keep the interface close to a spherical shape, the gravity effect should be negligible, which restricts these methods to the cases where the Bond number is small. As it is shown in Sec. 6, the presented computational framework can simulate many of the existing empirical methods that are commonly used to measure dynamic surface tension.
Dynamic measurement of the surface tension of the pulmonary surfactants goes back to the early work of Clements (1957), who used the Langmuir-Wilhelmy balance for this purpose. Since then, the method has been used for many years to present an isothermal process expressing the surface tension vs. surface area relationship of lung extracts (e.g. Hills, 1985). Although this tool is very suitable to spread a well-defined monolayer of lipids and proteins on the air-liquid interface and the surface area can be varied simply by moving a barrier on the interface, it is susceptible to liquid leakage and it is slow compared to the dynamics of the respiratory system (Veldhuizen and Haagsman, 2000). Being a successful alternative of the Langmuir-Wilhelmy balance, the pulsating bubble surfactometer (PBS) works based on the formation of a gas bubble inside a chamber of liquid. The surface tension is found by solving the Young-Laplace equation assuming that the bubble is quasi spherical. Although it has been prevalently used since its introduction by Enhorning (1977), it has technical problems like liquid leakage and inaccuracy at low surface tensions (Veldhuizen and Haagsman, 2000), which are improved in a modified version, namely captive bubble surfactometer introduced by Schürch et al. (1989). Recently, drop shape methods have become more popular for the dynamic measurement of pulmonary surfactants. As these methods are based on the inverse solution of the Young-Laplace equation (1) by means of numerical techniques, they are more flexible and can include external forces like gravity and different boundary conditions like specified contact angles (Saad and Neumann, 2016). Usually, the experiments are set up so that the drops are axially symmetric; thus, axisymmetric drop shape analysis (ADSA) techniques are used, which either use regression methods to fit the experimental drop profile to a Laplacian curve (e.g. Rotenberg et al., 1983) or solve an axisymmetric 2D variation of the Young-Laplace equation (1), which depends only on the arc length along the interface, through an iterative optimization procedure (e.g. Saad, 2011 andSaad andNeumann, 2016). The method is combined with the constrained sessile drop (CSD) test (Saad et al., 2010(Saad et al., , 2012 and it is also used for the pendant drop (PD) test by circular approximation of the droplet profile .

Theoretical models
Since the 1970s, several theoretical models are introduced to explain the diffusion-adsorption process (Horn and Davis, 1975;Otis et al., 1994;Morris et al., 2001;Krueger andGaver, 2000 andSaad et al., 2010), through which the molecules of surfactants are transferred from the bullk to the interface. All these models assume that the surface tension is determined only by the amount of surfactants at the interface. Surfactants are transfered to the interface in two phases: First, surfactants are conveyed from the liquid bulk to the layers adjacent to the interface through a diffusion process. Second, the surfactants are adsorbed from the layers to the interface. As diffusion and adsorption occur at different rates, most of the existing models consider either diffusion or adsorption as the key mechanism. Thus, in general, the existing theoretical models can be grouped into the diffusion-controlled (e.g. Loglio et al., 1991) and adsorption/desorptioncontrolled models (e.g. Otis et al., 1994;Morris et al., 2001;Krueger andGaver, 2000 andSaad et al., 2010). In the former approach, it is assumed that the diffusion of surfactants molecules through the bulk is the rate limiting process while, in the latter, the adsorption and desorption are assumed to be more dominant. In Sec. 4, the adsorption/desorption-controlled models of Otis et al. (1994) and Saad et al. (2010) are introduced in detail and it is shown how they can be incorporated into the presented dynamic finite element formulation.

Computational models
Although dynamic surface tension of pulmonary surfactants and lung extracts has been the subject of many empirical and theoretical studies in the last decades, there are few computational models that address the concentration-dependent surface tension of the pulmonary surfactants. In the earlier approaches, the surface tension vs. surface area curves at constant temperature obtained by experiments are used to define the surface tension as a function of surface area. For instance, Karakaplan et al. (1980) propose an alveolar model discretized by triangular planar finite elements that combines the elastic properties of the alveolar wall and the surface tension of the pulmonary surfactants. They assume that the behavior of the surfactants is defined by two experimentally-fitted functions, one for inflation and another for deflation, that depend on the changes of local area. In the same fashion, Kowe et al. (1986), use an exponential form to add the triangular planar surface tension elements to the earlier model of Dale et al. (1980) who use elastic pin-joined bar elements to build a 3D alveolus model approximated by a truncated octahedron consisting of fibers only. To distinguish between inflation and deflation, two similar exponential functions with different parameters to be fitted are used. The single alveolus model of Kowe et al. (1986) is extended by Schroter (1995, 1997) to an alveolar model constructed by an assembly of 36 truncated octahedra. Kojić et al. (2006Kojić et al. ( , 2009) introduce a scaling approach to consider the hysteresis effect of surface tension in their proposed material model of biological membranes. Similar to other approaches, they also find the surface tension through a quasi static function derived from a single surface tension vs. surface area hysteresis curve inferred from an experiment. However, as the experiments are restricted to a fixed range of area changes, they introduce a scaling method to be able to stretch their biological membranes to ranges that are different from a particular experiment. The computational models briefly introduced above are simple to implement but they are all time-independent; insensitive to the rate of changes in surface area; quasi-static rather dynamic and, excluding the models of Kojić et al. (2006Kojić et al. ( , 2009, they are restricted to the surface area range dictated by the experiment, from which the parameters are identified. In contrast to the time-independent approaches, there are time-dependent computational models that explicitly consider the exchange of surfactants between the liquid bulk and the liquid-gas interface. The very first model of this kind is the mathematical model of Archie (1973) that assumes a linear adsorption/desorption between the interface and the bulk. Accordingly, he assumes that the surfactants concentration is inversely proportional to the surface area and the surface tension is inversely proportional to the concentration. Schroter (2000, 2006) modify the earlier work of Schroter (1995, 1997) by using the dynamic surface tension model of Otis et al. (1994) instead of the quasi-static exponential model of surface tension. Wiechert et al. (2009) and Wiechert (2011) also use the dynamic model of Otis et al. (1994) to add the interfacial energy of the surfactant layer to the surface of the alveolar wall modeled by finite solid elements.

Liquid membrane theory
The liquid membrane formulation considered here is following the finite element formulation of Sauer et al. (2014) and Sauer (2014). This section briefly reviews the existing formulation and extends it to account for dynamic surface tension and constant line tension.

Membrane kinematics
In general, the membrane surface S can be characterized by a mapping from a parametric domain P with components ξ 1 and ξ 2 ∈ [−1, 1] as where t ∈ [0, T ] denotes time. From this follow the covariant tangent vectors The associated contra-variant tangent vectors a α = a αβ a β are obtained from the metric [a αβ ] = [a αβ ] −1 and a αβ = a α · a β . Accordingly, an infinitesimal surface area in the deformed configuration is related to the parametric domain as da = J a dξ 1 dξ 2 with J a = det a αβ and similarly dA = J A dξ 1 dξ 2 with J A = det A αβ in the reference configuration. Therefore, the surface stretch between the undeformed surface S 0 and the deformed surface S is Further the material time derivative of the surface stretch is (Sauer, 2018) where v :=ẋ is the material velocity. The surface divergence follows from div s v = tr ∇ s v, where ∇ s := (a α ⊗ a α ) · ∇ x is the surface gradient operator and ∇ x is the regular gradient operator in the current configuration. The surface normal is and the components of the curvature tensor b = b αβ a α ⊗ a β are given by where a α,β = ∂a α /∂ξ β . The mean and Gaussian curvature of surface are, respectively, and where

Constitution of liquid membranes
The stress state in liquid membranes has a hydrostatic component that is governed by the surface tension γ (Sauer, 2014) and a viscous component that is assumed here to follow a simple linear (i.e. Newtonian) viscosity model (Sahu et al., 2017). In this case, the surface stress tensor can be written as where are the contra-variant stress components and η is the kinematic surface viscosity. It is treated as constant here. γ on the other hand can change with time. This is different to the formulation of Sauer et al. (2014) that treats γ constant.
For many applications, such as pulmonary surfactants, surface tension shows hysteresis. Hence, to determine γ(t) at the current time step t, an evolution equation is needed to find the surface tension changes over time. This evolution law has the form of an ordinary differential equation (ODE) e.g.γ = f (γ, J, t). For such a formulation, one can use the dynamic CR and AR models introduced in Secs. 4.1 and 4.2, respectively. Later in Sec. 5.4, a time integration scheme of these models is introduced.

Strong form
According to Cauchy's formula, the membrane traction on a surface normal to a α is Furthermore, the balance of linear momentum provides us with the strong form of the equilibrium equation, which governs the membrane, as where the inertia effects are neglected since only problems at small length scale are considered here. Further, f is a distributed body force, which can be decomposed into in-plane and outof-plane components as f = f α a α + p n, where p is the net pressure exerted on the membrane surface. For the case that there is fluid only on one side, e.g. for a droplet, p equals to the difference of the fluid pressure p f and the external pressure p ext , acting on the membrane surface as The fluid pressure is composed of the hydrostatic pressure where ρ is the fluid density and where g = [0, 0, g] T is the gravity vector, and the capillary pressure, p v , across the liquid-gas interface S LG , that can be interpreted as the Lagrange multiplier for the volume constraint in case of incompressible droplets Sauer, 2014). Here, V 0 is the prescribed volume. For liquid droplets and menisci, where the internal fluid flow is not significant, it is computationally more efficient to model the internal behavior by the scalar Eq. (19), instead of discretizing the interior and solving a flow problem there. The solution of Eq. (15) is found by applying the Dirichlet and Neumann boundary conditions on the membrane boundary ∂S = ∂ u S ∪ ∂ t S. An example for the tractiont is the distributed contact line force described in Sec. 3.5.

Weak form
As shown by Sauer et al. (2014), the weak form of Eq. (15) is given by where are the internal and external virtual work contributions, respectively. Here, w ∈ W is a kinematically admissible variation of x ∈ S and W is a suitable space for w = w α a α + w n, where w α := w · a α and w := w · n.
Remark 3.1. For the extension of Eq. (21) to liquid films with bending resistance, e.g. for lipid bilayers, see .

Contact line
Any computational model for the simulation of droplets and liquid menisci should be able to handle the contact line. The numerical enforcement of equilibrium at the contact line in the context of nonlinear finite elements is explored in detail in the work of Sauer (2014), which also presents a general droplet contact model. The model is extended by Sauer (2016b) to allow for contact angle hysteresis, using a frictional sliding algorithm, arbitrary meniscus shapes and arbitrary substrate roughness, heterogeneity and compliance. As shown in Fig. 1, the droplet contact model consists of two bodies: The solid substrate B and the liquid droplet D. The droplet and solid substrate have separate interfaces with the surrounding gas, denoted by S SG and S LG , respectively and they share the S SL interface. The contact line C is the common intersection of the three mentioned interfaces. Along the contact line, the droplet surface meets the substrate at a contact angle θ c = θ c (x c ), which can depend on the position x c ∈ C if the substrate is heterogeneous or if contact angle hysteresis occurs. Here, we assume that the substrate is homogeneous and the contact angle is constant. For the behavior of droplets on chemically heterogeneous substrate surfaces see Luginsland and Sauer (2017).
In contrast to Sauer (2014Sauer ( , 2016a and Sauer and Luginsland (2018) who consider closed membrane models that explicitly account for the solid-liquid interface S SL , on which the contact constraints for non-penetration and sticking are enforced point-wise, here, S SL is not included in the model as the contact conditions at this interface are assumed to be homogeneous. Without S SL the model becomes an open membrane model. In this case, the contact impenetrability constraint is replaced by the following choice of boundary conditions at the contact line: 1) a pure Dirichlet boundary condition, where x c is fixed, 2) a pure Neumann boundary condition, wheret is given (see Fig. 2.a), and 3) a mixed Dirichlet-Neumann boundary condition. Here, the first approach is followed for pinned droplets and menisci and the third one is adopted for droplets and menisci that have fixed contact angles.
As shown in Fig. 2.b, the Neumann traction at the contact boundary is given bȳ with the componentst Further,t m = γ SG − γ SG . These component are defined in the basis {m c , a c , n c } illustrated in Fig. 1. The substrate normal n c is considered fixed here (e.g. n c = [0, 0, 1] T ); a c is the tangent to the contact line at x c and the outward unit normal of contact line C is given by In the following, we discuss the application of mixed Dirichlet-Neumann boundary conditions. By imposing a Dirichlet boundary condition in normal direction n c , the normal component of the contact force,t n , is found as the reaction force corresponding to this boundary condition.
Here, two approaches are proposed to define an appropriate value for the tangential component of the contact line force,t m , corresponding to a given contact angle θ c .

General shapes
According to Eq. (24), the contact line forcet m depends explicitly on the surface tension of the liquid-gas interface γ LG . Even though the surface tension is not constant, i.e. γ LG = γ(x c , t), Eq. (24) can still be used to enforce the contact angle at x c for droplets and menisci.

Droplets
Alternatively, Eq. (24) can be reformulated as such thatt m now depends on the unknown normal forcet n . In the contact-based models of Sauer (2014Sauer ( , 2016a,t n corresponds to the normal contact force at x c (and it is termed q n instead oft n ). Here, it corresponds to the reaction force at the Dirichlet boundary, which can be determined analytically for planar rigid substrates. As shown in Fig. 2.a, if the contact surface is horizontal and planar, the contact pressure p c is equal to the internal pressure of the fluid (17) As schematically shown in Fig. 3, considering an infinitesimal surface area da = 1 2 r c ds on the substrate,t n can be found by balancing the forces in the normal direction n c as Here, r c := x c − x 0 , where x 0 denotes the center of the circular contact surface. Plugging Eq. (29) into Eq. (26), the tangential force along the contact line is obtained as which is applied as a displacement-follower load distributed along the contact line to enforce the specified contact angle θ c . It should be noted that the latter expression fort m cannot be used for the particular case of liquid menisci that are modeled with two Dirichlet boundaries and accordinglyt n cannot be found analytically. To summarize, the droplet contact model introduced here is based on three assumptions: 1) As the solid substrate is rigid, horizontal and planar and accordingly the solid-liquid interface S SL remains planar, there is no need to explicitly consider this interface and the droplet can be modeled as an open membrane. 2) To model the contact between the droplet and the rigid substrate, a Dirichlet boundary condition is imposed on the contact line in the direction of the substrate normal n c . 3) To enforce a given contact angle θ c , a distributed forcet m =t m m c is applied along the contact line, wheret m is obtained according to Secs. 3.5.1 or 3.5.2.

Line tension
The contact angle introduced in Sec. 3.5 is called an intrinsic contact angle (Marmur, 1997) when the solid substrate has a smooth surface with homogeneous chemical and interfacial properties. As already mentioned, the Young equation is usually used to measure the contact angle. However, it has been experimentally observed for small systems that the contact angle predicted by Eq. (31) is different from the empirically measured contact angle, due to the effect of line tension in C. Therefore, Eq. (31) is modified as (Duncan et al., 1995 andMarmur, 1997) cos where λ (with the unit of force) is the line tension along the contact line C and κ g is the geodesic curvature of C at point x c ∈ C. For the special case of axially symmetric droplets and menisci, κ g = 1/r c . Here, r c = r c is the radius of the circular three-phase contact line in the plane of the substrate surface, described in Sec. 3.5.2. Further, θ ∞ c is the contact angle as r c → ∞, given by Eq. (31). In analogy to surface tension, line tension can be associated with the free line energy As shown in Appendix A, the spatial variation of the free energy W λ gives the contribution of the line tension to weak form (22.1) as 4 Material models for the surface tension In Sec. 2.2, the existing literature on the theoretical models for the surfactants-dependent dynamics of surface tension is briefly reviewed. Here, the adsorption-limited (AL) model of Otis et al. (1994) and the compression-relaxation (CR) model of Saad et al. (2010), which are in good agreement with experimental studies and can be efficiently implemented within a finite membrane element formulation, are discussed. The implementation is then discussed in Sec. 5. Otis et al. (1994) introduce a dynamic model for surfactant TA ®3 based on an experiment with a pulsating bubble surfactometer (PBS). For the adsorption and desorption processes, it is assumed that the behavior of surfactants concentration Γ at the liquid-gas interface can be described in three different regimes. These regimes are distinguished by two specific concentration values: The maximum equilibrium interfacial concentration Γ * and the maximum interfacial concentration Γ max , which are constants to be determined from experiments. Otis et al. (1994) assume that if the maximum surfactants concentration is reached and the surface area is reduced beyond the corresponding level of Γ max , surfactants are squeezed out of the interface so that the concentration of surfactants remains at Γ max . Furthermore, for the concentration values between Γ * and Γ max , it is assumed that there is neither adsorption nor desorption, which implies that the surfactant molecules remain at the interface. For the regime below Γ * , adsorption and desorption processes are governed by Langmuir kinetics (Miller et al., 1994). Therefore, three governing equations for the dynamics of the surface concentration can be postulated as (Otis et al., 1994)

Adsorption-limited (AL) model
where d(.)/dt is the material derivative, J is the local stretch and K 1 := k 1 C. The coefficients of adsorption onto and desorption from the surface are k 1 and k 2 , respectively. The bulk concentration C is assumed to be constant as diffusion effects, which occur at different rates, are neglected w.r.t. the adsorption/desorption process. Thus, k 1 and C are merged into one parameter K 1 . Here, for the sake of simplicity, the normalized concentration φ = Γ/Γ * and the normalized maximum concentration φ max = Γ max /Γ * are used.
Remark 4.1. Here, following continuum theory, the model of (Otis et al., 1994) is defined locally in terms of local stretch J rather than the total surface area of the membrane as it is done in Otis et al. (1994).
Remark 4.2. Morris et al. (2001) extend the formulation of Otis et al. (1994) by including diffusion effects. In this sense, C varies with time and position within the bulk. However, for thin liquid membranes, such sophisticated formulations are not efficient. Furthermore, Krueger and Gaver (2000) modify the model of Otis et al. (1994) by including a secondary collapse layer, which leads to a more complicated set of state equations.
From Eq. (35), one can derive the surface tension in form of an equation of state, which relates the surfactants interfacial concentration Γ to the surface tension γ. For this purpose, Otis et al. (1994) suggest two straight lines that meet at Γ = Γ * , which give where γ 0 is the temperature-dependent surface tension of water, e.g. ≈ 70 mN/m at 25 • C, and γ * is the minimum equilibrium surface tension, corresponding to the maximum equilibrium surfactants concentration Γ * . The parameters K 1 , k 2 and m 2 are the only unknowns that need to be identified by fitting model (35)-(36) to experimental results; on the other hand, m 1 can be found by knowing γ * and using Eq. (36.1). If there are no surfactants on the interface, i.e. φ = 0, the surface tension should be that of pure water, i.e. γ = γ 0 . Therefore, If normalized concentration values are used, Γ * is not needed; however, it can be assumed as 3 mg/m 2 (Morris et al., 2001). Two other parameters, namely γ * and γ min , need to be determined directly from experiments. The minimum equilibrium surface tension γ * can be determined by finding the lowest surface tension obtained for repeated measurements that are in equilibrium and have considerably high surfactants concentration in the bulk. For example, it can be 22.2 mN/m (Ingenito et al., 1999) or 25 mN/m (Morris et al., 2001). Similarly, the minimum surface tension γ min , can be obtained from cyclic experiments, which usually give 1 ∼ 2 mN/m (Morris et al., 2001). Finally, from Eq. (36.2), Thus, the AL model needs six independent parameters, i.e. K 1 , k 2 , m 2 , γ 0 , γ * , and γ min , among which the first three ones are found by curve fitting and the last three ones are directly determined from an experiment. Saad et al. (2010)

Compression-relaxation (CR) model
If the maximum equilibrium interfacial concentration Γ * is known, the instantaneous equilibrium interfacial concentration Γ eq , depends on the instantaneous bulk concentration C as (Otis et al., 1994) Γ eq Γ * = If the bulk concentration C is fixed, by plugging Eq. (40) into Eq. (39), one can derive where k := K 1 + k 2 is a coefficient for the dynamic adsorption and desorption. Assuming a linear relationship between surface tension and interfacial concentration (see Eq. (36)), the rate of changes in surface tension can be obtained as where the equilibrium surface tension γ eq corresponds to the specified equilibrium interfacial concentration Γ eq . For the adsorption or spreading process, where γ ≥ γ eq , the adsorption coefficient k = k a defines the adsorption rate. Similarly, for the desorption or relaxation, where γ < γ eq , the desorption rate is given by the desorption coefficient k = k r . Both k a and k r are parameters to be found by fitting the model to experimental results.
The changes of surface tension γ w.r.t. the surface stretch J can be characterized by the elasticity parameter = dγ d ln J that is assumed constant here. It corresponds to the elastic surface bulk modulus at J ≈ 1 . Rearranging and dividing by dt leads to where = c during compression and = e during expansion. The elasticity coefficients c > 0 and e > 0 should also be obtained through a parameter identification process. Furthermore, according to experiments (e.g. Otis et al., 1994), there is a maximum surfactants concentration, which is the limit where the layers of surfactants collapse. Therefore, the surface tension cannot be reduced below the minimum surface tension γ min by further compression.
Based on the four main processes that affect the dynamic response of surface tension -namely 1) adsorption or spreading, 2) desorption or relaxation, 3) elasticity during compression and 4) elasticity during expansion - Saad et al. (2010) postulate a unified formulation that combines all cases. In the framework of continuum mechanics, it can be formulated as Here, it is assumed that the adsorption/desorption and elasticity processes can happen at the same time. The CR model requires six different parameters: k a , k r , c and e are set by curve fitting of experimental results. γ min and γ eq are directly measured by experiments.
From a physical point of view, local variation of concentration results in a surface tension gradient, which in turn induces an interfacial flow known as Marangoni flow (Scriven andSternling, 1960 andVelarde et al., 2002). If the bulk viscosity is significant, the Marangoni flow will cause an internal flow in the bulk of the fluid through viscous forces. However, in the examples presented in this work, the bulk viscosity is neglected. Generally, Marangoni effect occurs due to temperature or concentration gradients and also at boundaries of different liquids. As the AL and CR models, introduced in Sec. 4, are explicitly or implicitly dependent on the surfactants concentration, Marangoni flow can be expected if the point-wise approach is followed. The Lagrangian description of the membrane presented here can capture such flows provided that they are not large, as shown in the example of Sec. 6.3.
If the surface flow becomes large, an arbitrary Lagrangian-Eulerian (ALE) formulation can be used (Yang andJames, 2007 andGanesan andTobiska, 2009), which is not considered here. The ALE approach requires expanding the material time derivatives into a spacial and a convective part (Stone, 1990) where the surface gradient operator, ∇ s , is given in Sec. 3.1. The rear term in Eq. (48) represents surface diffusivity assuming Fick's laws of diffusion, where D s is the diffusivity. f is a surfactant source term representing the surface adsorption/desorption process, which can be derived following the approach of Morris et al. (2001). The interfacial flows, governed by Eq. (48), are not considered in this study.

Finite element solution
The finite element method is used to solve the governing equations. Thus, the membrane surface S is approximated by the discretized surface S h , which is constructed from elements Ω e . In this section, first the FE interpolation is introduced and then the discretized weak form is summarized following Sauer et al. (2014) and Sauer (2014). Here also the contribution of line tension to the discretized weak form is derived. Third and foremost, the evolution equations of the surface tension for the AL and CR models are discretized within an implicit time integration scheme. As the BVP (21) is highly nonlinear, it should be linearized and solved by an iterative algorithm such as the Newton-Raphson method..

Finite element interpolation
Any point on the surface x ∈ S is approximated by x h ∈ S h as where x e is a vector that contains all control point positions belonging to Ω e and the shape functions are arranged as N := [N 1 1, N 2 1, ..., N ne 1], with n e being the number of control points or nodes of element Ω e . Similarly in the reference configuration, Due to Eq. (49), we have and Following a Galerkin approach, the variation w is approximated like the deformation, which gives w ≈ N w e .
Here, for the finite element interpolation, quadratic NURBS-based shape functions with n e = 9 are used. Here, {N e A } ne A=1 is the B-spline basis function expressed in terms of Bernstein polynomials using the Bézier extraction operator (Borden et al., 2011). However, as shown by Sauer et al. (2014), classical finite elements with quadratic Lagrange interpolation polynomials are also suitable for the modeling of liquid drops and menisci.

FE force vectors
Following the FE setting introduced in Sec. 5.1, the discretized version of weak form (21) is with the internal FE force vector f e int := f e intτ + f e intλ , where and the external FE force vector where and if the contact line model of Sec. 3.5.1 is used and if the contact line model of Sec. 3.5.2 is used. Here, it is assumed that the external distributed load per surface is f = f 0 /J + p n, where f 0 is a dead load per reference area and p is a live pressure as described in Sec. 3.3; however, in the examples of Sec. 6, dead loading is not considered. Besides, f e extt is the contribution of a prescribed distributed loadt on the Neumann boundaries Γ e t ⊂ ∂ t S other than the contact boundary; f e extc is the contribution of the tangential force q m along the contact line element Γ e c ⊂ C and N t is the array of shape functions corresponding to the contact line.

The Newton-Raphson iteration
As the discretized weak form (56) holds for any admissible w e , the governing equation can be reduced to the system of ODEs where f is the global residual force vector formed by an assembly of FE force vectors f e int and f e ext introduced in Sec. 5.2. x is a vector that collects all nodal positions (in a Cartesian coordinate system),ẋ := dx/dt, and p v is a single scalar Lagrange multiplier, see Sec. 3.3.ẋ is eliminated by a time discretization scheme (see Sec. 5.4) such that f = 0 becomes an algebraic equation for x and p v . Following the introduced FE setting, the volume V , enclosed by the discretized surface, is ) This equation is derived for closed droplets, but it also applies to open droplets, if the origin of the coordinate system lies on the contact surface. Thus, the complete system of equations to be solved by the finite element method is As F = F(x, p v ) is nonlinear, it can be solved by the Newton-Raphson method, which is an iterative algorithm. This requires the linearization of F w.r.t. x and p v as where K is the global tangent matrix and includes all degrees of freedom. The global tangent matrix can be arranged as where k := k int − k ext and The element tangent matrices l e ext and h e v can be found in Sauer et al. (2014). The internal stiffness matrices k e intτ and k e intλ and the tangent matrix k e extc for the contact line force vector are derived in Appendix B in detail.
Remark 5.1. For cases where the volume is not constrained, e.g. for free films, Eq. (66) reduces to ∆f = k ∆x.
Remark 5.2. It should be noted that, for the constitutive laws presented in Sec. 4, the internal stiffness tangent k int is influenced by the nonlinearity of the surface tension γ; therefore, it has extra terms compared to Sauer et al. (2014) Remark 5.3. All the tangents k e int , k e ext and l e ext are derived at the current time step, t n . Hence, the history quantities coming from the previous time steps (t n−1 , t n−2 , ...) are treated as constants in the Newton-Raphson iteration and thus do not need to be linearized.

Time integration
In the introduced dynamic setup, the surface tension and stress are time dependent. Thus, they should be updated with a discretized time integration scheme. Here, the backward Euler scheme is used, i.e.γ where ∆t n is the time step between t n and t n−1 . Equation (71) is implicit and requires solving an algebraic equation in terms of unknown γ n , which is solved iteratively by the Newton-Raphson method introduced in Sec. 5.3. The contribution of the dynamic concentration-dependent surface tension γ n , which varies by J n , to the internal tangent matrix is further explained in Secs. 5.4.1 and 5.4.2. The solution algorithm for the introduced dynamic problem is summarized in Table 1. 1) Retrieve relevant quantities from t n−1 , i.e. γ n−1 and J n−1 .
2) Update the current time step t n .

3) Newton-Raphson iteration
which gives the time-discretized version of the CR model as Then, the surface tension is updated according to the equation of state Eq. (36), which gives The linearization of γ = γ(x n , t n ) leads to where and The contribution of ∆ x γ to the material tangent matrix is shown in Appendix B.1.

CR model
Plugging Eq. (45) into Eq. (71), the time-discretized version of the CR model is obtained as The linearization of γ is also given by Eq. (76). From Eq. (80), we now have

Numerical examples
In Sec. 1, the typical approaches for dynamic surface tension measurement -namely the Langmuir-Wilhelmy balance (Clements, 1957), captive bubble surfactometer (CBS)  and constrained sessile drop (CSD) (Saad et al., 2010) -are introduced. Here, the capabilities of the presented membrane formulation to simulate similar cases are demonstrated. The main focus of the simulations is on the drop shape methods. First, the CSD test is examined following the experimental setup and the CR material model of Saad et al. (2010). Then, the behavior of a thin liquid film in cyclic compression/expansion and relaxation without external loads is analyzed. Next, the pendant drop (PD), sessile drop (SD) and liquid bridge (LB) tests are simulated. Finally, the influence of the line tension on the contact angle is studied. Fig. 4 shows the computational models with the corresponding boundary conditions of all the examples. As discussed by Alonso et al. (2004) and Hermans et al. (2015), in addition to the adsorption/desorption mechanisms, surface viscosity has also important contribution to the interfacial behavior of pulmonary surfactants. The corresponding values of the surface viscosity are not reported for the experiments from which the material parameters are taken here. Thus, the surface viscosity is taken as η = 1 mNs/m, which is within the ranges reported by Alonso et al. (2004) and Rudiger et al. (2005)

Constrained sessile drop (CSD) test
As the first example, the constrained sessile drop (CSD) test, which is one of the most popular methods for dynamic measurement of surface tension, is simulated. The method is also used by Saad et al. (2010) for axisymmetric drop shape analysis. As the droplet is constrained, it is modeled as a pinned droplet, where the contact angle is not fixed. Therefore, the contact line formulation presented in Sec. 3.5, is not needed for this problem. The droplet is modeled as a hemisphere with initial radius of 1.5 mm, whose opening radius is fixed on a rigid pedestal and its volume is decreased and increased cyclically. As shown in Fig. 4.a, only 1/4 of the droplet is modeled to exploit the symmetry of the problem in order to reduce the computational cost.
Following (Saad et al., 2010), for each test, two different concentrations of BLES 4 in humid and dry conditions are considered. As already mentioned, the CR model requires the elasticity of compression c , the elasticity of expansion e , the relaxation coefficient k r , the adsorption coefficient k a , the equilibrium surface tension γ eq and the minimum surface tension γ min . The corresponding parameters are adopted from Saad et al. (2010) as listed in Tab. 2. The droplet volume is controlled by the volume constraint (19). Fig. 5 show the droplet, modeled with the parameters of case 1 in Tab. 2, at different volumes colored by the surface tension γ [mN/m].    (Saad et al., 2010) as the volume of droplet changes, small deviations in the surface tension are expected. Thus, to compare the results, the average surface tension, is examined. Here, S is the total droplet surface area, in the deformed configuration, defined as Similarly, the initial droplet surface is defined as in the reference configuration. As shown in the example of Sec. 6.3, if enough time is given to the system to relax, it relaxes to a new configuration, where the surface tension is more homogeneous. Fig. 7 compares the standard deviation of the surface tension, σ(γ), and its mean valueγ. As the standard deviation of the surface tension is small w.r.t. the mean value for all the examples except the liquid film of Sec. 6.3, it is only reported here to keep the paper short. In Fig. 8, the corresponding isotherm is plotted for different materials. The simulation results are compared with the experimental results of Saad et al. (2010). The qualitative agreement between our FE simulations and the experiments of Saad et al. (2010) is very good.
However there are small differences. These differences can have different reasons: For example, all the required parameters to simulate the test exactly as the experiment are not provided in the literature. We have estimated some variables such as γ min and γ eq . Further, the loading pattern used here to control the droplet volume can be slightly different from the loading patter used in the experiment of Saad et al. (2010).

Liquid film expansion/compression
In the Langmuir-Wilhelmy balance (Butt et al., 2006) a movable barrier controls the total area available for a thin film of interfacial molecules, and the surface tension is measured from the wetting force acting over the Wilhelmy plate. Similarly, here a thin liquid film is modeled by a 2 L 0 × L 0 rectangular membrane, where L 0 is a length scale. The bulk effects are not considered and only the free surface is modeled. As shown in Fig. 4.b, three edges of the membrane are fixed as zero-displacement Dirichlet boundaries and the forth edge is controlled by a prescribed Dirichlet boundary conditionū = [ū x , 0, 0] T that changes with time as where T = 3 s and the loading amplitude is ∆L = 0.33 L 0 unless it is changed, e.g. in Fig. 14.
In this example, the behavior of the surfactants monolayer is described by the AL model of Sec. 4.1. The material constants are set as listed in Tab. 3 following Wiechert (2011).
0.016 48 140 10 70 Table 3: Parameters of the AL model (Wiechert, 2011) In Figs. 9 to 14, the changes of surface tension γ and the normalized interfacial surfactants concentration of surfactants φ are plotted against the surface area change J with different values for the model parameters. As is expected, our finite element results are in agreement with the finite element results of Wiechert (2011).

Relaxation of a liquid film
As already discussed in Remark 4.3, the point-wise approach to model the surface tension and interfacial concentration can result in stress gradients that in turn cause interfacial flows. Even though, the presented finite element formulation is Lagrangian and not Eulerian, it can represent such flows as long as they are not too large, which is demonstrated next. We therefore consider a L × L square liquid film, the edges of which are completely fixed. It is assumed that the film consists of BLES 0.5 mg/ml in humid conditions (see Tab. 2) and it is modeled by the CR model. Further, it is supposed that the surfactants concentration is not homogeneous, which causes the initial surface tension distribution shown in Fig. 15.a. As it can be seen in Fig. 15.b, if the system is allowed to relax, the surface deforms due to the local variation in the surface tension and finally reaches a new configuration, where the surface tension is not only homogeneous but also equal to the equilibrium surface tension γ eq = 25 mN/m. Fig. 15.c shows the corresponding deformation map colored by the surface stretch J in the final relaxed state and Fig. 15.d shows how the mean surface tensionγ, defined by Eq. (84), and the corresponding standard deviation σ(γ) evolve during the relaxation process. Hence, by employing the presented computational model, where the local stretch and the local surfactants concentration (and equivalently the local surface tension) are mutually related, the system will indeed find a new equilibrated configuration with a homogeneous surface tension.

Pendant drop (PD) test
In this section, the pendant drop test, introduced earlier, is simulated. In the reference configuration, the drop is modeled as 1/8 of a sphere with boundary conditions similar to the constrained droplet (see Fig. 4.a) with the difference that the drop is held downward in the direction of the gravity force. The initial radius of the drop is 1.0 mm. Here, the volume of the droplet is controlled following the volume constraint (19 show the finite element results for the cases where the volume is controlled in a monotone way. As a second example, the drop volume is altered cyclically as where T = 3 s and ∆V max = 0.3 V 0 here. Fig. 18 shows the dynamic changes of the mean surface tension γ, given by Eq. (84), and the surface area with cyclic alteration of the drop volume.
As it can be observed, for small changes in the volume, both models predict similar surface tension; however, as the drop volume increases, the models predict different surface tension and accordingly difference Bond numbers for an identical given volume. It should be noted that a detailed comparison of the two models is only possible when the parameters of each model is identified for the same experiment with the same type of surfactants. This is left for future work.
Here, the parameters are adopted from two different resources. Nevertheless, the difference in the behavior of these two models is well discussed in Saad et al. (2010) and Saad (2011).
Bo ≈ 4.87 Figure 16: The PD test: Drop shapes modeled by the AL model and colored by γ.
Bo ≈ 4.87  b) Figure 18: The PD test: Cyclic alteration of volume V , surface area S and the mean surface tensionγ for a) the AL and b) the CR models.

Unconstrained sessile drop (SD) test
In Sec. 6.1, a constrained sessile drop, that is pinned at the contact line, is studied. Here, an unconstrained sessile drop is simulated, where the contact line is free to move on the substrate and the contact angle is fixed. Particularly, it is assumed that the contact angle remains the same if the contact line is advancing or receding. One can use the proposed formulation of Sauer (2016b) to distinguish between the advancing and receding contact lines, which is not done here. The initial radius of the droplet is R = 1.0 mm and the computational model with the corresponding boundary conditions are shown in Fig. 4.c. Here, the same material parameters as in the example of Sec. 6.4 are used for both the AL and CR models. To enforce the contact line constraint, the formulation of Sec

Liquid bridge (LB) test
This section is devoted to show the capability of the presented computational model to simulate liquid bridges or menisci. As illustrated in Fig. 4 bridge is modeled as 1/4 of a cylinder to benefit from the symmetry of the problem. Accordingly, the symmetry boundary conditions are applied on the axial edges. On the top edge, the liquid bridge with length L = 2 mm is pinned to the end of a needle or a small flat circular holder with radius R = 2 mm. On the bottom edge, the contact angle constraint is enforced following the approach of Sec. 3.5.1. Here, only the results of the CR model are presented and the material parameters and the loading conditions are the same as in the example of Sec. 6.4. As it can be seen in Fig. 23, the cyclic deformation of a liquid bridge can be efficiently predicted by the presented formulation.

Sessile drop with line tension
In the examples shown so far, the influence of line tension was neglected. Here, it is investigated how the line tension can affect the apparent contact angle. For this purpose, the unconstrained sessile drop of Sec. 6.5 is re-examined. To study only the influence of the line tension, here the surface tension is assumed to be constant, namely γ = 22 mN/m and the line tension is set as λ = γ R/2. The contact angle θ ∞ c = 60 • is imposed following the droplet contact line model of Sec. 3.5.2. The volume of the droplet is increased from V = V 0 to V = 20 V 0 and gravity is taken into account similar to the previous drop examples. As shown in Fig. 24.a, here the contact angle observed during the FE simulation is compared with the contact angle predicted by Eq. (32). As it is expected, by increasing the radius of the contact line, the contribution of the line tension decreases and θ c approaches θ ∞ c . Here, the contact angle is measured at the same point along the contact line, denoted by a filled circle in Fig. 24.b, where θ ∞ c is also shown by thick solid lines. The small difference between the observed contact angle and the value predicted by Eq. (32) decreases as the mesh is refined. Therefore, depending on the values of the surface tension, line tension and contact line radius, the contribution of line tension can be efficiently incorporated into the model.

Conclusion
A new computational framework is presented to model dynamic concentration-dependent surface tension of liquids. In particular, the new formulation is adapted to two material models, namely the adsorption-limited (AL) model of Otis et al. (1994) and the compression-relaxation (CR) model of Saad et al. (2010), which have been developed to study dynamic surface tension of pulmonary surfactants. The models are adapted here to a general continuum membrane formulation in the framework of arbitrary large deformations. The presented model can also be used for other similar constitutive laws, where the surface tension depends on the local concentration of surfactant and the local interface stretch. As discussed, the dynamics of pulmonary surface tension plays an important role in the proper function of lungs and any deficiency or disorder in surfactant behavior can result in serious pulmonary diseases such as respiratory distress syndrome (RDS). Thus, the development of computational models, such as the presented framework, is of great importance for a better understanding of lung biomechanics and for developing new methods that can prevent pulmonary disorders or provide better treatments for them.
As it is shown through multiple numerical examples, the presented formulation can be used to simulate liquid films, drops (such as constrained and unconstrained sessile drops and pendant drops) and liquid menisci (or liquid bridges). An efficient formulation is presented to model the contact line for liquids of arbitrary shapes that are in contact with rigid planar smooth solid substrates. In this regard, a general formulation and a specific one for droplets is presented. Further, the contribution of the line tension to the contact line is investigated. The presented approach for modeling dynamic surface tension can be extended to liquid films with bending resistance  and can be coupled to fluid flow within the droplet (Sauer and Luginsland, 2018). The introduced liquid membrane formulation can be used through a parameter identification process to obtain the material parameters of the presented surfactant models. The application of such an inverse analysis is left for future work.

Acknowledgment
Financial support from the German Research Foundation (DFG) through grant GSC 111, is gratefully acknowledged. The authors also wish to thank Katharina Immel for her help in surveying the literature.

A Contribution of line tension to the weak form
Along the contact line C, ds = dx c = a c dξ , where ξ ∈ [−1, 1] is the parameterization of curve C according to x c = x c (ξ) such that and s denotes the arc length. It should be noted that this parameterization of curve C is associated to the parametrization of surface S introduced in Sec. 3.1. Strictly, ξ := ξ 1 with ξ 2 ± 1 or ξ := ξ 2 with ξ 1 ± 1. As shown by Sauer (2014), and from Eq. (89) δds = δ a c dξ .

B FE tangent matrices
In this appendix, the internal stiffness tangents that allow for the nonlinearities of surface tension γ = γ(x, t) and the tangent matrix of the contact line force vector (61) and (62)

B.1 Internal stiffness tangent
The spacial variation of the internal force vector f e int , given in Eq. (57), is where τ αβ := J σ αβ . Since for all the models considered here, γ depends on x through J only, the liquid stress tensor τ αβ yields ∆ x τ αβ = γ J ∆ x a αβ + a αβ γ + J ∂γ ∂J where the two front terms are the variations due the surface tension and the rear term is the variation due to the time-discretized viscous stress. Following Sauer and Duong (2017) ∆ x a αβ = a αβγδ ∆ x a γδ , with a αβγδ := − 1 2 a αγ a βδ + a αδ a βγ , and ∆ x a αβ = a α · N ,β + a β · N ,α ∆x e .
Thus, Eq. (95) can be reformulated as with c αβγδ := 2 ∂τ αβ ∂a γδ = 2 γ J a αβγδ + γ + J ∂γ ∂J J a αβ a γδ − η J ȧ αβ a γδ + 2 ∆t a αβγδ , where ∂γ/∂J is derived in Secs. 5.4.1 and 5.4.2 for the AL and CR models, respectively. To determine the element stiffness tangents, Eq. (94) can be arranged as  ∆f e int = k e mat + k e geo ∆x e , where the material stiffness matrix The linearization of Eq. (104) results in ∆f e extc = k e extc ∆x e + l e extc ∆p .
For a rigid planar substrate, where ∆n c = 0, (Sauer, 2014) ∆m c = − a c ⊗ m c ∆a c with a c to be the dual vector of a c , where a c ·a c = 0 as described in Sec. 3.1, and ∆a c = N t,ξ ∆x e . Furthermore, ∆ r c = r c −1 r c · ∆x c , where ∆x c = N t ∆x e . Thus, where the rear term explicitly depends on the nonlinearities of the surface tension γ = γ(J) = γ(x c ) in general and is vanished if the surface tension is constant as for the typical pure liquids.
If the effect of line tension is considered, following Eq. (58), the tangent matrix k e intλ := should be added to the internal tangent matrix.