关于非线性偏微分方程的数值求解

感谢贴吧大佬的帮助，在解决了上一次的问题后我对方程进行了一些改变，主要是加入了非线性项，代码如下
Clear["Global`*"];
pin = 100*(10^(-3)); g = 58*(10^-6); s = 100*(10^-6); w =
300*(10^-6); lamda0 = 46; d = 10*(10^-6); l = 200*(10^-6); L :=
d + l + lc; W := g + s/2 + w; T0 = 300; lc = 300*(10^-6); ka =
10*(10^7);
ds = 22*(10^-6);
\[CapitalPhi]0 := 2*pin/(ds*(g^2));
sigmod[x_, a_] := 1/(E^(-x*a) + 1);
x1[x_] := sigmod[x - (lc - g/4), ka];
x2[x_] := sigmod[x - lc, ka];
x3[x_] := sigmod[x, ka];
x4[x_] := sigmod[x - L, ka];
y1[y_] := sigmod[y - w, ka];
y2[y_] := sigmod[y - (w + g), ka];
y3[y_] := sigmod[y - W, ka];
y4[y_] := sigmod[y, ka];
ux[x_] := x1[x] - x2[x];
uy[y_] := y1[y] - y2[y];
uz[x_, y_] := ux[x]*uy[y]*\[CapitalPhi]0;
lamda[x_,
y_] := (1.5*lamda0*(x3[x] - x2[x]) +
lamda0*(x2[x] - x4[x]))*(y4[y] - y3[y]);
region = Rectangle[{0, 0}, {L, W}];
pdeSoln =
NDSolveValue[{\!\(TraditionalForm\`\((\(lamda(x, y)\)\ \*
FractionBox[
RowBox[{
SuperscriptBox["\[PartialD]", "2"],
RowBox[{"t", "(",
RowBox[{"x", ",", "y"}], ")"}]}],
RowBox[{
RowBox[{"\[PartialD]", "x"}], "\[ThinSpace]",
RowBox[{"\[PartialD]", "x"}]}],
MultilineFunction->None] + \*
FractionBox[
RowBox[{"\[PartialD]",
RowBox[{"lamda", "(",
RowBox[{"x", ",", "y"}], ")"}]}],
RowBox[{"\[PartialD]", "x"}],
MultilineFunction->None]\ \*
FractionBox[
RowBox[{"\[PartialD]",
RowBox[{"t", "(",
RowBox[{"x", ",", "y"}], ")"}]}],
RowBox[{"\[PartialD]", "x"}],
MultilineFunction->None])\) + \((\(lamda(x, y)\)\ \*
FractionBox[
RowBox[{
SuperscriptBox["\[PartialD]", "2"],
RowBox[{"t", "(",
RowBox[{"x", ",", "y"}], ")"}]}],
RowBox[{
RowBox[{"\[PartialD]", "y"}], "\[ThinSpace]",
RowBox[{"\[PartialD]", "y"}]}],
MultilineFunction->None] + \*
FractionBox[
RowBox[{"\[PartialD]",
RowBox[{"lamda", "(",
RowBox[{"x", ",", "y"}], ")"}]}],
RowBox[{"\[PartialD]", "y"}],
MultilineFunction->None]\ \*
FractionBox[
RowBox[{"\[PartialD]",
RowBox[{"t", "(",
RowBox[{"x", ",", "y"}], ")"}]}],
RowBox[{"\[PartialD]", "y"}],
MultilineFunction->None])\)\) + uz[x, y] ==
NeumannValue[0, y == W] + NeumannValue[0, x == 0],
{DirichletCondition[t[x, y] == T0, x == L],
DirichletCondition[t[x, y] == T0, y == 0]}},
t, {x, y} \[Element] region];
ContourPlot[pdeSoln[x, y], {x, y} \[Element] region,
Mesh -> None, ColorFunction -> "TemperatureMap",
Contours -> 30, AspectRatio -> Automatic,
FrameTicksStyle -> Directive[Black, 18]]
Plot3D[pdeSoln[x, y], {x, y} \[Element] region]
但是求解出现问题，如下
NDSolveValue::femdpop: The NDSolve`FEM`FEMStiffnessElements operator failed.
请问各位该如何解决