I am trying to implement Large eddy simulation with FEM for solving air flow simulation using Smagorinsky sub-grid stress model. In order to implement this there is a parameter (filter size) in the differential equation system that has to be calculated based on the grid size: $Delta=(Delta_xDelta_yDelta_z)^{1/3} $ where $Delta_x,Delta_y,Delta_z $ is the size of the mesh cell in the x, y and z direction, respectively. Is there any possibility to implement this with the NDsolve function? I did not find any solution in the documentation of the FEM.
Tag: equations
differential equations – How to solve NDSolve dependent on two NIntegrates?
Kab(n_, z_) := NIntegrate((1 + (t*n)/(f0(z))^2), {t, 0, 1});
S4(n_, z_) := NIntegrate((t/(n*(f0(z))^2)*Exp(-2*Kab(n, z)*z)), {t, 0, 1});
Com(n_, z_) := NDSolve({(f0'(z)) == n*Exp(-2*Kab(n, z)*z)*S4(n, z), f0(0) == 1, f0'(0) == 0}, f0, {z, 0, 3});
ca(n_, z_) = Com(2, z)
Plot(Evaluate(f0(z) /. ca(n, z)), {z, 0, 3})
My code is given above. I am unable to get the plot and getting some errors.
differential equations – how we can extract the value of the solution of a PDE in a point x? (NDSolve)
Please I need your help, I calculate the solution of heat equation using methode of line
This is my code:
n = 10
grid = 1/n Range[0, n];
d1 = NDSolve`FiniteDifferenceDerivative[Derivative[1], grid];
d2 = NDSolve`FiniteDifferenceDerivative[Derivative[2], grid];
M1 = d1["DifferentiationMatrix"];
M2 = d2["DifferentiationMatrix"];
y00[t_] := Sqrt[2] Sin[Pi t];
T= 0.02;
tab = Table[u[i]
tab1 = Table[u[i], {i, Length[grid]} ];
ux = M1.tab;
uxx = M2.tab;
solu1 = D[u[1]
solu2 = D[u[n + 1]
solu3 = Table[D[u[i]
solu4 = Table[u[i][0] == y00[grid[[i]]], {i, Length[grid]}];
sol1 = NDSolve[{solu1, solu2, solu3, solu4}, tab1, {t, 0, T},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MaxPoints" -> 25}}];
Now that I have calculated the solution “sol1”. I need to approximate the value of this solution in $t = T$ in the space $ grid$ something like
h = Table [sol1 [grid[[i]],T],{i,Length[grid]}]
but this is not working for me, can anyone help me with that please??
real analysis – How to get the reverse Hölder inequality for the weak solution of elliptic equations?
MathOverflow is a question and answer site for professional mathematicians. It only takes a minute to sign up.
Sign up to join this community
Anybody can ask a question
Anybody can answer
The best answers are voted up and rise to the top
Asked
Viewed
13 times
$begingroup$
Suppose that the elliptic operator is defined as
begin{eqnarray}
L_{epsilon}=-text{div}(A(x/epsilon)triangledown)=-frac{partial}{partial x_i}left(a_{ij}^{alphabeta}left(frac{x}{epsilon}right)right).epsilon>0,
end{eqnarray}
where $ A(y)=(a_{ij}^{alphabeta}(y)) $ with $ 1leq i,jleq d $ and $ 1leqalpha,betaleq m $ ia real, bounded measurable, and satisfies ellipticity condition as follows
begin{eqnarray}
&&left|Aright|_{infty}<mu^{-1}\
&&muint_{mathbb{R}^d}|triangledown u|^2 dxleqint_{mathbb{R}^d}Atriangledown ucdottriangledown u dxtext{ for any }uin C_{0}^{infty}(mathbb{R}^d;mathbb{R}^m)
end{eqnarray}
where $ mu>0 $ is a positive constant. Let $ B=B(0,1) $ is a ball with center $ 0 $ and radius $ 1 $, If $ u_{epsilon}in H^1(2B;mathbb{R}^m) $ is a weak solution of $ L_{epsilon}(u_{epsilon})=0 $ in $ 2B $. Then there exists some $ p>2 $, depending only on $ mu $(and $ d,m $) such that
begin{eqnarray}
left(frac{1}{|B|}int_B|triangledown u_{epsilon}|^pdxright)^{1/p}leq Cleft(frac{1}{|2B|}int_{2B}|triangledown u_{epsilon}|^2dxright)^{1/2}.
end{eqnarray}
I recently meet with the reverse inequality for the elliptic equation as above. However I am very confused about the proof of this statement. In the book, it proves that
begin{eqnarray}
left(frac{1}{|B|}int_B|triangledown u_{epsilon}|^2dxright)^{1/2}leq Cleft(frac{1}{|2B|}int_{2B}|triangledown u_{epsilon}|^qdxright)^{1/q}(1)
end{eqnarray}
where $ q=frac{2d}{d+2} $. How can I get the final estimates by (1)? Can you give me some hints or references?
asked 50 mins ago
$endgroup$
MathOverflow works best with JavaScript enabled
differential equations – Solving an ODE with paramters and taking the limit of the solution
I am very new to Mathematica and already spent a lot of time trying to do this but failed.
I am trying to solve an ODE:
solution = DSolve({-((m (1 + m) + 4/(9 (-2/3 + t) t)) y(t)) +
2 (-1/3 + t) y'(t) + (-2/3 + t) t y''(
t) == (-4 (1 + C/2))/(9 (-2/3 + t) t), y(1) == 1, y'(1) == C},
y, t)
where $m$ is a nonnegative integer and $C$ is a real number. (I am not sure how to incorporate that $m$ is a nonnegative integer in the code. )
I want to show that there exists a $C$ such that the solution is $0$ at infinity.
When I try that code:
Limit(y(t) /. solution((1)), t -> Infinity, m (Element) Integers)
it just spits out the same thing.
What should I do?
(Note that I don’t need to find that value of $C$; I just need to show that for every $m$, there is a number $C$ in which the solution vanishes at infinity. )
Any reference on solving functional and differenial equations involving power towers of a function?
I am curious if this is covered in any textbooks or papers. I am interested in finding out how to solve functional and differential (even pde) equations involving power towers in them. For example if tow(f(x),n)= g(x) where tow meant a diagonal power tower of f with n determining it’s height. Also, can this be solved if some of the f’s in the tower (not necessarily all of them) were derivatives , or even partial derivatives in a multivariable case? I’m also interested in what happens if n goes to infinity. Thanks.
differential equations – Numerically Solving Fourth Order PDE with Boundary Condition at Infinity
I am trying to numerically simulate (using NDSolve) the following differential equation with the following boundary conditions:
$$partial_x^2psi +cpartial_x^4psi – b(apsi + psi^4)= lambdapartial_tpsi$$
$$partial_tpsi(x,t=infty)= 0$$
$$psi(x,t=0)=f(x)$$
where $a,b,c,lambda$ are all constants. I have been encountering a couple of problems. For example, it seems like Mathematica cannot simulate fourth order partial differential equations(?) Also when working with simpler versions of the problem, I was running into trouble with the boundary condition at $tto infty$. The code I’ve been working with is something along the lines of this:
a = 1;
b = 1;
c = 1;
l = 1;
NDSolveValue({D((Psi)(x, t), {x, 2}) + D((Psi)(x, t), {x, 4}) -
b*(a*(Psi)(x, t) + (Psi)(x, t)^3) ==
l*D((Psi)(x, t), t), (Psi)(x, 0) ==
Tanh(x), (D((Psi)(x, t), t) /. t -> 1000) == 0}, (Psi), {x, -20,
20}, {t, 0, 1000}
)
I tried making the time boundary condition some large number, rather than putting in Infinity. But this entire scheme seems to be too simple.
How to solve this kind of ordinary differential equations $t y y^{prime prime}-2 tleft(y^{prime}right)^{2}+3 y y^{prime}=0$
Often, to solve an ODE, it if enough to classify it correctly and apply correspondent well-known method. But in this case I struggle
$t y y^{prime prime}-2 tleft(y^{prime}right)^{2}+3 y y^{prime}=0$
It is non-linear homogeneous second-order ODE with variable coefficients. In textbooks, people talk about similar equations where one variable is missing (x or y or y’). In that case, we use substitution like y = z(z) or y’ = p(x). But in my ODE, I do not see what to do.
differential equations – Modified Heat Transfer in Fluid Flow
I am trying to simulate Modified Heat Transfer in Fluid Flow in MMA 12.02 (based on Buoyancy-Driven Flow in a Square Cavity ).
The modified heat transfer takes the form:
with the solid volume fraction:
The whole model is defined as:
ClearAll("Global`*")
Needs("NDSolve`FEM`")
sizes = {length -> 1, hight -> 1};
(CapitalOmega) = Rectangle({0, 0}, {length, hight}) /. sizes;
em = ToElementMesh((CapitalOmega) , MaxCellMeasure -> .0001,
"MeshOrder" -> 2);
Pr = 50;
Ra = 2.27*10^5;
timetol = 0.00001;
Th = 1.0;
Tc = -0.01;
ClearAll((Nu))
op = {
!(*SuperscriptBox((u),
TagBox(
RowBox({"(",
RowBox({"1", ",", "0", ",", "0"}), ")"}),
Derivative),
MultilineFunction->None))(t, x, y) +
Inactive(Div)((- Inactive(Grad)(u(t, x, y), {x, y})), {x,
y}) + {u(t, x, y), v(t, x, y)} .
Inactive(Grad)(u(t, x, y), {x, y}) +
!(*SuperscriptBox((p),
TagBox(
RowBox({"(",
RowBox({"0", ",", "1", ",", "0"}), ")"}),
Derivative),
MultilineFunction->None))(t, x, y),
!(*SuperscriptBox((v),
TagBox(
RowBox({"(",
RowBox({"1", ",", "0", ",", "0"}), ")"}),
Derivative),
MultilineFunction->None))(t, x, y) +
Inactive(Div)((-Inactive(Grad)(v(t, x, y), {x, y})), {x,
y}) + {u(t, x, y), v(t, x, y)} .
Inactive(Grad)(v(t, x, y), {x, y}) +
!(*SuperscriptBox((p),
TagBox(
RowBox({"(",
RowBox({"0", ",", "0", ",", "1"}), ")"}),
Derivative),
MultilineFunction->None))(t, x, y) - Ra/Pr*T(t, x, y),
!(*SuperscriptBox((u),
TagBox(
RowBox({"(",
RowBox({"0", ",", "1", ",", "0"}), ")"}),
Derivative),
MultilineFunction->None))(t, x, y) +
!(*SuperscriptBox((v),
TagBox(
RowBox({"(",
RowBox({"0", ",", "0", ",", "1"}), ")"}),
Derivative),
MultilineFunction->None))(t, x, y),
!(*SuperscriptBox((T),
TagBox(
RowBox({"(",
RowBox({"1", ",", "0", ",", "0"}), ")"}),
Derivative),
MultilineFunction->None))(t, x,
y)*(1.0 + 1/0.045*20* Sech(40* (0.01 - T(t, x, y)))^2 ) +
Inactive(Div)((-Inactive(Grad)(T(t, x, y), {x, y})), {x,
y}) + {u(t, x, y), v(t, x, y)} .
Inactive(Grad)(T(t, x, y), {x, y})};
wall = DirichletCondition({u(t, x, y) == 0, v(t, x, y) == 0}, True);
reference = DirichletCondition(p(t, x, y) == 0, x == 0 && y == 0);
temperatures = {DirichletCondition(T(t, x, y) == Th, x == 0),
DirichletCondition(T(t, x, y) == 0, x == length)};
bcs = {wall, reference, temperatures} /. sizes;
ic = {u(0, x, y) == 0, v(0, x, y) == 0, p(0, x, y) == 0,
T(0, x, y) == 0};
Monitor(AbsoluteTiming({xVel, yVel, pressure, temperature} =
NDSolveValue({op == {0, 0, 0, 0}, bcs, ic}, {u, v, p,
T}, {x, y} (Element) em, {t, 0, timetol},
Method -> {"TimeIntegration" -> {"IDA",
"MaxDifferenceOrder" -> 2},
"PDEDiscretization" -> {"MethodOfLines",
"SpatialDiscretization" -> {"FiniteElement",
"InterpolationOrder" -> {u -> 2, v -> 2, p -> 1,
T -> 2}}}},
EvaluationMonitor :> (currentTime =
Row({"t = ", CForm(t)})));), currentTime)
However, this modified simulation model produces some numerical errors:
How to set the solver parameters that we can fix the bug “Matrix SparseArray(<1395840>, {101404, 101404}) is singular …”
differential equations – Different answer between HeatTransferPDEComponent and Matlab
I’ve been working on a 2-layer HeatTransferPDE Model.
The diagram of the model is
The two materials have an initial temperature of 37℃, and the outside is 75℃. The model is for heat insulation so it can simplify to a one-dimensional HeatTransferPDE.
Mathematica codes are:
outsideTemp = 273.15 + 75;
boundary1 = 0.6/1000;
boundary2 = (0.6 + 6)/1000;
Ω1 = Line({{0}, {boundary1}});
Ω2 = Line({{boundary1}, {boundary2}});
vars = {T(t, x), t, {x}};
pars1 = <|"ThermalConductivity" -> 0.082, "AmbientTemperature" -> 273.15 + 37, "Density" -> 300, "SpecificHeatCapacity" -> 1377|>;
B1 = HeatTransferValue(x == 0, vars, <|"HeatTransferCoefficient" -> 114, "AmbientTemperature" -> outsideTemp|>);
pde1 = {HeatTransferPDEComponent(vars, pars1) == B1, T(0, x) == 273.15 + 37};
Tfun1 = NDSolveValue(pde1, T, {t, 0, 100}, {x} ∈ Ω1);
pars2 = <|"ThermalConductivity" -> 0.37, "Density" -> 862, "SpecificHeatCapacity" -> 2100, "SurfaceTemperature" -> Tfun1(t, boundary1)|>;
B2 = HeatTemperatureCondition(x == boundary1, vars, pars2);
pde2 = {HeatTransferPDEComponent(vars, pars2) == 0, B2, T(0, x) == 273.15 + 37};
Tfun2 = NDSolveValue(pde2, T, {t, 0, 100}, {x}∈ Ω2);
Plot3D({Tfun1(t, x), Tfun2(t, x)}, {t, 0, 100}, {x, 0, boundary2}, PlotRange -> All)
The output is
However, I use finite difference with same parameters and the output is
The matlab code is too long to put in here(because it’s for a 4-layer model).
I checked the code for a day and still have no clue why there is a difference between them, the material properties is in the table below.
Thanks!
