Wafer Chucking Model

Question

I am trying to calculate the pressure required to flatten a bowed wafer. The simple model I am using for the study has a bottom layer (the bowed wafer) pushed against a harder flat surface as shown in Fig. 1.

Using the plane stress operator from the FEM tutorial I get a solution where the upper surface of the wafer protrudes from the bottom surface of the harder layer (when the pressure is high enough to overcome the gap). Fig. 2, where I combined the starting mesh boundary, in black, and the deformed one, in red, shows what I mean.

I am not a structural engineer, so I do not fully appreciate the subtleties of setting up a model for this case. Perhaps additional boundary conditions on the gap between the layers are needed. I could envision a condition involving inequalities (such as the difference between the top and the bottom deformation be less than the gap distance), but I am not sure how it could be accomplished in Mathematica's FEM framework. I am wondering how a solution that does not have this issue could be obtained using Mathematica FEM framework?

Here is the Mathematica code used for this example.

Needs["NDSolve`FEM`"];
pi = 3.1415926535897932384626433832795028845;
(* linear units *)
mm = N[1];
um = N[10^-3];
psi2GPa = 6.894757313303502*10^(-6);
(* structure geometry *)
h1 = 0.2 mm;
h2 = h1 + 1 mm;
s = 10 um;
w = 30 mm;
bmesh = ToBoundaryMesh[
   ImplicitRegion[(y >= -s/2 (1 + Cos[(2 pi)/w x]) && 
       y <= h1 - s/2 (1 + Cos[(2 pi)/w x])) || (y > h1 && 
       y <= h2), {{x, -w/2, w/2}, y}]];
mesh = ToElementMesh[bmesh, MaxCellMeasure -> 1*10^-3];
(* uncomment section below to visualize the mesh *)
(*
bmesh["Wireframe"[ImageSize\[Rule]"Large"]]
Show[bmesh["Wireframe"[ImageSize\[Rule]400,AspectRatio\[Rule]0.2]],\
PlotRange\[Rule]{All,{0-s-2 um,h1+5 um}}]
bmesh["Wireframe"["MeshElement"\[Rule]"BoundaryElements",\
"MeshElementMarkerStyle"\[Rule]Blue,ImageSize\[Rule]500,AspectRatio\
\[Rule]0.2]]
Show[bmesh["Wireframe"["MeshElement"\[Rule]"BoundaryElements",\
"MeshElementMarkerStyle"\[Rule]Blue,ImageSize\[Rule]500]],PlotRange\
\[Rule]{All,{h1-s-1 um,h1+5 um}},AspectRatio\[Rule]0.2]
Show[mesh["Wireframe"],ImageSize\[Rule]"Large",PlotRange\[Rule]{{-100 \
um,100 um},{h1-50 um,h1+50 um}}]
*)
(* Material parameters *)
Y1 = 64;
nu1 = 0.2;
Y2 = 345;
nu2 = 0.29;
With[{h1 = h1, h2 = h2,
   Y1 = Y1, Y2 = Y2,
   nu1 = nu1, nu2 = nu2},
  Y = If[y <= h1, Y1, Y2];
  nu = If[y < h1, nu1, nu2];
  ];
(* Boundary Conditions *)
Subscript[\[CapitalGamma], D] = 
  DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, y == h2];
psi = 0.015;  (* applied pressure in psi *)
p = psi*psi2GPa;
Subscript[\[CapitalGamma], N] = {0, 
  NeumannValue[p, 
   ElementMarker == 
    4]}; (* ElementMarker 4 is the bottom surface  *)
(* plane stress operator *)
pde = {Inactive[
     Div][({{0, -((Y nu)/(1 - 
             nu^2))}, {-((Y (1 - nu))/(2 (1 - nu^2))), 0}} . 
      Inactive[Grad][v[x, y], {x, y}]), {x, y}] + 
   Inactive[
     Div][({{-(Y/(1 - nu^2)), 
        0}, {0, -((Y (1 - nu))/(2 (1 - nu^2)))}} . 
      Inactive[Grad][u[x, y], {x, y}]), {x, y}], 
  Inactive[Div][({{0, -((Y (1 - nu))/(2 (1 - nu^2)))}, {-((Y nu)/(1 - 
             nu^2)), 0}} . Inactive[Grad][u[x, y], {x, y}]), {x, y}] +
    Inactive[
     Div][({{-((Y (1 - nu))/(2 (1 - nu^2))), 
        0}, {0, -(Y/(1 - nu^2))}} . 
      Inactive[Grad][v[x, y], {x, y}]), {x, y}]};
(* Solve  *)
{ufun, vfun} = 
  NDSolveValue[{pde == Subscript[\[CapitalGamma], N], 
    Subscript[\[CapitalGamma], D]}, {u, v}, {x, y} \[Element] mesh];
(* show the undeformed and deformed mesh boundary near the gap *)
Show[{mesh["Wireframe"["MeshElement" -> "BoundaryElements"]], 
  ElementMeshDeformation[mesh, {ufun, vfun}, "ScalingFactor" -> 1][
   "Wireframe"["MeshElement" -> "BoundaryElements", 
    "MeshElementStyle" -> Red]]}, ImageSize -> 500, 
 AspectRatio -> 0.5, PlotRange -> {All, {h1 - s - 2 um, h1 + 5 um}}]

Update: Given that this problem seems to involve internal boundary conditions, I am starting to think that it would be too difficult to solve using the current Mathematica finite element framework (I would love to be proven wrong). Although that remains my problem at hand, for the sake of making some progress, I made it simpler by eliminating the top (harder) layer (as also suggested by Oleksii's comment). The new geometry is shown in Fig. 3 (it is the same thin wafer, but in the graph the aspect ratio was modified to make the bow more evident).

There are only two points of contact and a Dirichlet condition can only be defined for those two points. To account for the presence of the hard constrain in between these two points, we can define a fictitious external pressure that changes from zero to a large value as soon as the displacement gets close to the gap value. This can be done by defining a non-linear Neumann value. The transition from zero to a large pressure value can be implemented using a smoothed step function. The attached code is the implementation of this idea. Figg. 4 and 5 show the solution without and with the fictitious pressure. Unfortunately I get a warning message ("The line search decreased the step size to within tolerance specified by AccuracyGoal..."), so I was wondering if anybody could help fix this numerical issue.

Needs["NDSolve`FEM`"];
SmoothedStepFunction[fmin_, fmax_, ts_, m_] :=
 Function[t, (fmin*Exp[ts*m] + fmax*Exp[m*t])/(Exp[ts*m] + Exp[m*t])]
steepnessFactor = 20000; 
\[Theta] = SmoothedStepFunction[1, 0, 0, steepnessFactor]; 
pi = 3.1415926535897932384626433832795028845;
(* linear units *)
mm = N[1];
um = N[10^-3];
psi2GPa = 6.894757313303502*10^(-6);
(* structure geometry *)
h1 = 0.2 mm;
s = 10 um;
w = 30 mm;
bmesh1 = ToBoundaryMesh[
   ImplicitRegion[(y >= -s/2 (1 + Cos[(2 pi)/w x]) && 
      y <= h1 - s/2 (1 + Cos[(2 pi)/w x])), {{x, -w/2, w/2}, y}]];
mesh1 = ToElementMesh[bmesh1, MaxCellMeasure -> 1*10^-3];
(* uncomment section below to visualize the mesh *)
(*
bmesh1["Wireframe"[ImageSize\[Rule]"Large"]]
Show[bmesh1["Wireframe"[ImageSize\[Rule]"Large",AspectRatio\[Rule]0.2]\
]]
bmesh1["Wireframe"["MeshElement"\[Rule]"BoundaryElements",\
"MeshElementMarkerStyle"\[Rule]Blue,ImageSize\[Rule]500,AspectRatio\
\[Rule]0.2]]
Show[mesh1["Wireframe"],ImageSize\[Rule]"Large"]
*)
(* Material parameters *)
Y1 = 64;
nu1 = 0.2;
Y2 = 345;
nu2 = 0.29;
With[{h1 = h1, h2 = h2,
   Y1 = Y1, Y2 = Y2,
   nu1 = nu1, nu2 = nu2},
  Y = If[y <= h1, Y1, Y2];
  nu = If[y < h1, nu1, nu2];
  ];
(* Boundary Conditions *)
Subscript[\[CapitalGamma], D] = 
  DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, 
   y == h1(*&&(x\[Equal]-w/2||x\[Equal]w/2)*)];
psi = 0.1;  (* applied pressure in psi *)
p = psi*psi2GPa;
Subscript[\[CapitalGamma], N1] = {0, 
  NeumannValue[p, 
   ElementMarker == 
    1]};  (* ElementMarker 1 is the bottom surface  *)Subscript[\
\[CapitalGamma], N2] = {0, 
  NeumannValue[-100 (1 - \[Theta][
       v[x, y] - s/2 (1 + Cos[(2 pi)/w x])]), (-w/2 < x < w/2) && 
    y == h1 - s/2 (1 + Cos[(2 pi)/w x])]}; 
(* plane stress operator *)
pde = {Inactive[
      Div][({{0, -((Y nu)/(1 - 
              nu^2))}, {-((Y (1 - nu))/(2 (1 - nu^2))), 0}} . 
       Inactive[Grad][v[x, y], {x, y}]), {x, y}] + 
    Inactive[
      Div][({{-(Y/(1 - nu^2)), 
         0}, {0, -((Y (1 - nu))/(2 (1 - nu^2)))}} . 
       Inactive[Grad][u[x, y], {x, y}]), {x, y}], 
   Inactive[
      Div][({{0, -((Y (1 - nu))/(2 (1 - nu^2)))}, {-((Y nu)/(1 - 
              nu^2)), 0}} . Inactive[Grad][u[x, y], {x, y}]), {x, 
      y}] + Inactive[
      Div][({{-((Y (1 - nu))/(2 (1 - nu^2))), 
         0}, {0, -(Y/(1 - nu^2))}} . 
       Inactive[Grad][v[x, y], {x, y}]), {x, y}]};
(* Solve  without the fictitious pressure*)
{u0, v0} = 
  NDSolveValue[{pde == Subscript[\[CapitalGamma], N1], 
    Subscript[\[CapitalGamma], D]}, {u, v}, {x, y} \[Element] mesh1];
(* show the undeformed and deformed mesh boundary *)
Print["Deformation without the fictitious pressure"];
Show[{mesh1["Wireframe"["MeshElement" -> "BoundaryElements"]], 
  ElementMeshDeformation[mesh1, {u0, v0}, "ScalingFactor" -> 1][
   "Wireframe"["MeshElement" -> "BoundaryElements", 
    "MeshElementStyle" -> Red]]}, ImageSize -> 500, AspectRatio -> 0.3]
(* Solve  with the fictitious pressure*)
{u1, v1} = 
  NDSolveValue[{pde == 
     Subscript[\[CapitalGamma], N1] + Subscript[\[CapitalGamma], N2], 
    Subscript[\[CapitalGamma], D]}, {u, v}, {x, y} \[Element] mesh1, 
   Method -> {"PDEDiscretization" -> {"FiniteElement", 
       "PDESolveOptions" -> {"FindRootOptions" -> {Method -> \
{"Newton"}}}}}];
Print["Deformation with the fictitious pressure"];
Show[{mesh1["Wireframe"["MeshElement" -> "BoundaryElements"]], 
  ElementMeshDeformation[mesh1, {u1, v1}, "ScalingFactor" -> 1][
   "Wireframe"["MeshElement" -> "BoundaryElements", 
    "MeshElementStyle" -> Red]]}, ImageSize -> 500, 
 AspectRatio -> 0.3]

Second Update: Looking at the user21 result, I also thought that the issue was the inability of the top surface to move in the x direction due to the Dirichlet condition u[x,y]==0 when y==h1. I think that, due to the symmetry, an appropriate replacement is u[x,y]==0 when (y<0 && x==0). Here are the modifications to the code.

delta = w/500;  (* a small fraction of w *)
(* Boundary Conditions *)
dbc = {DirichletCondition[v[x, y] == 0, y == h1], 
  DirichletCondition[u[x, y] == 0, 
   y <= 0 && -delta/2 <= x <= delta/2]};
nbc1 = {0, NeumannValue[p, y <= 0]};  nbc2 = {0, 
  NeumannValue[
   If[(y + v[x, y]) < h1, 0, 1/10^4*(0 - v[x, y])], (-w/2 < x < w/2) &&
     y >= h1 - s/2 (1 + Cos[(2 Pi)/w x])]}; 
{u1, v1} = 
  NDSolveValue[{pde == nbc1 + nbc2, dbc, Sequence[]}, {u, 
    v}, {x, y} \[Element] mesh1, Method -> {"PDEDiscretization" -> 
{"FiniteElement", 
    "PDESolveOptions" -> {"FindRootOptions" -> {Method -> {"Newton"}}}}
}];

Even though now the top surface can slide in the x direction, as shown in the plot detailing the top left corner, the solution still penetrates into the constraining surface. The FindRoot::lstol error message is generated during the computation. If the "Newton" method is not specified, then FindRoot fails.

Another interesting thing is that if I use the user21's code, i.e. with the problematic Dirichlet condition, and specify "Newton" in the FindRootOptions, then the solution is very similar to Fig. 7 (except for the top corners), i.e. the dip in the middle is below the blue line (a FindRoot::lstol error message is issued even in this case). The fact that we are getting two solutions (with different error messages) depending on the method used does not give me confidence that we are getting closer to the correct solution.

Answer 1

Suppose that upper body is ideally rigid and provides support of roller type in connection area. Also assume that contact area is very small compared with characteristic dimensions of the system. Thus we can suppose that connection takes place in two points (lines in 3D consideration). But these two points don't necessarily have to coincide with the corners of the wafer upper surface. The coordinates of the contact points have to be determined. The function describing upper surface of the wafer have to reach a local maximum in contact points. On the basis of this criteria we will find their coordinates. So, we will iterate over different pairs of points located on upper surface of the wafer. I do not know how to realize Dirichlet condition $v=0$ in single points when using NDSolveValue. But it is possible to implement it under the low level by means of finite element programming procedures MMA FE programming techniques. For this problem it is also convenient to make use of structured mesh which can be generated by function StructuredMesh from library FEMAddOns which in turn can be downloaded from Here.

Needs["NDSolve`FEM`"]
Needs["FEMAddOns`"]

Properties of the material

Y1 = 64;
nu1 = 0.2;
Y2 = 345;
nu2 = 0.29;             

Geometry of computational domain

h1 = 0.2 ;
h2 = h1 + 1 ;
s = 10 *10^-3;
w = 30 ;

Mesh generation

Clear[WafBottom, WafTop]
WafBottom[x_] := -s/2 (1 + Cos[(2 \[Pi])/w x]); (*bottom surface of wafer*)
WafTop[x_] := h1 - s/2 (1 + Cos[(2 \[Pi])/w x]);(*upper surface of wafer*)

Nx = 200;
NyW = 4;
hx = w/Nx;
raster1 = {
            Table[{x, WafBottom[x]}, {x, -w/2, w/2, w/Nx}],
            Table[{x, WafTop[x]}, {x, -w/2, w/2, w/Nx}]
          };

mesh1 = StructuredMesh[raster1, {Nx, NyW}]; (*mesh in wafer domain *)

NyH = 10;
raster2 = {{{-w/2, h1}, {w/2, h1}}, {{-w/2, h1 + h2}, {w/2, h1 + h2}}};
mesh2 = StructuredMesh[raster2, {Nx, NyH}];  (*mesh in rigid body*)

Show[
     mesh2["Wireframe"["MeshElementStyle" -> FaceForm[Blue]]],
     mesh1["Wireframe"["MeshElementStyle" -> FaceForm[Green]]],
     ImageSize -> 700
    ]

Extracting of important mesh characteristics

Vert = mesh1["Coordinates"];
Quads = mesh1["MeshElements"][[1]][[1]];
Bound = mesh1["BoundaryElements"][[1]][[1]];

TopEl = {};    (*elements on top surface*)

BottomEl = {}; (*elements on bottom surface*)
Do[
   {x, y} = Mean[Vert[[Bound[[i]]]]];
   If[Abs[y - WafTop[x]] <= 0.01*hx, 
      TopEl = Join[TopEl, {Bound[[i]]}]];
   If[Abs[y - WafBottom[x]] <= 0.01*hx, 
      BottomEl = Join[BottomEl, {Bound[[i]]}]];
   , {i, 1, Length[Bound]}
  ]
TopSurfaceNodes = DeleteDuplicates[Flatten[TopEl, 1]];
BottomSurfaceNodes = DeleteDuplicates[Flatten[BottomEl, 1]];

Applied load

psi = 0.015;
psi2GPa = 6.894757313303502*10^(-6);
NormPress = psi*psi2GPa;

Low level FE programming procedures

vd = NDSolve`VariableData[{"DependentVariables", 
     "Space"} -> {{u, v}, {x, y}}];
sd = NDSolve`SolutionData["Space" -> ToNumericalRegion[mesh1]];

methodData = InitializePDEMethodData[vd, sd];
diffusionCoefficients = 
  "DiffusionCoefficients" -> {{{{-(Y/(1 - \[Nu]^2)), 
        0}, {0, -((Y (1 - \[Nu]))/(2 (1 - \[Nu]^2)))}}, {{0, -((Y \
\[Nu])/(1 - \[Nu]^2))}, {-((Y (1 - \[Nu]))/(2 (1 - \[Nu]^2))), 
        0}}}, {{{0, -((Y (1 - \[Nu]))/(2 (1 - \[Nu]^2)))}, {-((Y \
\[Nu])/(1 - \[Nu]^2)), 0}}, {{-((Y (1 - \[Nu]))/(2 (1 - \[Nu]^2))), 
        0}, {0, -(Y/(1 - \[Nu]^2))}}}} /. {Y -> Y1, \[Nu] -> nu1};    

initCoeffs = 
  InitializePDECoefficients[vd, sd, {diffusionCoefficients}];
discretePDE = DiscretizePDE[initCoeffs, methodData, sd]

(*impose Neumann boundary conditions on bottom wafer surface *)

GammaNv = NeumannValue[NormPress, y == WafBottom[x]];
initBCs = InitializeBoundaryConditions[vd, sd, {{}, {GammaNv}}];
 discreteBCs = DiscretizeBoundaryConditions[initBCs, methodData, sd];

split = Span @@@ 
   Transpose[{Most[# + 1], Rest[#]} &[methodData["IncidentOffsets"]]];

(*extract stiffness matrix and load vector*)
{load, stiffness} = Take[discretePDE["SystemMatrices"], 2];
(*taking into accaunt of loading. stifness matrix does not changes*)
 DeployBoundaryConditions[{load, stiffness}, discreteBCs]; 

Searching for points satisfying the criterion

PlotArr = {}; 
IterNumber = 50;
UnknNumber = 
  Length[load];(*number of unknowns*)                                 \
            
ColumnInd = Table[i, {i, 1, UnknNumber}];          
     

Do[ 
     DirichletUnknownPosit = {};
     DirichletValues = {};
 
 (*x-coordinates of the node touching the rigid body *)
      
 xfix = w/2 - iter*hx;
 (*collection of info for realization of Dirichlet boundary \
conditions*)
      Do[
            k = TopSurfaceNodes[[i]];
            x = Vert[[k]][[1]]; 
            If[Abs[x - xfix] <= 0.01 hx || Abs[x + xfix] <= 0.01 hx,
                  NumV = Length[Vert] + k;
                  AppendTo[DirichletUnknownPosit, NumV];
                 AppendTo[DirichletValues, h1 - WafTop[x]]
               ]
           , {i, 1, Length[TopSurfaceNodes]} 
          ];
 
 (*realization of Dirichlet boundary conditions. such procedure makes \
an appropriate changes in stiffness matrix*)
     
 diagonal = ConstantArray[0, UnknNumber];
       diagonal[[DirichletUnknownPosit]] = 1.;
       SparseDiagon = 
  SparseArray[Table[{i, i} -> diagonal[[i]], {i, 1, UnknNumber}]];
        A = stiffness;
      A[[DirichletUnknownPosit, ColumnInd]] = 0.; 
        A = 
  A + SparseDiagon;                                                   \
       
       b = Flatten[load];
     b[[DirichletUnknownPosit]] = DirichletValues;
 (*solution of system of linear equations*)
     Clear[solution]; 
    solution = LinearSolve[A, b, Method -> "Pardiso"];
 (*displacements of mesh nodes*)
  
 displ = Transpose[{solution[[split[[1]]]], solution[[split[[2]]]]}];
 (*coordinates of nodes of deformed mesh*)
 
 VertDeformed = Vert + displ;
 
 pic = ListLinePlot[{VertDeformed[[TopSurfaceNodes]], 
    Vert[[TopSurfaceNodes]]}, 
   GridLinesStyle -> Directive[Thickness[0.005], Dashed], 
   GridLines -> {{-xfix, xfix}, {h1}}, 
   PlotRange -> {h1 - s, h1 + 2 s}, ImageSize -> 600, Frame -> True, 
   FrameLabel -> {"x", "y"}, FrameStyle -> RGBColor[0, 0, 0], 
   BaseStyle -> 20, PlotLegends -> {"under load", "initial shape"}, 
   PlotLabel -> "shape of the the wafer upper surface"];
 
 AppendTo[PlotArr, pic]
     , {iter, 0, IterNumber}
  ]

Manipulate[PlotArr[[i]], {i, 1, IterNumber, 1}]

By means of given iteration procedure we can find contact points and bending of the wafer that corresponds to such a position of supports.

Thus, if points with coordinates $x\approx \pm 11.25$ (for more precision finer FE mesh should be used) are chosen as a contacts, the wafer does not "penetrate" into rigid body. We can go further, considering deformation of upper body (substrate). Of course in such a case the upper body could no longer be considered as ideally rigid body. For this purpose the concentrated force have to be applied in contacts. Thereby boundary conditions for upper body will be the next:

$\sigma_{yy}(x,h_1)=0.5\cdot\delta(x-x_{c1})\cdot P\cdot w+0.5\cdot\delta(x-x_{c2})\cdot P\cdot w$, where $P$ is a loading applied to wafer, $x_{c1}, x_{c2}$ are the contacts coordinates which are already determined. $\sigma_{xx}(\pm 0.5w,y)=0, u(x,h_1+h_2)=v(x,h_1+h_2)=0$

Displacements $v$ ($y$ component of the displacement) on the bottom surface of the substrate under the different positions of the contact are the next

Corresponding 2D field $v(x,y)$ inside substrate in the vicinity of the contact looks like this

As we can see in substrate $v_{max} \approx 10^{-8}$ that is much smaller than displacements in wafer. Hereby our initial assumption of ideally rigid substrate works quite well.

Answer 2

This is not a solution but a somewhat crazy idea to try. Before I explain the idea, I'd like to show how one can hack a DirichletCondition with a NeumannValue. This is in general not a good idea to pursue, since the method is not consistent. Consider this:

u0 = 0;
if1 = NDSolveValue[{-Laplacian[u[x, y], {x, y}] == 1, 
   DirichletCondition[u[x, y] == u0, True]}, 
  u, {x, y} \[Element] Disk[]]

This solves a Poisson type equation over a Disk. Next we solve the same problem but replace the DirichletCondition with a NeumannValue.

if2 = NDSolveValue[{-Laplacian[u[x, y], {x, y}] == 
     1 + NeumannValue[1/\[Epsilon] (u0 - u[x, y]), 
       True] /. \[Epsilon] -> 10^-3}, u, {x, y} \[Element] Disk[]]

So we approximate the Dirichlet condition with 1/epsilon * (u0 - u). This method is the Ritz approximation and the basis for Nitsche's Method (not discussed in this post, but here)

Again, this should not be used in general but it works in this case:

Plot3D[if1[x, y] - if2[x, y], {x, y} \[Element] Disk[], 
 PlotRange -> All]

The idea is the following: Model just the domain of the waver. At the tilted boundary use {0, NeumannValue[1/f[v[x,y],{x,y}] * (0 - v[x,y]) , pred]}. f[v[x,y],{x,y}] computes a number that gets smaller and smaller the closer you get to h1 and continues to impose a constraint on v such that no movement beyond h1 is possible.

I have not tried this, it's just an idea that may or may not work.

Update: Here is a prototype.

Needs["NDSolve`FEM`"];
(*linear units*)
mm = 1;
um = 10^-3;
psi2GPa = 6.894757313303502*10^(-6);
(*structure geometry*)
h1 = 2/10 mm;
s = 10 um;
w = 30 mm;
region = ImplicitRegion[(y >= -s/2 (1 + Cos[(2 Pi)/w x]) && 
     y <= h1 - s/2 (1 + Cos[(2 Pi)/w x])), {{x, -w/2, w/2}, y}];
bmesh1 = ToBoundaryMesh[region];
bmesh1["Wireframe"[AspectRatio -> 0.2]]

mesh1 = ToElementMesh[bmesh1, MaxCellMeasure -> 1*10^-3];

Y1 = 64;
nu1 = 2/10;
Y = Y1;
nu = nu1;
psi = 0.1;(*applied pressure in psi*)
p = psi*psi2GPa;
pde = {Inactive[
      Div][({{0, -((Y nu)/(1 - 
              nu^2))}, {-((Y (1 - nu))/(2 (1 - nu^2))), 0}} . 
       Inactive[Grad][v[x, y], {x, y}]), {x, y}] + 
    Inactive[
      Div][({{-(Y/(1 - nu^2)), 
         0}, {0, -((Y (1 - nu))/(2 (1 - nu^2)))}} . 
       Inactive[Grad][u[x, y], {x, y}]), {x, y}], 
   Inactive[
      Div][({{0, -((Y (1 - nu))/(2 (1 - nu^2)))}, {-((Y nu)/(1 - 
              nu^2)), 0}} . Inactive[Grad][u[x, y], {x, y}]), {x, 
      y}] + Inactive[
      Div][({{-((Y (1 - nu))/(2 (1 - nu^2))), 
         0}, {0, -(Y/(1 - nu^2))}} . 
       Inactive[Grad][v[x, y], {x, y}]), {x, y}]};

Now, this is where things start to get interesting. The boundary conditions are as follows:

nbc = {0, NeumannValue[p, y <= 0]};

This DirichletCondition worries me a bit - this means that the top left and right hand corner are fixed. They can not move. Now, if you imagine a plate to be pressed in between them it's not clear where the body should move to, since it's fixated at these point. The waver, however, is not fixated in any manner like that. You'd need to think a bit about that.

dbc1 = DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, y == h1];
dbc2 = Sequence[];

I did experiment a bit with other Dirichlet BCs but could not come up with something satisfactory. And since this is not the main question, I left that as is.

Next, come the contract constraint. I have implemented this as a discontinuous contact. If y+v[x,y] < h1 then we use a Neumann 0 condition, else we impose a (hacked) 'weak form' of a Dirichlet condition, like shown above.

N2 = {0, NeumannValue[
    If[(y + v[x, y]) < h1, 0, 
     1/10^4*(0 - v[x, y])], (-w/2 < x < w/2) && 
     y >= (h1 - s/2 (1 + Cos[(2 Pi)/w x]))]};

The 1/10^4 coefficient I found by trial an error. Setting it too high or to low ruins the weak Dirichlet condition effect.

Solve the PDE:

{u1, v1} = 
  NDSolveValue[{pde == nbc + N2, dbc1, dbc2}, {u, 
    v}, {x, y} \[Element] mesh1];

This gives a message:

FindRoot::cvmit: Failed to converge to the requested accuracy or precision within 100 iterations.

My take on this is that's because the Dirichlet constraints are problematic. You can lower the "MaxIterations" (e.g. 50) and you will see that the contact does not quite match, or increase (e.g. 150) and not much improvement is seen.

Show[{
  Graphics[{Blue, Line[{{-w/2, h1}, {w/2, h1}}]}],
  mesh1["Wireframe"["MeshElement" -> "BoundaryElements"]], 
  ElementMeshDeformation[mesh1, {u1, v1}, "ScalingFactor" -> 1][
   "Wireframe"["MeshElement" -> "BoundaryElements", 
    "MeshElementStyle" -> Red]]}, ImageSize -> 500, 
 AspectRatio -> 0.3]

In the visualization I added a blue line where the contact is. A better visualization is, however,

Show[{
  Graphics[{Blue, Line[{{-w/2, h1}, {w/2, h1}}]}],
  ElementMeshDeformation[mesh1, {u1, v1}, "ScalingFactor" -> 1][
   "Wireframe"["MeshElement" -> "BoundaryElements", 
    "MeshElementStyle" -> Red]]}, ImageSize -> 500, 
 AspectRatio -> 0.3, PlotRange -> {All, {h1 - 0.01, h1 + 0.01}}, 
 Axes -> True]

We still have penetration, but that's not unheard of: see for example here. If you continue to experiment with this an find a better solution, I'd be keen to hear about it.