网页资讯视频图片知道文库贴吧地图采购
进入贴吧全吧搜索
mathematica吧 关注：19,734贴子：73,657
  • 1回复贴，共1页

求助：用有限元方法NeumannValue边界条件无法生效

\[CapitalOmega] = ImplicitRegion[True, {{z, 0, 200}}];
\[CapitalGamma]2 = NeumannValue[0., x == 0 || x == 200]
s = NDSolveValue[{D[n[t, z], t] -
D[n[t, z], {z, 2}] - \[Tau]/
nb \[Alpha] (1 -
R) Pdensity/(Sqrt[2 \[Pi] ] \[Tau]L hv ) gauss[t,
0, \[Tau]L/\[Tau]] Exp[-\[Alpha] Ldi z] == -Ldi e^2 nb Ldi/(k \
T \[Epsilon] \[Epsilon]b) (0 + eF0O[z]) D[n[t, z] , z] -
Ldi e^2 nb Ldi/(k T \[Epsilon] \[Epsilon]b) (p0[z] - n0[z] 0 +
1 - E^(-2 u)) n[t, z] + \[CapitalGamma]2,
n[-10 \[Tau]L/\[Tau], z] == 0
},
n,
{t, -10 \[Tau]L/\[Tau], 5},
{z} \[Element] \[CapitalOmega],
PrecisionGoal -> 35,
Method -> {"PDEDiscretization" -> {Automatic,
"SpatialDiscretization" -> {"FiniteElement",
"MeshOptions" -> {MaxCellMeasure -> 0.001}}}}]
Plot[QuantityMagnitude[nb, "Centimeters^-3"] s[t0/n\[Tau],
z/nLdi], {z, 0, 200 nLdi }, PlotRange -> All], {t0, -10 n\[Tau]L,
5 n\[Tau]}]
Manipulate[
Plot[QuantityMagnitude[nb, "Centimeters^-3"] ND[
s[t0/n\[Tau], z/nLdi], z, z0], {z0, 0, 200 nLdi },
PlotRange -> All], {t0, -10 n\[Tau]L, 5 n\[Tau]}]
代码如上所示，一维问题，NeumannValue边界条件设置了两侧边界导数值为0，但实际求解出来的两侧导数值并不为0，似乎是方法自己生成的某种边界条件。报错的提示显示NeumannValue函数无法正确识别边界。第一张图是报错的提示，第二张图是解的导数的图，可以看到两侧边界导数并不为0，第三章是代码图。求大佬帮忙
Manipulate[



二楼补充一下，代码中出现的一些函数和常数引用自之前的计算结果，本楼把相关代码一并贴出，需要先运行本楼的代码，再运行上述代码
<< NumericalCalculus`
Needs["NDSolve`FEM`"] ;
gauss[x_, x0_, \[Sigma]_] := Exp[-(x - x0)^2/(2 \[Sigma]^2)];
gauss[x_, x0_, \[Sigma]_] :=
Piecewise[{{Exp[-(x - x0)^2/(2 \[Sigma]^2)],
Abs[x - x0] <= 10 \[Sigma]}, {0, Abs[x - x0] > 10 \[Sigma]}}];
\[Epsilon] = Quantity[1, "ElectricConstant"];
e = Quantity[1, "ElementaryCharge"];
k = Quantity[1, "BoltzmannConstant"];
T = Quantity[300, "Kelvins"];
\[Epsilon]b = 7.1;
ni = Quantity[2 10^11, "Centimeters^-3"];
Wfi = Quantity[0.33, "eV"];
u = Wfi/(k T);
nb = ni E^u;
nnb = QuantityMagnitude[nb, "Centimeters^-3"]
Ldi = Sqrt[(\[Epsilon] \[Epsilon]b k T)/(
e^2 nb) ] (*先点一下下一行再点会本行，然后光标放在右侧方框处等待，才会跳出建议栏目*)
nLdi = QuantityMagnitude[Ldi, "Nanometers"];
nVs = -0.42;
Vs = Quantity[nVs, "Volts"](*单位是电压不是电子伏特*);
vs = Vs e/(k T)
F = Function[{x}, Sqrt[2 (Cosh[u + x]/Cosh[u] - x Tanh[u] - 1)]];
Pdensity = Quantity[2.09, "Millijoules/Centimeters^2"];
\[Lambda] = Quantity[618 3/2, "Nanometers"];
\[Alpha] = 1/\[Lambda];
R = 0.426;
\[Tau]L = Quantity[300, "Femtoseconds"];
n\[Tau]L = QuantityMagnitude[\[Tau]L, "Picoseconds"];
\[Mu] = Quantity[0.202, "Centimeters^2 /Volts Seconds"];
dif = k T \[Mu]/e;
\[Tau] = Ldi^2/dif;
n\[Tau] = QuantityMagnitude[\[Tau], "Picoseconds"]
hv = Quantity[1.55, "Electronvolts"];
耗尽区方程
vO = NDSolveValue[{v'[z] == F[v[z]], v[0] == vs}, v, {z, 0, 200},
WorkingPrecision -> 35, MaxSteps -> 10^6, AccuracyGoal -> 30]
n0[z_] := E^vO[z]
Plot[QuantityMagnitude[ni, "Centimeters^-3"] E^u n0[z/nLdi], {z, 0,
20 nLdi}, AxesLabel -> {"z(nm)", "n(cm^-3)"}, PlotRange -> All]
p0[z_] := E^(-2 u) E^(-vO[z])
Plot[QuantityMagnitude[ni, "Centimeters^-3"] E^u p0[z/nLdi], {z, 0,
20 nLdi}, AxesLabel -> {"z(nm)", "n(cm^-3)"}, PlotRange -> All]
eF0O = NDSolveValue[{D[eF0[z], z] == -n0[z] + p0[z] + 1 - E^(-2 u),
eF0[200] == 0}, eF0, {z, 0, 200}, WorkingPrecision -> 35,
MaxSteps -> 10^6, AccuracyGoal -> 30]
Plot[QuantityMagnitude[e nb Ldi/(\[Epsilon] \[Epsilon]b), "V/m"] eF0O[
z/nLdi], {z, 0, 20 nLdi}, AxesLabel -> {"z(nm)", "E(V/m)"},
PlotRange -> All]


扫二维码下载贴吧客户端

下载贴吧APP
看高清直播、视频！
  • 1回复贴，共1页
分享到: