From: wolf
-
@@ -162,15 +162,135 @@ in the velocity variables:
O(h) O(h2 O(h3)
+ O(h) O(h2) O(h3)
-
The result concerning the convergence order is the same here.
O(h) O(h2 O(h3)
+ O(h) O(h2) O(h3)
+Realistic flow computations for ground water or oil reservoir simulations will +not use a constant permeability. Here's a first, rather simple way to change +this situation: we use a permeability that decays very rapidly away from a +central flowline until it hits a background value of 0.001. This is to mimick +the behavior of fluids in sandstone: in most of the domain, the sandstone is +homogenous and, while permeably to fluids, not overly so; on the other stone, +the stone has cracked, or faulted, along one line, and the fluids flow much +easier along this large crask. Here is how we could implement something like +this: +
+template <int dim>
+void
+KInverse<dim>::value_list (const std::vector<Point<dim> > &points,
+ std::vector<Tensor<2,dim> > &values) const
+{
+ Assert (points.size() == values.size(),
+ ExcDimensionMismatch (points.size(), values.size()));
+
+ for (unsigned int p=0; p<points.size(); ++p)
+ {
+ values[p].clear ();
+
+ const double distance_to_flowline
+ = std::fabs(points[p][1]-0.2*std::sin(10*points[p][0]));
+
+ const double permeability = std::max(std::exp(-(distance_to_flowline*
+ distance_to_flowline)
+ / (0.1 * 0.1)),
+ 0.001);
+
+ for (unsigned int d=0; d<dim; ++d)
+ values[p][d][d] = 1./permeability;
+ }
+}
+
+Remember that the function returns the inverse of the permeability tensor.
+
-different quadrature formula for errors
++With a significantly higher mesh resolution, we can visualize this, here with +x- and y-velocity: +
+ + + +It is obvious how fluids flow essentially only along the middle line, and not +anywhere else. + + ++Another possibility would be to use a random permeability field. A simple way +to achieve this would be to scatter a number of centers around the domain and +then use a permeability field that is the sum of (negative) exponentials for +each of these centers. Flow would then try to hop from one center of high +permeability to the next one. This is an entirely unscientific attempt at +describing a random medium, but one possibility to implement this behavior +would look like this: +
+template <int dim>
+class KInverse : public TensorFunction<2,dim>
+{
+ public:
+ KInverse ();
+
+ virtual void value_list (const std::vector<Point<dim> > &points,
+ std::vector<Tensor<2,dim> > &values) const;
+
+ private:
+ std::vector<Point<dim> > centers;
+};
+
+
+template <int dim>
+KInverse<dim>::KInverse ()
+{
+ const unsigned int N = 40;
+ centers.resize (N);
+ for (unsigned int i=0; i<N; ++i)
+ for (unsigned int d=0; d<dim; ++d)
+ centers[i][d] = 2.*rand()/RAND_MAX-1;
+}
+
+
+template <int dim>
+void
+KInverse<dim>::value_list (const std::vector<Point<dim> > &points,
+ std::vector<Tensor<2,dim> > &values) const
+{
+ Assert (points.size() == values.size(),
+ ExcDimensionMismatch (points.size(), values.size()));
+
+ for (unsigned int p=0; p<points.size(); ++p)
+ {
+ values[p].clear ();
+
+ double permeability = 0;
+ for (unsigned int i=0; i<centers.size(); ++i)
+ permeability += std::exp(-(points[p]-centers[i]).square()
+ / (0.1 * 0.1));
+
+ const double normalized_permeability
+ = std::max(permeability, 0.005);
+
+ for (unsigned int d=0; d<dim; ++d)
+ values[p][d][d] = 1./normalized_permeability;
+ }
+}
+
+
++With a permeability field like this, we would get x-velocities and pressures as +follows: +
+ +