When parametric functions are used to blend 3D surfaces, geometric continuity of displacements and derivatives until to the surface boundary must be satisfied. By the traditional blending techniques, however, arbitrariness of the solutions arises to cause a difficulty in choosing a suitable blending surface. Hence to explore new blending techniques is necessary to construct good surfaces so as to satisfy engineering requirements. In this paper, a blending surface is described as a flexibly elastic plate both in partial differential equations and in their variational equations, thus to lead to a unique solution in a sense of the minimal global surface curvature. Boundary penalty finite element methods (BP-FEMs) with and without approximate integration are proposed to handle the complicated constraints along the blending boundary. Not only have the optimal convergence rate $O (h^2)$ of second order generalized derivatives of the solutions in the solution domain been obtained, but also the high convergence rate $O (h^4)$ of the tangent boundary condition of the solutions can be achieved, where $h$ is the maximal boundary length of rectangular elements used. Moreover, useful guidance in computation is discovered to deal with interpolation and approximation in the boundary penalty integrals. A numerical example is also provided to verify perfectly the main theoretical analysis made. This paper yields a framework of mathematical modelling, numerical techniques and error analysis to the general and complicated blending problems.