In this document we demonstrate the adaptive solution of the Young Laplace equation with contact angle boundary conditions. We start by reviewing the physical background in the context of a representative model problem, and then discuss the spine-based representation of free contact lines and the implementation of the contact angle boundary condition along such lines.
Acknowledgement: This tutorial and the associated driver codes were developed jointly with Cedric Ody (Ecole Polytechnique, Paris; now Rennes).
|
A model problem
The figure below shows a sketch of a T-junction in a microchannel with a rectangular cross-section. (The front wall has been removed for clarity). Fluid is being pushed quasi-steadily along the (vertical) main channel and is in the process of entering the T-junction. We assume that the air-liquid interface (shown in red) remains pinned at the two sharp edges (at ) where the channels meet, while the meniscus forms a quasi-static contact angle, , with the smooth front and back walls.
A typical problem: Fluid propagates quasi-steadily through a T-junction that connects two channels of rectangular cross-section.
It is of interest to determine the maximum pressure that the meniscus can withstand: if the driving pressure is less than that value, the fluid will not be able to propagate past the T-junction.
Theory and implementation
Spine-based representation of the meniscus
Recall that we parametrised the meniscus by two intrinsic coordinates as , where . Furthermore, we parametrised the domain boundary, , by a scalar coordinate so that,
The normal to the meniscus is then given by
where a comma denotes partial differentiation with respect to one of the intrinsic coordinates,
Along the contact line we define two unit vectors, and , that are tangential to the meniscus. is tangent to the contact line while is normal to it and points away from the meniscus, as shown in the sketch below.
We split the domain boundary so that and assume that along the meniscus is pinned,
where is given. On the meniscus meets the wall at a prescribed contact angle so that
where is the outer unit normal to the wall as shown in this sketch:
Sketch of the meniscus, the contact line along which it meets the wall, and the spine-based representation of the meniscus.
The figure also illustrates the spine-based representation of the meniscus in the form
where the spine basis and spines are pre-determined vector fields, chosen such that
Computation of the contact-angle term in the variational principle
Recall that the variational principle that determines the shape of the meniscus contained the line term
Along the line integral vanishes because . The line integral can therefore be written as
or, using the spine-based representation of the meniscus, (2),
We shall now demonstrate that the integrand in this expression can be expressed in terms of the contact angle boundary condition (1). We start with several observations:
- is tangential to the wall.
- Since is normal to the wall, is tangential to the wall and orthogonal to .
- is tangential to the wall and can therefore be decomposed into its components parallel to and as
for some values of and . In fact,
- During the computation it is most convenient to perform all calculations in terms of quantities that are easily obtained from the parametrisation of the meniscus as this avoids having to specify explicitly. For this purpose we exploit that and are tangential to the wall and not parallel to each other (unless the parametrisation of the meniscus by (2) is no longer one-to-one). Therefore can be obtained from quantities that are intrinsic to the meniscus representation via
and thus
- Given (3) and the fact that , we have
and with (1):
Hence, the line integral may be written as
where is given by (4).
Equation (5) is easily discretised by finite elements. Within oomph-lib
, the line integral is decomposed into FaceElements
that are attached to the "bulk" Young-Laplace elements that are adjacent to the contact line. The imposition of the contact angle boundary condition for the Young Laplace equation is therefore as easy as the application of Neumann boundary conditions for a Poisson equation, say.
Results
The animation below illustrates the variation in the quasi-steady meniscus shape as the fluid enters the T-junction.
Animation of the meniscus shapes.
The computation was performed with full spatial adaptivity. The plot below illustrates how the automatic mesh adaptation has strongly refined the mesh towards the corners of the domain where the meniscus shape has a singularity. (The singularity develops because in the corners of the domain the contact angle boundary condition along the side walls is inconsistent with the contact angle enforced by the pinned boundary condition along the sharp edges.)
Illustration of the adaptive mesh refinement.
Finally, here is a plot of the "load-displacement diagram", i.e. a plot of the meniscus deflection as a function of its curvature (i.e. the applied pressure drop). The limit point indicates the maximum pressure that can be withstood by the static meniscus.
The load-displacement diagram for the meniscus.
The driver code
The modifications to the driver code required to impose the contact angle boundary conditions are very similar to those used in other driver codes for problems with Neumann-type boundary conditions. We attach FaceElements
to the appropriate faces of the "bulk" Young-Laplace elements detach/re-attach them before and after any spatial adaptation of the "bulk" mesh.
The global namespace
The namespace that defines the problem parameters is very similar to that used in the previous example without contact angle boundary conditions. We provide storage for the cosine of the contact angle, and the prescribed meniscus height that is used by the displacement control method.
{
double Cos_gamma=cos(MathematicalConstants::Pi/6.0);
Namespace for "global" problem parameters.
double Controlled_height
Height control value.
double Cos_gamma
Cos of contact angle.
As before, we use the spine basis to establish a reference configuration in which the flat meniscus is located in the plane and occupies the domain
Vector<double>& spine_B,
Vector< Vector<double> >& dspine_B)
{
spine_B[0] = x[0];
spine_B[1] = x[1];
spine_B[2] = 0.0 ;
dspine_B[0][0] = 1.0 ;
dspine_B[1][0] = 0.0 ;
dspine_B[0][1] = 0.0 ;
dspine_B[1][1] = 1.0 ;
dspine_B[0][2] = 0.0 ;
dspine_B[1][2] = 0.0 ;
}
double L_x
Length and width of the domain.
void spine_base_function(const Vector< double > &x, Vector< double > &spine_B, Vector< Vector< double > > &dspine_B)
Spine basis: The position vector to the basis of the spine as a function of the two coordinates x_1 a...
double L_y
Width of domain.
As in the previous example, we rotate the spines in the -direction to allow the representation of meniscus shapes that cannot be projected onto the -plane.
double Alpha_min = MathematicalConstants::Pi/2.0*1.5;
double Alpha_max = MathematicalConstants::Pi/2.0*0.5;
Vector<double>& spine,
Vector< Vector<double> >& dspine)
{
spine[0]=0.0;
dspine[0][0]=0.0;
dspine[1][0]=0.0;
dspine[0][1]=0.0;
dspine[0][2]=0.0;
}
}
double Alpha_max
Max. spine angle against horizontal plane.
void spine_function(const Vector< double > &x, Vector< double > &spine, Vector< Vector< double > > &dspine)
Spine: The spine vector field as a function of the two coordinates x_1 and x_2, and its derivatives w...
double Alpha_min
Min. spine angle against horizontal plane.
The driver code
We start by defining the output directory and open a trace file to record the load-displacement curve.
{
DocInfo doc_info;
ofstream trace_file;
doc_info.set_directory("RESLT");
char filename[100];
sprintf(filename,"%s/trace.dat",doc_info.directory().c_str());
trace_file.open(filename);
trace_file << "VARIABLES=\"<GREEK>k</GREEK>\",\"h\"" << std::endl;
trace_file << "ZONE" << std::endl;
Next, we create the problem object, refine the mesh uniformly and output the initial guess for the solution: a flat interface which, unlike the previous case, is not a solution of the problem because it does not satisfy the contact-angle boundary condition; see the section Comments and Exercises for a more detailed discussion of this issue.
problem.refine_uniformly();
doc_info.number()++;
2D RefineableYoungLaplace problem on rectangular domain, discretised with 2D QRefineableYoungLaplace ...
void doc_solution(DocInfo &doc_info, ofstream &trace_file)
Doc the solution. DocInfo object stores flags/labels for where the output gets written to and the tra...
Finally, we perform a parameter study by slowly incrementing the control displacement and recomputing the meniscus shape.
double increment=0.1;
unsigned nstep=2;
for (unsigned istep=0;istep<nstep;istep++)
{
unsigned max_adapt=1;
problem.newton_solve(max_adapt);
doc_info.number()++;
}
trace_file.close();
}
The problem class
The problem class contains the usual member functions. The functions actions_before_adapt()
and actions_after_adapt()
are used to detach and re-attach (and rebuild) the contact angle elements on the appropriate boundaries of the "bulk" mesh.
template<class ELEMENT>
{
public:
{
rebuild_global_mesh();
}
{
for (unsigned e=0;e<nel;e++)
{
YoungLaplaceContactAngleElement<ELEMENT> *el_pt =
dynamic_cast<YoungLaplaceContactAngleElement<ELEMENT>*>(
}
rebuild_global_mesh();
}
void create_contact_angle_elements(const unsigned &b)
Create YoungLaplace contact angle elements on the b-th boundary of the bulk mesh and add them to cont...
void actions_after_newton_solve()
Update the problem after solve: Empty.
~RefineableYoungLaplaceProblem()
Destructor (empty)
void actions_before_newton_solve()
Update the problem specs before solve: Empty.
Mesh * Contact_angle_mesh_pt
Pointer to the contact angle mesh.
RefineableYoungLaplaceProblem()
Constructor:
void actions_after_adapt()
Actions after adapt: Rebuild the mesh of contact angle elements.
void delete_contact_angle_elements()
Delete contact angle elements.
void actions_before_adapt()
Actions before adapt: Wipe the mesh of contact angle elements.
Two private helper functions are provided to create and delete the contact angle elements. The class also provides storage for the pointers to the various meshes, to the node at which the meniscus displacement is prescribed by the displacement control method, and to the Data
object whose one-and-only value stores the (unknown) meniscus curvature.
private:
void create_contact_angle_elements(const unsigned& b);
void delete_contact_angle_elements();
RefineableRectangularQuadMesh<ELEMENT>* Bulk_mesh_pt;
Mesh* Contact_angle_mesh_pt;
Mesh* Height_control_mesh_pt;
Node* Control_node_pt;
};
Data * Kappa_pt
Pointer to Data object that stores the prescribed curvature.
The problem constructor
We start by creating the "bulk" mesh of refineable Young Laplace elements and specify the error estimator.
template<class ELEMENT>
{
unsigned n_x=8;
unsigned n_y=8;
Bulk_mesh_pt=new RefineableRectangularQuadMesh<ELEMENT>(n_x,n_y,l_x,l_y);
Bulk_mesh_pt->spatial_error_estimator_pt()=new Z2ErrorEstimator;
Bulk_mesh_pt->max_permitted_error()=1.0e-4;
Bulk_mesh_pt->min_permitted_error()=1.0e-6;
We identify the node (in the centre of the mesh) at which we apply displacement control. We pass a pointer to this node to the constructor of the displacement control element and store that element in its own mesh.
if ((n_x%2!=0)||(n_y%2!=0))
{
cout << "n_x n_y should be even" << endl;
abort();
}
ELEMENT* prescribed_height_element_pt= dynamic_cast<ELEMENT*>(
Bulk_mesh_pt->element_pt(n_y*n_x/2+n_x/2));
Control_node_pt= static_cast<Node*>(prescribed_height_element_pt->node_pt(0));
std::cout << "Controlling height at (x,y) : (" << Control_node_pt->x(0)
<< "," << Control_node_pt->x(1) << ")" << "\n" << endl;
HeightControlElement* height_control_element_pt=new HeightControlElement(
Height_control_mesh_pt = new Mesh;
Height_control_mesh_pt->add_element_pt(height_control_element_pt);
Kappa_pt=height_control_element_pt->kappa_pt();
Next we create the mesh that stores the contact-angle elements. We attach these elements to boundaries 1 and 3 of the "bulk" mesh.
Contact_angle_mesh_pt=new Mesh;
create_contact_angle_elements(1);
create_contact_angle_elements(3);
The various sub-meshes are now added to the problem and the global mesh is built.
add_sub_mesh(Bulk_mesh_pt);
add_sub_mesh(Height_control_mesh_pt);
add_sub_mesh(Contact_angle_mesh_pt);
build_global_mesh();
As usual, we enforce only the essential boundary conditions directly by pinning the meniscus displacement along mesh boundaries 0 and 2:
unsigned n_bound = Bulk_mesh_pt->nboundary();
for(unsigned b=0;b<n_bound;b++)
{
if ((b==0)||(b==2))
{
unsigned n_node = Bulk_mesh_pt->nboundary_node(b);
for (unsigned n=0;n<n_node;n++)
{
Bulk_mesh_pt->boundary_node_pt(b,n)->pin(0);
}
}
}
The build of the "bulk" Young Laplace elements is completed by specifying the function pointers to the spine functions and the pointer to the Data
object that stores the curvature.
unsigned n_bulk=Bulk_mesh_pt->nelement();
for(unsigned i=0;i<n_bulk;i++)
{
ELEMENT *el_pt = dynamic_cast<ELEMENT*>(Bulk_mesh_pt->element_pt(i));
}
Finally, we complete the build of the contact line elements by passing the pointer to the double that stores the cosine of the contact angle.
unsigned nel=Contact_angle_mesh_pt->nelement();
for (unsigned e=0;e<nel;e++)
{
YoungLaplaceContactAngleElement<ELEMENT> *el_pt =
dynamic_cast<YoungLaplaceContactAngleElement<ELEMENT>*>(
Contact_angle_mesh_pt->element_pt(e));
}
All that's now left to do is to assign the equation numbers:
cout <<"\nNumber of equations: " << assign_eqn_numbers() << endl;
cout << "\n********************************************\n" << endl;
}
Creating the contact angle elements
The function create_contact_angle_elements()
attaches the FaceElements
that apply the contact angle boundary condition to the specified boundary of the "bulk" mesh. Pointers to the newly-created FaceElements
are stored in a separate mesh.
template<class ELEMENT>
const unsigned &b)
{
unsigned n_element = Bulk_mesh_pt->nboundary_element(b);
for(unsigned e=0;e<n_element;e++)
{
ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
Bulk_mesh_pt->boundary_element_pt(b,e));
int face_index = Bulk_mesh_pt->face_index_at_boundary(b,e);
YoungLaplaceContactAngleElement<ELEMENT>* contact_angle_element_pt = new
YoungLaplaceContactAngleElement<ELEMENT>(bulk_elem_pt,face_index);
Contact_angle_mesh_pt->add_element_pt(contact_angle_element_pt);
}
}
Deleting the contact angle elements
The function delete_contact_angle_elements()
deletes the contact angle elements and flushes the associated mesh.
template<class ELEMENT>
{
unsigned n_element = Contact_angle_mesh_pt->nelement();
for(unsigned e=0;e<n_element;e++)
{
delete Contact_angle_mesh_pt->element_pt(e);
}
Contact_angle_mesh_pt->flush_element_and_node_storage();
}
Post-processing
We output the load-displacement data, the meniscus shape, and various contact line quantities.
template<class ELEMENT>
ofstream& trace_file)
{
trace_file << -1.0*
Kappa_pt->value(0) <<
" ";
trace_file << Control_node_pt->value(0) ;
trace_file << endl;
unsigned npts=5;
ofstream some_file;
char filename[100];
sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
doc_info.number());
some_file.open(filename);
Bulk_mesh_pt->output(some_file,npts);
some_file.close();
ofstream tangent_file;
sprintf(filename,"%s/tangent_to_contact_line%i.dat",
doc_info.directory().c_str(),
doc_info.number());
tangent_file.open(filename);
ofstream normal_file;
sprintf(filename,"%s/normal_to_contact_line%i.dat",
doc_info.directory().c_str(),
doc_info.number());
normal_file.open(filename);
ofstream contact_angle_file;
sprintf(filename,"%s/contact_angle%i.dat",
doc_info.directory().c_str(),
doc_info.number());
contact_angle_file.open(filename);
Vector<double> tangent(3);
Vector<double> normal(3);
Vector<double> r_contact(3);
unsigned n_element = Contact_angle_mesh_pt->nelement();
for(unsigned e=0;e<n_element;e++)
{
tangent_file << "ZONE" << std::endl;
normal_file << "ZONE" << std::endl;
contact_angle_file << "ZONE" << std::endl;
YoungLaplaceContactAngleElement<ELEMENT>* el_pt =
dynamic_cast<YoungLaplaceContactAngleElement<ELEMENT>*>(
Contact_angle_mesh_pt->element_pt(e));
Vector<double> s(1);
for (unsigned i=0;i<npts;i++)
{
s[0]=-1.0+2.0*double(i)/double(npts-1);
dynamic_cast<ELEMENT*>(el_pt->bulk_element_pt())->
position(el_pt->local_coordinate_in_bulk(s),r_contact);
el_pt->contact_line_vectors(s,tangent,normal);
tangent_file << r_contact[0] << " "
<< r_contact[1] << " "
<< r_contact[2] << " "
<< tangent[0] << " "
<< tangent[1] << " "
<< tangent[2] << " " << std::endl;
normal_file << r_contact[0] << " "
<< r_contact[1] << " "
<< r_contact[2] << " "
<< normal[0] << " "
<< normal[1] << " "
<< normal[2] << " " << std::endl;
contact_angle_file << r_contact[1] << " "
<< el_pt->actual_cos_contact_angle(s)
<< std::endl;
}
}
tangent_file.close();
normal_file.close();
contact_angle_file.close();
cout << "\n********************************************" << endl << endl;
}
Comments and Exercises
How to generate a good initial guess for the solution
We already commented on the need to provide a "good" initial guess for the solution in order to ensure the convergence of the Newton iteration. In the previous example this was easy because the flat meniscus (clearly a solution of the Young-Laplace equations for zero curvature) also satisfied the boundary conditions. In the present example, and in many others, this is not the case. In such problems it may be difficult to generate initial guesses for the meniscus shape that are sufficiently close to actual solution.
In such cases it may be necessary to compute the initial solution to the problem whose behaviour we wish to investigate during the actual parameter study via a preliminary auxiliary continuation procedure that transforms an easier-solve-problem (for which a good initial guess can be found) into the actual problem.
Explore this approach in the present problem by implementing the following steps:
- Set the contact angle to and solve the problem, using the "flat" meniscus as the initial guess. The "flat" meniscus is, of course, the exact solution for zero control displacement and/or zero curvature.
- Now start a preliminary continuation procedure in which the contact angle is adjusted in small steps until it reaches the desired value. Keep the prescribed control displacement (or the meniscus curvature) constant during this procedure.
- The solution for the desired contact angle may now be used as the initial guess for the actual parameter study in which the control displacement (or the meniscus curvature) are increased while the contact angle is kept fixed.
Limitations of the current approach – suggestions for improvement
One of the main disadvantages of the approach adopted here is that the spine vector fields and must be specified a priori. For sufficiently complicated meniscus shapes (or for menisci that undergo large changes in shape as their curvature is varied) the choice of suitable spines may be very difficult.
One (possible) solution to this problem could be (we haven't tried it!) to occasionally update the spine representation. For instance, assume that we have computed a meniscus shape in the form
with an associated normal vector . We can reparametrise this shape by setting
and
before continuing the computation. Provided this is done sufficiently frequently, i.e. long before the displacement along the spines has become so large that the mapping from to is about to become non-one-to-one, this should allow the computation of arbitrarily large meniscus deflections. Try it out and let us know how it works!
Zero contact angles
Our problem formulation suffers from an additional, more fundamental problem: it cannot be used to solve problems with zero contact angle. This is because for zero contact angles the equilibrium solution is no longer a minimiser of the variational principle: given a solution at which the meniscus meets the wall at zero contact angle, it is always possible to extend the meniscus with an arbitrary-length "collar" along the wall without changing the overall energy of the system. As a result, the position of the contact line becomes increasingly ill-defined as the contact angle is reduced, causing the Newton method to converge very slowly (and ultimately not at all) as
Source files for this tutorial
PDF file
A pdf version of this document is available.