0

我尝试在 MATLAB 中解决具有节点热源的四面体有限元上的热扩散问题，这取决于解向量。非线性方程组如下所示：

B U' + A U = q(T)

其中 B 是热容量矩阵，A 是传导率矩阵，q 是源项，U 是温度。我使用带有 Picard 迭代和时间步控制的 Adams-Bashforth/Trapezoid Rule 预测器-校正器方案。源项的温度在最后一个时间步的温度和预测变量的温度之间精确评估。这是预测器-校正器代码的简化版本。源的计算是一个函数。

    % predictor
    K0 = t(n)-t(n-1);              
    Upre(dirichlet) = u_d_t(coordinates(dirichlet,:));
    Upre(FreeNodes) = U(FreeNodes,n) + (dt/2)*((2+dt/K0)*U_dt(FreeNodes,3) - (dt/K0)*U_dt(FreeNodes,1));      % predictor step
    Uguess = Upre;     % save as initial guess for Picard iteration

    % corrector with picard iteration
    while res >= picard_tolerance 
        T_theta = Uguess*theta + (1-theta)*U(:,n);
        b = q(T_theta);
        % Building right-hand side vector (without Dirichlet boundary conditions yet)
        rhs = ((2/dt)*B*U(:,n) + B*U_dt(:,1))+b;

        % Applying Dirichlet Boundary Conditions to the Solution Vector
        Ucor(dirichlet) = u_d_t(coordinates(dirichlet,:));
        rhs = rhs - ((2/dt)*B+A)*Ucor;

        % Solving the linearized system using the backslash operator
        % P*U(n+1) = f(Un) => U(n+1) = P\f(Un)
        Ucor(FreeNodes) = ((2/dt)*B(FreeNodes,FreeNodes)+A(FreeNodes,FreeNodes))\rhs(FreeNodes);
        res = norm(Uguess-Ucor);
        Uguess = Ucor;
        U(:,n+1) = Ucor;
    end

如您所见，我使用反斜杠运算符来解决系统问题。系统的非线性应该不会太差。然而，随着时间步长的增加，picard 方法收敛得更慢，最终完全停止收敛。不过，我需要更大的时间步长，所以我将整个校正步骤放入一个函数中，并尝试用 fsolve 来解决它，而不是看看我是否实现了更快的收敛。不幸的是 fsolve 似乎从来没有完成第一个时间步骤。我想我没有正确配置 fsolve 的选项。谁能告诉我，如何为大型稀疏非线性系统配置 fsolve（我们正在谈论数千到上万个方程）。或者对于这个问题，是否有比 fsolve 更好的解决方案？帮助和 - 因为我不是专家或计算工程师 - 非常感谢明确的建议！

4

1 回答 1

0

根据我的经验，非线性方程是通过线性化来求解温度增量并使用 Newton Raphson 求解器等迭代收敛来求解的。因此，如果您使用的是隐式积分模式，则您有一个外部时间步长解决方案，其中一个内部非线性解决方案用于时间步长上的温度步长。

于 2014-12-08T20:49:17.590 回答