2D inhomogeneous biharmonic equation with wedged edge

Question

I'm solving bending of rectangular plate while, boundary conditions are

I have found similar problem solved: datavoreconsulting.com/programming-tips/numerically-solving-pdes-mathematica-finite-differences/ and used this method to solve my problem

Define Functions Make the domain into a grid

\[Nu] = 0.3;
ee = 210*10^9;
p = -19000
h = 7;
DD = (ee h^3)/(12 (1 - \[Nu]^2));
a = 1000;
b = 2 a;
xmin = -(1/2) a;
xmax = 1/2 a;
ymin = -(1/2) b;
ymax = 1/2 b;
xdivisions = 10;
ydivisions = 10;
dx = (xmax - xmin)/xdivisions
dy = (ymax - ymin)/ydivisions
q0 = p;
xgrid = Range[xmin, xmax, dx];
ygrid = Range[ymin, ymax, dy];
grid = Outer[{#1, #2} &, xgrid, ygrid];

Make an array of values for the solution. Each entry in the array corresponds to the value of the unknown function at a point in space

W = Array[w, {xdivisions + 1, ydivisions + 1}, {{xmin, xmax}, {ymin, ymax}}];

-19000

100

200

Specify boundary condition W=0

wl[y_] = 0;
wr[y_] = 0;
wb[x_] = 0;
wt[x_] = 0;

Find finite difference appoximations for different derivatives of u

dwdx = NDSolve`FiniteDifferenceDerivative[{1, 0}, {xgrid, ygrid}, W];
dwdx2 = NDSolve`FiniteDifferenceDerivative[{2, 0}, {xgrid, ygrid}, W];
dwdy = NDSolve`FiniteDifferenceDerivative[{0, 1}, {xgrid, ygrid}, W];
dwdy2 = NDSolve`FiniteDifferenceDerivative[{0, 2}, {xgrid, ygrid}, W];
dwdxdy = NDSolve`FiniteDifferenceDerivative[{1, 1}, {xgrid, ygrid}, W];
dwdx4 = NDSolve`FiniteDifferenceDerivative[{4, 0}, {xgrid, ygrid}, W];
dwdy4 = NDSolve`FiniteDifferenceDerivative[{0, 4}, {xgrid, ygrid}, W];
dwdx2dy2 = NDSolve`FiniteDifferenceDerivative[{2, 2}, {xgrid, ygrid}, W];

Boundary Conditions

(*Consider the left-side bc *)
leftbc = Table[dwdx[[1, k + 1]], {k, 0, Length[ygrid] - 1}];
leftbcw = Table[W[[1, k + 1]], {k, 0, Length[ygrid] - 1}];

(*Do the same for the right-side bc*)
rightbc = Table[dwdx[[-1, k + 1]], {k, 0, Length[ygrid] - 1}];
rightbcw = Table[W[[-1, k + 1]], {k, 0, Length[ygrid] - 1}];

(*And for the bottom side boundary condition*)
bottombc = Table[dwdy[[k + 1, 1]], {k, 0, Length[xgrid] - 1}];
bottombcw = Table[W[[k + 1, 1]], {k, 0, Length[xgrid] - 1}];

(*And for the top side bc*)
topbc = Table[dwdy2[[k + 1, -1]], {k, 0, Length[xgrid] - 1}];
topbcw = Table[W[[k + 1, -1]], {k, 0, Length[xgrid] - 1}];

Takes the left boundary conditions, and solves them to yield the values of \ u along the left side

   wleft = NSolve[Map[# == 0 &, leftbc], Table[W[[1, k]], {k, 1, Length[ygrid]}]];
    wleftw = NSolve[Map[# == 0 &, leftbcw], 
       Table[W[[1, k]], {k, 1, Length[ygrid]}]];

takes the right boundary conditions, and solves them to yield the values of \ u along the right side

wright = NSolve[Map[# == 0 &, rightbc], 
   Table[W[[-1, k]], {k, 1, Length[ygrid]}]];
wrightw = NSolve[Map[# == 0 &, rightbcw], 
   Table[W[[-1, k]], {k, 1, Length[ygrid]}]];

takes the bottom boundary conditions, and solves them to yield the values \ of u along the bottom side

wbottom = NSolve[Map[# == 0 &, bottombc], 
   Table[W[[k, 1]], {k, 1, Length[xgrid]}]];
wbottomw = 
  NSolve[Map[# == 0 &, bottombcw], Table[W[[k, 1]], {k, 1, Length[xgrid]}]];

takes the top boundary conditions, and solves them to yield the values of u along the top side

 wtop = NSolve[Map[# == 0 &, topbc], Table[W[[k, -1]], {k, 1, Length[xgrid]}]];
    wtopw = NSolve[Map[# == 0 &, topbcw], 
       Table[W[[k, -1]], {k, 1, Length[xgrid]}]];

To make sure boundary conditions are consitent, we only want one u value \ per point

wleft1 = wleft[[1, 1 ;; Length[ygrid]]];
wtop1 = wtop[[1, 2 ;; Length[xgrid] - 1]];
wright1 = wright[[1, 1 ;; Length[ygrid]]];
wbottom1 = wbottom[[1, 2 ;; Length[xgrid] - 1]];
wleft2 = wleftw[[1, 1 ;; Length[ygrid]]];
wtop2 = wtopw[[1, 2 ;; Length[xgrid] - 1]];
wright2 = wrightw[[1, 1 ;; Length[ygrid]]];
wbottom2 = wbottomw[[1, 2 ;; Length[xgrid] - 1]];
(*This is a list of all the boundary values of u*)
boundary1 = Join[wleft1, wright1, wtop1, wbottom1];
boundary2 = Join[wleft2, wright2, wtop2, wbottom2];


boundary = boundary2 /. boundary1;

Solving

We now create a set of equations using the PDE. We make a table, with each entry corresponding to an interior grid point. Each entry in the table becomes an equation, from the discretized PDE. We use our knowledge of the boundary conditions to eliminate the the values of u on the boundary

equations = 
  Map[(q0/DD == #) &, 
   Flatten[Table[
       dwdx4[[i, j]] + 2*dwdx2dy2[[i, j]] + dwdy4[[i, j]], {i, 2, 
        Length[xgrid] - 1}, {j, 2, Length[ygrid] - 1}], 1] /. boundary /. 
    boundary2];

(*We solve these equations, for the interior values of u*)
intSol = NSolve[equations, 
    Flatten[W[[2 ;; Length[xgrid] - 1, 2 ;; Length[ygrid] - 1]]]][[1]];

We substitute the interior values of u into the equations that determine the boundary values of u

boundarySol = boundary2 /. intSol;

solutionArray is an array of the values of u. I'm not sure why I had to do the last replacement rule twice!

solutionArray = ((W /. intSol) /. boundarySol) /. boundarySol;

Use this to make a table of 3d coordinates

dataPoints = 
  Table[{xmin + i*dx, ymin + j*dy, solutionArray[[i + 1, j + 1]]}, {i, 0, 
    Length[xgrid] - 1}, {j, 0, Length[ygrid] - 1}];
Min[solutionArray]

-90.0146

Plot the solution

ListPlot3D[Flatten[dataPoints, 1], 
 AxesLabel -> {Style[x, Medium, Blue], Style[y, Medium, Blue], 
   Style[u, Medium, Blue]}, PlotRange -> All]

but solution differs from FEM softvare solution and some other mathematica solution, my code gives max deflection 90.0146 but it should be around 1.5.

Am I applying boundary conditions wrong way or does mathematica method give huge error?

I have solved the problem by programming my own finite difference parameters.

new finite difference approximation functions:
dwdx = Table[(W[[i + 1, j]] - W[[i - 1, j]])/(2*dx), {i, 3, 
    Length[xgrid] - 2}, {j, 3, Length[ygrid] - 2}];
dwdy = Table[(W[[i, j + 1]] - W[[i, j - 1]])/(2*dy), {i, 3, 
    Length[xgrid] - 2}, {j, 3, Length[ygrid] - 2}];
dwdx2 = Table[(W[[i + 1, j]] - 2 W[[i, j]] + 
      W[[i - 1, j]])/(dx^2), {i, 3, Length[xgrid] - 2}, {j, 3, 
    Length[ygrid] - 2}];
dwdy2 = Table[(W[[i, j + 1]] - 2 W[[i, j]] + 
      W[[i, j - 1]])/(dy^2), {i, 3, Length[xgrid] - 2}, {j, 3, 
    Length[ygrid] - 2}];
dwdx4 = Table[(W[[i + 2, j]] - 4 W[[i + 1, j]] + 6 W[[i, j]] - 
      4 W[[i - 1, j]] + W[[i - 2, j]])/(dx^4), {i, 3, 
    Length[xgrid] - 2}, {j, 3, Length[ygrid] - 2}];
dwdy4 = Table[(W[[i, j + 2]] - 4 W[[i, j + 1]] + 6 W[[i, j]] - 
      4 W[[i, j - 1]] + W[[i, j - 2]])/(dy^4), {i, 3, 
    Length[xgrid] - 2}, {j, 3, Length[ygrid] - 2}];
dwdx2dy2 = 
  Table[(W[[i + 1, j + 1]] - 2 W[[i + 1, j]] + W[[i + 1, j - 1]] - 
      2 W[[i, j + 1]] + 4 W[[i, j]] - 2 W[[i, j - 1]] + 
      W[[i - 1, j + 1]] - 2 W[[i - 1, j]] + W[[i - 1, j - 1]])/(dx^2*
      dy^2), {i, 3, Length[xgrid] - 2}, {j, 3, Length[ygrid] - 2}];

Now I have another problem to solve, bending of a plate with hole

Answer 1

I haven't check your code, but here's my solution based on my pdetoae. You may feel it confusing at first glance, but I believe you'll be able to understand it if you settle down and read through it. The idea is simple: generate the difference equations automatically with ptoafunc, then remove those redundant equations with del, extract the coefficients of the linear equations system with CoefficientArrays and solve it with LinearSolve (you can also solve the equation system with Solve directly, but it'll be much slower), and built a InterpolatingFunction with ListInterpolation.

ν = 3/10;
Ε = 210*10^9;
qz = -19000;
h = 7;
\[ScriptCapitalD] = (Ε h^3)/(12 (1 - ν^2));
a = 1000;
b = 2 a;

lap = Laplacian[#, {x, y}] &;

With[{w = w[x, y]},
 eq = lap@lap@w == qz/\[ScriptCapitalD];
 (*bc={
 {w == 0,D[w,x,x] == 0}/.{{x -> 0},{x -> a}},
 {{w == 0,D[w,y,y] == 0}/.y -> 0,{w == 0,D[w,y] == 0}/.y -> b}}*)
 bc = {
   {w == 0, D[w, x] == 0} /. {{x -> 0}, {x -> a}},
   {{w == 0, D[w, y] == 0} /. y -> 0, {w == 0, D[w, y, y] == 0} /. y -> b}};]

points = 50;
{grid@x, grid@y} = Array[# &, points, #] & /@ {{0, a}, {0, b}};

difforder = 4;
ptoafunc = pdetoae[w[x, y], grid /@ {x, y}, difforder];
del = #[[3 ;; -3]] &;

ae = del /@ del@ptoafunc@eq;
aebc = MapAt[del, ptoafunc@bc, {1, All, All}];

var = Outer[w, grid@x, grid@y] // Flatten;
{blst, mat} = CoefficientArrays[{ae, aebc} // Flatten, var];
sollst = LinearSolve[mat, -N@blst];    
sol = ListInterpolation[Partition[sollst, points], {grid@x, grid@y}];

Plot3D[sol[x, y], {x, 0, a}, {y, 0, b}]

Please notice the b.c.s in your first and last picture are different, and it's not clear to me which boundary is wedged in the last case so I chose one arbitrarily. Anyway the result agrees well with the one in the picture:

NMinimize[{sol[x, y], 0 < x < a, 0 < y < b}, {x, y}]
(* {-7.43643, {x -> 500., y -> 1086.09}} *)

---

Update

If you feel confused about del, the following is an alternative method without it:

fullsys = ptoafunc@{eq, bc};
{blst, mat} = CoefficientArrays[fullsys // Flatten, var];
sollst = LeastSquares[mat, -N@blst, Method -> Direct]; // AbsoluteTiming
sol = ListInterpolation[Partition[sollst, points], {grid@x, grid@y}];