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

求助：MMA有限元方法双变量耦合方程组怎么设置诺伊曼边界条件

\[CapitalOmega] = ImplicitRegion[True, {{z, 0, 200}}];
\[CapitalGamma]1 = NeumannValue[0 , z == 0 || z == 200];
TraditionalForm[\[CapitalGamma]1]
{nO, pO} = 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] + 1 -
E^(-2 u) + p[t, z] - n[t, z]) n[t, z] + \[CapitalGamma]1,
D[p[t, z], t] ==
0.1 D[p[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]
- 0.1 Ldi e^2 nb Ldi/(k T \[Epsilon] \[Epsilon]b) (0 +
eF0O[z]) D[p[t, z] , z] -
0.1 Ldi e^2 nb Ldi/(k T \[Epsilon] \[Epsilon]b) (p0[z] -
n0[z] + 1 - E^(-2 u) + p[t, z] - n[t, z]) p[t,
z] + \[CapitalGamma]1,
n[-10 \[Tau]L/\[Tau], z] == 0,
p[-10 \[Tau]L/\[Tau], z] == 0},
{n, p},
{t, -10 \[Tau]L/\[Tau], 5},
{z} \[Element] \[CapitalOmega],
Method -> {"PDEDiscretization" -> {"MethodOfLines",
"SpatialDiscretization" -> {"FiniteElement",
"MeshOptions" -> {MaxCellMeasure -> 0.1}}}}]
Manipulate[
Plot[QuantityMagnitude[nb, "Centimeters^-3"] nO[t0/n\[Tau],
z0/nLdi], {z0, 0, 200 nLdi}, PlotRange -> All], {t0, -10 n\[Tau]L,
5 n\[Tau]}]

以上是源代码，我想求解的方程包含两个偏微分方程的耦合， 分别是n[t,z]和p[t,z]，想要分别设置诺伊曼边界条件
dn/dx==0,x==0||x==200,dp/dx==0,x==0||x==200，请问格式要怎么写？我按照图中的写法，运算一直卡住算不出来。

另外如果希望设置参数是随时间变化的话要怎么写？比如把eF0[z]从一个已知的函数换成一个由n[t,z]和p[t,z]决定的函数，比如eF[t,z]满足
d eF[t,z]/dz=p[t,z]-n[t,z]
eF[t,0]==0
我试过直接把这些加紧NDSolve函数里面，但是会是别成三个偏微分方程。


<< 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页
分享到: