Discussion Forum

Hi,

The most probable reason is that you see an interpolation or extrapolation artefact. The external material is evaluated in the Gauss points during solution, so there the results are in some sense 'exact'.

In a surface plot, the results must cover the whole element, so the question arises exactly how values are determined at other locations in the element. This depends on a lot of factors: Plot type; Refinement; Averaging between elements; Whether or not Gauss point results are requested etc.

A vector that has unit length but different orientation at different locations will typically not have unit length when interpolated (How do you interpolate a vector?)

In order to be more specific, you need to provide more information: 2D/3D; exactly which expressions are you using, etc.

Hi, Henrik:

Thank you for your reply. The model is 2D. And the principle vectors used to plot the direction are extracted from the exernal material model by using the state variable. i.e. extmat1.state.vecs1 and extmat1.state.vecs2. Actually I don't know how these two variables are interpolated in the surface plot. What I was doing is just put the norm of this vector, i.e. extmat1.state.vecs1^2+extmat1.state.vecs2^2, into the Expression. I found that if I use the directon provided by Comsol, i.e. solid.sp1x and solid.sp1y, the norm of this vector is exactly one at each position. Is it possible the function used to calculate the principle direction, i.e. csext_eig(sTmatrix, vals, vecs), is not accurate enough to determine the value?

In addition, I found that the other physics, like Darcy's module, cannot recognize the state variable from the external material model, i.e. extmat1.state.vecs1. In particular, I want to use the time derivative of the state variable. I input the d(extmat1.state.S1, TIME) to the Darcy's module, but it does not play any roles to the flow process. Do I need to calculate the time derivative in the C code? If possible, how to do that?

Thanks a lot for your time and help.

Hi,

There are two distinct questions.

Time Derivatives

A state variable defined in an external material should be possible to access everywhere. However, you cannot use the built-in time derivative operator on a state variable. There are a few tricks that you can use, however.

One is to compute the time derivative yourself in the external material function, and store as an extra state. This is probably easiest.

Also, you can create another field outside of the external material, to which you map your state. For mapping, you can use a weak expression like

(my_new_S1-nojac(extmat1.state.S1))*test(my_new_S1)

As the new field is an ordinary degree of freedom, it can be differentiated with respect to time.

Plotting the Principal Vectors

Here it becomes tricky.

States as such 'live' at the Gauss points. When accessed from other locations, then it is possible that your external material is called and evaluated in that other point.

I suggest that you experiment with switching on and off averaging between elements. (Smoothing), and also to add vector plots to figure out what is going on. Try plotting individual components too.

One possible explanation is that the vector is pointing out of the plane somewhere. Then, the norm of the in-plane components is no longer 1 as assumed. Remember that the function csext_eig assumes full 3D.

Hi, Henrik:

Thank you for your reply. Now I have three questions.

1. If I calculate the time derivative in C code, how should I determine the time or time step used in the current step? There are some available input quantities in the external material model, but I donot know how to use them. There is no information about them in the help documents.
2. I try to input the weak expression you gave in the ODE, but the it did not work. Should I input the weak form in other places?
3. In terms of the time derivative of a historical maximum value, I try to use the previous solution node to record the maximum value,i.e. H_max, and then use d(H_max,TIME) to determine the time derivative. But it does not work. Should I still use the weak expresion you gave, or there are any other ways to achieve it?

Thanks a lot for your help.

1. You have to pass down and store times yourself.
2. You need to define appropriate shape functions for the fields you are mapping to.
3. It is difficult to understand how a historical max can have a time derivative. If you are interested in the time derivative at the time when the max is reached, then you can just store the time derivative at the same time as the max is reached in a new DOF.