The classic Burmester problem is concerned with computing dimensions of planar four-bar linkages consisting of all revolute joints for five-pose problems. We define extended Burmester problem as the one where all types of planar four-bars consisting of dyads of type RR, PR, RP, or PP (R: revolute, P: prismatic) and their dimensions need to be computed for n-geometric constraints, where a geometric constraint is an algebraically expressed constraint on the pose, pivots, or something equivalent. In addition, we extend it to linear, nonlinear, exact, and approximate constraints. This extension also includes the problems when there is no solution to the classic Burmester problem, but designers would still like to design a four-bar that may come closest to capturing their intent. Machine designers often grapple with such problems while designing linkage systems where the constraints are of different varieties and usually imprecise. In this paper, we present (1) a unified approach for solving the extended Burmester problem by showing that all linear and nonlinear constraints can be handled in a unified way without resorting to special cases, (2) in the event of no or unsatisfactory solutions to the synthesis problem, certain constraints can be relaxed, and (3) such constraints can be approximately satisfied by minimizing the algebraic fitting error using Lagrange multiplier method. We present a new algorithm, which solves new problems including optimal approximate synthesis of Burmester problem with no exact solutions.
Introduction
According to an article by Professor McCarthy [1], over the past 40 yr, a total number of 2, 688 U.S. patents have been awarded that involve four-bar linkages, while in the same period, the next natural extension to more complex linkages, viz., six-bar linkages have been awarded a mere 84 patents. This is a clear demonstration of the popularity and ubiquity of planar four-bar linkages. Therefore, it is not surprising that the four-bar linkage synthesis and analysis problem still receives a considerable attention of researchers. Several text books, such as McCarthy and Soh [2], Sandor and Erdman [3], Hunt [4], Hartenberg and Denavit [5], Suh and Radcliffe [6], and Lohse [7], cover the science and art of planar four-bar and higher-order linkages. Kerle et al. [8] have given a historical overview of the development of mechanisms for motion generation. Despite a glut of literature on this topic, simultaneous computation of type and dimensions and accommodating practical geometric constraint for the motion generation problem has not been explored much. Erdman and Sandor in their seminal mechanism design and analysis text [9] clearly mention that assuming the wrong type (linkage topology and type of joints) to compute the dimensions of a linkage system may result in either none or suboptimal solutions.
This paper is a continuation of our research [10–12], wherein we have presented a task-driven approach to simultaneous type and dimensional synthesis of planar dyads for the motion generation problem. A four-bar linkage is constructed as combination of any two of the synthesized dyads. This dyadic construction simplifies the synthesis process and renders the method as modular building block for synthesis of mechanisms with more links such as six-bar mechanisms [10]. By using the concepts of kinematic mapping [13,14] and planar quaternions [15,16], we obtained a unified form of kinematic constraints of the planar dyads and created an algorithm for unified type and dimensional synthesis of planar four-bar linkages. This is accomplished via a two-step process. The first step is algebraic fitting of image points on a pencil of G-manifolds using singular value decomposition (SVD). This pencil of G-manifolds forms a candidate solution space for constraints accounted for in the SVD process. In the second step, we impose two fundamental quadratic conditions on the candidate solutions to extract the dyad types and their dimensions.
Ravani and Roth [17,18] were the first to use kinematic mapping approach for mechanism synthesis. Thereafter, Bodduluri and McCarthy [19], Bodduluri [20], Ge and Larochelle [21], Larochelle [22,23], Husty et al. [24], Hayes et al. [25,26], Wu et al. [27], Purwar and Gupta [28], and Schröcker et al. [29] have used this approach for the motion generation problem.
The original contribution of this paper is in the reformulation of our framework [10–12] in a general way to extend the Burmester problem by accommodating a variety of geometric constraints. In addition, the new formulation solves problems which our previous approach could not solve. These problems are (1) finding optimal approximate solutions for Burmester problem with no or suboptimal solutions, and (2) finding optimal linkages that minimize the algebraic fitting error of nonlinear geometric constraints. We note that in this paper, we do not restrict Burmester problem to only five poses. Burmester [30] showed that only a finite set of four-bars can be synthesized for five precision pose motion generation problems. In this paper, we show that other than poses, other practical geometric constraints can also be accommodated. For more than five geometric constraints, typically, only an approximate motion synthesis can be performed. Holte et al. [31], Sabada and coworkers [32], and Venkataraman et al. [33] have presented some techniques for mixed exact-approximate synthesis of planar mechanisms, albeit only for poses. Larochelle [34] has presented dimensional synthesis technique for solving the mixed exact and approximate motion synthesis problem for planar RR kinematic chains. Lin et al. [35] have presented pole curve transformation-based approach for motion synthesis of planar mechanisms. Al-Widyan et al. [36] have presented a numerically robust algorithm to solve the classic Burmester problem. Bourrelle et al. [37] presented a graphical user interface that uses the algorithm developed in Ref. [36] to solve the classic Burmester problem.
In general, there are two basic approaches to linkage synthesis for motion generation problem: (1) precision pose approach and (2) error minimization approach. In this paper, we use a combination of both methods in our framework.
Optimal synthesis of linkages is often a constrained nonlinear and multimodal problem in a multidimensional design space. In case of optimal linkage synthesis, relatively more work has been reported for path and function generation problem than for motion generation. Mariappan and Krishnamurty [38] used a generalized exact gradient method for planar mechanism synthesis. Vallejo et al. [39] developed an optimization method for planar mechanisms with lower pairs of any type, which uses error function as deformation of dimensions that mechanism has to undergo to perform the task. Yao and Angeles [40] solved path-generation problem using least-squared error function. They employ contour method to find out all stationery points of the problem. Kramer and Sandor [41] have presented a method of optimal design of planar mechanisms called selective precision synthesis. Selective precision synthesis can incorporate design requirements such as link ratios, fixed pivot locations, and transmission angle range. Gogate and Matekar [42] used method of differential evolution to find optimal linkages based on three error functions, which ensure that the mechanism is crank rocker, branch, and circuit defect free with minimized positional and angular deviations. Venkataraman et al. [33] used tolerance-based approach for four-position problem. They searched for optimal design parameters along suitable ranges of center point curve. Cabrera et al. [43] have presented genetic algorithm on planar 4R mechanism synthesis for its simple implementation and fast convergence. Hegedüs et al. [44] recently presented synthesis of spatial 6R linkages for interpolating four given poses using factorization of motion polynomials. All of the above said methods first select the type of linkage to be optimized, while our approach takes unified representation for all types of planar dyads into account. In contrast to other methods, which directly use linkage parameters, we formulate the optimization problem in terms of intermediate parameters obtained by geometric fitting. This allows us to use a two-step process wherein we first perform the least-squares fitting of the geometric constraints and then use a Lagrange multiplier method with additional linear and nonlinear constraints to extract the dyad types and dimensions. This approach allows formulation of an objective function consisting of squared error for the constraints to be minimized while keeping certain other constraints exactly satisfied. Ultimately, the minimization procedure leads to a system of quadratic equations, which can be solved using computer algebra software like mathematica. Each real solution of this system of equation corresponds to one optimum dyad; all such dyads are computed and then combined pairwise to obtain a set of four-bar linkage solutions.
The organization of the paper is as follows: Section 2 reviews the concept of kinematic mapping and planar quaternions as a special case of dual quaternions. Section 3 reviews unified form of kinematic constraints of the planar dyads as G-manifolds in the image space, while Sec. 4 presents various geometric constraints in the image space. Section 5 presents how we algebraically fit various linear geometric constraints to the G-manifolds using SVD. Section 6 presents the Lagrange multiplier method to minimize error functions while keeping certain linear or nonlinear constraints satisfied exactly. We finally present two practical examples in Sec. 7.
Kinematic Mapping
where and X = (X1, X2, X3) and L = (L1, L2, L3) are corresponding point and line coordinates in F.
Generalized (G-) Constraint Manifold
where (a0, a1, a2, a3) are the homogeneous coordinates of the constraint circle expressed in fixed reference frame, and (x, y) are the coordinates of the circle point expressed in moving reference frame. Here, a0 is the homogenizing factor. Equation (4) represents a generalized quadric (G-manifold) in the image space whose actual form would depend on the dyad type. For RR dyad, the quadric is a hyperboloid of one sheet, while for other dyads, it is a hyperbolic paraboloid [16].
Equations (4) and (5) are said to define the constraint manifold of all types of dyads (for details, see Ref. [10]). For a PR dyad, we have a0 = 0 and therefore, p1 = p2 = p3 = 0; for the RP dyad, we have p1 = p4 = p5 = 0; and for the PP dyad, we have p1 = p2 = p3 = p4 = p5 = 0. In all of these cases, the quadratic conditions in Eq. (5) are clearly satisfied. Thus, all planar dyads can be represented in the same form by Eqs. (4) and (5), and we can determine the type of a planar dyad by looking at the zeros in the coefficients pi (called signature of a dyad). This leads to a unified algorithm for simultaneous type and dimensional synthesis of planar dyads. In our approach, we first obtain the homogeneous coordinates pi, determine the dyad type from the signature of coefficient array pi, and then compute the dyad parameters using inverse relationships in Eq. (6). The coefficient array pi forms an eight-dimensional vector; henceforth, it is called a dyad vector p.
Task-Driven Geometric Constraints
In this section, we present various geometric constraints for the motion generation problem and show their representation in the image space. Some such constraints can be classified in following ways: (1) pose (position and orientation) constraint on the coupler, (2) constraint on moving or fixed pivot locations, namely, point, line, or general quadratic curve.
Pose Constraint.
For the classic Burmester problem, five poses would be specified, each of which would give one such equation.
Point Constraint.
where (xm, ym) are the coordinates of moving pivot location in moving reference frame.
Line Constraint.
Quadratic Curve Constraint.
The above is a homogeneous equation of a quadric in the image space and its degenerate forms reduce to the point and line constraints given in Eqs. (9) and (11). A similar equation can be obtained when the moving pivot is constrained to the quadratic curve of the form in Eq. (13). Without any loss of generalization, for the demonstration purposes, let us look at an elliptical-curve constraint now.
All PR dyads (p1 = p2 = p3 = 0) trivially satisfy Eq. (18), so extraneous PR dyads need to be filtered out.
Algebraic Fitting of Linear Constraints
where αi are (8 − n) homogeneous parameters; thus, without loss of generality, we can assume α1 = 1. For n = 5, the nullity of the matrix [A] is three; therefore, three singular vectors corresponding to the smallest singular-values are selected to form the dyad vector p. Substituting for p in the quadratic conditions given by Eq. (5), one can solve for two unknowns α2 and α3. Examining the signature of p, the dyad type is determined and by inverse computation in Eq. (7), the mechanism parameters are obtained. This fitting process also works for n > 5 cases; however, only an approximate solution may be obtained in this case. Typically, for large values of n, such solutions mostly produce unsatisfactory results.
Now, consider the possibility that no real or unacceptable solutions emerge from the above fitting process for n = 5 case. This is not an impractical scenario—very often, for five pose problems, the fitting process yields none, non-Grashof, or defective linkages. Another possibility is that the designer does not care about the five poses to be interpolated exactly. Manufacturing errors in links and play at joints would never lead to an exact interpolation of constraints anyways. For a pick and place operation, the first and last poses may be critical, while in between poses may be desired to be reached in a minimum-error sense. Moreover, the designer may want to introduce more important constraints in the problem, such as adding a location for pivots. The SVD also does not admit a nonlinear constraint of the form (13). In that case, the above formulation breaks down. We now need a method to relax certain constraints while simultaneously satisfying exact constraints. In Sec. 6, we show how such problems are solved.
Error Functions and Constraints
where ki are defined in terms of the parameters of dyad vector p and the known constraint parameters; e.g., Xf and Yf are the known constraint parameters for the fixed-pivot constraint.
where k0, ki, and kij are also defined in terms of the parameters of dyad vector p and the constraint parameters.
where k1ij and k2ij are defined in terms of dyad vector parameters pi.
Inequality Constraints.
Minimization Using Lagrange Multipliers.
In Sec. 5, we created algebraic fitting function for various linear and nonlinear constraints. Now, we present the Lagrange multiplier method to solve the optimization problem.
The system of polynomial equations comprised of Eqs. (30) and (31) is solved by computation of Grobner basis followed by eigen-system methods to extract numerical roots. This is achieved using Wolfram Mathematica's NSolve [46] routine. Solutions obtained are subjected to feasibility test using Eq. (32) to obtain feasible solutions, which give rise to a set of dyad vectors by substituting them into Eq. (21). Examining the signature of the dyad vector p, we can determine the dyad type and using the inverse relations in Eq. (7), we can obtain the mechanism parameters. A pool of mechanical dyads is obtained from the set of dyad vectors using inverse relations given in Eq. (7). Any two of these mechanical dyads can be combined to form a four-bar linkage as the solution for motion generation problem. A complete method using the Lagrange multiplier is given in Algorithm 1.
Algorithm 1: Algorithm for optimal four-bar linkage synthesis
1: procedureLagrange Multiplier Method
2: linear constraints (n) ≤ 5
3: for each linear constraint do
4: add constraint eq. to [A]
5: perform SVD and pick n − 8 Singular Vectors
6: f = 0
7: error function for each of the relaxed constraint → fi
8: for each fi
9:
10: form h1 and h2 using Eq. (24)
11: if nonlinear constraints exist then
12: number of nonlinear constraints = m
13: for each mdo
14: form h(2+m) using Eq. (28)
15: Minimizef subjected to h1 = 0,..., h(2+m) = 0, g ≤ 0
16: Minimize F by Lagrange Multipliers
17:
18: Take partial derivatives of F with respect to α2,..., α(8−n), λ1,..., λ(2+m) to form a system of Equations along with Eq. (31),
19: Solve (9 − n + m) equations
20: Compute dyad parameters
Examples
Now we present two examples, which demonstrate efficacy of our framework. We do not presume linkage types and determine best types and dimensions from the task requirements.
Optimal Solution for Five Positions With No Exact Solution.
Table 1 contains five precision poses as the input to the classic Burmester problem. Unfortunately for this problem, no exact planar four-bar solutions of any type exist. Thus, only approximate solutions can be computed by relaxing some of the constraints. Our previous approach in Ref. [10] fails to find exact or approximate solution for this problem because SVD cannot find approximate solutions for fully constrained problem. The new algorithm presented in this paper can deal with such problems. If the first and last poses are critical, algorithm can treat in-between poses to be approximate. However, when one exact constraint is relaxed, solution space increases by ∞1 along with number of optimization variables αi. Thus, the recommended approach is to relax as few constraints as possible to keep the optimization computationally cheaper. Hence, third precision pose is relaxed.
Poses | X | Y | ϕ (deg) |
---|---|---|---|
Pose 1 | −5.74803 | −0.00787402 | 88.5679 |
Pose 2 | −4.12598 | 0.795276 | 2.16642 |
Pose 3 | −2.72441 | 1.67717 | 356.968 |
Pose 4 | −1.54331 | 0.433071 | 1.03102 |
Pose 5 | 1.22835 | −0.590551 | 345.624 |
Poses | X | Y | ϕ (deg) |
---|---|---|---|
Pose 1 | −5.74803 | −0.00787402 | 88.5679 |
Pose 2 | −4.12598 | 0.795276 | 2.16642 |
Pose 3 | −2.72441 | 1.67717 | 356.968 |
Pose 4 | −1.54331 | 0.433071 | 1.03102 |
Pose 5 | 1.22835 | −0.590551 | 345.624 |
Now redefined task is to find a four-bar that interpolates four precision poses exactly and approximates third pose as closely as possible. The sense of approximation here is termed as algebraic fitting error of pose and constraint manifold of each dyad. First step is the algebraic fitting of geometric constraints, which in this case are first two and last two poses. Therefore, matrix [A] of size 4 × 8 is formed using Eq. (19) and four singular vectors corresponding to near-zero singular values are obtained by SVD, which are tabulated in Table 2. None of these singular vectors correspond to a mechanical dyad as they do not satisfy conditions in Eq. (5).
Dyad vector | p1 | p2 | p3 | p4 | p5 | p6 | p7 | p8 |
---|---|---|---|---|---|---|---|---|
p1 | −0.01874 | −0.2516 | 0.5174 | −0.3985 | 0.02237 | 0.03616 | 0.4980 | 0.5097 |
p2 | 0.007703 | 0.4284 | −0.3233 | 0.5047 | −0.02227 | −0.06264 | 0.4823 | 0.4688 |
p3 | 0.02871 | 0.08073 | 0.1114 | 0.1484 | −0.0120 | 0.9778 | −0.03979 | 0.01384 |
p4 | 0.07561 | −0.1628 | 0.02429 | 0.1901 | 0.9647 | −0.009039 | −0.007750 | 0.01223 |
Dyad vector | p1 | p2 | p3 | p4 | p5 | p6 | p7 | p8 |
---|---|---|---|---|---|---|---|---|
p1 | −0.01874 | −0.2516 | 0.5174 | −0.3985 | 0.02237 | 0.03616 | 0.4980 | 0.5097 |
p2 | 0.007703 | 0.4284 | −0.3233 | 0.5047 | −0.02227 | −0.06264 | 0.4823 | 0.4688 |
p3 | 0.02871 | 0.08073 | 0.1114 | 0.1484 | −0.0120 | 0.9778 | −0.03979 | 0.01384 |
p4 | 0.07561 | −0.1628 | 0.02429 | 0.1901 | 0.9647 | −0.009039 | −0.007750 | 0.01223 |
Dyad vectors are computed by substituting these solutions in Eq. (21). Table 3 contains these two dyad vectors. The four-bar linkages obtained by assembling these two dyads with coupler are shown in Fig. 1, where we can clearly see that coupler approximates third pose while interpolating remaining poses.
Optimal Linkage for Four Precision Poses With Region Constraint.
Figure 2 shows five positions of a landing gear moving from the landing position to the retracted position. Table 4 contains position and orientation data for five poses. It is desirable that fixed pivots should lie inside the circle of radius 2.3 with center located at (3.33, 2.04). The task is to synthesize a mechanism, which interpolates through precision poses (1, 2, 4, 5) and minimizes the algebraic error for the third pose while keeping fixed pivot locations inside the allowed region as shown in the figure.
Poses | X | Y | ϕ (deg) |
---|---|---|---|
Pose 1 | −0.0125 | −0.0374 | 66.3 |
Pose 2 | 0.303 | 0.634 | 35.5 |
Pose 3 | 0.599 | 1.83 | 352 |
Pose 4 | 0.268 | 2.30 | 331 |
Pose 5 | 0.606 | 1.31 | 22.2 |
Poses | X | Y | ϕ (deg) |
---|---|---|---|
Pose 1 | −0.0125 | −0.0374 | 66.3 |
Pose 2 | 0.303 | 0.634 | 35.5 |
Pose 3 | 0.599 | 1.83 | 352 |
Pose 4 | 0.268 | 2.30 | 331 |
Pose 5 | 0.606 | 1.31 | 22.2 |
First step is to extract all four geometric constraints, i.e., four precision poses and form matrix [A] using Eq. (19). Here, n = 4, which means solution space consists of four singular vectors, which are obtained using SVD and tabulated in Table 5. Once this linear algebraic fitting is done, optimization problem can be formulated.
Vector | p1 | p2 | p3 | p4 | p5 | p6 | p7 | p8 |
---|---|---|---|---|---|---|---|---|
p1 | −0.585 | −0.0211 | −0.506 | 0.165 | 0.134 | −0.350 | 0.368 | 0.313 |
p2 | 0.0640 | −0.330 | 0.190 | −0.804 | 0.236 | −0.264 | 0.194 | 0.204 |
p3 | −0.280 | 0.0484 | −0.466 | −0.398 | 0.232 | 0.603 | −0.136 | −0.329 |
p4 | −0.137 | 0.145 | 0.469 | 0.250 | 0.747 | 0.229 | 0.259 | 0.0111 |
Vector | p1 | p2 | p3 | p4 | p5 | p6 | p7 | p8 |
---|---|---|---|---|---|---|---|---|
p1 | −0.585 | −0.0211 | −0.506 | 0.165 | 0.134 | −0.350 | 0.368 | 0.313 |
p2 | 0.0640 | −0.330 | 0.190 | −0.804 | 0.236 | −0.264 | 0.194 | 0.204 |
p3 | −0.280 | 0.0484 | −0.466 | −0.398 | 0.232 | 0.603 | −0.136 | −0.329 |
p4 | −0.137 | 0.145 | 0.469 | 0.250 | 0.747 | 0.229 | 0.259 | 0.0111 |
Solving these equations followed by filtering on the basis of feasibility using Eq. (32) yields four unique and feasible solutions tabulated in Table 6. All of these solutions satisfy Karush–Kuhn–Tucker condition for optimality given by Eq. (33). Dyad vectors are calculated by substituting these solutions into Eq. (21) and are given in Table 7. Any of these four dyad vectors when substituted in Eq. (4) forms a quartic equation, which when projected on hyperplane Z4 = 1 represents a quadric surface. Figure 3 shows intersection of hyperboloid and hyperbolic paraboloid formed from first and second dyad vectors. The intersection curve represents workspace of the corresponding four-bar linkage. Table 8 contains the minimized algebraic fitting error of objective function. From this table, we can see that dyad s4 has least pose fitting error. All dyads except s1 are of RR type dyads while s1 is an RP dyad. Figure 4 shows a branch defect free four-bar mechanism formed by combining s1 and s2.
Dyad | α1 | α2 | α3 | α4 | λ1 | λ2 | μ |
---|---|---|---|---|---|---|---|
s1 | 1 | −0.40 | −0.027 | 1.09 | 0.0 | 0.0 | 0.0 |
s2 | 1 | −2.04 | −1.81 | 2.06 | 0.0 | 0.0 | 0.0 |
s3 | 1 | −3.64 | −4.24 | 5.37 | 0.0 | 0.0 | 0.0 |
s4 | 1 | 5.08 | −2.58 | −2.94 | 0.0 | 0.0 | 0.0 |
Dyad | α1 | α2 | α3 | α4 | λ1 | λ2 | μ |
---|---|---|---|---|---|---|---|
s1 | 1 | −0.40 | −0.027 | 1.09 | 0.0 | 0.0 | 0.0 |
s2 | 1 | −2.04 | −1.81 | 2.06 | 0.0 | 0.0 | 0.0 |
s3 | 1 | −3.64 | −4.24 | 5.37 | 0.0 | 0.0 | 0.0 |
s4 | 1 | 5.08 | −2.58 | −2.94 | 0.0 | 0.0 | 0.0 |
Vector | p1 | p2 | p3 | p4 | p5 | p6 | p7 | p8 |
---|---|---|---|---|---|---|---|---|
s1 | 0.0254 | 0.192 | −0.202 | −0.0965 | 0.00562 | 0.723 | −0.387 | −0.490 |
s2 | 0.276 | −0.406 | 0.211 | −0.425 | −0.161 | −0.561 | 0.251 | 0.360 |
s3 | 0.315 | −0.356 | 0.189 | −0.397 | −0.318 | −0.598 | 0.129 | 0.324 |
s4 | 0.208 | −0.302 | 0.141 | −0.335 | −0.321 | −0.691 | 0.134 | 0.368 |
Vector | p1 | p2 | p3 | p4 | p5 | p6 | p7 | p8 |
---|---|---|---|---|---|---|---|---|
s1 | 0.0254 | 0.192 | −0.202 | −0.0965 | 0.00562 | 0.723 | −0.387 | −0.490 |
s2 | 0.276 | −0.406 | 0.211 | −0.425 | −0.161 | −0.561 | 0.251 | 0.360 |
s3 | 0.315 | −0.356 | 0.189 | −0.397 | −0.318 | −0.598 | 0.129 | 0.324 |
s4 | 0.208 | −0.302 | 0.141 | −0.335 | −0.321 | −0.691 | 0.134 | 0.368 |
Conclusion
In this paper, we presented a task-driven approach to unified and optimal synthesis of planar four-bar linkages for extended Burmester problem. In this formulation, various geometric constraints are treated equivalently, which in turn leads to a much simpler two-step based algorithm for computing planar dyads of four-bar linkages. Original contributions of this paper have been into reforming a mixed exact-approximate algebraic fitting problem into problem of task oriented optimal fitting of algebraic manifold. The framework presented here can accommodate linear as well as nonlinear equality and inequality geometric constraints and minimize objective functions that can be expressed in terms of dyadic parameters. Although adding nonlinear geometric constraints increase computational complexity, computer algebra software like mathematica could be used to compute solutions of quadratic system of equations in a reasonable amount of time. Experimentations show that mathematica takes less than 3 s on a MacBook Pro with 2.4 GHz Intel core i5 processor and 8 GB RAM for computing solutions for the system of seven quadratic equations. The framework also preserves previously achieved real-time solutions for linear geometric constraints with no optimality criterion. Two examples demonstrating computation of optimal type and dimensions of dyads that minimize task oriented objective function are presented.
Acknowledgment
This work has been financially supported by National Science Foundation under a research grant to Stony Brook University (A. Purwar and Q. J. Ge). All findings and results presented in this paper are those of the authors and do not represent those of the funding agencies.
Funding Data
Directorate for Engineering (Grant No. CMMI-1563413).